矩阵计算程序

#ifndef _MATRIX_H
#define  _MATRIX_H

class  Matrix
{
     private :
        int      row;     //  矩阵的行数
         int      col;     //  矩阵的列数
         int      n;       //  矩阵元素个数
         double *     mtx;  //  动态分配用来存放数组的空间
     public :
        Matrix( int  row = 1 ,  int  col = 1 );              //  带默认参数的构造函数
        Matrix( int  row,  int  col,  double  mtx[]);     //  用数组创建一个矩阵
        Matrix( const  Matrix  & obj);                  //  copy构造函数
        ~ Matrix() { delete[]  this -> mtx; }        

         void     print() const ;                      //  格式化输出矩阵
         int      getRow() const  {  return   this -> row; } //  访问矩阵行数
         int      getCol() const  {  return   this -> col; } //  访问矩阵列数
         int      getN() const    {  return   this -> n;   } //  访问矩阵元素个数
         double *  getMtx() const  {  return   this -> mtx; } //  获取该矩阵的数组
        
         //  用下标访问矩阵元素
         double       get ( const   int  i,  const   int  j) const ; 
         //  用下标修改矩阵元素值
         void       set ( const   int  i,  const   int  j,  const   double  e); 

         //  重载了一些常用操作符，包括 +，-，x，=，负号，正号，
         //  A = B
        Matrix  & operator =  ( const  Matrix  & obj);
         //  +A
        Matrix   operator +  () const  {  return   * this ; }
         //  -A
        Matrix   operator -  () const ;
         //  A + B
        friend  Matrix   operator +  ( const  Matrix  & A,  const  Matrix  & B);
         //  A - B
        friend  Matrix   operator -  ( const  Matrix  & A,  const  Matrix  & B);
         //  A * B 两矩阵相乘
        friend  Matrix   operator *  ( const  Matrix  & A,  const  Matrix  & B);
         //  a * B 实数与矩阵相乘
        friend  Matrix   operator *  ( const   double   & a,  const  Matrix  & B);
         //  A 的转置
        friend  Matrix  trv( const  Matrix  & A);
         //  A 的行列式值，采用列主元消去法
         //  求行列式须将矩阵化为三角阵,此处为了防止修改原矩阵，采用传值调用
        friend   double   det(Matrix A);
         //  A 的逆矩阵，采用高斯-若当列主元消去法
        friend  Matrix  inv(Matrix A);
};
#endif         

实现
#include < iostream.h >
#include < math.h >                 
#include < stdlib.h >                 
#include < iomanip.h >              
#include " matrix.h "

//  带默认参数值的构造函数
//  构造一个row行col列的零矩阵
Matrix::Matrix( int  row,  int  col) {                
    this->row = row;
        this->col = col;                                                         
    this->n   = row * col;
        this->mtx = new double[n];
                                               
    for(int i=0; i<n; i++)
           this->mtx[i] = 0.0;
        }                                                                                     
//  用一个数组初始化矩阵                                                        
Matrix::Matrix( int  row,  int  col,  double  mtx[])                                  
{
    this->row = row;
        this->col = col;                                                         
    this->n   = row * col;    
        this->mtx = new double[n];                                               
    
        for(int i=0; i<n; i++)                                                     
       this->mtx[i] = mtx[i];                                                
}                                                                                                   
//  拷贝构造函数，因为成员变量含有动态空间，防止传递参数 
//  等操作发生错误
Matrix::Matrix( const  Matrix  & obj) {            
    this->row = obj.getRow();                                         
    this->col = obj.getCol();                                          
    this->n      = obj.getN();                                            
    this->mtx = new double[n];
                                            
    for(int i=0; i<n; i++)
       this->mtx[i] = obj.getMtx()[i];
}                                                                                                    
//  格式化输出矩阵所有元素
void  Matrix::print() const {                
    for(int i=0; i<this->row; i++){                                         
       for(int j=0; j<this->col; j++)
          if(fabs(this->get(i,j)) <= 1.0e-10)
                cout << setiosflags(ios::left) << setw(12) << 0.0 << ' ';
              else
                cout << setiosflags(ios::left) 
                     << setw(12) << this->get(i,j) << ' ';
           cout <<endl;      
    }          
}

//  获取矩阵元素
//  注意这里矩阵下标从(0,0)开始
double  Matrix:: get ( const   int  i,  const   int  j) const {                
    return this->mtx[i*this->col + j];                           
}

//  修改矩阵元素
void  Matrix:: set ( const   int  i,  const   int  j,  const   double  e) {                
    this->mtx[i*this->col + j] = e;                               
}                                                             
                
//  重载赋值操作符，由于成员变量中含有动态分配
Matrix  & Matrix:: operator =  ( const  Matrix  & obj) {                
    if(this == &obj)    // 将一个矩阵赋给它自己时简单做返回即可        
       return *this;
    delete[] this->mtx; // 首先删除目的操作数的动态空间
    this->row = obj.getRow();
    this->col = obj.getCol();
    this->n   = obj.getN();
    this->mtx = new double[n]; // 重新分配空间保存obj的数组

    for(int i=0; i<n; i++)
       this->mtx[i] = obj.getMtx()[i];                              
    return *this;
}    
           
//  负号操作符，返回值为该矩阵的负矩阵，原矩阵不变
Matrix Matrix:: operator -  () const {                
    // 为了不改变原来的矩阵，此处从新构造一个矩阵
         Matrix _A(this->row, this->col);

    for(int i=0; i<_A.n; i++)
       _A.mtx[i] = -(this->mtx[i]);
    return _A;       
}                                                                                                   
//  矩阵求和，对应位置元素相加
Matrix  operator +  ( const  Matrix  & A,  const  Matrix  & B) {                
    Matrix AB(A.row, A.col);

    if(A.row!=B.row || A.col!=B.col){
       cout << "Can't do A+B\n"; // 如果矩阵A和B行列数不一致则不可相加
             exit(0);        
    }

    for(int i=0; i<AB.n; i++)
       AB.mtx[i] = A.mtx[i] + B.mtx[i];
    
        return AB;       
}                                                                                                   
//  矩阵减法，用加上一个负矩阵来实现
Matrix  operator -  ( const  Matrix  & A,  const  Matrix  & B) {                
    return (A + (-B));
}                                                                                                   
//  矩阵乘法
Matrix  operator *  ( const  Matrix  & A,  const  Matrix  & B) {                
    if(A.col != B.row){    // A的列数必须和B的行数一致
       cout << "Can't multiply\n";
       exit(0);          
    }                                                                          

    Matrix AB(A.row, B.col); // AB用于保存乘积

    for(int i=0; i<AB.row; i++)
       for(int j=0; j<AB.col; j++)
          for(int k=0; k<A.col; k++)
        AB.set(i, j, AB.get(i,j) + A.get(i,k)*B.get(k,j));
          return AB;       
}

//  矩阵与实数相乘
Matrix  operator *  ( const   double   & a,  const  Matrix  & B) {
    Matrix aB(B);
    for(int i=0; i<aB.row; i++)
       for(int j=0; j<aB.col; j++)
          aB.set(i,j, a*B.get(i,j));    
    return aB;
}

//  矩阵的转置 将(i,j)与(j,i)互换
//  此函数返回一个矩阵的转置矩阵，并不改变原来的矩阵
Matrix trv( const  Matrix  & A) {
    Matrix AT(A.col, A.row);
    for(int i=0; i<AT.row; i++)
       for(int j=0; j<AT.col; j++)
          AT.set(i, j, A.get(j,i));
    return AT;       
}
                                                                                  
//  矩阵行列式值，采用列主元消去法
double  det(Matrix A) {
    if(A.row != A.col) {    // 矩阵必须为n*n的才可进行行列式求值       
        cout << "error" << endl;
        return 0.0;     //如果不满足行列数相等返回0.0    }

    double detValue = 1.0;  // 用于保存行列式值
    for(int i=0; i<A.getRow()-1; i++){ // 需要n-1步列化零操作

         //------------------ 选主元 ---------------------------------    
       double max = fabs(A.get(i,i));    // 主元初始默认为右下方矩阵首个元素     
        int    ind = i;  
          // 主元行号默认为右下方矩阵首行      
    
       for(int j=i+1; j<A.getRow(); j++){  // 选择列主元
        if(fabs(A.get(j,i)) > max){     // 遇到绝对值更大的元素
           max = fabs(A.get(j,i));      // 更新主元值
           ind = j;                       // 更新主元行号
        }                  
      }//loop j    
    
       //------------------- 移动主元行 -----------------------------
       if(max <= 1.0e-10) return 0.0;    // 右下方矩阵首行为零，显然行列式值为零      
        if(ind != i){// 主元行非右下方矩阵首行                                  
         for(int k=i; k<A.getRow(); k++){     // 将主元行与右下方矩阵首行互换          
         double temp = A.get(i,k);
        A.set(i,k,A.get(ind,k));
        A.set(ind,k,temp);
        }
          detValue = -detValue;             // 互换行列式两行，行列式值反号    
    }                            
    
       //------------------- 消元 ----------------------------------                                  
        for(j=i+1; j<A.getRow(); j++){     // 遍历行  
            double temp = A.get(j,i)/A.get(i,i);        
           for(int k=i; k<A.getRow(); k++) // 遍历行中每个元素，行首置0
              A.set(j, k, A.get(j,k)-A.get(i,k)*temp);        }                    
          detValue *= A.get(i,i); // 每步消元都会产生一个对角线上元素，将其累乘
       }// loop i     
      
      // 注意矩阵最后一个元素在消元的过程中没有被累乘到
      return detValue * A.get(A.getRow()-1, A.getRow()-1);}//det()        
                                  
      // A的逆矩阵 高斯-若当消去法，按列选主元
      Matrix inv(Matrix A){                
    if(A.row != A.col){ // 只可求狭义逆矩阵，即行列数相同
       cout << "Matrix should be N x N\n";
       exit(0);       
    }                                                                                                
    // 构造一个与A行列相同的单位阵B
    Matrix B(A.row,A.col);
      
    for(int r=0; r<A.row; r++)
       for(int c=0; c<A.col; c++)
          if(r == c) B.set(r,c,1.0); 
     
    // 对矩阵A进行A.row次消元运算，每次保证第K列只有对角线上非零
    // 同时以同样的操作施与矩阵B，结果A变为单位阵B为所求逆阵
    for(int k=0; k<A.row; k++){
    //------------------ 选主元 --------------------------------------    
      double max = fabs(A.get(k,k));    // 主元初始默认为右下方矩阵首个元素     
     int    ind = k;            // 主元行号默认为右下方矩阵首行
           // 结果第ind行为列主元行
      for(int n=k+1; n<A.getRow(); n++){
         if(fabs(A.get(n,k)) > max){// 遇到绝对值更大的元素             
        max = fabs(A.get(n,k));    // 更新主元值
                  ind = n;// 更新主元行号
      }
       }                   
       //------------------- 移动主元行 --------------------------------
       if(ind != k){// 主元行不是右下方矩阵首行
           for(int m=k; m<A.row; m++){// 将主元行与右下方矩阵首行互换         
         double tempa = A.get(k,m);
         A.set(k, m, A.get(ind,m));
             A.set(ind, m, tempa);            
          }
      for(m=0; m<B.row; m++){
             double tempb = B.get(k,m); // 对矩阵B施以相同操作                  
   B.set(k, m, B.get(ind,m)); // B与A阶数相同，可在一个循环中               
  B.set(ind, m, tempb);    
          }    
       }        
      //--------------------- 消元 -----------------------------------
      // 第k次消元操作，以第k行作为主元行，将其上下各行的第k列元素化为零
       // 同时以同样的参数对B施以同样的操作，此时可以将B看作A矩阵的一部分          
      for(int i=0; i<A.col; i++){                                 
         if(i != k){
          double Mik = -A.get(i,k)/A.get(k,k);                         
          for(int j=k+1; j<A.row; j++)
                A.set(i, j, A.get(i,j) + Mik*A.get(k,j));                     
          for(j=0; j<B.row; j++)
                B.set(i, j, B.get(i,j) + Mik*B.get(k,j));                   
        }//end if 
      }//loop i  
      double Mkk = 1.0/A.get(k,k);
      
      for(int j=0; j<A.row; j++)
         A.set(k, j, A.get(k,j) * Mkk);
      for(j=0; j<B.row; j++)
         B.set(k, j, B.get(k,j) * Mkk);
     }//loop k      
     return B;
}//inv()

转载于:https://www.cnblogs.com/JoeDZ/archive/2008/06/30/1232741.html