19

我是 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，任何帮助将不胜感激。

4

3 回答 3

21

您需要先分解矩阵（通过调用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);
}
于 2012-04-11T18:58:15.893 回答
10

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

于 2013-09-08T21:22:50.403 回答
2

如果你想使用 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驱动程序函数的列表。

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

于 2012-10-02T18:01:36.990 回答