// This may look like C code, but it is really -*- C++ -*-
// RTSS --- Railway Total System Simulator
// (c) TAKAGI Ryo (at Kogakuin University)
// matrix.hh --- class matrix, class diagmatrix
// -----
// ChangeLog:
// 2008. 7. 7
//  Included many features for speed-ups.
// 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(...).
// -----

#ifndef ___MATRIX_HH
#define ___MATRIX_HH	// こうしておくと2度インクルードしても大丈夫！

#include <cstdio>
#include <cstdlib>
#include <vector>
#include <iostream>

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


class diagmatrix ;


// ----
// class matrix: matrix calculations
// class matrix: 行列とその演算クラス
// ----
class matrix
  : public vector < double >
{
  // ----
  // Basic informatio for the matrix this object holds
  // このオブジェクトが保持する行列に関する基本情報
  // ----
  size_type m ;			// rows / 行
  size_type n ;			// columns / 列

  // ----
  // Row order information for LU decomposition
  // LU分解における行入れ替え情報
  // ----
  vector < size_type > ip ;

  bool LU_decomposed ;		// true if LU decomposition is done

  // ----
  // Private member functions
  // プライベートメンバ関数
  // ----
  void swapf ( double & , double & ) ;	// Swapping variables / 変数値入換
  void matrixinst ( size_type , size_type ) ; // Initialising / 初期化
    
public :

  // ----
  // Constructors.
  // コンストラクタ。
  // ----
  matrix ( size_type mm , size_type nn )
    : ip () , LU_decomposed ( false ) { matrixinst ( mm , nn ) ; }

  explicit matrix ( size_type mm )
    : ip () , LU_decomposed ( false ) { matrixinst ( mm , 1 ) ; }

  matrix ()
    : vector < double > () , ip () , LU_decomposed ( false )
  { m = n = 0 ; }

  ~matrix () {}		// destructor

  matrix ( const matrix & ) ;	// 参照による初期化
    
  double & operator() ( size_type i , size_type j = 0 ) {	// 要素の参照
    if ( i < 0 || j < 0 || i >= m || j >= n ) {
      cerr << "Error: subscript out of range in matrix element call" << endl ;
      exit ( 1 ) ;
    }
    return ( * this ) [ i * n + j ] ;
  }
    
  double element ( size_type i , size_type j = 0 ) const {	// 要素の参照
    if ( i < 0 || j < 0 || i >= m || j >= n ) {
      cerr << "Error: subscript out of range in matrix element call" << endl ;
      exit ( 2 ) ;
    }
    return ( * this ) [ i * n + j ] ;
  }
    
  // ----
  // Minus sign.
  // マイナス符号
  // ----
  matrix operator- () ;
    
  // ----
  // Substitution operators.
  // 代入等を行う演算子群。
  // ----
  matrix & operator= ( const matrix & ) ;	// Substitution / 代入演算子
  matrix & operator+= ( const matrix & ) ;
  matrix & operator+= ( const diagmatrix & ) ;	// Addition / 加える
  matrix & operator-= ( const matrix & ) ;
  matrix & operator-= ( const diagmatrix & ) ;	// Subtraction / 差し引く
  matrix & operator*= ( const matrix & ) ;
  matrix & operator*= ( const diagmatrix & ) ;	// Multiplication / 掛ける
  matrix & operator*= ( double ) ;
    
  // ----
  // Old LU decomposition implementation (1992 by RT...)
  // LU分解による一次方程式求解の古いインプリメンテーション (1992 by RT...)。
  // ----
//   int eqnsolve () ;
//   int inverse () ;	// Inverse matrix / 逆行列

  // ----
  // New LU decomposition implementation. Reference: Numerical Recipe.
  // 新しいLU分解による一次方程式求解ルーチン。
  // ----
  int decompose_LU () ;
  int eqnsolve ( const matrix & , matrix & ) ;
  int inverse ( matrix & ) ;
  void resetLUDecompFlag () { LU_decomposed = false ; }

  // ----
  // Transposition. Result will be stored in the object given as argument.
  // 転置。結果は引数として与えたオブジェクトに格納される。
  // ----
  void transpose ( matrix & ) ;
    
  void renewmn ( size_type , size_type = 1 ) ;	// 行／列数の変更
  void renew ( double = 0.0 ) ;
  size_type row () const { return m ; }	// 行数
  size_type col () const { return n ; }	// 列数
    
  // Dividing the matrix into four / 行列を４つに分割
  void quaddiv
  ( size_type , size_type , matrix & , matrix & , matrix & , matrix & ) ;

  // Concatenating four matrices to one / 分割された4行列を1つに結合
  void quadjoin ( matrix & , matrix & , matrix & , matrix & ) ;

  friend matrix operator* ( const diagmatrix & , const matrix & ) ;
  friend matrix operator* ( const matrix & , const diagmatrix & ) ;
  friend matrix operator* ( const matrix & , const matrix & ) ;
  friend matrix operator+ ( const matrix & , const matrix & ) ;
  friend matrix operator- ( const matrix & , const matrix & ) ;
  friend matrix operator* ( double , const matrix & ) ;

} ;




// ----
// class diagmatrix: diagonal matrix calculations
// class diagmatrix: 対角行列とその演算クラス
// ----
class diagmatrix
  : public vector < double >
{
  size_type m ;	// 行 = 列
  double x ;    // = 0.0

public :

  diagmatrix ( size_type ) ;	// constructor: 引数: m
  diagmatrix () : m ( 0 ) , x ( 0.0 ) {}	// empty constructor
  ~diagmatrix () {}	// destructor

  diagmatrix ( const diagmatrix & ) ;	// 参照による初期化
    
  double & operator() ( size_type i , size_type j ) {	// 要素の参照
    x = 0.0 ;
    if ( i < 0 || j < 0 || i >= m || j >= m ) {
      cerr << "Error: subscript out of range in diagmatrix element call" 
           << endl ; exit ( 1 ) ;
    } else if ( i != j ) return x ;
    else return ( * this ) [ i ] ;
  }

  double & operator() ( size_type i ) {		// 要素の参照
    if ( i < 0 || i >= m ) {
      cerr << "Error: subscript out of range in diagmatrix element call"
           << endl ; exit ( 2 ) ;
    }
    return ( * this ) [ i ] ;
  }

  double element ( size_type i , size_type j ) const {	// 要素の参照
    if ( i < 0 || j < 0 || i >= m || j >= m ) {
      cerr << "Error: subscript out of range in diagmatrix element call"
           << endl ; exit ( 3 ) ;
    } else if ( i != j ) return 0.0 ;
    else return ( * this ) [ i ] ;
  }

  double element ( size_type i ) const {	// 要素の参照
    if ( i < 0 || i >= m ) {
      cerr << "Error: subscript out of range in diagmatrix element call"
           << endl ; exit ( 4 ) ;
    } else return ( * this ) [ i ] ;
  }

  diagmatrix & operator= ( const diagmatrix & ) ;	// 代入演算子

  // ----
  // Minus sign.
  // マイナス符号
  // ----
  diagmatrix operator- () ;

  // ----
  // Substitution operators.
  // 代入等を行う演算子群。
  // ----
  diagmatrix & operator+= ( const diagmatrix & ) ;
  diagmatrix & operator-= ( const diagmatrix & ) ;
  diagmatrix & operator*= ( const diagmatrix & ) ;
  diagmatrix & operator*= ( double ) ;

  // ----
  // Calculation of inverse matrix.
  // 逆行列
  // ----
  int inverse () ;

  void renewmn ( size_type ) ;		// 行／列数の変更

  size_type row () const { return m ; }	// 行数
  size_type col () const { return m ; }	// 列数

  friend diagmatrix operator* ( const diagmatrix & , const diagmatrix & ) ;
  friend matrix operator* ( const diagmatrix & , const matrix & ) ;
  friend matrix operator* ( const matrix & , const diagmatrix & ) ;

  friend diagmatrix operator* ( double , const diagmatrix & ) ;

} ;




// ----
// Operators
// 演算子群
// ----
matrix operator+ ( const matrix & , const matrix & ) ;
matrix operator- ( const matrix & , const matrix & ) ;
matrix operator* ( const matrix & , const matrix & ) ;
matrix operator* ( const diagmatrix & , const matrix & ) ;
matrix operator* ( const matrix & , const diagmatrix & ) ;
diagmatrix operator* ( const diagmatrix & , const diagmatrix & ) ;
matrix operator* ( double , const matrix & ) ;
diagmatrix operator* ( double , const diagmatrix & ) ;


#endif /* ! ___MATRIX_HH */
