用一个简单的例子来理解C++中的LAPACK调用

Question

我是初学者，有LAPACK和C++/Fortran接口。我需要在Mac OS-X Lion上使用LAPACK/BLAS解决线性方程和特征值问题。 OS-X Lion提供优化的BLAS和LAPACK库（在/ usr/lib中），我将这些库链接起来，而不是从netlib下载它们。

我的程序（转载如下）编译和运行良好，但它给了我错误的答案。我已经在Web和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，任何帮助将不胜感激。

Answer 1

您需要因子矩阵（通过调用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); 
} 

Answer 2

如果你想使用LAPACK从C++，你可能希望有一个看起来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); 

这里一个list of FLENS-LAPACK驱动程序的功能。

声明：是的，FLENS是我的宝贝！这意味着我编码了大约95％，每行代码都值得。

Answer 3

对于那些不想与加速框架打扰谁，我提供斯蒂芬佳能的代码（感谢他，当然），什么也没有，但纯库链接

// 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