tmatrix.h

是一个matrix类

       vector<double>::iterator  maxloc = max_element(vt.begin(),vt.end());//Find the max element in the column:
       if((*maxloc + 1.0) == 1.0){
          CERR("Sigular matrix!Any key to quit!");
       }
       s = maxloc - vt.begin();
       if(s)                                    //swap the line which have the max elem;
           {Mt.SwapRows(i,s+i); I.SwapRows(i,s+i);}
       for(int j = 1;j+i < rowCount ;j++){      //Elementary operation:
            k = Mt.ReadCells(j+i,i)/Mt.ReadCells(i,i);
            for(int u = 0;u < colCount;u++){
                vj[u] = Mt.ReadCells(i+j,u) - k * Mt.ReadCells(i,u);
                vi[u] = I.ReadCells(i+j,u) - k * I.ReadCells(i,u);
            }
            Mt.SetRows(i+j,vj);
            I.SetRows(i+j,vi);
       }
   } //The Mt matrix now is a  upper triangular matrix

   for(int i = rowCount - 1;i >= 0 ;i--){
      for(int j = 1;i-j >= 0;j++){
         k = Mt.ReadCells(i - j,i)/Mt.ReadCells(i,i);
         for(int u = 0;u < colCount;u++)
            vi[u] =I.ReadCells(i-j,u) - k*I.ReadCells(i,u);
         I.SetRows(i-j,vi);
      }
   }
   for(int i = 0;i < colCount;i++){
      k = Mt.ReadCells(i,i);
      vi = I.ReadRows(i);
      transform(vi.begin(),vi.end(),vi.begin(),bind2nd(divides<complex_double_type>(),k));
      I.SetRows(i,vi);
   }
   return I;       //Inverse matrix
}

//---------------------------------------------------------------------------
//"~" operator to get the turn of a matrix:
template<class type> matrix<type>
matrix<type>::operator~(){
    matrix<type> MT(colCount,rowCount);
    for(int i = 0;i< MT.RowCount() ;i++)
        for(int j = 0;j< MT.ColCount() ;j++){
            MT.SetCells(i,j,this->ReadCells(j,i));
        }
    return MT;
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
// Member Binary operators:
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//+= opreator:
template<class type> matrix<type>&
matrix<type>::operator +=(const matrix<type>& M){
    *this = operator+(*this,M);
    return *this;
}
//---------------------------------------------------------------------------
//-= operator:
template<class type> matrix<type>&
matrix<type>:: operator -=(const matrix<type>& M){
    *this = operator-(*this,M);
    return *this;
}
//---------------------------------------------------------------------------
//*= operator:Matrix multply
template<class type> matrix<type>&
matrix<type>::operator *=(const matrix<type>& M){
    *this = operator*(*this,M);
    return *this;
}
//*= operator:Matrix multply with const type:
template<class type> matrix<type>&
matrix<type>::operator *=(const type& T){
    *this = operator*(*this,T);
    return *this;
}
//---------------------------------------------------------------------------
//"^" matrix pow:
template<class type> matrix<type>
matrix<type>::operator ^(const int P){
    matrix<type> Mt = *this;
    for(int i = 1;i < P ;i++){
        Mt *=(*this);
    }
    return Mt;
}
//---------------------------------------------------------------------------
template<class type>  bool
matrix<type>::operator ==(const matrix<type>& M)
{
    if(M.RowCount != rowCount || M.ColCount != colCount)
        return false;
    else
        for(int i = 0;i< rowCount ;i++)
            for(int j = 0;j< colCount ;j++){
                if(Matrix[i][j] != M.ReadCells(i,j))
                    return false;
            }
     return true;
}
//---------------------------------------------------------------------------
template<class type>  bool
matrix<type>::operator !=(const matrix<type>& M )
{
    if(M.RowCount != rowCount || M.ColCount != colCount)
        return true;
    else
        for(int i = 0;i< rowCount ;i++)
            for(int j = 0;j< colCount ;j++){
                if(Matrix[i][j] != M.ReadCells(i,j))
                    return true;
            }
     return false;
}
//---------------------------------------------------------------------------
template<class type> const vector<type>
matrix<type>:: operator[ ](const int& RowIndex) const{
    if(RowIndex >= rowCount && RowIndex < 0)
        CERR("Invalid row index!Any key to quit");
    return Matrix[RowIndex];
};
//---------------------------------------------------------------------------
//***************************************************************************
//---------------------------------------------------------------------------
//Non-member functions;Binary operators:
//---------------------------------------------------------------------------
//***************************************************************************
//---------------------------------------------------------------------------

//"+" operator: Plus
//---------------------------------------------------------------------------

template<class type> matrix<type>
operator+(const matrix<type>& MLeft,const matrix<type>& MRight)
{
    vector<type> v1,v2,v3;
    //Ensure MLeft and MRight have the same dims:
    if(MLeft.RowCount()!= MRight.RowCount() || MLeft.ColCount() != MRight.ColCount() ){
       CERR("Matrixs are not conformable for add.Any key to quit!")
    }
    //"+" operation:
    matrix<type> Mt(MLeft.RowCount(),MLeft.ColCount());
    for(int i = 0;i <MLeft.RowCount();i++){
        v1 = MLeft.ReadRows(i);
        v2 = MRight.ReadRows(i);
        v3 = Mt.ReadRows(i);
        transform(v1.begin(),v1.end(),v2.begin(),v3.begin(),plus<type>());
        Mt.SetRows(i,v3);
        }
    return Mt;
}

//---------------------------------------------------------------------------
//"-" operator: Minus
//---------------------------------------------------------------------------
template<class type> matrix<type>
operator-(const matrix<type>& MLeft,const matrix<type>& MRight)
{
     vector<type> v1,v2,v3;
     //Ensure MLeft and MRight have the same dims:
    if(MLeft.RowCount()!= MRight.RowCount() || MLeft.ColCount() != MRight.ColCount() ){
       CERR("Matrixs are not conformable for minus.Any key to quit!")
    }
    //"-" operation:
    matrix<type> Mt(MLeft.RowCount(),MLeft.ColCount());
    for(int i = 0;i <MLeft.RowCount();i++){
        v1 = MLeft.ReadRows(i);
        v2 = MRight.ReadRows(i);
        v3 = Mt.ReadRows(i);
        transform(v1.begin(),v1.end(),v2.begin(),v3.begin(),minus<type>());
        Mt.SetRows(i,v3);
        }
    return Mt;
}
//---------------------------------------------------------------------------
//"*" operator: Multiplay
//---------------------------------------------------------------------------
//A Matrix multiplays another.
template<class type> matrix<type>
operator*(const matrix<type>& MLeft,const matrix<type>& MRight)
{
    vector<type> v1,v2,v3(MLeft.ColCount()),v4(1,0);
    if(MLeft.ColCount() != MRight.RowCount()){
         CERR("Matrixs are not conformable for multiplication.Any key to quit!");
    }
    matrix<type> Mt(MLeft.RowCount(),MRight.ColCount());
    for(int i = 0;i< MLeft.RowCount() ;i++){
        v1 = MLeft.ReadRows(i);
        for(int j = 0;j< MRight.ColCount() ;j++,v4[0] = 0){
            v2 = MRight.ReadCols(j);
            transform(v1.begin(),v1.end(),v2.begin(),v3.begin(),multiplies<type>());
            for(int i = 0;i <v3.size();i++){
                v4[0] += v3[i];
            }
            Mt.SetCells(i,j,v4[0]);
        }
    }
    return Mt;
}

//const cell (prefix)multiplay Matrix
template<class type> matrix<type>
operator*(const type& T,const matrix<type>& M)
{
    matrix<type> Mt(M);
    vector<type> v1(Mt.ColCount());
    for(int i = 0;i< Mt.RowCount() ;i++)
        for(int j = 0;j< Mt.ColCount() ;j++){
            v1[j] = M.ReadCells(i,j) * T;
            Mt.SetCells(i,j,v1[j]);
        }
    return Mt;
}

//const cell (postfix)multiplay Matrix
template<class type> matrix<type>
operator*(const matrix<type>& M,const type& T)
{
    return operator*(T,M);

}
//---------------------------------------------------------------------------
//"<<" and ">>" operators: output/input matrix's cells:
template<class type> void
operator<<(ostream& os,const matrix<type>& M){
    for(int i = 0;i < M.RowCount() ;i++){
        for(int j = 0;j < M.ColCount() ;j++){
 //           os.precision(3);
//            os.setf(ios::scientific,ios::floatfield);
            os<<M.ReadCells(i,j)<<"\t";
        }
        os<<"\n";
    }
}

template<class type> void
operator>>(istream& is,const matrix<type>& M){
    type t;
    for(int i = 0;i < M.RowCount() ;i++){
        for(int j = 0;j < M.ColCount() ;j++){

            is>>t;
            M.SetCells(i,j,t);
       }
    }
}
//---------------------------------------------------------------------------//
//***************************************************************************
//Non-member functions:
//***************************************************************************
//---------------------------------------------------------------------------
template<class type> matrix<type>
IdentityMatrix(const int RowCount)
{
    matrix<type> M(RowCount,RowCount);
    for(int i = 0;i < RowCount ;i++){
        M.SetCells(i,i,1);
    }
    return M;
}
//---------------------------------------------------------------------------
//Kronecker: (直积)
template<class type> matrix<type>
kroneck(const matrix<type>& M1,const matrix<type>& M2)
{
    int m,n;
    m = M1.RowCount()* M2.RowCount();
    n = M1.ColCount()* M2.ColCount();
    matrix<type> Mt(M2.RowCount(),M2.ColCount());
    matrix<type> M(m,n);
    M2.RowCount();
    vector<type> vt,v1;
    for(int i = 0;i< M1.RowCount() ;i++){
        for(int j = 0;j< M1.ColCount() ;j++){
            Mt = M1.ReadCells(i,j) * M2;
            for(int k = 0;k < M2.RowCount();k++){
                for(int p = 0;p < M2.ColCount() ;p++){
                    M.SetCells(i*M2.RowCount()+k,j*M2.ColCount()+p,Mt.ReadCells(k,p));
                }
            }
        }
    }
    return M;
}

#endif  //TMATRIX_H