使用Lapack with 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,SIZE).
任何帮助,将不胜感激.
作为旁注,我也是C的初学者,所以我感谢您对我的代码提出的任何批评/建议.
最佳答案 你是对的,问题在于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(没有临时向量),你将用y1替换x1,从而产生错误的答案.