用一个简单的例子理解 C++ 中的 LAPACK 调用
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/10112135/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Understanding LAPACK calls in C++ with a simple example
提问by RDK
I am a beginner with LAPACK and C++/Fortran interfacing. I need to solve linear equations and eigenvalues problems using LAPACK/BLAS on Mac OS-X Lion. OS-X Lion provides optimized BLAS and LAPACK libraries (in /usr/lib) and I am linking these libraries instead of downloading them from netlib.
我是 LAPACK 和 C++/Fortran 接口的初学者。我需要在 Mac OS-X Lion 上使用 LAPACK/BLAS 解决线性方程和特征值问题。OS-X Lion 提供优化的 BLAS 和 LAPACK 库（在 /usr/lib 中），我链接这些库而不是从 netlib 下载它们。
My program (reproduced below) is compiling and running fine, but it is giving me wrong answers. I have researched in the web and Stackoverflow and the issue may have to deal with how C++ and Fortran store arrays in differing formats (row major vs Column major). However, as you will see in my example, the simple array for my example should look identical in C++ and fortran. Anyway here goes.
我的程序（转载如下）正在编译和运行良好，但它给了我错误的答案。我已经在网络和 Stackoverflow 中进行了研究，这个问题可能必须处理 C++ 和 Fortran 如何以不同的格式（行专业与列专业）存储数组。但是，正如您将在我的示例中看到的，我示例的简单数组在 C++ 和 fortran 中看起来应该相同。无论如何，这里是。
Lets say we want to solve the following linear system:
假设我们要解决以下线性系统：
x + y = 2
x + y = 2
x - y = 0
x - y = 0
The solution is (x,y) = (1,1). Now I tried to solve this using Lapack as follows
解是 (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);
}
This code was compiled as follows:
这段代码编译如下：
g++ -Wall -llapack -lblas lapacktest.cpp
g++ -Wall -llapack -lblas lapacktest.cpp
On running this, however, I get the solution as (-2,2) which is obviously wrong. I have tried all combination of row/column re-arrangement of my matrix a. Also observe the matrix representation of ashould be identical in row and column formats. I think getting this simple example to work will get me started with LAPACK and any help will be appreciated.
但是，在运行它时，我得到的解决方案为 (-2,2)，这显然是错误的。我已经尝试了矩阵的所有行/列重新排列组合a。还要注意a行和列格式的矩阵表示应该是相同的。我认为让这个简单的例子工作将使我开始使用 LAPACK，任何帮助将不胜感激。
回答by Stephen Canon
You need to factorthe matrix (by calling dgetrf) before you can solve the system using dgetrs. Alternatively, you can use the dgesvroutine, which does both steps for you.
您需要因子矩阵（通过调用dgetrf），然后才能使用解决系统dgetrs。或者，您可以使用dgesv例程，它为您完成两个步骤。
By the way, you don't need to declare the interfaces yourself, they are in the Accelerate headers:
顺便说一下，您不需要自己声明接口，它们在 Accelerate 标头中：
// 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);
}
回答by The Quantum Physicist
For those who don't want bother with the Accelerate Framework, I provide the code of Stephen Canon (thanks to him, of course) with nothing but pure library linking
对于那些不想打扰 Accelerate Framework 的人，我提供了 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);
}
And about the manual, there's a full PDF version available at Intel's website. Here's a sample of their HTML documentation.
关于手册，英特尔网站上有完整的 PDF 版本。这是他们的 HTML 文档示例。
回答by Michael Lehn
If you want to use LAPACK from C++ you might want to have a look a FLENS. It defines low- and high-level interfaces to LAPACK but also re-implements some LAPACK functions.
如果您想从 C++ 使用 LAPACK，您可能想看看FLENS。它定义了 LAPACK 的低级和高级接口，但也重新实现了一些 LAPACK 功能。
With the low-level FLENS-LAPACK interface you can use your own matrix/vector types (if they have a LAPACK conform memory layout). Your call of dgetrfwould look like that:
通过低级 FLENS-LAPACK 接口，您可以使用自己的矩阵/向量类型（如果它们具有符合 LAPACK 的内存布局）。你的电话dgetrf看起来像这样：
info = lapack::getrf(NoTrans, dim, nrhs, a.begin(), LDA, ipiv);
and for dgetrs
并为 dgetrs
lapack::getrs(NoTrans, dim, nrhs, a.begin(), LDA, ipiv, b.begin(), LDB);
So the low-level FLENS-LAPACK functions are overloaded with respect to the element types. Consequently LAPACK function sgetrs, dgetrs, cgetrs, zgetrsare in the low-level interface of FLENS-LAPACK lapack::getrs. You also pass parameters by value/reference and not as pointer (e.g. LDAinstead of &LDA).
因此，低级 FLENS-LAPACK 函数在元素类型方面被重载。因此LAPACK功能sgetrs，dgetrs，cgetrs，zgetrs是在FLENS-LAPACK的低层接口lapack::getrs。您还可以通过值/引用而不是作为指针（例如LDA代替&LDA）传递参数。
If you use the FLENS matrix-types you can code it as
如果您使用 FLENS 矩阵类型，您可以将其编码为
info = lapack::trf(NoTrans, A, ipiv);
if (info==0) {
lapack::trs(NoTrans, A, ipiv, b);
}
Or you just use the LAPACK driver function dgesv
或者你只是使用 LAPACK 驱动程序功能 dgesv
lapack::sv(NoTrans, A, ipiv, b);
Here a list of FLENS-LAPACKdriver functions.
这是FLENS-LAPACK驱动程序函数的列表。
Disclaimer: Yes, FLENS is my baby! That means I coded about 95% of it and every line of code was worth it.
免责声明：是的，FLENS 是我的宝贝！这意味着我编写了大约 95% 的代码，每一行代码都是值得的。
