19

私はLAPACKとC++/Fortranインターフェースの初心者です。MacOS-XLionで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を使い始めることができ、助けていただければ幸いです。

4

3 に答える 3

21

を使用して系を解く前に、( を呼び出して) 行列を因数分解する必要があります。または、両方の手順を実行するルーチンを使用することもできます。dgetrfdgetrsdgesv

ところで、インターフェイスを自分で宣言する必要はありません。これらは 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 に答える
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はzgetrsFLENS-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 に答える