Producto exterior que usa CBLAS

Question

Tengo problemas para utilizar CBLA para realizar un producto exterior. Mi código es el siguiente:

//===SET UP===//
double x1[] = {1,2,3,4};
double x2[] = {1,2,3};
int dx1 = 4;
int dx2 = 3;
double X[dx1 * dx2];
for (int i = 0; i < (dx1*dx2); i++) {X[i] = 0.0;}

//===DO THE OUTER PRODUCT===//
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, dx1, dx2, 1, 1.0, x1, dx1, x2, 1, 0.0, X, dx1);

//===PRINT THE RESULTS===//
printf("\nMatrix X (%d x %d) = x1 (*) x2 is:\n", dx1, dx2);
for (i=0; i<4; i++) {
    for (j=0; j<3; j++) {
        printf ("%lf ", X[j+i*3]);
    }
    printf ("\n");
}

Yo obtengo:

Matrix X (4 x 3) = x1 (*) x2 is:
1.000000 2.000000 3.000000 
0.000000 -1.000000 -2.000000 
-3.000000 0.000000 7.000000 
14.000000 21.000000 0.000000 

Pero la respuesta correcta se encuentra aquí:https://www.sharcnet.ca/help/index.php/blas_and_cblas_usage_and_examples

He visto: Cálculo eficiente de productos de Kronecker en C

Pero, no me ayuda porque en realidad no dicen cómo utilizar DGEMM para hacer esto ...

¿Alguna ayuda? ¿Qué estoy haciendo mal aquí?

Answer 1

Puede hacerlo con DGEMM, pero sería más estilísticamente correcto usar DGER, que es una implementación dedicada de productos externos. Como tal, es algo más fácil de usar correctamente:

cblas_dger(CblasRowMajor, /* you’re using row-major storage */
           dx1,           /* the matrix X has dx1 rows ...  */
           dx2,           /*  ... and dx2 columns.          */
           1.0,           /* scale factor to apply to x1x2' */
           x1,
           1,             /* stride between elements of x1. */
           x2,
           1,             /* stride between elements of x2. */
           X,
           dx2);          /* leading dimension of matrix X. */

dgemm lo hace tener la buena característica que pase \beta = 0 Inicializa la matriz de resultados para usted, lo que le ahorra de la necesidad de cero explícitamente antes de la llamada. La respuesta de @artem shinkarov proporciona una buena descripción de cómo usar DGEMM.

Answer 2

Las interfaces no son muy convenientes en BLAS, sin embargo, intentemos resolverlo. En primer lugar, digamos que todas nuestras matrices están en Rowmajor. Ahora tenemos la siguiente configuración

     row  col
x1:  dx1   1   (A)
x2:   1   dx2  (B)
 X:  dx1  dx2  (C)

Ahora, solo necesitamos completar la llamada de acuerdo con la documentación, que se especifica en términos de

C = \alpha A*B + \beta C

Entonces tenemos:

cblas_dgemm (CblasRowMajor, CblasNoTrans, CblasNoTrans,
             (int)dx1, /* rows in A         */
             (int)dx2, /* columns in B      */
             (int)1,   /* columns in A      */
             1.0, x1,  /* \alpha, A itself  */
             (int)1,   /* Colums in A       */
             x2,       /* B itself          */
             (int)dx2, /* Columns in B      */
             0.0, X,   /* \beta, C itself   */
             (int)dx2  /* Columns in C  */);

Entonces eso debería hacer el trabajo que espero. Aquí hay una descripción de los parámetros de DGEMM: Enlace