matrix.cpp

贝叶斯算法实现

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()