Матрица-векторный продукт с dgemm/dgemv

Question

Использование Lapack с C++ дает мне небольшую головную боль. Я нашел функции, определенные для fortran немного эксцентричным, поэтому я попытался сделать несколько функций на C++, чтобы облегчить мне чтение того, что происходит.

Во всяком случае, я не получаю продукт матрицы-вектора, работающий по своему усмотрению. Вот небольшой пример программы.

smallmatlib.cpp:

#include <cstdio> 
#include <stdlib.h> 


extern "C"{ 
    // product C= alphaA.B + betaC            
    void dgemm_(char* TRANSA, char* TRANSB, const int* M, 
       const int* N, const int* K, double* alpha, double* A, 
       const int* LDA, double* B, const int* LDB, double* beta, 
       double* C, const int* LDC); 
    // product Y= alphaA.X + betaY            
    void dgemv_(char* TRANS, const int* M, const int* N, 
       double* alpha, double* A, const int* LDA, double* X, 
       const int* INCX, double* beta, double* C, const int* INCY); 
    } 

void initvec(double* v, int N){ 
    for(int i= 0; i<N; ++i){ 
     v[i]= 0.0; 
     } 
    } 

void matvecprod(double* A, double* v, double* u, int N){ 
    double alpha= 1.0, beta= 0.0; 
    char no= 'N', tr= 'T'; 
    int m= N, n= N, lda= N, incx= N, incy= N; 
    double* tmp= new double[N]; 
    initvec(tmp, N); 
    dgemv_(&no,&m,&n,&alpha,A,&lda,v,&incx,&beta,tmp,&incy); 
    for(int i= 0; i<N; ++i){ 
     u[i]= tmp[i]; 
     } 
    delete [] tmp; 
    } 

void vecmatprod(double* v, double* A, double* u, int N){ 
    double alpha= 1.0, beta= 0.0; 
    char no= 'N', tr= 'T'; 
    int m= N, n= 1, k= N, lda= N, ldb= N, ldc= N; 
    double* tmp= new double[N]; 
    initvec(tmp, N); 
    dgemm_(&no,&no,&m,&n,&k,&alpha,A,&lda,v,&ldb,&beta,tmp,&ldc); 
    for(int i= 0; i<N; ++i){ 
     u[i]= tmp[i]; 
     } 
    delete [] tmp; 
    } 

smallmatlib.h:

#ifndef SMALLMATLIB_H 
#define SMALLMATLIB_H 

void initvec(double* v, int N); 

void matvecprod(double* A, double* v, double* u, int N); 

void vecmatprod(double* v, double* A, double* u, int N); 

#endif 

smallmatlab.cpp:

#include "smallmatlib.h" 
#include <cstdio> 
#include <stdlib.h> 
#define SIZE 2 

int main(){ 
    double A[SIZE*SIZE]= 
    { 1,2, 
     3,4 }; 
    double v[SIZE]= {2,5.2}; 
    double u[SIZE]= {0,0}; 
    matvecprod(A,v,u,SIZE); 
    printf("%f %f\n",u[0],u[1]); 
    vecmatprod(v,A,u,SIZE); 
    printf("%f %f\n",u[0],u[1]); 
    return 0; 
    } 

компилирование:
г ++ -c smallmatlib.cpp
г ++ smallmatlab.cpp smallm atlib.o -L/usr/local/lib -lclapack -lcblas

Теперь проблема с функцией matvecprod. С примером матрицы А и пример вектора V, он должен производить вывод, как

12.4.. 26.8.. 

, но вместо этого он печатает

2.00.. 0.00.. 

Я пытался произвести правильный результат как с DGEMM и dgemv, но не смог. У меня есть подозрение, что мои переменные incx и incy не имеют правильных значений, но трудно найти объяснение для них, которое я бы понял.

Меньшая проблема заключается в том, что в данный момент я не могу использовать их как vecmatprod (v, A, v, SIZE) - то есть, я всегда должен определить вектор и, который будет содержать результат, отдельно, и вызовите vecmatprod (v, A, u, SIZE). Любая помощь будет оценена по достоинству.

Как примечание, я также начинаю на C++, поэтому я ценю любую критику/предложение, которое может иметь о моем коде.

Answer 1

Вы правы, проблема в значении incx - это должно быть 1, взгляните на reference.

INCX is INTEGER 
    On entry, INCX specifies the increment for the elements of 
    X. INCX must not be zero. 

Таким образом, это значение должно быть использован, когда значение вектора x не помещается один за другим (вектором комплексных значений, например, когда вы хотите использовать только вещественную часть).

Также вы не можете использовать vecmatprod(v,A,v,SIZE) с v как x, так и y. Это связано с тем, как работает умножение матричных векторов (например, wikipedia). Вам нужны значения оригинала x все время для получения правильного результата.Небольшой пример:

y = A * x 

где

y = [ y1 y2 ] 
A = [ [a11 a12] [a21 a22] ] 
x = [ x1 x2 ] 

И мы рассчитываем y как этот

y1 = a11 * x1 + a12 * x2 
y2 = a21 * x1 + a22 * x2 

Вы можете увидеть, что при вычислении y2 нам нужно x1 и x2, но если вы используете x = A * x (без временного вектора) вы замените x1 на y1 t Муха произносит неправильный ответ.