matrix6.h

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

//------------------------------------------------------------------------
template<class T> void
Matrix<T>::AddCol(T value)
{
    for(int i=0;i<matrix.size();++i)
        matrix[i].push_back(value);
    ++ColSize;
}
//------------------------------------------------------------------------
template<class T> void
Matrix<T>::DeleteRow(const int& row)
{
    if(row < 0 || row >= RowSize )
        throw OutOfRange(row,-1,"at DeleteRow");
    matrix.erase(matrix.begin() + row);
    --RowSize;
}
//------------------------------------------------------------------------
template<class T> void
Matrix<T>::DeleteCol(const int& col)
{
    if(col < 0 || col >= ColSize )
        throw OutOfRange(-1,col,"at DeleteCol");
    for(int i=0;i<RowSize;++i)
        matrix[i].erase(matrix[i].begin() + col);
    --ColSize;
}
//------------------------------------------------------------------------
template<class T> void
Matrix<T>::InsertRow(const int& row,T value)
{
    if(row < 0 || row >= RowSize )
        throw OutOfRange(row,-1,"at InsertRow");
    vector<T> vec(ColSize,value);
    matrix.insert(matrix.begin()+row,vec);
    ++RowSize;
}
//------------------------------------------------------------------------
template<class T> void  
Matrix<T>::InsertCol(const int& col,T value)
{
    if(col < 0 || col >= ColSize )
        throw OutOfRange(-1,col,"at InsertCol");
    vector<T> vec(RowSize,value);
    for(int i=0;i<RowSize;++i)
        matrix[i].insert(matrix[i].begin()+col,value);
    ++ColSize;   
}
//------------------------------------------------------------------------
template<class T> void  
Matrix<T>::SetRowVector(const int& row,const vector<T>& vec)
{ 
    vector<T> tempvec(vec);
    if(row >= RowSize || row < 0 )
        throw OutOfRange(row,-1,"at SetRowVector");
    if(vec.size()!=(unsigned int)matrix[0].size())
        tempvec.resize((unsigned int)matrix[0].size());
    swap(matrix[row],tempvec);
}
//------------------------------------------------------------------------   
template<class T> void  
Matrix<T>::SetColVector(const int& col,const vector<T>& vec)
{
    if(col >= ColSize || col < 0 )
        throw OutOfRange(-1,col,"at SetColVector");
    vector<T> tempvec(vec);   
    if(vec.size()!=(unsigned int)matrix.size())
        tempvec.resize((unsigned int)matrix.size());
    for(int i = 0 ;i < RowSize ;++i)
        matrix[i][col] = tempvec[i];   
}
//------------------------------------------------------------------------   
template<class T> void  
Matrix<T>::SwapRow(const int& row1,const int& row2)
{
    if(row1 >= RowSize || row1 < 0 )
        throw OutOfRange(row1,-1,"at SwapRow");
    if(row2 >= RowSize || row2 < 0 )
        throw OutOfRange(row2,-1,"at SwapRow");
    if(row1 != row2)
        swap(matrix[row1],matrix[row2]);
}
//------------------------------------------------------------------------   
template<class T> void  
Matrix<T>::SwapCol(const int& col1,const int& col2)
{
    if(col1 >= ColSize || col1 < 0)
        throw OutOfRange(-1,col1,"at SwapCol");
    if(col2 >= ColSize || col2 < 0)
        throw OutOfRange(-1,col2,"at SwapCol");   
    if(col1 != col2)
    {
        vector<T> vec1 = GetColVector(col1);
        vector<T> vec2 = GetColVector(col2);
        swap(vec1,vec2);
        SetColVector( vec1 , col1);
        SetColVector( vec2 , col2);
    }
}
//-----------------------------------------------------------------------
template<class T> Matrix<T>
Matrix<T>::Turn()
{
    Matrix<T> temp(ColSize,RowSize);
    for(int i=0;i<temp.GetRowSize();++i)
       for(int j=0;j<temp.GetColSize();++j)
            temp.SetValue(i,j,GetValue(j,i));           
    return temp;
}
//--------------------------------------------------------------------
// Inverse function
template<class T> Matrix<T>
Matrix<T>::Inverse()
{
    T t = Det();
    if( t == T(0) )
       throw MisMatch("Error in Inverse(): The det of the matrix is zero .");  
    Matrix<T> temp(RowSize,ColSize);
    for( int i=0;i<RowSize;++i)
        for(int j=0;j<ColSize;++j)
        {
            T mi(1),n(-1);
            for(int k=0;k<i+j;++k)
                mi *= n;
            temp.SetValue(j,i,GetAij(i,j).Det()*mi);
        }   
    temp=(T(1)/t)*temp;
    return temp;
}
//--------------------------------------------------------------------
template<class T> Matrix<T>
Matrix<T>::Unit(const int& row)
{
    if(row<=0||row>SafeSize)
        throw OutOfRange(row,-1,"at Unit");
    Matrix<T> M(row,row);
    for(int i=0;i<row;++i)
        M.SetValue(i,i,T(1));
    return M;           
}
//------------------------------------------------------------------------
template<class T> Matrix<T>
Matrix<T>::Diag()
{
    if(RowSize!=ColSize)
        throw MisMatch("Error in Diag(): Row not equal to Col");
    Matrix<T> M(RowSize,1);
    for(int i=0;i<RowSize;++i)
        M.SetValue(i,0,GetValue(i,i));
    return M;
}
//--------------------------------------------------------------------
// Det function [determinant]
template<class T>
T Matrix<T>::Det()
{
    if(RowSize != ColSize)
        throw MisMatch("Error in Det(): The matrix is not square");
    if(RowSize == 1)
        return GetValue(0,0);
    if(RowSize == 2)
        return GetValue(1,1)*GetValue(0,0)-GetValue(1,0)*GetValue(0,1);
    T t(0);
    for(int i=0;i<RowSize;++i)
    {
        T mi(1),n(-1);
        for(int j=0;j<i;++j)
            mi *= n;
        t += GetValue(0,i)*GetAij(0,i).Det()*mi;
    }
    return t;
}
//--------------------------------------------------------------------
// Trace function
template<class T>
T Matrix<T>::Trace()
{  
    if(RowSize != ColSize)
        throw MisMatch("Error in Trace(): The matrix is not square");
     T t=T(0);
     for(int i=0;i<RowSize;++i)
         t += GetValue(i,i);
     return t;
}
//--------------------------------------------------------------------
// Norm function [only float,double long ,double is support]
template<class T> T
Matrix<T>::Norm(const int& p)
{
    return sqrt((((*this).Turn())*(*this)).Trace());
}
//------------------------------------------------------------------------
template<class T> T
Matrix<T>::CondNumber()
{
    Matrix<T> temp;
    temp = Inverse();
    return temp.Norm()*Norm();
}
//------------------------------------------------------------------------
template<class T> T
Matrix<T>::Radius()
{
    if(RowSize!=ColSize)
        throw MisMatch("Error in Radius() : RowSize not equal to Colszie");
    Matrix<T> A(RowSize,1,1),C,B;
    T U(0),RMD(1);
    while(1)
     {
        C=(*this)*A;
        B=C*(T(1)/C.Norm());
        RMD=(B.Turn()*C).GetValue(0,0);
        if(abs(RMD-U)<1e-6)
            break;
        U=RMD;
        A=B;
     }
     return RMD;
}
//------------------------------------------------------------------------
template<class T> int
Matrix<T>::Rank()
{
    int rank=1;
    Matrix<T> M(*this),Zero(RowSize,ColSize);
    if(M==Zero)
        return 0;
    if(RowSize>ColSize)
        M=M.Turn();
    if(M.GetRowSize()==1)
        return 1;
    bool allzero=true;   
    for(int i=0;i<M.GetRowSize();++i)
        if(M.GetValue(i,0)!=T(0))
           {
            allzero=false;
            break;
           }
    if(allzero)
       {
        M.DeleteCol(0);
        return M.Rank();
       }    
    if(M.GetValue(0,0)==T(0))
    {
         for(int i=1;i<M.GetRowSize();++i)
         {
          if(M.GetValue(i,0)!=T(0))
            {
             M.SwapRow(0,i);
             break;
            }
         }
    }
 for( i=1;i<M.GetRowSize();++i)
 {
  T first=M.GetValue(i,0);
  for( int j=0;j<M.GetColSize();++j)
   {
  T value=(M.GetValue(i,j)-first*M.GetValue(0,j)/M.GetValue(0,0));
  M.SetValue(i,j,value);
   }
 }
 rank+=M.GetAij(0,0).Rank();
 return rank;
}
//------------------------------------------------------------------------
// GetAij  function
template<class T> Matrix<T>
Matrix<T>::GetAij(const int& row ,const int& col)
{
    if(RowSize==1&&ColSize==1)
        throw MisMatch("Error in GetAij(): It should be a matrix");
    if(row < 0 || row >= RowSize)
        throw OutOfRange(row,-1,"at GetAij");
    if(col < 0 || col >= ColSize)
        throw OutOfRange(-1,col,"at GetAij");
    Matrix<T> N(*this);
    N.DeleteRow(row);
    N.DeleteCol(col);
    return N ;
}
//------------------------------------------------------------------------
template<class T> CMatrix< complex<double> >
Matrix<T>::ToCDouble()
{
    CMatrix< complex<double> > C(RowSize,ColSize);
    for(int i=0;i<RowSize;++i)
        for(int j=0;j<ColSize;++j)
            {
             complex<double> cd(GetValue(i,j));
             C.SetValue(i,j,cd);
            }
    return C;
}
//------------------------------------------------------------------------
// binary operator "*"
template<class T>  Matrix<T>
operator * (const Matrix<T>& M,T t)
{
    Matrix<T> temp(M.GetRowSize(),M.GetColSize());
    T d;
    for(int i=0;i<M.GetRowSize();++i)
        for(int j=0;j<M.GetColSize();++j)
        {
         d = M.GetValue(i,j)*t;
         temp.SetValue(i,j,d);
        }
    return   temp;
}
//--------------------------------------------------------------------
// binary operator "*"
template<class T> Matrix<T>
operator * ( T t,const Matrix<T>& M)
{
    return M*t;
}
//-----------------------------------------------------------------------
template<class T> ostream&
operator <<(ostream& o,const Matrix<T>& M)
{
    for(int i=0;i<M.GetRowSize();++i)
       {
        for(int j=0;j<M.GetColSize();++j)
                o<<M.GetValue(i,j)<<"  ";
        o<<endl;               
       }
    return o;   
}
//--------------------------------------------------------------------
template<class T> Matrix<T>
Kronecker(const Matrix<T>& M1,const Matrix<T>& M2)
{
    Matrix<T> M(M1.GetRowSize()*M2.GetRowSize(),M1.GetColSize()*M2.GetColSize());
    for(int i=0;i<M1.GetRowSize();++i)
        for(int j=0;j<M1.GetColSize();++j)
        {
         Matrix<T> temp(M2);
         temp=M1.GetValue(i,j)*temp;
         for(int k=0;k<M2.GetRowSize();++k)
            for(int l=0;l<M2.GetColSize();++l)
                M.SetValue(i*M2.GetRowSize()+k,j*M2.GetColSize()+l,temp.GetValue(k,l));
        }
    return M;   
}   
//-----------------------------------------------------------------------
    /**********RMatrix**********/
//-----------------------------------------------------------------------
template<class T >
class RMatrix : public Matrix<T>{
public:
    RMatrix():Matrix<T>(){};
    RMatrix(int row,int col,T value=T(0)):Matrix<T>(row,col,value){}
    RMatrix(RMatrix<T>& N):Matrix<T>(N){}
// operator override
    RMatrix<T>& operator = (const Matrix<T>&);
// convert function
    RMatrix<double> ToDouble();
// JACOBI