2016-11-18 290 views
1

我有兩個GSL矩陣,ATA矩陣乘法使用GSL

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

我想乘這些矩陣並將結果存儲在矩陣ATA,但我無法理解的GSL BLAS文檔。

我初始化了ATA矩陣:

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

我看到,我可以使用gsl_blas_zgemm繁衍複雜的矩陣,但是這些矩陣並不複雜。我會怎麼做呢?

UPDATE

我曾嘗試:

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

它引起以下錯誤:

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

回答

1

我的矩陣相乘倒退。正確的行是:

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