CODEWITHSUNDEEP
clinear-algebrascientific-computingnumerical
The Dude
The Dude

Reputation: 681

Outer product using CBLAS

I am having trouble utilizing CBLAS to perform an Outer Product. My code is as follows:

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

I get:

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

But the correct answer is found here: https://www.sharcnet.ca/help/index.php/BLAS_and_CBLAS_Usage_and_Examples

I have seen: Efficient computation of kronecker products in C

But, it doesn't help me because they don't actually say how to utilize dgemm to actually do this...

Any help? What am I doing wrong here?

Upvotes: 4

Views: 3667

Answers (2)

Stephen Canon
Stephen Canon

Reputation: 106167

You can do it with dgemm, but it would be more stylistically correct to use dger, which is a dedicated outer-product implementation. As such it's somewhat easier to use correctly:

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 does have the nice feature that passing \beta = 0 initializes the result matrix for you, which saves you from needing to explicitly zero it out yourself before the call. @Artem Shinkarov’s answer provides a nice description of how to use dgemm.

Upvotes: 11

Artem Shinkarov
Artem Shinkarov

Reputation: 379

The interfaces are not very convenient in BLAS, however, let's try to figure it out. First of all, let's say that all our matrices are in RowMajor. Now we have the following set-up

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

Now, we just need to fill the call according to the documentation, which is specified in terms of 

C = \alpha A*B + \beta C

So we get:

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

So that should do the job I would hope. Here is a description of the parameters of dgemm: Link

Upvotes: 6

Related Questions