Hiểu các cuộc gọi LAPACK trong C++ với một ví dụ đơn giản

Question

Tôi là người mới bắt đầu với LAPACK và C++/Fortran giao tiếp. Tôi cần phải giải các phương trình tuyến tính và các vấn đề riêng biệt bằng cách sử dụng LAPACK/BLAS trên Mac OS-X Lion. OS-X Lion cung cấp các thư viện BLAS và LAPACK được tối ưu hóa (trong/usr/lib) và tôi liên kết các thư viện này thay vì tải chúng xuống từ netlib.

Chương trình của tôi (được sao chép bên dưới) đang biên soạn và chạy tốt, nhưng nó cho tôi câu trả lời sai. Tôi đã nghiên cứu trong web và Stackoverflow và vấn đề có thể phải đối phó với cách mảng C++ và Fortran trong các định dạng khác nhau (hàng chính so với cột chính). Tuy nhiên, như bạn sẽ thấy trong ví dụ của tôi, mảng đơn giản cho ví dụ của tôi sẽ trông giống hệt trong C++ và fortran. Dù sao đi nữa.

phép nói rằng chúng ta muốn giải quyết các hệ thống tuyến tính sau:

x + y = 2

x - y = 0

Giải pháp là (x, y) = (1,1). Bây giờ tôi đã cố gắng để giải quyết việc này bằng LAPACK như sau

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

Mã này đã được biên soạn như sau:

g++ -Wall -llapack -lblas lapacktest.cpp

: Khi chạy phần này, tuy nhiên, tôi nhận được các giải pháp như (-2,2) rõ ràng là sai. Tôi đã thử tất cả kết hợp sắp xếp lại hàng/cột của ma trận a. Đồng thời quan sát biểu diễn ma trận của a phải giống hệt nhau ở định dạng hàng và cột. Tôi nghĩ rằng việc làm ví dụ đơn giản này sẽ giúp tôi bắt đầu với LAPACK và mọi trợ giúp sẽ được đánh giá cao.

Answer 1

Bạn cần yếu tố ma trận (bằng cách gọi dgetrf) trước khi bạn có thể giải quyết hệ thống bằng cách sử dụng dgetrs. Ngoài ra, bạn có thể sử dụng thói quen dgesv, thực hiện cả hai bước cho bạn.

Bằng cách này, bạn không cần phải khai báo giao diện chính mình, họ đang có trong các tiêu đề Tăng tốc:

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

Nếu bạn muốn sử dụng LAPACK từ C++ bạn có thể muốn có một cái nhìn một FLENS . Nó định nghĩa các giao diện cấp thấp và cao cho LAPACK nhưng cũng thực hiện lại một số hàm LAPACK.

Với giao diện FLENS-LAPACK cấp thấp, bạn có thể sử dụng các loại ma trận/vector của riêng mình (nếu chúng có bố cục bộ nhớ phù hợp LAPACK). Cuộc gọi của bạn của dgetrf sẽ trông như thế:

info = lapack::getrf(NoTrans, dim, nrhs, a.begin(), LDA, ipiv); 

và cho dgetrs

lapack::getrs(NoTrans, dim, nrhs, a.begin(), LDA, ipiv, b.begin(), LDB); 

Vì vậy, các cấp thấp FLENS-LAPACK chức năng đang quá tải đối với các loại phần tử với. Do đó chức năng LAPACK sgetrs, dgetrs, cgetrs, zgetrs nằm trong giao diện cấp thấp của FLENS-LAPACK lapack::getrs. Bạn cũng chuyển các tham số theo giá trị/tham chiếu và không phải là con trỏ (ví dụ: LDA thay vì &LDA).

Nếu bạn sử dụng FLENS matrix-loại bạn có thể mã hóa nó như

info = lapack::trf(NoTrans, A, ipiv); 
if (info==0) { 
    lapack::trs(NoTrans, A, ipiv, b); 
} 

Hoặc bạn chỉ cần sử dụng chức năng điều khiển LAPACK dgesv

lapack::sv(NoTrans, A, ipiv, b); 

Dưới đây là một list of FLENS-LAPACK chức năng điều khiển.

Tuyên bố từ chối trách nhiệm: Có, FLENS là con của tôi! Điều đó có nghĩa là tôi đã mã hóa khoảng 95% của nó và mọi dòng mã đều đáng giá.

Answer 3

Đối với những người không muốn bận tâm với Tăng tốc Framework, tôi cung cấp mã của Stephen Canon (nhờ anh ấy, tất nhiên) không có gì nhưng thư viện liên kết tinh khiết

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

Và về hướng dẫn, có một phiên bản PDF đầy đủ có sẵn tại trang web của Intel. Đây là một mẫu tài liệu HTML của họ.

http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-A02DB70F-9704-42A4-9071-D409D783D911.htm