行列のできるプログラミング言語(2)

C++編

■boost::numeric::ublas編

uBLAS_Util.hpp

#ifndef __UBLAS_UTIL_HPP__
#define __UBLAS_UTIL_HPP__

#include <iostream>
#include <cmath>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/triangular.hpp>

using namespace std;
using namespace boost;
using namespace boost::numeric;
using namespace boost::numeric::ublas;

typedef ublas::matrix<double> dmatrix;
typedef ublas::matrix<complex<double> > cdmatrix;
typedef ublas::vector<double> dvector;
typedef ublas::vector<complex<double> > cdvector;

/*
 matrixの表示
 */
template<class Matrix>
void m_print(Matrix x) {
 int n1 = x.size1();
 int n2 = x.size2();
 for (int i = 0; i < n1; i++) {
  for (int j = 0; j < n2; j++) {
   cout << x(i, j) << "\t";
  }
  cout << endl;
 }
 return;
}

/*
 vectorの表示
 */
template<class Vector>
void v_print(Vector x) {
 int n = x.size();
 for (int i = 0; i < n; i++) {
  cout << x(i) << endl;
 }
 return;
}

#endif

main.cpp

#include <iostream>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include "uBlas_Util.hpp"
using namespace std;
using namespace boost::numeric::ublas;

int main() {
 // make matrix A
 dmatrix A(3, 3);
 double m[][3] = { { 3, -4, 2 }, { 2, 5, 3 }, { -2, 3, -2 } };
 for (int i = 0; i < 3; i++) {
  for (int j = 0; j < 3; j++) {
   A(i, j) = m[i][j];
  }
 }
 cout << "matrix A:" << endl;
 m_print(A);
 // make A^-1
 dmatrix I = identity_matrix<double>(3);// 単位行列
 dmatrix copyA(A);// Aのコピー
 // LU分解時の行交換情報が入る行列
 permutation_matrix<> pm(A.size1());
 // LU分解
 lu_factorize(A,pm);
 // 行交換を考慮した前進消去と後退代入
 lu_substitute(A,pm,I);

 // 解の表示
 cout << "A^-1" << endl;
 m_print(I);
 // A*A^-1
 cout << "A*A^-1" << endl;
 m_print(prod(copyA,I));

 return 0;
}