CODEWITHSUNDEEP
cmatrixmatrix-multiplicationblasgsl
Sara Fuerst
Sara Fuerst

Reputation: 6098

Matrix Multiplication Using GSL

I have two GSL matrices, AT and A: 

gsl_matrix *A;    /* coefficient matrix A     */
gsl_matrix *AT;   /* coefficient matrix A'    */

AT   = gsl_matrix_alloc(nc, nr); /* Data matrix */
A    = gsl_matrix_alloc(nr, nc); /* Data matrix */

/* Initialize A */
for(i = 0; i < nr; i++){
  gsl_matrix_set(A, i, 0, 1.0); 
}

for(i = 0; i < nr; i++){
  for(j = 1; j < nc; j++){
    gsl_matrix_set(A, i, j, 1.0 / (double)(i + j + 1)); 
  }
}

gsl_matrix_transpose_memcpy(AT, A);

I would like to multiply these matrices and store the result in the matrix ATA, but I am having trouble understanding the GSL BLAS documentation.

I have initialized the ATA matrix: 

gsl_matrix *ATA;  /* coefficient matrix A'A   */
ATA  = gsl_matrix_alloc(nc, nc); /* Data matrix */

I see that I can use gsl_blas_zgemm to multiply complex matrices, but these matrices are not complex. How would I go about this?

UPDATE

I have tried: 

gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A, AT, 0.0, ATA);

It's causing the following error:

gsl: blas.c:1354: ERROR: invalid length Default GSL error handler invoked. Aborted (core dumped)

Upvotes: 3

Views: 6452

Answers (1)

Sara Fuerst
Sara Fuerst

Reputation: 6098

I was multiplying the matrices backwards. The proper line is: 

gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, AT, A, 0.0, ATA);

Upvotes: 5

Related Questions