cholesky.cpp

// ////////////////////////////////////////////////////////////////
// 対称行列のコレスキー分解を行う
// ////////////////////////////////////////////////////////////////
#include <iostream>

#include <tnt_array2d.h>
#include <tnt_array2d_utils.h>
#include <jama_cholesky.h>

template<class T>
TNT::Array2D<T>
transpose(const TNT::Array2D<T>& X)
{
  int m = X.dim1();
  int n = X.dim2();
  TNT::Array2D<T> Y(n, m);
  int i, j;
  for (i = 0; i < m; ++i) {
    for (j = 0; j < n; ++j) {
      Y[j][i] = X[i][j];
    }
  }
  return Y;
}

int
main(int, char *[])
{
  TNT::Array2D<double> L(5, 5);
  L[0][0] = 4.0; L[0][1] = 0.0; L[0][2] = 0.0; L[0][3] = 0.0; L[0][4] = 0.0;
  L[1][0] = 3.0; L[1][1] = 2.0; L[1][2] = 0.0; L[1][3] = 0.0; L[1][4] = 0.0;
  L[2][0] = 2.0; L[2][1] = 6.0; L[2][2] = 6.0; L[2][3] = 0.0; L[2][4] = 0.0;
  L[3][0] = 2.0; L[3][1] = 8.0; L[3][2] = 6.0; L[3][3] = 5.0; L[3][4] = 0.0;
  L[4][0] = 5.0; L[4][1] = 8.0; L[4][2] = 3.0; L[4][3] = 2.0; L[4][4] = 9.0;
  TNT::Array2D<double> Ld = transpose(L);
  TNT::Array2D<double> A = TNT::matmult(L, Ld);
  std::cout << "A: " << std::endl;
  std::cout << A << std::endl;

  JAMA::Cholesky<double> chol(A);
  std::cout << "getL(): " << std::endl;
  std::cout << chol.getL() << std::endl;

  return 0;
}