Матрично-векторное произведение с 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;
  }

Компиляция:
g ++ -c smallmatlib.cpp
g ++ smallmatlab.cpp smallmatlib.o -L / usr / local / lib -lclapack -lcblas

Теперь функция matvecprod является проблемой. С примерной матрицей A и примерным вектором v это должно привести к выводу типа

12.4..    26.8..

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

2.00..    0.00..

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

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

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

Answer 1

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

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, Это связано с тем, как работает умножение матрицы на вектор (см., Например, википедию). Вам нужны значения оригинала 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 таким образом, дать неправильный ответ.