Умножение матрицы на вектор BLAS против умножения вектора на матрицу. Один работает; другой терпит неудачу

Мне удалось умножить матрицу на вектор при работе с cgemv функцией BLAS lvl 2 в Lapack, но когда я пытаюсь транспонировать, я получаю неправильный ответ. Можете ли вы указать мне на мою ошибку? (На самом деле я использую оболочку 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. Я понятия не имею, что может быть не так.


c linear-algebra lapack blas
person JellicleCat    schedule 13.06.2017    source источник
comment
Последний элемент 'a' установлен в (complex){r:6, i:0}, но отличается от первого уравнения в вопросе..?   -  person roygvib    schedule 14.06.2017

Ответы (1)


arrow_upward
1
arrow_downward

Поскольку вектор 2 умножается на матрицу 2x2, выполнение операции с использованием пера и куска бумага не слишком сложна. Если используется транспонированная матрица:

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

Любопытно:

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

Таким образом, вывод MATLAB, скорее всего, является выводом, полученным с использованием комплексно-сопряженное транспонирование. Тот же вывод может быть получен BLAS, если параметр TRANS cgemv_ установлено значение 'C'.

Вот пример кода, основанный на нашем, показывающий, что на самом деле вычисляет 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;
}
person francis    schedule 20.06.2017