通过一个简单的示例了解C ++中的LAPACK调用

问题描述

我是LAPACK和C ++ / Fortran接口的初学者。我需要解决线性方程和特征值问题使用LAPACK / BLAS在Mac OS X狮子。 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存储不同格式的数组（行主要vs列主要）。但是，正如你在我的例子中看到的，我的例子的简单数组应该在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 = 2

x - y = 0

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

这篇关于通过一个简单的示例了解C ++中的LAPACK调用的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！