Multiplication matrice-vecteur BLAS vs multiplication de matrice-vecteur. On travaille; l'autre échoue

Question

Je réussis la multiplication matricielle-vectorielle en travaillant avec cgemv une fonction BLAS lvl 2 dans Lapack, mais quand j'essaye la transposition, j'obtiens une mauvaise réponse. Pouvez-vous m'instruire dans mon erreur? (J'utilise en fait l'emballage C, non FORTRAN.)

Je tentais

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

Pour être clair, le premier réussit. Le second donne une sortie incorrecte.

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

Après le premier appel cgemv_, y détient 16.0000+11.0000i 48.0000+37.0000i, qui confirme Matlab être correcte.

Mais après le deuxième appel cgemv_, y détient 38.0000+17.0000i 21.0000+6.0000i, alors que MATLAB dit qu'il devrait être 42.0000-1.0000i 21.0000+6.0000i. Je n'ai aucune idée de ce qui pourrait être faux.

Answer 1

Puisqu'un vecteur 2 est multiplié par une matrice 2x2, effectuer l'opération en utilisant un pen et un morceau de paper n'est pas trop complexe. Si la matrice transposée est utilisée:

(4 + i) * (3 + 2i) + (14 + 3i) * 2 = 38 + 17i

Curieusement:

(4-i) * (3 + 2i) + (14-3i) * 2 = 42-i

La sortie de MATLAB est donc probablement la sortie obtenue en utilisant le complex conjugate transpose. La même sortie peut être produite par BLAS si le paramètre TRANS de cgemv_ est défini sur 'C'.

Voici un exemple de code basé sur le nôtre montrant ce que BLAS calcule réellement pour les différentes valeurs de TRANS. Il peut être calculé par 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; 
}