用一個簡單的例子來理解C++中的LAPACK調用

Question

我是初學者，有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，任何幫助將不勝感激。

Answer 1

您需要因子矩陣（通過調用dgetrf）纔可以使用dgetrs來解決系統問題。或者，您可以使用dgesv例程，該例程爲您執行兩個步驟。

順便說一句，你不需要接口自行申報的，他們都在加速標題：

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

Answer 2

如果你想使用LAPACK從C++，你可能希望有一個看起來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％，每行代碼都值得。

Answer 3

對於那些不想與加速框架打擾誰，我提供斯蒂芬佳能的代碼（感謝他，當然），什麼也沒有，但純庫鏈接

// 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