用一个简单的例子来理解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例程,它为您执行这两个步骤.

顺便说一句,您不需要自己声明接口,它们位于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);
}

Answer 2

对于那些不想使用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);
}

关于手册,英特尔网站上提供了完整的PDF版本.这是他们的HTML文档示例.

http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-A02DB70F-9704-42A4-9071-D409D783D911.htm