矩阵求值

typedef struct 
{
  double *datap;
  int    lie;//列
  int    hang;//
}Dmatrix_t;

__inline double *DMatrix_ij_P(Dmatrix_t *X,int i,int j )
{
    return ( X->datap+ i*X->lie+ j);    
}
int DMatrix_huanhang(Dmatrix_t *X,int d,int s )
{
    double temp;
    int i;
    for(i=0;i<X->lie;i++)
    {
        temp= *DMatrix_ij_P(X,d,i);
        *DMatrix_ij_P(X,d,i)=*DMatrix_ij_P(X,s,i);
        *DMatrix_ij_P(X,s,i)=temp;    
    }
    return 1;

}
/********************
思路,化为
 a a a a a
 0 a a a a
 0 0 a a a
 0 0 0 a a  
 0 0 0 0 a  
 对角线再相乘得到 

原矩阵被破坏，调用时注意备份
******************/
int DMatrix_Det( Dmatrix_t *X,double *zhi )
{
    
    
    if(X->hang!= X->lie ) return 0;
    if(X->hang==0) return 0;
    if(X->hang==1) {*zhi= *X->datap;return 1;    }
    if(X->hang==2) {*zhi=   X->datap[0]*X->datap[3]- X->datap[1]*X->datap[2];      return 1;    }
      double temp;
     int i,j,k,n;
     k=1;
    for(i=0;i<X->lie;i++ )//列 循环一次把这一列除头全部清零 
    {      
      if( double_cmp_equar(*DMatrix_ij_P(X,i,i ) ,0))//本列不能为 0，从后面换一行过来 ,如果没有 ，直接得0 
      {
          int huanghang=0; 
       for(j=i+1;j<X->hang;j++  )
       {
           if(double_cmp_equar_no(*DMatrix_ij_P(X,j,i ) ,0))
           {
            DMatrix_huanhang(X,i,j  );                 
            k*=-1;       
         huanghang=1;  

     break;    
        }
           
       }
       if(huanghang==0)
       {
           *zhi=0; return 1;
       }
        }  
        
    
     for(j=i+1;j<X->hang;j++)//行 
     {
       temp =   *DMatrix_ij_P(X,j,i  ) / *DMatrix_ij_P(X,i,i  );
         for( n=i;n<X->lie;n++ )//又变列 ，把本行后面的列全处理一遍 
         {
               *DMatrix_ij_P(X,j, n) -=         *DMatrix_ij_P(X,i, n) *temp; 
              if(double_cmp_equar( *DMatrix_ij_P(X,j, n),0))
              {
                   *DMatrix_ij_P(X,j, n)=0.0;
              }
         }
      }
              
    }
    temp=(double)k;
    for(i=0;i< X->lie ;i++)
    {
    temp *=     *DMatrix_ij_P(X,i, i);    
    }

    *zhi=  temp;
  return 1;
    
    
}