fork download
copy 
  1. #include <iostream>
  2. #include <vector>
  3. #include <cstdlib>
  4. #include <iomanip>
  5.  
  6. const int N1 = 3;
  7. const int N2 = 4;
  8.  
  9. template <class T>
  10. class LUDecompose {
  11.   typedef std::vector<std::vector<T> > VV;
  12.   typedef std::vector<T> V;
  13.   VV a, l, u, x, y;
  14.   V b;
  15.   int n;
  16.  
  17.   void MakeSquareMatrix(VV& a) {
  18.     a.resize(n);
  19.     for (int i = 0; i < n; i++)
  20.       a[i].resize(n);
  21.   }
  22.   void MakeVector(V& b) {
  23.     b.resize(n);
  24.   }
  25.   // LU分解(ドゥーリトル法)
  26.   void LUDecomp(VV& a, VV& l, VV& u) {
  27.     int i, j, k;
  28.  
  29.     for (i = 0; i < n; i++)
  30.       l[i][i] = 1.0;
  31.     for (i = 0; i < n; i++)
  32.       u[0][i] = a[0][i];
  33.  
  34.     for (i = 1; i < n; i++) {
  35.       if (u[i - 1][i - 1] == 0.0) {
  36.         std::cerr << "行列Uは正則ではありません" << std::endl;
  37.         std::exit(1);
  38.       }
  39.       for (j = i; j < n; j++) {
  40.         // L成分の算出
  41.         l[j][i - 1] = a[j][i - 1];
  42.         for (k = 0; k < i - 1; k++)
  43.           l[j][i - 1] -= l[j][k] * u[k][i - 1];
  44.         l[j][i - 1] /= u[i - 1][i - 1];
  45.         // U成分の算出
  46.         u[i][j] = a[i][j];
  47.         for (k = 0; k < i; k++)
  48.           u[i][j] -= l[i][k] * u[k][j];
  49.       }
  50.     }
  51.   }
  52.   void PrintMatrix(VV& a, const char* s) {
  53.     std::cout << "*** " << s << "***\n";
  54.  
  55.     for (int i = 0; i < n; i++) {
  56.       for (int j = 0; j < n; j++) {
  57.         std::cout << std::setw(11) << std::setprecision(7) << std::fixed
  58.           << std::showpoint << std::right << a[i][j];
  59.         std::cout << "  ";
  60.       }
  61.       std::cout << '\n';
  62.     }
  63.   }
  64.   // LU分解を利用して逆行列を求める
  65.   void LUInv(VV& l, VV& u, VV& x) {
  66.     int i, j, k;
  67.     double r;
  68.     VV lu;
  69.  
  70.     // 上三角行列と下三角行列を合成する
  71.     MakeSquareMatrix(lu);
  72.  
  73.     for (i = 0; i < n; i++)
  74.       for (j = 0; j < n; j++)
  75.         lu[i][j] = (i <= j) ? u[i][j] : l[i][j];
  76.  
  77.     // 逆行列を求める
  78.     for (k = 0; k < n; k++) {
  79.       for (i = 0; i < n; i++) {
  80.         r = (i == k) ? 1.0 : 0.0;
  81.         for (j = 0; j < n; j++)
  82.           r -= lu[i][j] * x[j][k];
  83.         x[i][k] = r;
  84.       }
  85.       for (i = n - 1; i >= 0; i--) {
  86.         r = x[i][k];
  87.         for (j = i + 1; j < n; j++)
  88.           r -= lu[i][j] * x[j][k];
  89.         x[i][k] = r / lu[i][i];
  90.       }
  91.     }
  92.   }
  93.   // 行列の積
  94.   void MatMul(VV& z, VV& x, VV& y) {
  95.     int i, j, k;
  96.     double s;
  97.  
  98.     for (i = 0; i < n; i++)
  99.       for (j = 0; j < n; j++) {
  100.         s = 0.0;
  101.       for (k = 0; k < n; k++)
  102.         s += x[i][k] * y[k][j];
  103.       z[i][j] = s;
  104.     }
  105.   }
  106.   /* 三角行列の行列式は対角成分の積 */
  107.   double DetTriangularMatrix(VV& a) {
  108.     double det = 1.0;
  109.  
  110.     for (int i = 0; i < n; i++)
  111.       det *= a[i][i];
  112.  
  113.     return det;
  114.   }
  115.  
  116. public:
  117.   LUDecompose(int nn, T* init) : n(nn) {
  118.     MakeSquareMatrix(a);
  119.     for (int i = 0; i < n; i++)
  120.       for (int j = 0; j < n; j++)
  121.         a[i][j] = init[i * n + j];
  122.     MakeSquareMatrix(l);
  123.     MakeSquareMatrix(u);
  124.     MakeSquareMatrix(x);
  125.     MakeSquareMatrix(y);
  126.   }
  127.   void Calculate() {
  128.     LUDecomp(a, l, u);
  129.     LUInv(l, u, x);
  130.     MatMul(y, a, x);
  131.   }
  132.   void PrintMatrix() {
  133.     PrintMatrix(a, "行列 A");
  134.     PrintMatrix(l, "行列 L");
  135.     PrintMatrix(u, "行列 U");
  136.     PrintMatrix(x, "逆行列 X(A^(-1))");
  137.     PrintMatrix(y, "行列 Y(A * X)");
  138.     std::cout << "行列式 |A| = " << DetTriangularMatrix(l) * DetTriangularMatrix(u) << std::endl;
  139.   }
  140. };
  141.  
  142. int main(void)
  143. {
  144.   double a1[] = {2.0, 3.0, -1.0, 4.0, 4.0, -3.0, 2.0, -3.0, 1.0};
  145.   double a2[] = {1.0, 1.0, 0.0, 3.0, 2.0, 1.0, -1.0, 1.0, 3.0, -1.0, -1.0, 2.0, -1.0, 2.0, 3.0, -2.0};
  146.  
  147.   LUDecompose<double> lu1(N1, a1);
  148.   lu1.Calculate();
  149.   lu1.PrintMatrix();
  150.  
  151.   std::cout << '\n';
  152.  
  153.   LUDecompose<double> lu2(N2, a2);
  154.   lu2.Calculate();
  155.   lu2.PrintMatrix();
  156. }
  157.  
Not running #stdin #stdout 0s 0KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/qeepAY
language:
C++ (gcc 8.3)
created:
11 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!