LAPACK求解系統（矩形全包裝（RFP）格式）

Question

我正在使用RFP格式存儲矩陣並嘗試解決方程組。但結果是錯誤的。 =（可能有人請幫忙。

我得到 B = {5.5，10，8.5}

而應該 B = {2.875，4.75，3.5}

我不明白， 。在那裏我已經犯了一個錯誤只需簡單地使用標準功能：分解，然後解因式分解矩陣

謝謝

#include "mkl.h" 
#include "mkl_lapacke.h" 
#include <math.h> 
#include <iostream> 
using namespace std; 
#define NI 3 
#define NJ 1 

int main(int argc, char* argv[]) 
{ 

    double a[NI][NI]; 
    double b[NI][NJ]; 

    a[0][0] = 2; a[0][1] = -1; a[0][2] = 0; 
    a[1][0] = 0; a[1][1] = 2; a[1][2] = -1; 
    a[2][0] = 0; a[2][1] = 0; a[2][2] = 2; 

    b[0][0] = 1; 
    b[0][1] = 6; 
    b[0][2] = 7; 

    cout << "A1 = \n"; 
    for(int i = 0; i < NI; i++) { 
    for(int j = 0; j < NI; j++) { 
     cout << a[i][j] << "\t"; 
    } 
    cout << "\n"; 
    } 
    cout << "\n"; 



    cout << "B1 = \n"; 
    for(int i = 0; i < NI; i++) { 
    for(int j = 0; j < NJ; j++) { 
     cout << b[i][j] << "\t"; 
    } 
    cout << "\n"; 
    } 
    cout << "\n"; 

    char transr = 'N'; 
    char uplo = 'U'; 
    lapack_int n = NI; 
    lapack_int lda = NI; //LDA is used to define the distance in memory between elements of two consecutive columns which have the same row index. 
    double * arf = new double[ n * (n + 1)/2 ]; 
    lapack_int info = -1; 

CONVER。牛逼一般矩陣RFP格式

info = LAPACKE_dtrttf(LAPACK_ROW_MAJOR, transr, uplo, n, *a, lda, arf); 
    //lapack_int LAPACKE_<?>trttf(int matrix_order, char transr, char uplo, lapack_int n, const <datatype>* a, lapack_int lda, <datatype>* arf); 
    cout << "LAPACKE_dtrttf = " << info << "\n"; 

    cout << "Rectangular full packed: \n"; 
    //cout.setf(std::ios::scientific); 
    for(int i = 0; i < NI; i++) { 
    for(int j = 0; j < (NI+1)/2; j++) { 
     cout << arf[i * (NI+1)/2 + j] << "\t"; 
     //cout << arf[i] << "\t"; 
    } 
    cout << "\n"; 
    } 
    cout << "\n"; 

比化矩陣

int matrix_order = LAPACK_ROW_MAJOR; 
    transr = 'N'; 
    uplo = 'U'; 
    n = NI; 

    info = LAPACKE_dpftrf(matrix_order, transr, uplo, n, arf); 
    //lapack_int LAPACKE_<?>pftrf(int matrix_order, char transr, char uplo, lapack_int n, <datatype>* a); 
    cout << "INFO LAPACKE_dpftrf = " << info << "\n"; 
    cout << "Factorized matrix: " << endl; 

    for(int i = 0; i < NI; i++) { 
    for(int j = 0; j < (NI+1)/2; j++) { 
     cout << arf[i*(NI+1)/2+j] << "\t"; 
    } 
    cout << "\n"; 
    } 
    cout << "\n"; 

解決系統

lapack_int nrhs = NJ; 
    lapack_int ldb = NJ; 
    info = LAPACKE_dpftrs(matrix_order, transr, uplo, n, nrhs, arf, &b[0][0], ldb); 
    //lapack_int LAPACKE_<?>pftrs(int matrix_order, char transr, char uplo, lapack_int n, lapack_int nrhs, const <datatype>* a, <datatype>* b, lapack_int ldb); 
    cout << "INFO LAPACKE_dpftrs = " << info << "\n"; 

結果

cout << "Solved = \n"; 
    for(int i = 0; i < NI; i++) { 
    for(int j = 0; j < NJ; j++) { 
     cout << b[i][j] << "\t"; 
    } 
    cout << "\n"; 
    } 
    cout << "\n"; 
    delete [] arf; 

    char ch; 
    cin.get(ch); 

    return 0; 
} 

Answer 1

你肯定NJ的值是正確的如下：

b[NI][NJ]; 
b[0][0] = 1; 
b[0][1] = 6; 
b[0][2] = 7; 

其中NI = 3，NJ = 1 ...

我想你已經與柱放錯位置的行的值。如果是這個問題，那麼你可能想知道它爲什麼沒有給出任何錯誤，這也是一個已經在這裏的問題： Accessing an array out of bounds gives no error, why?