0

行列の格納に 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;

一般的なマトリックスを 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;
}
4

1 に答える 1

1

NJ の値が次の場合に正しいと確信していますか?

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

ここで、NI = 3 および NJ = 1...

行と列の値を間違えたと思います。それが問題だった場合は、エラーが発生しなかった理由を知りたいと思うかもしれません。これは、すでにここにある質問です 。範囲外の配列にアクセスしてもエラーが発生しないのはなぜですか?

于 2013-03-20T16:01:25.993 に答える