BLAS矩陣矢量乘法與矢量矩陣乘法。一個工作;另一個失敗

Question

在使用Lapack中的BLAS lvl 2函數時，我成功實現了矩陣向量乘法，但是當我嘗試轉置時，我得到了錯誤的答案。你能指導我解決我的錯誤嗎？ （實際上，我使用C包裝，不FORTRAN）

我試圖

| 4+i 3 | | 3+2i |   | 4+i 3 |^T | 3+2i | 
| 14+3i 2 | * | 2 | (AND) | 14+3i 2 | * | 2 | 

需要明確的是，第一個成功。第二個輸出不正確。

/* config variables */ 
char normal = 'N'; 
char transpose = 'T'; 
integer m = 2; 
complex alpha = {r:1,i:0}; 
complex beta = {r:0,i:0}; 
integer one = 1; 
/* data buffers */ 
complex a[4] = {(complex){r:4, i:1},(complex){r:14, i:3},(complex){r:3, i:0},(complex){r:6, i:0}}; 
complex x[2] = {(complex){r:3, i:2},(complex){r:2, i:0}}; 
complex y[2]; 
/* execution */ 
cgemv_(&normal, &m, &m, &alpha, &a[0], &m, &x[0], &one, &beta, &y[0], &one); 
cgemv_(&transpose, &m, &m, &alpha, &a[0], &m, &x[0], &one, &beta, &y[0], &one); 

第一cgemv_電話後，y持有16.0000+11.0000i 48.0000+37.0000i，其中MATLAB證實是正確的。

但第二個cgemv_調用後，y保存38.0000+17.0000i 21.0000+6.0000i，而MATLAB說它應該是42.0000-1.0000i 21.0000+6.0000i。我不知道可能會出現什麼問題。

Answer 1

由於2向量乘以2x2矩陣，因此使用pen和paper執行操作並不太複雜。如果轉置矩陣被用於：

（4 + I）*（3 + 2I）+（14 + 3I）* 2 = 38 + 17I

奇怪：

（4-I）* （3 + 2i）+（14-3i）* 2 = 42 -i

因此，MATLAB的輸出可能是使用complex conjugate transpose獲得的輸出。如果參數的TRANS設置爲'C'，則BLAS可以產生相同的輸出。

這裏是一個基於我們的示例代碼，顯示了BLAS真正爲TRANS的不同值計算的代碼。它可以通過gcc main.c -o main -lblas -Wall來計算。

#include <stdlib.h> 
#include <stdio.h> 

#include <complex.h> 


extern int cgemv_(char* trans, int * m, int * n, float complex* alpha, float complex* A, int * lda,float complex * x, int* incx, float complex * beta, float complex * y,int* incy); 

int main(void) { 

    /* config variables */ 
    char normal = 'N'; 
    char transpose = 'T'; 
    char ctranspose = 'C'; 
    int m = 2; 
    float complex alpha = 1.0+0.*I; 
    float complex beta = 0.0+0.*I; 
    int one = 1; 
    /* data buffers */ 
    float complex a[4]= {4+1.*I,14+3.*I,3.+0.*I,6.+0.*I}; 
    float complex x[2] = {3.+2.*I,2+0.*I}; 
    float complex y[2]; 
    /* execution */ 

    float complex ye[2]; 
    ye[0]=a[0]*x[0]+a[2]*x[1]; 
    ye[1]=a[1]*x[0]+a[3]*x[1]; 

    cgemv_(&normal, &m, &m, &alpha, &a[0], &m, &x[0], &one, &beta, &y[0], &one); 

    printf("N\n"); 
    printf("y[0]=%2.6f + %2.6f I expected %6f + %6f I\n",creal(y[0]),cimag(y[0]),creal(ye[0]),cimag(ye[0])); 
    printf("y[1]=%2.6f + %2.6f I expected %6f + %6f I\n",creal(y[1]),cimag(y[1]),creal(ye[1]),cimag(ye[1])); 

    //float complex ye[2]; 
    ye[0]=a[0]*x[0]+a[1]*x[1]; 
    ye[1]=a[2]*x[0]+a[3]*x[1]; 

    cgemv_(&transpose, &m, &m, &alpha, &a[0], &m, &x[0], &one, &beta, &y[0], &one); 

    printf("T\n"); 
    printf("y[0]=%2.6f + %2.6f I expected %6f + %6f I\n",creal(y[0]),cimag(y[0]),creal(ye[0]),cimag(ye[0])); 
    printf("y[1]=%2.6f + %2.6f I expected %6f + %6f I\n",creal(y[1]),cimag(y[1]),creal(ye[1]),cimag(ye[1])); 

    ye[0]=conj(a[0])*x[0]+conj(a[1])*x[1]; 
    ye[1]=conj(a[2])*x[0]+conj(a[3])*x[1]; 

    cgemv_(&ctranspose, &m, &m, &alpha, &a[0], &m, &x[0], &one, &beta, &y[0], &one); 

    printf("C\n"); 
    printf("y[0]=%2.6f + %2.6f I expected %6f + %6f I\n",creal(y[0]),cimag(y[0]),creal(ye[0]),cimag(ye[0])); 
    printf("y[1]=%2.6f + %2.6f I expected %6f + %6f I\n",creal(y[1]),cimag(y[1]),creal(ye[1]),cimag(ye[1])); 


    return 0; 
}