行列のLU分解を行いたい

Tips

解説

正方行列のピボット付LU分解を行います。

サンプルプログラム

LU.cpp

#include <iostream>
#include <octave/config.h>
#include <octave/CmplxLU.h>

using namespace std; 

int main()
{ 

  cout.precision(3);
  cout.setf(ios::fixed); 

  int const N = 4;
  Complex const IU(0.0, 1.0);

  ComplexMatrix A(N, N);
  A(0,0) = 3.0 + IU; A(0,1) = 2.0       ; A(0,2) = IU     ; A(0,3) = -IU;
  A(1,0) = 4.0     ; A(1,1) = 1.0 + IU  ; A(1,2) = 6.0    ; A(1,3) = 2.0;
  A(2,0) = 5.0*IU  ; A(2,1) = 3.0 - IU  ; A(2,2) = -4.0*IU; A(2,3) = -3.0;
  A(3,0) = -2.0    ; A(3,1) = 1.0-5.0*IU; A(3,2) = 3.0    ; A(3,3) = -1.0+IU;
  cout << "A = " << endl << A << endl; 

  ComplexLU lu(A); 

  ComplexMatrix L = lu.L();
  ComplexMatrix U = lu.U();
  Matrix P = lu.P(); 

  cout << "L = " << endl << L << endl;
  cout << "U = " << endl << U << endl;
  cout << "P = " << endl << P << endl;
  
  cout << "P A = " << endl << P*A << endl;
  cout << "L U = " << endl << L*U << endl; 

  return 0;  
}

出力

A = 
 (3.000,1.000) (2.000,0.000) (0.000,1.000) (-0.000,-1.000)
 (4.000,0.000) (1.000,1.000) (6.000,0.000) (2.000,0.000)
 (0.000,5.000) (3.000,-1.000) (-0.000,-4.000) (-3.000,0.000)
 (-2.000,0.000) (1.000,-5.000) (3.000,0.000) (-1.000,1.000)

L = 
 (1.000,0.000) (0.000,0.000) (0.000,0.000) (0.000,0.000)
 (0.000,0.400) (1.000,0.000) (0.000,0.000) (0.000,0.000)
 (0.000,-0.800) (-0.515,0.340) (1.000,0.000) (0.000,0.000)
 (0.200,-0.600) (-0.289,0.351) (0.276,0.145) (1.000,0.000)

U = 
 (0.000,5.000) (3.000,-1.000) (-0.000,-4.000) (-3.000,0.000)
 (0.000,0.000) (0.600,-6.200) (1.400,0.000) (-1.000,2.200)
 (0.000,0.000) (0.000,0.000) (9.922,-0.476) (2.233,-0.926)
 (0.000,0.000) (0.000,0.000) (0.000,0.000) (0.333,-1.883) 

P = 
 0.000 0.000 1.000 0.000
 0.000 0.000 0.000 1.000
 0.000 1.000 0.000 0.000
 1.000 0.000 0.000 0.000 

P A = 
 (0.000,5.000) (3.000,-1.000) (0.000,-4.000) (-3.000,0.000)
 (-2.000,0.000) (1.000,-5.000) (3.000,0.000) (-1.000,1.000)
 (4.000,0.000) (1.000,1.000) (6.000,0.000) (2.000,0.000)
 (3.000,1.000) (2.000,0.000) (0.000,1.000) (0.000,-1.000) 

L U = 
 (0.000,5.000) (3.000,-1.000) (0.000,-4.000) (-3.000,0.000)
 (-2.000,0.000) (1.000,-5.000) (3.000,0.000) (-1.000,1.000)
 (4.000,0.000) (1.000,1.000) (6.000,0.000) (2.000,0.000)
 (3.000,1.000) (2.000,-0.000) (0.000,1.000) (0.000,-1.000)

コメント