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

2022-01-14 00:00:00 fortran lapack c++

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

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.

我的程序(转载如下)正在编译和运行良好，但它给了我错误的答案.我在网络和 Stackoverflow 中进行了研究，问题可能必须处理 C++ 和 Fortran 如何以不同格式存储数组(行专业与列专业).但是，正如您将在我的示例中看到的那样，我示例的简单数组在 C++ 和 fortran 中应该看起来相同.反正就这样吧.

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.

假设我们要解决以下线性系统:

Lets say we want to solve the following linear system:

x + y = 2

x - y = 0

解是 (x,y) = (1,1).现在我尝试使用 Lapack 解决这个问题，如下所示

The solution is (x,y) = (1,1). Now I tried to solve this using Lapack as follows

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

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 a should 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.

推荐答案

你需要分解矩阵(通过调用dgetrf)才能使用<代码>dgetrs.或者，您可以使用 dgesv 例程，它会为您完成这两个步骤.

You need to factor the matrix (by calling dgetrf) before you can solve the system using dgetrs. Alternatively, you can use the dgesv routine, which does both steps for you.

顺便说一句，你不需要自己声明接口，它们在 Accelerate 头文件中:

By the way, you don't need to declare the interfaces yourself, they are in the Accelerate headers:

// 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);
}

相关文章