Außenprodukt unter Verwendung von CBLAs

Question

Ich habe Probleme, CBLAs zu verwenden, um ein äußeres Produkt auszuführen. Mein Code ist wie folgt:

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

Ich bekomme:

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 

Aber die richtige Antwort finden Sie hier:https://www.sharcnet.ca/help/index.php/blas_and_cblas_usage_and_examples

Ich habe gesehen: Effiziente Berechnung von Kronecker -Produkten in C.

Aber es hilft mir nicht, weil sie nicht wirklich sagen, wie man DGEMM nutzt, um dies tatsächlich zu tun ...

Irgendeine Hilfe? Was mache ich hier falsch?

Answer 1

Sie können es mit DGEMM tun, aber es wäre stilistisch korrekter, DGER zu verwenden, was eine dedizierte Implementierung im Außenprodukt ist. Als solches ist es etwas einfacher, richtig zu verwenden:

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 tut Habe das schöne Feature, dass das Passieren \beta = 0 Initialisiert die Ergebnismatrix für Sie, die Sie vor dem Anruf vor dem Anruf explizit auf Null ausziehen muss. @Artem Shinkarovs Antwort enthält eine schöne Beschreibung, wie man DGEMM verwendet.

Answer 2

Die Schnittstellen sind in BLAS nicht sehr bequem, versuchen wir es jedoch herauszufinden. Nehmen wir zunächst an, dass alle unsere Matrizen in Rowmajor sind. Jetzt haben wir das folgende Setup

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

Jetzt müssen wir nur den Anruf entsprechend der Dokumentation, die in Bezug auf angegeben ist

C = \alpha A*B + \beta C

Also bekommen wir:

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

Das sollte also den Job machen, den ich hoffen würde. Hier finden Sie eine Beschreibung der Parameter von DGEMM: Verknüpfung