首页 文章

通过一个简单的例子来理解C语言中的LAPACK调用

提问于
浏览 次
17

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

c++ fortran lapack

3 回答

  • 7

    在使用 dgetrs 求解系统之前，需要对矩阵进行分解（通过调用 dgetrf ） . 或者，您可以使用 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);
}
    回复于
  • 21

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

    回复于
  • 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);

    这里有一个list of FLENS-LAPACK驱动程序 .

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

    回复于

相关问题