Matrix Template Library

最新推荐文章于 2026-06-23 08:00:30 发布
原创 最新推荐文章于 2026-06-23 08:00:30 发布 · 1.4k 阅读 · 0 ·
本内容遵循CC 4.0 BY-SA版权协议
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
标签
#matrix #library #function #class #string
#vector

/*
  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 dbflying@163.com
  Author     :   Dbflying
  Bug report :   dbflying@163.com      
  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 = 100;
//------------------------------------------------------------------------
// 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;
}
//------------------------------------------------------------------------
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 µü´ú·¨Çó½âʵ¶Ô³Æ¾ØÕóÌØÕ÷Öµ¡¢ÌØÕ÷ÏòÁ¿
    RMatrix<T> Jacobi(RMatrix<T>& value=RMatrix<T>());
};
//-----------------------------------------------------------------------
template<class T> RMatrix<T>&
RMatrix<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> RMatrix<double>
RMatrix<T>::ToDouble()
{
    Matrix< double > temp(RowSize,ColSize);
    double d;
    for(int i=0;i<RowSize;++i)
         for(int j=0;j<ColSize;++j)
         {
          d = matrix[i][j];
          temp.SetValue(i,j,d);
         }
    return temp;    
}
//-----------------------------------------------------------------------
template<class T> RMatrix<T>
RMatrix<T>::Jacobi(RMatrix<T>& value)
{
    if(RowSize!=ColSize)
        throw MisMatch("Error in Jacobi(): Jacobi method requres row=col");
    if(*this!=(*this).Turn())
        throw MisMatch("Error in Jacobi() :matrix should equal to turn of matrix.");
    RMatrix<T> M(*this);
    RMatrix<T> R(RowSize,ColSize),R1(R);
    int p,q;   
    R1=R1.Unit(R1.GetRowSize());
    for( int step=0;step<50;++step)
    {   
        p=1;
        q=0;
        T MaxValue=M.GetValue(1,0);
        for(int i=1;i<M.GetRowSize();++i)
            for(int j=0;j<i;++j)
            {
                if(abs(MaxValue)<abs(M.GetValue(i,j)))
                {
                    MaxValue=M.GetValue(i,j);
                    p=i;
                    q=j;
                }
            }
        R=R.Unit(R.GetRowSize());
        T V(M.GetValue(p,q)/(M.GetValue(q,q)-M.GetValue(p,p)));
        T C(1/sqrt(V*V+1));
        T S(C*V);
        R[p][p]=C;
        R[q][q]=C;
        R[p][q]=S;
        R[q][p]=-S;          
        M=R.Turn()*M*R;
        R1=R1*R;
        if(abs(MaxValue)<=1e-6)
            {
                value = R1;
                M=M.Diag();
                return M;
            } 
    }
    value = R1;
    M=M.Diag();
    return M;
}
//-----------------------------------------------------------------------
    /**********CMatrix**********/
//-----------------------------------------------------------------------
template<class T >
class CMatrix : public Matrix<T>{
public:
    CMatrix():Matrix<T>(){};
    CMatrix(int row,int col,T value=T(0)):Matrix<T>(row,col,value){}
    CMatrix(RMatrix<T>& N):Matrix<T>(N){}
// operator override   
    CMatrix<T>& operator = (const Matrix<T>&);   
// special function
    CMatrix<T>  Turn();
    RMatrix<T>  Real();
    RMatrix<T>  Imag();
};
//-----------------------------------------------------------------------
template<class T> CMatrix<T>&
CMatrix<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> CMatrix<T>
CMatrix<T>::Turn()
{
    CMatrix<T> temp(ColSize,RowSize);
    for(int i=0;i<temp.GetRowSize();++i)
       for(int j=0;j<temp.GetColSize();++j)
            temp.SetValue(i,j,conj(GetValue(j,i)));
    return temp;
}
//-----------------------------------------------------------------------
template<class T> RMatrix<T>
CMatrix<T>::Real()
{
    RMatrix<T> R(RowSize,ColSize);
    for(int i=0;i<RowSize;++i)
        for(int j=0;j<ColSize;++j)
            R.SetValue(i,j,GetValue(i,j).real());
    return R;       
}
//-----------------------------------------------------------------------
template<class T> RMatrix<T>
CMatrix<T>::Imag()
{
    RMatrix<T> C(RowSize,ColSize);
    for(int i=0;i<RowSize;++i)
        for(int j=0;j<ColSize;++j)
            C.SetValue(i,j,GetValue(i,j).imag());
    return C;
}
//------------------------------------------------------------------------
}//**namespace soubam**//
//------------------------------------------------------------------------
#endif
dblfly
关注
科学计算库mtl-2.1.2-20.zip
02-21
科学计算库 mtl-2.1.2-20 相当的好啊 不喜欢写数学计算函数的用它
VC matrix template library matrix头文件
09-30
Openvc 用到的matix类 VC matrix template library matrix头文件
C++各类库
HideInTime的博客
09-23 2069
1、C++各大有名库的介绍——C++标准库 2、C++各大有名库的介绍——准标准库Boost 3、C++各大有名库的介绍——GUI 4、C++各大有名库的介绍——网络通信 5、C++各大有名库的介绍——XML 6、C++各大有名库的介绍——科学计算 7、C++各大有名库的介绍——游戏开发 8、C++各大有名库的介绍——线程 9、C++各大有名库的介绍——序列化 10、C++各大有名库的介绍——字符...
常用数值计算库
menjiawan的专栏
08-02 4521
转载自：http://www.xuebuyuan.com/1706852.html 最近在使用数值库进行稀疏矩阵线性方程组的求解，接触了一些数值计算库，主要包括：MKL, GMM++等库。 MKL库是intel的一个运算库，大名鼎鼎的matlab中都使用MKL进行数值运算，但唯一的缺点是它是需要收费的。MKL有30天的试用版，其功能与完整的版本一样。 GMM++库是一个GNU开源库，其无需安
mtl库在GCC编译器下的使用
weixin_34297704的博客
03-26 156
最近一直在改造算法库，将其移植到Linux平台下。使用GCC编译器是发现MTL库中出现一大堆的问题。使用的MTL库下载地址为：http://osl.iu.edu/research/mtl/download.php3。在上面的网页中，需要下载下面三个压缩包，其中第一个是MTL库的源代码，第二个是在Visual Studio平台下编译的补丁包，第三个是gcc编译器的补丁包。http://osl....
MTL矩阵运算库精讲：模板与性能优化实战
weixin_42601134的博客
08-19 1016
MTL（Matrix Template Library）是一个高度优化的C++库，专门用于进行矩阵计算与操作。它提供了一个灵活且效率极高的框架，使得开发者可以轻松地执行各种矩阵相关的任务，从基本的线性代数运算到复杂的数值分析，都能够在MTL中找到相应的实现。稀疏矩阵是只包含少量非零元素的矩阵。在许多应用中，这些矩阵可能大部分元素为零，其中非零元素的数量与矩阵中所有元素的总数相比很少。稀疏矩阵的概念对于任何涉及大规模数据集的问题至关重要，因为它们极大地减少了存储需求和计算量。
Matrix Template Library Help
Soubam Software Workshop
07-30 1397
部分函数功能说明// element access functions    int GetRowSize()const    int GetColSize()const    T   GetValue(const int& row,const int& col)const    const vector GetRowVector(const int& row)const    const vec
与矩阵相关的C++，Java库
sulliy的专栏
10-27 2355
A generic, reusable, and extendable matrix class ： http://www.codeproject.com/KB/architecture/ymatrix.aspx 提供了一个可扩展的框架，比较适合学习，C++ 模板   JAMA， Java Matrix Package,  http://math.nist.gov/javanumeric
MTL--c++科学计算库--介绍
wsh6759的CGIS专栏
04-05 1689
MTL--c++科学计算库--介绍   本文来源是：http://www.osl.iu.edu/research/mtl/ 矩阵模板库(Matrix Template Library,缩写为MTL)是一个提供了针对矩阵格式的复杂的线性代数泛函的可以广泛使用的高性能的泛型组件库。 像标准模板库(STL)一样，MTL使用了五种元素，包括泛型函数，容器，迭代器，适配器以及函数对象。这些都是为了满
950基础矩阵乘法TLA示例
gitblog_01123的博客
06-23 1055
**Note: The community package does not currently support 950 capabilities. Stay tuned for a future supported version.** ## Code Organization ``` ├── 43_ascend950_basic_matmul │   ├── CMakeLists.txt
CANN/catlass分组矩阵乘反量化示例
gitblog_00468的博客
06-23 343
``` ├── 10_grouped_matmul_slice_m_per_token_dequant │ ├── CMakeLists.txt # CMake build file │ ├── README.md │ └── 10_grouped_matmul_slice_m_per_token_dequant.cpp # Main file ``` ## Function Th
matrix-template:强大JavaScript模板引擎
04-06
矩阵模板 强大JavaScript模板引擎 快速开始 安装matrix-template npm i matrix-template 在nodejs中使用它 const matrixTemplate = require('matrix-template'); console.log(matrixTemplate.process('hello ${this.name}', { name: 'zhangchi' })); //hello zhangchi 执照 张池（c）2013年至今
在C++中使用matrix<lib>
01-20
本文讲述了在C++中使用matrix<lib>的详细过程，有助于需要的朋友！
CANN/catlass GEMV示例
最新发布
gitblog_00702的博客
06-23 484
Note: The community package does not support the 950 capability currently. Stay tuned for future versions. ## Code Organization ``` ├── 50_ascend950_basic_matmul_gemv │ ├── CMakeLists.txt # CM
cann/catlass动态优化矩阵乘法示例
gitblog_00891的博客
06-23 873
In the field of high-performance computing and deep learning, matrix multiplication (Matmul) serves as a core operator, making its computational efficiency and generalization capability highly critica
CANN/catlass Tile组件开发详解
gitblog_00809的博客
06-23 846
Tile components are the lowest-level computation and data operation units in the CATLASS template library. They reside at the end of the entire hierarchy and interact directly with the hardware. Their
CANN/catlass示例设计文档
gitblog_00119的博客
06-23 705
This document summarizes the design ideas and code breakdown of certain samples. You can refer to the specific content as needed. - [102_DynamicOptimizedMatmul](https://link.gitcode.com/i/1243b027186
CANN/catlass动态优化量化矩阵乘法示例
gitblog_00041的博客
06-23 430
Based on generalization requirements, this example adds a generalized engineering implementation for quantized Matmul in addition to example 102. This example currently supports the basic template for
C++模板用多了编译报错？手把手教你用CMake跨平台解决MSVC/GCC的bigobj问题
weixin_30631587的博客
06-06 426
本文详细解析C++模板工程中常见的bigobj编译错误问题，提供基于CMake的跨平台解决方案，涵盖MSVC的/bigobj选项和GCC的-Wa,-mbig-obj方案。通过工业级验证的代码示例和最佳实践，帮助开发者高效解决模板实例化导致的文件过大问题，提升项目编译效率。
CANN/catlass量化矩阵乘全载A示例
gitblog_01034的博客
06-23 874
``` ├── 44_quant_matmul_full_loadA_tla │ ├── CMakeLists.txt # CMake build file │ ├── README.md │ ├── 44_quant_matmul_full_loadA_tla.md # Design document │ └── quant_matmul_full_loadA_tla.cpp #
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值