问题描述
在 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 的矩阵向量乘积的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!