我是LAPACK和C++/Fortran接口的初学者。我需要在Mac OS-X Lion上使用LAPACK/BLAS解决线性方程和特征值问题。OS-X Lion 提供了优化的 BLAS 和 LAPACK 库(在/usr/lib 中),我正在链接这些库,而不是从 netlib 下载它们。
我的程序(下面转载)编译和运行良好,但它给了我错误的答案。我已经在网络和Stackoverflow中进行了研究,这个问题可能必须处理C++和Fortran如何以不同的格式存储数组(行主要与列主要)。但是,正如您将在我的示例中看到的那样,我示例中的简单数组在 C++ 和 fortran 中应该看起来相同。反正来了。
假设我们要求解以下线性系统:
x + y = 2
x - y = 0
解是 (x,y) = (1,1)。现在我尝试使用 Lapack 解决此问题,如下所示
// LAPACK test code
#include<iostream>
#include<vector>
using namespace std;
extern "C" void dgetrs(char *TRANS, int *N, int *NRHS, double *A,
int *LDA, int *IPIV, double *B, int *LDB, int *INFO );
int main()
{
char trans = 'N';
int dim = 2;
int nrhs = 1;
int LDA = dim;
int LDB = dim;
int info;
vector<double> a, b;
a.push_back(1);
a.push_back(1);
a.push_back(1);
a.push_back(-1);
b.push_back(2);
b.push_back(0);
int ipiv[3];
dgetrs(&trans, &dim, &nrhs, & *a.begin(), &LDA, ipiv, & *b.begin(), &LDB, &info);
std::cout << "solution is:";
std::cout << "[" << b[0] << ", " << b[1] << ", " << "]" << std::endl;
std::cout << "Info = " << info << std::endl;
return(0);
}
此代码编译如下:
g++ -Wall -llapack -lblas lapacktest.cpp
但是,在运行此功能时,我得到的解决方案为 (-2,2),这显然是错误的。我已经尝试了矩阵a
的行/列重新排列的所有组合。还要注意a
的矩阵表示在行和列格式上应该相同。我认为让这个简单的例子工作将使我开始使用 LAPACK,任何帮助都将不胜感激。
您需要先分解矩阵(通过调用 dgetrf
),然后才能使用 dgetrs
求解系统。 或者,您可以使用 dgesv
例程,它为您执行这两个步骤。
顺便说一下,你不需要自己声明接口,它们在加速头中:
// LAPACK test code
#include <iostream>
#include <vector>
#include <Accelerate/Accelerate.h>
using namespace std;
int main()
{
char trans = 'N';
int dim = 2;
int nrhs = 1;
int LDA = dim;
int LDB = dim;
int info;
vector<double> a, b;
a.push_back(1);
a.push_back(1);
a.push_back(1);
a.push_back(-1);
b.push_back(2);
b.push_back(0);
int ipiv[3];
dgetrf_(&dim, &dim, &*a.begin(), &LDA, ipiv, &info);
dgetrs_(&trans, &dim, &nrhs, & *a.begin(), &LDA, ipiv, & *b.begin(), &LDB, &info);
std::cout << "solution is:";
std::cout << "[" << b[0] << ", " << b[1] << ", " << "]" << std::endl;
std::cout << "Info = " << info << std::endl;
return(0);
}
对于那些不想打扰加速框架的人,我提供了Stephen Canon的代码(当然,感谢他),除了纯粹的库链接之外什么都没有
// LAPACK test code
//compile with: g++ main.cpp -llapack -lblas -o testprog
#include <iostream>
#include <vector>
using namespace std;
extern "C" void dgetrf_(int* dim1, int* dim2, double* a, int* lda, int* ipiv, int* info);
extern "C" void dgetrs_(char *TRANS, int *N, int *NRHS, double *A, int *LDA, int *IPIV, double *B, int *LDB, int *INFO );
int main()
{
char trans = 'N';
int dim = 2;
int nrhs = 1;
int LDA = dim;
int LDB = dim;
int info;
vector<double> a, b;
a.push_back(1);
a.push_back(1);
a.push_back(1);
a.push_back(-1);
b.push_back(2);
b.push_back(0);
int ipiv[3];
dgetrf_(&dim, &dim, &*a.begin(), &LDA, ipiv, &info);
dgetrs_(&trans, &dim, &nrhs, & *a.begin(), &LDA, ipiv, & *b.begin(), &LDB, &info);
std::cout << "solution is:";
std::cout << "[" << b[0] << ", " << b[1] << ", " << "]" << std::endl;
std::cout << "Info = " << info << std::endl;
return(0);
}
关于手册,英特尔网站上有一个完整的PDF版本。这是他们的 HTML 文档示例。
http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-A02DB70F-9704-42A4-9071-D409D783D911.htm
从C++使用LAPACK,你可能想看看FLENS。 它定义了 LAPACK 的低级和高级接口,但也重新实现了 LAPACK 的一些函数。
使用低级 FLENS-LAPACK 接口,您可以使用自己的矩阵/矢量类型(如果它们具有符合 LAPACK 的内存布局)。 您的dgetrf
调用如下所示:
info = lapack::getrf(NoTrans, dim, nrhs, a.begin(), LDA, ipiv);
和dgetrs
lapack::getrs(NoTrans, dim, nrhs, a.begin(), LDA, ipiv, b.begin(), LDB);
因此,低级 FLENS-LAPACK 函数相对于元素类型是重载的。 因此LAPACK函数sgetrs
、dgetrs
、cgetrs
、zgetrs
都位于FLENS-LAPACK lapack::getrs
的低级接口中。 您还可以按值/引用而不是指针传递参数(例如 LDA
而不是&LDA
)。
如果使用 FLENS 矩阵类型,则可以将其编码为
info = lapack::trf(NoTrans, A, ipiv);
if (info==0) {
lapack::trs(NoTrans, A, ipiv, b);
}
或者你只是使用LAPACK驱动程序功能dgesv
lapack::sv(NoTrans, A, ipiv, b);
这里是FLENS-LAPACK驱动程序功能的列表。
免责声明:是的,弗伦斯是我的宝贝!这意味着我编码了大约 95% 的代码,每一行代码都是值得的。