matrix6.h

This is an matrix class which has been adapted for complex number computations. It ll be of help for

/*
  Name       :   Matrix Template Library(MTL) 6.0
  Copyright  :   All rights reserved
                 You can copy and use it freely but if you want to changed it
                 you should inform the author  please mail to [e-mail]
  Author     :   Dbflying
  Bug report :   [e-mail]      
  Date       :   03-09-03 14:28
  Description:   Many Thanks to my friend Warren Wo(HangZhou ZheJiang province)
*/
#ifndef MATRIXH
#define MATRIXH
//------------------------------------------------------------------------
#include <complex>
#include <vector>
#include <string>
#include <cmath>
//------------------------------------------------------------------------
using namespace std;
//**namesapce soubam**//
namespace soubam{
//------------------------------------------------------------------------
const int SafeSize = 130;
//------------------------------------------------------------------------
// forward declearation
template<class T>
class CMatrix;
template<class T>
class RMatrix;
// out of range error class
class OutOfRange
{
public:
    OutOfRange(int row=-1,int col=-1,string s="Error")
    :RowIndex(row),ColIndex(col),ErrString(s){}
    int InvalidRow(){return RowIndex;}
    int InvalidCol(){return ColIndex;}
    string GetError(){return ErrString;}
private:
    int RowIndex;
    int ColIndex;
    string ErrString;
};
//------------------------------------------------------------------------
// mismatch error class
class MisMatch
{
public:
    MisMatch(string s="Error"):ErrString(s){}
    string GetError(){return ErrString;}
private:
    string ErrString;   
};
//------------------------------------------------------------------------
template<class T>
class Matrix{
public:
// constructors   
    Matrix();
    Matrix(int row,int col,T value = T(0));
    Matrix(const Matrix<T>& M);
//operator override   
    Matrix<T>& operator = (const Matrix<T>& M);
    Matrix<T>  operator + (const Matrix<T>& M);
    Matrix<T>  operator - (const Matrix<T>& M);
    Matrix<T>  operator * (const Matrix<T>& M);
    Matrix<T>  operator ^ (const int&   times);
    Matrix<T>& operator +=(const Matrix<T>& M);
    Matrix<T>& operator -=(const Matrix<T>& M);
    Matrix<T>& operator *=(const Matrix<T>& M);
    Matrix<T>& operator ^=(const int&   times);
    vector<T>& operator [](const int&     row);
    bool       operator ==(const Matrix<T>& M);
    bool       operator !=(const Matrix<T>& M);

// element access functions
    int GetRowSize()const{return RowSize;}
    int GetColSize()const{return ColSize;}
    T   GetValue(const int& row,const int& col)const
        {return matrix[row][col];}
    const vector<T> GetRowVector(const int& row)const;
    const vector<T> GetColVector(const int& col)const;
    const Matrix<T> GetRow(const int& row)const;
    const Matrix<T> GetCol(const int& col)const;
   
// modify    
    void Resize(const int& row,const int& col,T value=T(0));
    void SetValue(const int& row,const int& col,T value = T(0));    
    void AddRow(T value = T(0));
    void AddCol(T value = T(0));
    void DeleteRow(const int& row);
    void DeleteCol(const int& col);
    void InsertRow(const int& row,T value = T(0));
    void InsertCol(const int& col,T value = T(0));
    void SetRowVector(const int& row,const vector<T>& vec = vector<T>(0));
    void SetColVector(const int& col,const vector<T>& vec = vector<T>(0));
    void SwapRow(const int& row1,const int& row2);
    void SwapCol(const int& col1,const int& col2);
// special functions
    Matrix<T> Turn();
    Matrix<T> Inverse();
    Matrix<T> Unit(const int& row);
    Matrix<T> Diag();
    T         Det();
    T         Trace();
    T         Norm(const int& p=2);
    T         CondNumber();
    T         Radius();// spectral  radius
    int       Rank();
    Matrix<T> GetAij(const int& row,const int& col);
// transform
    CMatrix< complex<double> > ToCDouble(); 
protected:
    typedef vector<T>           VectorType;
    typedef vector<VectorType>  MatrixType;
    VectorType    RowVector;
    MatrixType    matrix; 
    int        RowSize;
    int        ColSize;  
 
};           
//------------------------------------------------------------------------
template<class T>
inline Matrix<T>::Matrix():RowSize(1),ColSize(1)
{
    Resize(1,1);
}
//------------------------------------------------------------------------
template<class T>
inline Matrix<T>::Matrix(int row,int col,T value):RowSize(row),ColSize(col)
{
    Resize(row,col,value);
}
//------------------------------------------------------------------------
template<class T>
inline Matrix<T>::Matrix(const Matrix<T>& M)
{
    *this = M;
}
//------------------------------------------------------------------------
template<class T> Matrix<T>&
Matrix<T>::operator = (const Matrix<T>& M)
{
    if(this!=&M)
    {
        Resize(M.GetRowSize(),M.GetColSize());       
        for(int i = 0;i< M.GetRowSize();++i)
            matrix[i] = M.GetRowVector(i);
    }              
    return *this;
}
//------------------------------------------------------------------------   
template<class T> Matrix<T>
Matrix<T>::operator + (const Matrix<T>& M)
{
    if(RowSize != M.GetRowSize() || ColSize != M.GetColSize())
        throw MisMatch("Error at + : The tow matrixs are not match for '+' .");
    Matrix<T> temp(RowSize,ColSize);
    for(int i=0;i<RowSize;++i)
        for(int j=0;j<ColSize;++j)
            temp.SetValue(i,j,(GetValue(i,j)+M.GetValue(i,j)));
    return temp;           
}
//------------------------------------------------------------------------
template<class T> Matrix<T>
Matrix<T>::operator - (const Matrix<T>& M)
{
    if(RowSize != M.GetRowSize() || ColSize != M.GetColSize())
        throw MisMatch("Error at - : The tow matrixs are not match for '-' .");
    Matrix<T> temp(RowSize,ColSize);
    for(int i=0;i<RowSize;++i)
        for(int j=0;j<ColSize;++j)
            temp.SetValue(i,j,(GetValue(i,j)-M.GetValue(i,j)));
    return temp;           
}
//------------------------------------------------------------------------
template<class T> Matrix<T>
Matrix<T>::operator * ( const Matrix<T>& M)
{
    if(ColSize != M.GetRowSize())
        throw MisMatch("Error ar * : The two matrix are not match for '*'.");
    Matrix<T> temp(RowSize,M.GetColSize());
    vector<T> vec1,vec2;
    T t;
    for(int i=0;i<RowSize;++i)
    {
         for(int j=0;j<M.GetColSize();++j)
         {
              vec1 = GetRowVector(i);
              vec2 = M.GetColVector(j);
              t=T(0);
              for(size_t k=0;k<vec1.size();++k)
                  t += vec1[k]*vec2[k];
              temp.SetValue(i,j,t);
         }
    }
    return temp;
}
//------------------------------------------------------------------------
template<class T>  Matrix<T>
Matrix<T>::operator ^(const int& times)
{
    if(times < 0)
        throw MisMatch("Error at ^ : In a^t where t<0.");
    if(RowSize != ColSize)
       throw MisMatch("Error at ^ : The matrix is not square");
    Matrix<T> M;
    M=*this;
    if(times == 0)
    {
        for(int i=0;i<RowSize;++i)
           for(int j=0;j<ColSize;++j)
               M.SetValue(i,j,T(0));
        for(int i=0;i<RowSize;++i)
           M.SetValue(i,i,T(1));
    }
    else
    {
        for(int i=0;i<times;++i)
            M = M*M ;
    }
    return M;
}
//------------------------------------------------------------------------
template<class T> Matrix<T>&
Matrix<T>::operator += ( const Matrix<T>& M)
{
    if(RowSize != M.GetRowSize()|| ColSize != M.GetColSize())
        throw  MisMatch("Error at +=: The two matrixs are not match for '+'.");
    *this = *this + M;
     return *this;
}
//------------------------------------------------------------------------
template<class T> Matrix<T>&
Matrix<T>::operator -= ( const Matrix<T>& M)
{
    if(RowSize != M.GetRowSize()|| ColSize != M.GetColSize())
        throw MisMatch("Error at -= : The two matrix are not match for '-' .");
    *this = *this - M;
    return *this;
}
//------------------------------------------------------------------------
template<class T>  Matrix<T>&
Matrix<T>::operator *= ( const Matrix<T>& M)
{
    if(ColSize != M.GetRowSize())
        throw MisMatch("Error at *= : The two matrix are not match for '*'.");
    *this = *this * M;
    return *this;
}
//------------------------------------------------------------------------
template<class T>  Matrix<T>&
Matrix<T>::operator ^=(const int& times)
{
    if( times < 0)
        throw MisMatch("Error at ^= : In a^t where t<0");       
    if(RowSize != ColSize )
        throw MisMatch("Error at ^= : The matrix is not square");
    *this=*this^times;
    return  *this;
}
//------------------------------------------------------------------------
template<class T> vector<T>&
Matrix<T>::operator [](const int& row)
{
  if(row >= RowSize || row < 0)
      throw OutOfRange(row,-1,"at []");
  return  matrix[row];   
}
//------------------------------------------------------------------------
template<class T>
bool Matrix<T>::operator ==(const Matrix<T>& M)
{
    if( RowSize != M.GetRowSize() || ColSize !=M.GetColSize())
        return false;
    for(int i=0;i<RowSize;++i)  
       for(int j=0;j<ColSize;++j)
         if(GetValue(i,j) != M.GetValue(i,j))
            return false;
    return true;  
}
//------------------------------------------------------------------------
template<class T>
bool Matrix<T>::operator !=(const Matrix<T>& M)
{
    return !(*this==M);
}
//------------------------------------------------------------------------
template<class T> const vector<T>
Matrix<T>::GetRowVector(const int& row)const
{
    if(row < 0 || row >= RowSize)
        throw OutOfRange(row,-1,"at GetRowVector");
    return matrix[row];   
}
//------------------------------------------------------------------------
template<class T> const vector<T>
Matrix<T>::GetColVector(const int& col)const
{
    if(col < 0 || col >= ColSize)
        throw OutOfRange(-1,col,"at GetColVector");
    vector<T> vec(RowSize);
    for(int i=0;i<RowSize;++i)
        vec[i]=matrix[i][col];
    return vec; 
}
//------------------------------------------------------------------------
template<class T> const Matrix<T>
Matrix<T>::GetRow(const int& row)const
{
    if(row < 0 || row >= RowSize)
      throw OutOfRange(row,-1,"at GetRow");
    Matrix<T> M(1,ColSize);
    for(int i=0;i<ColSize;++i)
        M.SetValue(0,i,GetValue(row,i));
    return M;   
}
//------------------------------------------------------------------------
template<class T> const Matrix<T>
Matrix<T>::GetCol(const int& col)const
{
    if(col < 0 || col >= ColSize)
        throw OutOfRange(-1,col,"at GetCol");
    Matrix<T> M(RowSize,1);
    for(int i=0;i<RowSize;++i)
        M.SetValue(i,0,GetValue(i,col));
    return M;         
}
//------------------------------------------------------------------------
template<class T> void
Matrix<T>::Resize(const int& row,const int& col,T value)
{
    if(row <= 0 || row > SafeSize )
        throw OutOfRange(row,-1,"at Resize");
    if(col <= 0 || col > SafeSize )
        throw OutOfRange(-1,col,"at Resize");
    RowVector.resize(col,value);
    matrix.resize(row,RowVector);
    RowSize = row;
    ColSize = col;   
}
//------------------------------------------------------------------------
template<class T> void
Matrix<T>::SetValue(const int& row,const int& col,T value)
{
    if(row < 0 || row >= RowSize)
        throw OutOfRange(row,-1,"at SetValue");
    if(col < 0 || col >= ColSize)
        throw OutOfRange(-1,col,"at SetValue");
    matrix[row][col]=value;   
}
//------------------------------------------------------------------------
template<class T> void
Matrix<T>::AddRow(T value)
{
    vector<T> vec(matrix[0].size(),value);
    matrix.push_back(vec);
    ++RowSize;
}