// This may look like C code, but it is really -*- C++ -*-
// RTSS --- Railway Total System Simulator
// (c) TAKAGI Ryo (at Kogakuin University)
// matrix.cc --- class matrix, class diagmatrix functions
// -----
// ChangeLog:
// 2008. 7. 6
//  Improved speed of matrix::eqnsolve(). Re-designed both matrix and
//  diagmatrix classes so that they are derived from vector < double > .
// 2007. 11. 22
//  Removed fprintf(...).
// -----

#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <vector>
#include <iostream>
#include "matrix.hh"

#define MATRIX_INVLIM 1e-30

using std :: cerr ;
using std :: endl ;


// ----
// LU decomposition. The original matrix will be lost. The resulting data
// can then be used by the eqnsolve() routine.
// LU分解。元の行列は失われる。結果得られるデータはその後 eqnsolve() 関数で
// 使うことができる。
// ----
int
matrix :: decompose_LU
()
{
  if ( m != n )
  {
    cerr << "Error: matrix must be square if you want to do LU decomposition"
         << endl ;
    exit ( 1 ) ;
  }
  if ( LU_decomposed )
  {
    cerr << "Warning: this matrix seems to have been LU-decomposed before"
         << endl ;
  }
  else
  {
    LU_decomposed = true ;
  }
  vector < double > weight ;
  weight . resize ( n , 0.0 ) ;
  ip . resize ( n , 0 ) ;
  for ( size_type k = 0 ; k < n ; ++ k )
  {
    ip [ k ] = k ;
    double u = 0 ;
    for ( size_type j = 0 ; j < n ; ++ j )
    {
      double t = fabs ( element ( k , j ) ) ;
      if ( t > u )
      {
        u = t ;
      }
    }
    if ( u < MATRIX_INVLIM )
    {
      cerr << "Warning: no inverse matrix" << endl ;
      return 0 ;
    }
    weight [ k ] = 1 / u ;
  }
  for ( size_type k = 0 ; k < n ; ++ k )
  {
    double u = - 1 ;
    size_type j = n ;
    for ( size_type i = k ; i < n ; ++ i )
    {
      size_type ii = ip [ i ] ;
      double t = fabs ( element ( ii , k ) ) * weight [ ii ] ;
      if ( t > u )
      {
        u = t ;
        j = i ;
      }
    }
    size_type ik = ip [ j ] ;
    if ( j != k )
    {
      ip [ j ] = ip [ k ] ;
      ip [ k ] = ik ;
    }
    u = element ( ik , k ) ;
    if ( fabs ( u ) < MATRIX_INVLIM )
    {
      cerr << "Warning: no inverse matrix" << endl ;
      return 0;
    }
    for ( size_type i = k + 1 ; i < n ; ++ i )
    {
      size_type ii = ip [ i ] ;
      ( * this ) ( ii , k ) /= u ;
      double t = element ( ii , k ) ;
      for ( size_type j = k + 1 ; j < n ; ++ j )
      {
        ( * this ) ( ii , j ) -= t * element ( ik , j ) ;
      }
    }
  }
  return 1 ;
}



// ----
// Solve x of A x = b using the data obtained from the LU decomposition above.
// A is "this" matrix, and b, x will be given as arguments. Returns 1 if
// successful.
// 上のLU分解で得られたデータを用い，A x = b の x を求める。A がこの行列，b
// および x は引数として与える。成功すると1を返す。
// ----
int
matrix :: eqnsolve
( const matrix & b ,
  matrix & x )
{
  if ( ! LU_decomposed )
  {
    cerr << "Error: this matrix has never been LU_decomposed. Do that first "
         << "before doing eqnsolve" << endl ;
    exit ( 1 ) ;
  }

  x . renewmn ( b . row () ) ;

  for ( size_type i = 0 ; i < n ; ++ i )
  {
    size_type ii = ip [ i ] ;
    double t = b . element ( ii ) ;
    for ( size_type j = 0 ; j < i ; ++ j )
    {
      t -= element ( ii , j ) * x ( j ) ;
    }
    x ( i ) = t ;
  }
  for ( size_type ix = n ; ix > 0 ; -- ix )
  {
    size_type i = ix ;
    -- i ;
    double t = x ( i ) ;
    size_type ii = ip [ i ] ;
    for ( size_type j = i + 1 ; j < n ; ++ j )
    {
      t -= element ( ii , j ) * x ( j ) ;
    }
    x ( i ) = t / element ( ii , i ) ;
  }
  return 1 ;
}


// ----
// Old style eqnsolve.
// 古い eqnsolve() 関数。遅いので使わない方がよい。
// ----
// int matrix::eqnsolve() {
//     if ( m > n ) {
//       cerr << "Error: matrix::eqnsolve(): more rows than columns" << endl ;
//       exit ( 22 ) ;		// error 22: 行列演算エラー（単純なミス）
//     }
//     size_type i , j , k ;
//     size_type pivr , indexk ;
//     double w , piv ;
//     static size_type * index = 0 ;
//     static size_type indexsize = 0 ;
//     if ( m > indexsize )
//     {
//       delete index ;
//       index = new size_type [ indexsize = m ] ;
//     }

//     // ＬＵ分解法 --- 前進部分
//     for (k = 0; k < m - 1; k++) {
//       // ピボットを探す
// 	piv = 0.0;
// 	pivr = 0;
// 	for (i = k; i < m; i++) {
// 	    w = fabs(operator()(i, k));
// 	    if (piv < w) {
// 		piv = w;
// 		pivr = i;
// 	    }
// 	}
// 	index[k] = pivr;
// 	if (fabs(piv) < MATRIX_INVLIM) goto END;
//       // 行の入換: pivot を対角要素に持ってくる
// 	if (pivr != k) {
// 	    for (j = 0; j < n; j++)
// 		swapf(operator()(pivr, j), operator()(k, j));
// 	}
// 	w = 1.0 / operator()(k, k);
// 	for (j = k; j < n; j++) {
// 	    operator()(k, j) *= w;
// 	    for (i = k; i < m; i++)
// 		operator()(i, j) -= operator()(i, k) * operator()(k, j);
// 	}
//     }

//     // 前進部分 --- 最後の行の分
//     piv = fabs( operator() ( m - 1 , m - 1 ) ) ;
//     index [ k ] = m - 1 ;
//     if ( fabs ( piv ) < MATRIX_INVLIM ) goto END ;
//     w = 1.0 / operator() ( m - 1 , m - 1 ) ;
//     for ( j = m ; j < n ; ++ j )
//       operator() ( m - 1 , j ) *= w ;

//     // 後退部分
//     for ( size_type ix = m - 1 ; ix > 0 ; -- ix )
//     {
//       i = ix ;
//       -- i ;
//       for ( j = m ; j < n ; ++ j )
//       {
//         for ( k = i + 1 ; k < m ; ++ k )
//           operator() ( i , j ) -= operator() ( i , k ) * operator() ( k , j ) ;
//       }
//     }

//     // 列の交換
//     for (size_type kx = m ; kx > 0 ; -- kx )
//     {
//       k = kx ;
//       -- k ;
//       indexk = index [ k ] ;
//       if ( indexk != k )
//       {
//         for ( i = 0 ; i < m ; ++ i )
//           swapf ( operator() ( i , indexk ) , operator() ( i , k ) ) ;
//       }
//     }
//     return 1 ;

//  END :
// //    cerr << "Error: no inverse matrix" << endl ;
//     cerr << "Warning: no inverse matrix" << endl ;
//     return 0 ;
// //    exit(21);		// error 21: 逆行列がなかった
// }



// ----
// Resetting all elements to a value (default = 0).
// 行列のすべての要素をある値にリセットする。(デフォルト = 0)
// ----
void
matrix :: renew
( double x /* = 0.0 */ )
{
  assign ( size () , x ) ;
}



// // ----
// // Transform "this" object to an inverse matrix. Returns 1 if successful.
// // ##### TODO: rather than transforming this to an inverse matrix, I prefer
// // ##### creating an inverse matrix inside the object.
// // このオブジェクトを逆行列に変化させる。成功すれば1を返す。
// // ★★★今後行うこと: このオブジェクトを変化させる代わりに，逆行列を別に生成
// // ★★★してメンバとして保持する方がよさそう。
// // ----
// int
// matrix :: inverse
// ()
// {
//   if ( m != n )
//   {
//     cerr << "Error: numbers of rows and columns differ" << endl ;
//     exit ( 22 ) ;	// error 22: data error in matrix calculations
//   }

//   static matrix z ;
//   z . renewmn ( m , m * 2 ) ;
//   for ( size_type i = 0 ; i < m ; ++ i )
//   {
//     for ( size_type j = 0 ; j < m ; ++ j )
//       z ( i , j ) = operator() ( i , j ) ;
//     z ( i , i + m ) = 1.0 ;
//   }

//   if ( z . eqnsolve () == 0 ) return 0 ;

//   for ( size_type i = 0 ; i < m ; ++ i )
//   {
//     for ( size_type j = 0 ; j < m ; ++ j )
//       operator() ( i , j ) = z ( i , j + m ) ;
//   }
//   return 1 ;
// }



// ----
// Getting the inverse of a diagonal matrix. Returns 1 if successful.
// 対角行列の逆行列の計算（要素の逆数）。成功すれば1を返す。
// ----
int
diagmatrix :: inverse
()
{
  double det = 1.0 ;
  for ( size_type i = 0 ; i < m ; ++ i )
  {
    if ( det < MATRIX_INVLIM ) goto END ;
    ( * this ) [ i ] = 1.0 / ( * this ) [ i ] ;
  }
  return 1 ;

 END :
  cerr << "Warning: no inverse matrix" << endl ;
  return 0 ;	// error 21: 逆行列がなかった
}



// ----
// Subprogram of constructors. Arguments: m & n (rows and columns)
// コンストラクタ, 引数: m, n
// ----
void
matrix :: matrixinst
( matrix :: size_type mm ,
  matrix :: size_type nn )
{
  if ( mm < 0 || nn < 0 )
  {
    cerr << "Error: matrix::matrixinst(...): "
         << "illegal number of rows and columns "
         << "required for the matrix constructor" << endl ;
    exit ( 22 ) ;		// error 22: 行列演算エラー（単純なミス）
  }
  m = mm ;
  n = nn ;
  resize ( m * n , 0.0 ) ;
}



// ----
// Constructor. Arguments: m (rows and columns)
// コンストラクタ, 引数: m
// ----
diagmatrix :: diagmatrix
( size_type mm )
  : x ( 0.0 )
{
  if ( mm < 0 )
  {
    cerr << "Error: diagmatrix::diagmatrix(int): "
         << "illegal number of rows and columns "
         << "required for the matrix constructor" << endl ;
    exit ( 22 ) ;		// error 22: 行列演算エラー（単純なミス）
  }
  m = mm ;
  resize ( m , 0.0 ) ;
}



// ----
// Copy constructor.
// 参照による初期化
// ----
matrix :: matrix
( const matrix & x )
  : ip () ,
    LU_decomposed ( false )
{
  m = x . row () ;
  n = x . col () ;
  resize ( m * n , 0.0 ) ;
  for ( size_type i = 0 ; i < x . size () ; ++ i )
  {
    ( * this ) [ i ] = x [ i ] ;
  }
}



// ----
// Copy constructor.
// 参照による初期化
// ----
diagmatrix :: diagmatrix
( const diagmatrix & x )
  : x ( 0.0 )
{
  m = x . row () ;
  resize ( m , 0.0 ) ;
  for ( size_type i = 0 ; i < x . size () ; ++ i )
  {
    ( * this ) [ i ] = x [ i ] ;
  }
}



// ----
// Substitution operator.
// 代入演算子
// ----
matrix &
matrix :: operator=
( const matrix & x )
{
  if ( & x == this ) return * this ;
  m = x . row () ;
  n = x . col () ;
  assign ( x . begin () , x . end () ) ;
  LU_decomposed = false ;
  return * this ;
}



// ----
// Substitution operator.
// 代入演算子
// ----
diagmatrix &
diagmatrix :: operator=
( const diagmatrix & x )
{
  if ( & x == this ) return * this ;
  m = x . row () ;
  resize ( m , 0.0 ) ;
  for ( size_type i = 0 ; i < x . size() ; ++ i )
  {
    ( * this ) [ i ] = x [ i ] ;
  }
  return * this ;
}




// ----
// Substitution operator.
// 代入演算子
// ----
matrix &
matrix :: operator+=
( const matrix & x )
{
  if ( ( m != x . row () ) || ( n != x . col () ) )
  {
    cerr <<  "Error: can't add these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  if ( m == 0 || n == 0 ) return * this ;
  for ( size_type i = 0 ; i < m ; ++ i )
  {
    for ( size_type j = 0 ; j < n ; ++ j )
      operator() ( i , j ) += x . element ( i , j ) ;
  }
  return * this ;
}




// ----
// Substitution operator.
// 代入演算子
// ----
matrix &
matrix :: operator-=
( const matrix & x )
{
  if ( ( m != x . row () ) || ( n != x . col () ) )
  {
    cerr << "Error: can't add these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  if ( m == 0 || n == 0 ) return * this ;
  for ( size_type i = 0 ; i < m ; ++ i )
  {
    for ( size_type j = 0 ; j < n ; ++ j )
      operator() ( i , j ) -= x . element ( i , j ) ;
  }
  return * this ;
}




// ----
// Substitution operator.
// 代入演算子
// ----
matrix &
matrix :: operator+=
( const diagmatrix & x )
{
  if ( ( m != x . row () ) || ( n != x . col () ) )
  {
    cerr <<  "Error: can't add these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  if ( m == 0 || n == 0 ) return * this ;
  for ( size_type i = 0 ; i < m ; ++ i )
  {
    for ( size_type j = 0 ; j < n ; ++ j )
      operator() ( i , j ) += x . element ( i , j ) ;
  }
  return * this ;
}




// ----
// Substitution operator.
// 代入演算子
// ----
matrix &
matrix :: operator-=
( const diagmatrix & x )
{
  if ( ( m != x . row () ) || ( n != x . col () ) )
  {
    cerr << "Error: can't add these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  if ( m == 0 || n == 0 ) return * this ;
  for ( size_type i = 0 ; i < m ; ++ i )
    operator() ( i , i ) -= x . element ( i ) ;
  return * this ;
}




// ----
// Substitution operator.
// 代入演算子
// ----
diagmatrix &
diagmatrix :: operator+=
( const diagmatrix & x )
{
  if ( m != x . row () )
  {
    cerr << "Error: can't add these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  if ( m == 0 ) return * this ;
  for ( size_type i = 0 ; i < m ; ++ i )
    operator() ( i ) += x . element ( i ) ;
  return * this ;
}




// ----
// Substitution operator.
// 代入演算子
// ----
diagmatrix &
diagmatrix :: operator-=
( const diagmatrix & x )
{
  if ( m != x . row () )
  {
    cerr <<  "Error: can't add these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  if ( m == 0 ) return * this ;
  for ( size_type i = 0 ; i < m ; ++ i )
    operator() ( i ) -= x . element ( i ) ;
  return * this ;
}




// ----
// Substitution operator.
// 代入演算子
// ----
matrix &
matrix :: operator*=
( const matrix & x )
{
  if ( ( n != x . row () ) )
  {
    cerr <<  "Error: can't multiply these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  matrix z = * this ;
  renewmn ( z . m , x . n ) ;
  for ( size_type i = 0 ; i < z . m ; ++ i )
  {
    for ( size_type k = 0 ; k < x . n ; ++ k )
    {
      for ( size_type j = 0 ; j < z . n ; ++ j )
        operator() ( i , k ) += z ( i , j ) * x . element ( j , k ) ;
    }
  }
  return * this ;
}



// ----
// Substitution operator.
// 代入演算子
// ----
matrix &
matrix :: operator*=
( const diagmatrix & x )
{
  if ( ( n != x . row () ) )
  {
    cerr << "Error: can't multiply these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  matrix z = * this ;
  renewmn ( z . m , x . col () ) ;
  for ( size_type i = 0 ; i < z . m ; ++ i )
  {
    for ( size_type k = 0 ; k < x . col () ; ++ k )
    {
      for ( size_type j = 0 ; j < z . n ; ++ j )
        operator() ( i , k ) += z ( i , j ) * x . element ( j , k ) ;
    }
  }
  return * this ;
}




// ----
// Multiplying by a constant.
// 定数倍
// ----
matrix &
matrix :: operator*=
( double x )
{
  for ( size_type i = 0 ; i < m ; ++ i )
    for ( size_type j = 0 ; j < n ; ++ j )
      operator() ( i , j ) *= x ;
  return * this ;
}




// ----
// Multiplying by a constant.
// 定数倍
// ----
diagmatrix &
diagmatrix :: operator*=
( double x )
{
  for ( size_type i = 0 ; i < m ; ++ i )
    operator() ( i ) *= x ;
  return * this ;
}




// ----
// Substitution operator.
// 代入演算子
// ----
diagmatrix &
diagmatrix :: operator*=
( const diagmatrix & x )
{
  if ( m != x . m )
  {
    cerr <<  "Error: can't multiply these matrices" << endl ;
    exit ( 22 ) ;		// error 22: 行列演算エラー（単純なミス）
  }
  for ( size_type i = 0 ; i < m ; ++ i )
    operator() ( i ) *= x . element ( i ) ;
  return * this ;
}




// ----
// Multiplicaiton operator.
// 行列の積
// ----
matrix operator*
( const matrix & x ,
  const matrix & y )
{
  if ( x . n != y . m )
  {
    cerr << "Error: can't multiply these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  matrix z = x ;
  z *= y ;
  return z ;
}




// ----
// Multiplicaiton operator.
// 行列の積
// ----
matrix
operator*
( const diagmatrix & x ,
  const matrix & y )
{
  if ( x . m != y . m )
  {
    cerr <<  "Error: can't multiply these matrices" << endl ;
    exit ( 22 ) ;		// error 22: 行列演算エラー（単純なミス）
  }
  matrix z ( x . m , y . n ) ;
  for ( diagmatrix :: size_type i = 0 ; i < x . m ; ++ i )
  {
    for ( matrix :: size_type j = 0 ; j < y . n ; ++ j )
      z ( i , j ) = x . element ( i ) * y . element ( i , j ) ;
  }
  return z ;
}




// ----
// Multiplicaiton operator.
// 行列の積
// ----
matrix operator*
( const matrix & x ,
  const diagmatrix & y )
{
  if ( x . n != y . m )
  {
    cerr << "Error: can't multiply these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  matrix z = x ;
  z *= y ;
  return z ;
}




// ----
// Multiplicaiton operator.
// 行列の積
// ----
diagmatrix
operator*
( const diagmatrix & x ,
  const diagmatrix & y )
{
  if ( x . m != y . m )
  {
    cerr << "Error: can't multiply these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  diagmatrix z = x ;
  z *= y ;
  return z ;
}



// ----
// Addition operator.
// 行列の和
// ----
matrix
operator+
( const matrix & x ,
  const matrix & y )
{
  if ( x . m != y . m || x . n != y . n )
  {
    cerr << "Error: can't add these matrices" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  matrix z = x ;
  z += y ;
  return z ;
}



// ----
// Subtraction operator.
// 行列の差
// ----
matrix
operator-
( const matrix & x ,
  const matrix & y )
{
  if ( x . m != y . m || x . n != y . n )
  {
    cerr << "Error: can't subtract the former matrix from the latter"
         << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  matrix z = x ;
  z -= y ;
  return z ;
}




// ----
// Minus sign.
// マイナス符号
// ----
matrix
matrix :: operator-
()
{
  matrix z ;
  z . renewmn ( m , n ) ;
  for ( size_type i = 0 ; i < m ; ++ i )
  {
    for ( size_type j = 0 ; j < n ; ++ j )
      z ( i , j ) = - operator() ( i , j ) ;
  }
  return z ;
}



// ----
// Minus sign.
// マイナス符号
// ----
diagmatrix
diagmatrix :: operator-
()
{
  diagmatrix z ;
  z . renewmn ( m ) ;
  for ( size_type i = 0 ; i < m ; ++ i )
    z ( i ) = - operator() ( i ) ;
  return z ;
}




// ----
// Multiplying by a constant.
// 定数倍
// ----
matrix
operator*
( double x ,
  const matrix & y )
{
  matrix z = y ;
  z *= x ;
  return z ;
}




// ----
// Multiplying by a constant.
// 定数倍
// ----
diagmatrix
operator*
( double x ,
  const diagmatrix & y )
{
  diagmatrix z = y ;
  z *= x ;
  return z ;
}



// ----
// Renewing the number of rows and columns. Arguments are new m and n.
// 行,列数の更新  引数: m, n
// ----
void
matrix :: renewmn
( size_type mm ,
  size_type nn /* = 1 */ )
{
  if ( mm < 0 || nn < 0 )
  {
    cerr << "Error: matrix::renewmn(...): illegal number of rows and "
         << "columns required for the matrix renewer" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  m = mm ;
  n = nn ;
  resize ( m * n , 0.0 ) ;
  renew () ;
  LU_decomposed = false ;
}



// ----
// Renewing the number of rows and columns. Argument is new m.
// 行,列数の更新  引数: m
// ----
void
diagmatrix :: renewmn
( size_type mm )
{
  if ( mm < 0 )
  {
    cerr << "Error: diagmatrix::renewmn(...): illegal number of rows and "
         << "columns required for the matrix renewer" << endl ;
    exit ( 22 ) ;	// error 22: 行列演算エラー（単純なミス）
  }
  m = mm ;
  resize ( m , 0.0 ) ;
  assign ( size () , 0.0 ) ;
}




// ----
// Transposing the matrix.
// 行列の転置
// ----
void
matrix :: transpose
( matrix & z )
{
  z . renewmn ( n , m ) ;
  for ( size_type i = 0 ; i < m ; ++ i )
  {
    for ( size_type j = 0 ; j < n ; ++ j )
    {
      z ( j , i ) = operator() ( i , j ) ;
    }
  }
}



// ----
// Swapping the element ... This need not be a member function, but
// there it is... (by RT '08.7.6)
// 要素の交換…これはメンバ関数でなくてもいいようにも思うけど… (RT '08.7.6)
// ----
void
matrix :: swapf
( double & a ,
  double & b )
{
  double temp ;
  temp = a ;
  a = b ;
  b = temp ;
}



// ----
// Dividing a matrix into four.
// 行列の4分割
// ----
void
matrix :: quaddiv
( size_type mm ,
  size_type nn ,
  matrix & aa ,
  matrix & ab ,
  matrix & ba ,
  matrix & bb )
{
  if ( mm < 0 || nn < 0 || mm > m || nn > n )
  {
    cerr << "Error: subscript out of range in matrix division" << endl ;
    exit ( 22 ) ;  // error 22: data error in matrix calculations
  }
  aa . renewmn ( mm , nn ) ;
  ab . renewmn ( mm , n - nn ) ;
  ba . renewmn ( m - mm , nn ) ;
  bb . renewmn ( m - mm , n - nn ) ;

  for ( size_type i = 0 ; i < mm ; ++ i )
  {
    for ( size_type j = 0 ; j < nn ; ++ j )
    {
      aa ( i , j ) = operator() ( i , j ) ;
    }
    for ( size_type j = nn ; j < n ; ++ j )
    {
      ab ( i , j - nn ) = operator() ( i , j ) ;
    }
  }
  for ( size_type i = mm ; i < m ; ++ i )
  {
    for ( size_type j = 0 ; j < nn ; ++ j )
    {
      ba ( i - mm , j ) = operator() ( i , j ) ;
    }
    for ( size_type j = nn ; j < n ; ++ j )
    {
      bb ( i - mm , j - nn ) = operator() ( i , j );
    }
  }
}



// ----
// Join four matrices and create one.
// 4分割した行列を再結合
// ----
void
matrix :: quadjoin
( matrix & aa ,
  matrix & ab ,
  matrix & ba ,
  matrix & bb )
{
  if ( aa . m != ab . m || ba . m != bb . m
       || aa . n != ba . n || ab . n != bb . n )
  {
    cerr << "Error: cannot join these four matrices" << endl ;
    exit ( 22 ) ;  // error 22: data error in matrix calculations
  }

  renewmn ( aa . m + ba . m , aa . n + ab . n ) ;

  for ( size_type i = 0 ; i < aa . m ; ++ i )
  {
    for ( size_type j = 0 ; j < aa . n ; ++ j )
    {
      operator() ( i , j ) = aa ( i , j ) ;
    }
    for ( size_type j = aa . n ; j < n ; ++ j )
    {
      operator() ( i , j ) = ab ( i , j - aa . n ) ;
    }
  }
  for ( size_type i = aa . m ; i < m ; ++ i )
  {
    for ( size_type j = 0 ; j < aa . n ; ++ j )
    {
      operator() ( i , j ) = ba ( i - aa . m , j );
    }
    for ( size_type j = aa . n ; j < n ; ++ j )
    {
      operator() ( i , j ) = bb ( i - aa . m , j - aa . n ) ;
    }
  }
}



// ----
// Inverse matrix.
// 逆行列.
// ----
int
matrix :: inverse
( matrix & ret )
{
  if ( ! LU_decomposed )
  {
    cerr << "Error: matrix::inverse(matrix&) must be run after LU decomposition"
         << endl ;
    exit ( 1 ) ;
  }
  ret . renewmn ( row () , col () ) ;
  for ( matrix :: size_type j = 0 ; j < col () ; ++ j )
  {
    matrix tmpcol ;
    tmpcol . renewmn ( row () ) ;
    matrix rtncol ;
    for ( matrix :: size_type i = 0 ; i < row () ; ++ i )
    {
      tmpcol ( i ) = 0.0 ;
      tmpcol ( j ) = 1.0 ;
    }
    eqnsolve ( tmpcol , rtncol ) ;
    for ( matrix :: size_type i = 0 ; i < row () ; ++ i )
    {
      ret ( i , j ) = rtncol ( i ) ;
    }
  }
  return 1 ;
}
