#ifndef __MATRIX_H #define __MATRIX_H #include <string> #include <exception> /** * Matrixクラスの例外クラス。 * 例として、演算が成立しない場合など。 * */ class MatrixException: public exception{ private: string what_str; public: MatrixException(const string &what_arg): what_str(what_arg){} ~MatrixException() throw() {} /** * エラー内容を取得します。 * * @return (chsr *) エラー内容 */ const char *what() const throw(){ return what_str.c_str(); } }; #include <string.h> #include <ostream> #include "complex.h" template <class T> class Matrix{ private: T *m_Values; int *ref; //参照カウンタ unsigned int m_Rows; unsigned int m_Columns; public: /** * Matrixクラスのコンストラクタ。 * 指定の行数、指定の列数で行列を生成します。 * また成分はすべてT(0)で初期化されます。 * * @param rows 行数 * @param columns 列数 */ Matrix(const unsigned int &rows, const unsigned int &columns) throw(MatrixException) : m_Rows(rows), m_Columns(columns){ m_Values = new T[rows * columns]; ref = new int(0); if((!m_Values) || (!ref)){throw MatrixException("Lack of memory!!");} (*ref)++; for(int i = 0; i < rows * columns; i++){*(m_Values + i) = T(0);} } /** * コピーコンストラクタ。 * シャローコピーを生成します。 * * @param matrix コピー元 */ Matrix(const Matrix &matrix){ if(m_Values = matrix.m_Values){ (*(ref = matrix.ref))++; m_Rows = matrix.m_Rows; m_Columns = matrix.m_Columns; } } /** * デストラクタ。 */ ~Matrix(){ if(m_Values && ((--(*ref)) <= 0)){delete [] m_Values; delete ref;} } /** * 指定の単位行列を生成します。 * * @parma size 指定の行数(列数) */ static Matrix<T> getI(const int &size){ Matrix<T> result(size, size); for(int i = 0; i < size; i++){result(i, i) = T(1);} return result; } /** * 行数を返します。 * * @return (int) 行数 */ unsigned int rows() const{return m_Rows;} /** * 列数を返します。 * * @return (int) 列数 */ unsigned int columns() const{return m_Columns;} /** * 指定した行の行ベクトルを生成します。 * * @param row 行インデックス * @return (Matrix<T>) 行ベクトル */ Matrix<T> rowVector(const unsigned int &row) const throw(MatrixException){ if(row >= m_Rows){throw MatrixException("Index Incorrect!!");} Matrix<T> vector(1, m_Columns); for(int j = 0; j < m_Columns; j++){vector(0, j) = (*const_cast<Matrix<T> *>(this))(row, j);} return vector; } /** * 指定した列の列ベクトルを生成します。 * * @param column 列インデックス * @return (Matrix<T>) 列ベクトル */ Matrix<T> columnVector(const unsigned int &column) const throw(MatrixException){ if(column >= m_Columns){throw MatrixException("Index Incorrect!!");} Matrix<T> vector(m_Rows, 1); for(int i = 0; i < m_Rows; i++){vector(i, 0) = (*const_cast<Matrix<T> *>(this))(i, column);} return vector; } /** * 指定した行列成分を返します。 * * @param row 行インデックス(開始番号は0〜) * @param column 列インデックス(開始番号は0〜) * @return (T) 成分 */ T &operator()(const int &row, const int &column) throw(MatrixException){ if((row < 0 || row >= m_Rows) || (column < 0 || column >= m_Columns)){ throw MatrixException("Index incorrect"); } return *(m_Values + (row * m_Columns) + column); } /** * 指定した行列成分を返します。 * * @param row 行番号(開始番号は1〜) * @param column 列番号(開始番号は1〜) * @return (T) 成分 */ T &index(const int &row, const int &column) throw(MatrixException){ return (*this)(row -1, column -1); } /** * 行列を複製(ディープコピー)します。 * * @return (Matrix<T>) コピー */ Matrix<T> copy() const throw(MatrixException){ Matrix<T> matrix(m_Rows, m_Columns); /* for(int i = 0; i < m_Rows; i++){ for(int j = 0; j < m_Columns; j++){ matrix(i, j) = (*const_cast<Matrix<T> *>(this))(i, j); } } */ memcpy(matrix.m_Values, m_Values, sizeof(T) * m_Rows * m_Columns); return matrix; } /** * 行を入れ替えます。破壊的メソッドです。 * * @param row1 行インデックス1 * @param row2 行インデックス2 * @return (Matrix<T>) 自分自身 */ Matrix<T> &exchangeRows(const unsigned int &row1, const unsigned int &row2) throw(MatrixException){ if(row1 >= m_Rows || row2 >= m_Rows){throw MatrixException("Index incorrect");} T temp; for(int j = 0; j < m_Columns; j++){ temp = (*this)(row1, j); (*this)(row1, j) = (*this)(row2, j); (*this)(row2, j) = temp; } return *this; } /** * 列を入れ替えます。破壊的メソッドです。 * * @param column1 列インデックス1 * @param column2 列インデックス2 * @return (Matrix<T>) 自分自身 */ Matrix<T> &exchangeColumns(const unsigned int &column1, const unsigned int &column2){ if(column1 >= m_Columns || column2 >= m_Columns){throw MatrixException("Index incorrect");} T temp; for(int i = 0; i < m_Rows; i++){ temp = (*this)(i, column1); (*this)(i, column1) = (*this)(i, column2); (*this)(i, column2) = temp; } return *this; } /** * 代入演算子。 * * @return (Matrix<T>) 自分自身 */ Matrix &operator=(const Matrix &matrix){ if((this == &matrix) || (m_Values == matrix.m_Values)){return *this;} if(m_Values && ((--(*ref)) <= 0)){delete [] m_Values; delete ref;} if(m_Values = matrix.m_Values){ (*(ref = matrix.ref))++; m_Rows = matrix.rows(); m_Columns = matrix.columns(); } return *this; } /** * 正方行列かどうか調べます。 * * @return (bool) 正方行列である場合true、それ以外の場合false */ bool isSquare() const{return m_Rows == m_Columns;} /** * 対称行列かどうか調べます。 * * @return (bool) 対称行列である場合true、それ以外の場合false */ bool isSymmetric() const{ if(isSquare()){ for(int i = 0; i < m_Rows; i++){ for(int j = i + 1; j < m_Columns; j++){ if((*const_cast<Matrix *>(this))(i, j) != (*const_cast<Matrix *>(this))(j, i)){return false;} } } return true; }else{return false;} } /** * 行列の成分全てを指定倍します。破壊的メソッドです。 * * @param scalar 倍数 * @return (Matrix<T>) 自分自身 */ Matrix<T> &operator*=(const T &scalar){ for(int i = 0; i < m_Rows; i++){ for(int j = 0; j < m_Columns; j++){ (*this)(i, j) *= scalar; } } return *this; } /** * 行列の成分全てを指定倍します。 * * @param scalar 倍数 * @return (Matrix<T>) 結果 */ Matrix<T> operator*(const T &scalar) const{return (copy() *= scalar);} /** * 行列の成分全てを指定倍します。 * * @param scalar 倍数 * @return (Matrix<T>) 結果 */ friend Matrix<T> operator*(const T &scalar, const Matrix<T> &matrix){return matrix * scalar;} /** * 行列の成分全てを除算します。破壊的メソッドです。 * * @param scalar 倍数 * @return (Matrix<T>) 自分自身 */ Matrix<T> &operator/=(const T &scalar){return (*this) *= (1 / scalar);} /** * 行列の成分全てを除算します。 * * @param scalar 倍数 * @return (Matrix<T>) 結果 */ Matrix<T> operator/(const T &scalar) const{return (copy() /= scalar);} /** * 行列の成分全てを除算します。 * * @param scalar 倍数 * @return (Matrix<T>) 結果 */ friend Matrix<T> operator/(const T &scalar, const Matrix<T> &matrix){return matrix / scalar;} /** * 行列を転置します。 * * @return (Matrix<T>) 転置行列 */ Matrix<T> transpose() const{ Matrix<T> result(m_Columns, m_Rows); for(int i = 0; i < result.rows(); i++){ for(int j = 0; j < result.columns(); j++){ result(i, j) = (*const_cast<Matrix<T> *>(this))(j, i); } } return result; } /** * 行列を成分ごとに加算します。破壊的メソッドです。 * * @param matrix 加える行列 * @return (Matrix<T>) 自分自身 */ Matrix<T> &operator+=(const Matrix<T> &matrix){ if(m_Rows != matrix.rows() || m_Columns != matrix.columns()){throw MatrixException("Operatiorn void!!");} for(int i = 0; i < m_Rows; i++){ for(int j = 0; j < m_Columns; j++){ (*this)(i, j) += (*const_cast<Matrix<T> *>(&matrix))(i, j); } } return *this; } /** * 行列を成分ごとに加算します。 * * @param matrix 加える行列 * @return (Matrix<T>) 結果 */ Matrix<T> operator+(const Matrix<T> &matrix) const{return (copy() += matrix);} /** * 行列を成分ごとに減算します。 * * @param matrix 引く行列 * @return (Matrix<T>) 自分自身 */ Matrix<T> &operator-=(const Matrix<T> &matrix){ if(m_Rows != matrix.rows() || m_Columns != matrix.columns()){throw MatrixException("Operatiorn void!!");} for(int i = 0; i < m_Rows; i++){ for(int j = 0; j < m_Columns; j++){ (*this)(i, j) -= (*const_cast<Matrix<T> *>(&matrix))(i, j); } } return *this; } /** * 行列を成分ごとに減算します。 * * @param matrix 引く行列 * @return (Matrix<T>) 結果 */ Matrix<T> operator-(const Matrix<T> &matrix) const{return (copy() -= matrix);} /** * 行列を乗算します。 * * @param matrix かける行列 * @return (Matrix<T>) 結果 */ Matrix<T> operator*(const Matrix<T> &matrix) const throw(MatrixException){ if(m_Columns != matrix.rows()){throw MatrixException("Operatiorn void!!");} Matrix<T> result(m_Rows, matrix.columns()); for(int i = 0; i < result.rows(); i++){ for(int j = 0; j < result.columns(); j++){ for(int k = 0; k < m_Columns; k++){ result(i, j) += (*const_cast<Matrix<T> *>(this))(i, k) * (*const_cast<Matrix<T> *>(&matrix))(k, j); } } } return result; } /** * 行列を乗算します。破壊的メソッドです。 * * @param matrix かける行列 * @return (Matrix<T>) 自分自身 */ Matrix<T> &operator*=(const Matrix<T> &matrix) throw(MatrixException){ return (*this = (*this * matrix)); } /** * 単項演算子-。 * 効果は matrix * -1と同じです。 * * @return (Matrix<T>) -matrix */ Matrix<T> operator-() const{return (copy() *= -1);} /** * 補行列(余因子行列)を求めます。 * * @param row 行インデックス * @param column 列インデックス * @return (Matrix<T>) 補行列 */ Matrix<T> coMatrix(const int &row, const int &column) const throw(MatrixException){ if(row < 0 && row >= m_Rows && column < 0 && column >= m_Rows){ throw MatrixException("Index incorrect"); } Matrix<T> result(m_Rows -1, m_Columns -1); for(int i = 0; i < m_Rows -1; i++){ for(int j = 0; j < m_Columns -1; j++){ result(i, j) = (*const_cast<Matrix<T> *>(this))((i < row ? i : i + 1), (j < column ? j : j + 1)); } } return result; } /** * 行列式を計算します。 * * @return (T) 結果 */ T determinant() const throw(MatrixException){ if(!isSquare()){throw MatrixException("Operatiorn void!!");} if(m_Rows == 1){ return (*const_cast<Matrix<T> *>(this))(0, 0); }else{ T sum(0); for(int i = 0; i < m_Rows; i++){ if((*const_cast<Matrix<T> *>(this))(i, 0) != 0){ sum += (*const_cast<Matrix<T> *>(this))(i, 0) * (coMatrix(i, 0).determinant()) * (i % 2 == 0 ? 1 : -1); } } return sum; } } /** * LU分解をします。 * (0, 0)〜(n-1, n-1): L行列 * (0, n)〜(n-1, 2n-1):U行列 * * @return (Matrix<T>) LU分解 */ Matrix<T> divideLU() const throw(MatrixException){ if(!isSquare()){throw MatrixException("Operatiorn void!!");} Matrix<T> LU(m_Rows, m_Columns * 2); for(int i = 0; i < m_Rows; i++){ for(int j = 0; j < m_Columns; j++){ if(i >= j){ LU(i, j) = (*const_cast<Matrix *>(this))(i, j); for(int k = 0; k < j; k++){ LU(i, j) -= (LU(i, k) * LU(k, j + m_Columns)); } }else{ LU(i, j + m_Columns) = (*const_cast<Matrix *>(this))(i, j); for(int k = 0; k < i; k++){ LU(i, j + m_Columns) -= (LU(i, k) * LU(k, j + m_Columns)); } LU(i, j + m_Columns) /= LU(i, i); } } LU(i, i + m_Columns) = 1; } return LU; } /** * 逆行列を求めます。 * * @return (Matrix<T>) 逆行列 */ Matrix<T> inverse() const throw(MatrixException){ if(!isSquare()){throw MatrixException("Operatiorn void!!");} //クラメール(遅い) /* Matrix<T> result(m_Rows, m_Columns); T det; if((det = determinant()) == 0){throw MatrixException("Operatiorn void!!");} for(int i = 0; i < m_Rows; i++){ for(int j = 0; j < m_Columns; j++){ result(i, j) = coMatrix(i, j).determinant() * ((i + j) % 2 == 0 ? 1 : -1); } } return result.transpose() / det; */ //ガウス消去法 Matrix<T> left = copy(); Matrix<T> right = Matrix<T>::getI(m_Rows); if(left(0, 0) == T(0)){ //(0, 0)が存在するように並べ替え for(int i = 1; i <= m_Rows; i++){ if(i == m_Rows){throw MatrixException("Operatiorn void!! ; Invert matrix not exist");} if(left(i, 0) != T(0)){ left.exchangeRows(0, i); right.exchangeRows(0, i); break; } } } for(int i = 0; i < m_Rows; i++){ if(left(i, i) == T(0)){throw MatrixException("Operatiorn void!! ; Invert matrix not exist");} if(left(i, i) != T(1)){ for(int j = 0; j < m_Columns; j++){right(i, j) /= left(i, i);} for(int j = i+1; j < m_Columns; j++){left(i, j) /= left(i, i);} left(i, i) = T(1); } for(int k = 0; k < m_Rows; k++){ if(k == i){continue;} if(left(k, i) != T(0)){ for(int j = 0; j < m_Columns; j++){right(k, j) -= right(i, j) * left(k, i);} for(int j = i+1; j < m_Columns; j++){left(k, j) -= left(i, j) * left(k, i);} left(k, i) = T(0); } } } return right; //LU分解 /* */ } /** * 逆行列をかけます。破壊的メソッドです。 * * @param matrix 行列 * @return (Matrix<T>) 自分自身 */ Matrix<T> &operator/=(const Matrix<T> &matrix) throw(MatrixException){return (*this) *= matrix.inverse();} /** * 逆行列をかけます。 * * @param matrix 行列 * @return (Matrix<T>) 結果 */ Matrix<T> operator/(const Matrix<T> &matrix) const throw(MatrixException){return (copy() /= matrix);} /** * ピボットを指定して、加算します。 * 破壊的です。 * * @param row 行インデックス * @param column 列インデックス * @param matrix 足す行列 */ Matrix<T> pivotMerge(const int &row, const int &column, const Matrix &matrix){ for(int i = 0; i < matrix.rows(); i++){ if(row + i < 0){continue;} else if(row + i >= m_Rows){break;} for(int j = 0; j < matrix.columns(); j++){ if(column + j < 0){continue;} else if(column + j >= m_Columns){break;} (*this)(row + i, column + j) += (*const_cast<Matrix *>(&matrix))(i, j); } } return *this; } /** * ピボットを指定して、加算します。 * * @param row 行インデックス * @param column 列インデックス * @param matrix 足す行列 */ Matrix<T> pivotAdd(const int &row, const int &column, const Matrix &matrix) const{ return copy().pivotMerge(row, column, matrix); } /** * ハウスホルダー変換をしてヘッセンベルク行列を得ます。 * * @param transform 変換に用いた行列の積を格納するポインタ * @return (Matrix<T>) ヘッセンベルク行列 */ Matrix<T> hessenberg(Matrix *transform) const throw(MatrixException){ if(!isSquare()){throw MatrixException("Operation void!!");} Matrix<T> result = copy(); for(int j = 0; j < m_Columns -2; j++){ T t(0); for(int i = j + 1; i < m_Rows; i++){ t += pow(result(i, j), 2); } T s = sqrt(t); if(result(j + 1, j) < 0){s *= -1;} Matrix<T> omega(m_Rows -(j+1), 1); { for(int i = 0; i < omega.rows(); i++){ omega(i, 0) = result(j+i+1, j); } omega(0, 0) += s; } Matrix<T> P = Matrix<T>::getI(m_Rows).pivotMerge(j+1, j+1, -(omega * omega.transpose() / (t + result(j + 1, j) * s))); result = P * result * P; if(transform){(*transform) *= P;} } //ゼロ処理 bool sym = isSymmetric(); for(int j = 0; j < m_Columns -2; j++){ for(int i = j + 2; i < m_Rows; i++){ result(i, j) = T(0); if(sym){result(j, i) = T(0);} } } return result; } /** * ハウスホルダー変換をしてヘッセンベルク行列を得ます。 * * @return (Matrix<T>) ヘッセンベルク行列 */ Matrix<T> hessenberg() const throw(MatrixException){return hessenberg(NULL);} /** * 2次小行列の固有値を求めます。 * * @param row 2次小行列の左上項の行インデックス * @param column 2次小行列の左上項の列インデックス * @param upper 結果(固有値1) * @param lower 結果(固有値2) */ void eigen22(const int &row, const int &column, Complex<T> * &upper, Complex<T> * &lower) const throw(MatrixException){ T a = (*const_cast<Matrix *>(this))(row, column), b = (*const_cast<Matrix *>(this))(row, column + 1), c = (*const_cast<Matrix *>(this))(row + 1, column), d = (*const_cast<Matrix *>(this))(row + 1, column + 1); T root = pow((a -d), 2) + b * c * 4; if(root >= T(0)){ root = sqrt(root); upper = new Complex<T>((a + d + root) / 2); lower = new Complex<T>((a + d -root) / 2); }else{ root = sqrt(root * -1); upper = new Complex<T>((a + d) / 2, root / 2); lower = new Complex<T>((a + d) / 2, root / 2 * -1); } } /** * 固有値、固有ベクトルを求めます。 * 返却値はMatrix<Complex<T> >型で、 * (0,0)〜(n-1,n-1)要素が固有ベクトルの行列 * (0,n)〜(n-1,n)要素が対応する固有値の列ベクトル * になっています。 * (固有ベクトル1,固有ベクトル2,…,固有ベクトルn,固有値)のn×(n+1)行列 * * @param threshold 収束判定に用いる閾値 * @return (Matrix<T>) 固有値、固有ベクトル */ Matrix<Complex<T> > eigen(const T &threshold) const throw(MatrixException){ if(!isSquare()){throw MatrixException("Operatiorn void!!");} //パワー法(べき乗法) /*Matrix<T> result(m_Rows, m_Rows + 1); Matrix<T> source = copy(); for(int i = 0; i < m_Columns; i++){result(0, i) = T(1);} for(int i = 0; i < m_Columns; i++){ while(true){ Matrix<T> approxVec = source * result.columnVector(i); T approxVal(0); for(int j = 0; j < approxVec.rows(); j++){approxVal += pow(approxVec(j, 0), 2);} approxVal = sqrt(approxVal); for(int j = 0; j < approxVec.rows(); j++){result(j, i) = approxVec(j, 0) / approxVal;} T before = result(i, m_Rows); if(abs(before - (result(i, m_Rows) = approxVal)) < threshold){break;} } for(int j = 0; (i < m_Rows - 1) && (j < m_Rows); j++){ for(int k = 0; k < m_Rows; k++){ source(j, k) -= result(i, m_Rows) * result(j, i) * result(k, i); } } } return result;*/ //ダブルQR法 /* <手順> * ハウスホルダー法を適用して、上ヘッセンベルク行列に置換後、 * ダブルQR法を適用。 * 結果、固有値が得られるので、固有ベクトルを計算。 */ //固有値の計算 Matrix<Complex<T> > lambda(m_Rows, 1); //固有値 T mu_sum = T(0), mu_multi = T(0); Complex<T> p1, p2; int m = m_Rows; bool first = true; Matrix transform = getI(m_Rows); Matrix A = hessenberg(&transform); Matrix A_ = A; while(true){ //m = 1 or m = 2 if(m == 1){ lambda(0, 0) = A(0, 0); break; }else if(m == 2){ Complex<T> *lambda1, *lambda2; A.eigen22(0, 0, lambda1, lambda2); lambda(0, 0) = *lambda1; lambda(1, 0) = *lambda2; break; } //μ、μ*の更新(4.143) { Complex<T> *p1_new, *p2_new; A.eigen22(m-2, m-2, p1_new, p2_new); if(first ? (first = false) : true){ if((*p1_new -p1).abs() > p1_new->abs() / 2){ if((*p2_new -p2).abs() > p2_new->abs() / 2){ mu_sum = (p1 + p2).real(); mu_multi = (p1 * p2).real(); }else{ mu_sum = p2_new->real() * 2; mu_multi = pow(p2_new->real(), 2); } }else{ if((*p2_new -p2).abs() > p2_new->abs() / 2){ mu_sum = p1_new->real() * 2; mu_multi = pow(p1_new->real(), 2); }else{ mu_sum = (*p1_new + *p2_new).real(); mu_multi = (*p1_new * *p2_new).real(); } } } p1 = *p1_new, p2 = *p2_new; } //ハウスホルダー変換を繰り返す T b1, b2, b3, r; for(int i = 0; i < m -1; i++){ if(i == 0){ b1 = pow(A(0, 0), 2) -mu_sum * A(0, 0) + mu_multi + A(0, 1) * A(1, 0); b2 = A(1, 0) * (A(0, 0) + A(1, 1) -mu_sum); b3 = A(2, 1) * A(1, 0); }else{ b1 = A(i, i -1); b2 = A(i + 1, i -1); b3 = (i == m -2 ? T(0) : A(i + 2, i -1)); } r = sqrt(pow(b1, 2) + pow(b2, 2) + pow(b3, 2)); Matrix omega(3, 1); { omega(0, 0) = b1 + r * (b1 >= T(0) ? 1 : -1); omega(1, 0) = b2; if(b3 != T(0)){omega(2, 0) = b3;} } Matrix P = getI(m_Rows).pivotMerge(i, i, omega * omega.transpose() * -2 / (omega.transpose() * omega)(0, 0)); A = P * A * P; } //収束判定 T epsilon = threshold + threshold * abs(abs(A(m-2, m-2)) < abs(A(m-1, m-1)) ? A(m-2, m-2) : A(m-1, m-1)); if(abs(A(m-1, m-2)) < epsilon){ lambda(--m, 0) = A(m, m); }else if(abs(A(m-2, m-3)) < epsilon){ Complex<T> *lambda1, *lambda2; A.eigen22(m-2, m-2, lambda1, lambda2); lambda(--m, 0) = *lambda1; lambda(--m, 0) = *lambda2; } } //固有ベクトルの計算 Matrix<Complex<T> > x(m_Rows, m_Rows); //固有ベクトル A = A_; for(int j = 0; j < m_Rows; j++){ int n = m_Rows; for(int i = 0; i < j; i++){if(lambda(j, 0) == lambda(i, 0)){--n;}} x(--n, j) = 1; while(n--> 0){ x(n, j) = x(n+1, j) * (lambda(j, 0) -A(n+1, n+1)); for(int i = n+2; i < m_Rows; i++){ x(n, j) -= x(i, j) * A(n+1, i); } x(n, j) /= A(n+1, n); } } /*Matrix<Complex<T> > lambda2(m_Rows, m_Rows); for(int i = 0; i < m_Rows; i++){ lambda2(i, i) = lambda(i, 0); } cout << "A:" << A << endl; //cout << "x * x^-1" << x * x.inverse() << endl; cout << "x * lambda * x^-1:" << x * lambda2 * x.inverse() << endl;*/ //結果の格納 Matrix<Complex<T> > result(m_Rows, m_Rows + 1); for(int j = 0; j < x.columns(); j++){ result(j, m_Rows) = lambda(j, 0); for(int i = 0; i < x.rows(); i++){ for(int k = 0; k < transform.columns(); k++){ result(i, j) += transform(i, k) * x(k, j); } } //正規化 Complex<T> _norm; for(int i = 0; i < m_Rows; i++){ _norm += result(i, j).abs2(); } T norm = sqrt(_norm.real()); for(int i = 0; i < m_Rows; i++){ result(i, j) /= norm; } } return result; } /** * 固有値、固有ベクトルを求めます。 * 返却値はMatrix<T>型で、 * (0,0)〜(n-1,n-1)要素が固有ベクトルの行列 * (0,n)〜(n-1,n)要素が対応する固有値の列ベクトル * になっています。 * (固有ベクトル1,固有ベクトル2,…,固有ベクトルn,固有値)のn×(n+1)行列 * * @return (Matrix<T>) 固有値、固有ベクトル */ Matrix<Complex<T> > eigen() const throw(MatrixException){return eigen(10E-10);} /** * 行列を見やすい形で出力します。 * */ friend ostream &operator<<(ostream &out, const Matrix<T> &matrix){ if(matrix.m_Values){ out << "{"; for(int i = 0; i < matrix.rows(); i++){ out << (i == 0 ? "" : ",") << endl << "{"; for(int j = 0; j < matrix.columns(); j++){ out << (j == 0 ? "" : ",") << (*const_cast<Matrix<T> *>(&matrix))(i, j); } out << "}"; } out << endl << "}"; } return out; } }; #endif /* __MATRIX_H */