Produto externo usando CBLAs

Question

Estou tendo problemas para utilizar CBLAs para realizar um produto externo. Meu código é o seguinte:

//===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");
}

Eu recebo:

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 

Mas a resposta correta é encontrada aqui:https://www.sharcnet.ca/help/index.php/blas_and_cblas_usage_and_examples

Eu tenho visto: Computação eficiente de produtos Kronecker em C

Mas, isso não me ajuda, porque eles realmente não dizem como utilizar o DGEMM para realmente fazer isso ...

Qualquer ajuda? O que eu estou fazendo errado aqui?

Answer 1

Você pode fazer isso com o DGEMM, mas seria mais estilisticamente correto usar o DGER, que é uma implementação dedicada de produtos externos. Como tal, é um pouco mais fácil de usar corretamente:

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 faz tem o bom recurso que passa \beta = 0 Inicializa a matriz de resultados para você, o que evita que você precise se concentrar explicitamente antes da chamada. A resposta de @artem shinkarov fornece uma boa descrição de como usar o DGEMM.

Answer 2

As interfaces não são muito convenientes nas blas, no entanto, vamos tentar descobrir. Primeiro de tudo, digamos que todas as nossas matrizes estão em Rowmajor. Agora temos a seguinte configuração

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

Agora, só precisamos preencher a chamada de acordo com a documentação, que é especificada em termos de

C = \alpha A*B + \beta C

Então nós conseguimos:

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  */);

Então isso deve fazer o trabalho que eu espero. Aqui está uma descrição dos parâmetros do DGEMM: Link