带有 dgemm/dgemv 的矩阵向量乘积 [英] Matrix-vector product with dgemm/dgemv

问题描述

在 C++ 中使用 Lapack 让我有点头疼.我发现为 fortran 定义的函数有点古怪，所以我尝试在 C++ 上创建一些函数，以便我更容易阅读正在发生的事情.

Using Lapack with C++ is giving me a small headache. I found the functions defined for fortran a bit eccentric, so I tried to make a few functions on C++ to make it easier for me to read what's going on.

无论如何，我没有按照我的意愿得到矩阵向量乘积.这是该程序的一个小示例.

Anyway, I'm not getting th matrix-vector product working as I wish. Here is a small sample of the program.

smallmatlib.cpp:

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:

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:

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

Compiling:
g++ -c smallmatlib.cpp
g++ smallmatlab.cpp smallmatlib.o -L/usr/local/lib -lclapack -lcblas

现在函数 matvecprod 是问题所在.使用示例矩阵 A 和示例向量 v，它应该产生类似

Now the function matvecprod is the problem. With the example matrix A and example vector v, it should produce an output like

12.4..    26.8..

而是打印出来

2.00..    0.00..

我尝试使用 dgemm 和 dgemv 产生正确的结果，但未能成功.我有一种预感，我的变量 incx 和 incy 没有正确的值，但很难找到我能理解的解释.

I tried to produce the correct result with both dgemm and dgemv, but wasn't able to. I have a hunch that my variables incx and incy do not have correct values, but it's difficult to find an explanation for them that I'd understand.

一个较小的问题是目前我不能像vecmatprod(v,A,v,SIZE)- 也就是说，我总是必须定义向量 u，它将单独保存结果，并调用 vecmatprod(v,A,u,SIZE).任何帮助将不胜感激.

A smaller problem is that at the moment I can't use them like vecmatprod(v,A,v,SIZE) - that is, I always have to define the vector u, that will hold the result, separately, and call vecmatprod(v,A,u,SIZE). Any help would be appreciated.

顺便说一句，我也是 C++ 的初学者，所以我感谢您对我的代码提出的任何批评/建议.

As a side note, I'm also a beginner at C++, so I appreciate any criticism/suggestion you might have about my code.

推荐答案

你说得对，问题出在incx值——应该是1，看看参考.

You are right, the problem is in incx value - it should be 1, take a look at reference.

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

所以当向量 x 的值不是一一放置时(例如，当您只想使用实部时的复数值向量)时，应使用此值.

So this value should be used when values of vector x is not placed one by one (vector of complex values for example, when you want to use only real part).

你也不能使用 vecmatprod(v,A,v,SIZE) 和 v 作为 x 和 y.这是因为矩阵向量乘法的工作原理(例如，参见 wikipedia).您始终需要原始 x 的值才能产生正确的结果.小例子:

Also you can not use vecmatprod(v,A,v,SIZE) with v as both x and y. This is because how matrix-vector multiplication works (see wikipedia for example). You need values of original x the whole time to produce correct result. Small example:

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(没有临时向量)你将用 y1 替换 x1 从而产生错误的答案.

You can see, that when we calculate y2 we need x1 and x2, but if you use x = A * x(without temporary vector) you will replace x1 with y1 thus produce wrong answer.

这篇关于带有 dgemm/dgemv 的矩阵向量乘积的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！