src/sysmat.cpp
/******************************************                                     
*  SOLVING A LINEAR MATRIX SYSTEM AX = B  *
*  with Gauss-Jordan method using full    *
*  pivoting at each step. During the pro- *
*  cess, original AA and BB matrices are  *
*  destroyed to spare storage location.   *
* --------------------------------------- *
* INPUTS:    AA   MATRIX N*N              *                                     
*            BB   MATRIX N*M              *                                     
* --------------------------------------- *                                     
* OUTPUS:    AA   INVERSE OF AA N*N       *                                     
*            DET  DETERMINANT OF AA       *                                     
*            BB   SOLUTION MATRIX N*M     *                                     
* --------------------------------------- *                                     
* NOTA - If M=0 inversion of AA matrix    *
*        only.                            *                                     
******************************************/
SOLVING A LINEAR MATRIX SYSTEM AX = B
with Gauss-Jordan method using full
pivoting at each step. During the pro
cess, original AA and BB matrices are
destroyed to spare storage location.


SOLVING A LINEAR MATRIX SYSTEM AX = B with Gauss - Jordan method using full pivoting at each step . During the pro cess original AA and BB matrices are destroyed to spare storage location .
void MATINV(int N,int M,MAT AA,MAT1 BB,REAL *DET)   { 	
      REAL PC[NMAX],PL[NMAX],CS[NMAX];                                                  
      REAL PV,PAV,temp,TT;   
      int I,IK,J,JK,K;
      *DET= 1.0; 
      for (I=0; I<N; I++)  {                                                                
        PC[I]= 0.0;                                                                
        PL[I]= 0.0;                                                                
        CS[I]= 0.0;                                                                
      }              
      for (K=0; K<N; K++) {                                                                  
        PV=AA[K][K];                                                              
        IK=K;                                                                    
        JK=K;                                                                    
        PAV= ABS(PV);                                                            
        for (I=K; I<N; I++)                                                                
          for (J=K; J<N; J++)  {     
            temp = ABS(AA[I][J]);                                                        
            if (temp > PAV) {                                      
              PV=AA[I][J];                                                        
              PAV= ABS(PV);                                                      
              IK=I;                                                              
              JK=J;                                                              
            }                                                               
          }                                                                 
	PC[K]=JK;                                                                
        PL[K]=IK;                                                                
        if (IK!=K) *DET=-*DET;                                                   
        if (JK!=K) *DET=-*DET;                                                   
        *DET=*DET*PV;  
        temp= ABS(*DET);                                                            
        if (temp < MACH_EPS) {                                          
          return;                                                                  
        }                                                                   
        if(IK!=K)                                                         
          for (I=0; I<N; I++) {                                                              
            TT=AA[IK][I];                                                         
            AA[IK][I]=AA[K][I];                                                    
            AA[K][I]=TT;                                                          
          }                                                                 
        if(M!=0)                                                           
          for (I=0; I<M; I++) {                                                               
            TT=BB[IK][I];                                                           
            BB[IK][I]=BB[K][I];                                                      
            BB[K][I]=TT;                                                            
          }                                                                   
        if(JK!=K)                                                         
          for (I=0; I<N; I++) {                                                              
            TT=AA[I][JK];                                                         
            AA[I][JK]=AA[I][K];                                                    
            AA[I][K]=TT;                                                          
          }                                                                 
        for (I=0; I<N; I++) {                                                                
          CS[I]=AA[I][K];                                                         
          AA[I][K]= 0.0;                                                          
        }                                                                   
        CS[K]= 0.0;                                                                
        AA[K][K]= 1.0;                                                              
        temp= ABS(PV);                                          
        if(temp < MACH_EPS) {                                            
          return;                                                                  
        }                                                                   
        for (I=0; I<N; I++)                                                                
          AA[K][I]=AA[K][I]/PV;                                                    
        if (M!=0)                                                         
          for (I=0; I<M; I++)                                                             
            BB[K][I]=BB[K][I]/PV;                                                  
        for (J=0; J<N; J++) {                                                                
          if (J==K) J++;                                                  
          for (I=0; I<N; I++)                                                              
            AA[J][I]=AA[J][I]-CS[J]*AA[K][I];
          if (M!=0)                                                       
            for (I=0; I<M; I++)                                                          
              BB[J][I]=BB[J][I]-CS[J]*BB[K][I];                                     
        }                                                                   
      } 
      for (I=N-1; I>=0; I--) {                                                               
        IK=(int) PC[I];                                                                
        if (IK==I) goto fin1;                                                   
        for (J=0; J<N; J++) {                                                                
          TT=AA[I][J];                                                            
          AA[I][J]=AA[IK][J];                                                      
          AA[IK][J]=TT;                                                           
        }
        if (M!=0)                                                         
          for (J=0; J<M; J++) {                                                             
            TT=BB[I][J];                                                          
            BB[I][J]=BB[IK][J];                                                    
            BB[IK][J]=TT;                                                         
          }                                                                 
        fin1: ;     			  
      } 
      for (J=N-1; J>=0; J--) {                                                                         
        JK=(int) PL[J];                                                                
        if (JK==J) goto fin2;                                                   
        for (I=0; I<N; I++) {                                                                
          TT=AA[I][J];                                                            
          AA[I][J]=AA[I][JK];                                                      
          AA[I][JK]=TT;                                                           
        } 
        fin2: ;                                                                  
      }                                                                     
      return;  
  }
void MATINV ( int N , int M , MAT AA , MAT1 BB , REAL * DET ) { REAL PC [ NMAX ] , PL [ NMAX ] , CS [ NMAX ] ; REAL PV , PAV , temp , TT ; int I , IK , J , JK , K ; * DET = 1.0 ; for ( I = 0 ; I < N ; I ++ ) { PC [ I ] = 0.0 ; PL [ I ] = 0.0 ; CS [ I ] = 0.0 ; } for ( K = 0 ; K < N ; K ++ ) { PV = AA [ K ] [ K ] ; IK = K ; JK = K ; PAV = ABS ( PV ) ; for ( I = K ; I < N ; I ++ ) for ( J = K ; J < N ; J ++ ) { temp = ABS ( AA [ I ] [ J ] ) ; if ( temp > PAV ) { PV = AA [ I ] [ J ] ; PAV = ABS ( PV ) ; IK = I ; JK = J ; } } PC [ K ] = JK ; PL [ K ] = IK ; if ( IK != K ) * DET = - * DET ; if ( JK != K ) * DET = - * DET ; * DET = * DET * PV ; temp = ABS ( * DET ) ; if ( temp < MACH_EPS ) { return ; } if ( IK != K ) for ( I = 0 ; I < N ; I ++ ) { TT = AA [ IK ] [ I ] ; AA [ IK ] [ I ] = AA [ K ] [ I ] ; AA [ K ] [ I ] = TT ; } if ( M != 0 ) for ( I = 0 ; I < M ; I ++ ) { TT = BB [ IK ] [ I ] ; BB [ IK ] [ I ] = BB [ K ] [ I ] ; BB [ K ] [ I ] = TT ; } if ( JK != K ) for ( I = 0 ; I < N ; I ++ ) { TT = AA [ I ] [ JK ] ; AA [ I ] [ JK ] = AA [ I ] [ K ] ; AA [ I ] [ K ] = TT ; } for ( I = 0 ; I < N ; I ++ ) { CS [ I ] = AA [ I ] [ K ] ; AA [ I ] [ K ] = 0.0 ; } CS [ K ] = 0.0 ; AA [ K ] [ K ] = 1.0 ; temp = ABS ( PV ) ; if ( temp < MACH_EPS ) { return ; } for ( I = 0 ; I < N ; I ++ ) AA [ K ] [ I ] = AA [ K ] [ I ] / PV ; if ( M != 0 ) for ( I = 0 ; I < M ; I ++ ) BB [ K ] [ I ] = BB [ K ] [ I ] / PV ; for ( J = 0 ; J < N ; J ++ ) { if ( J == K ) J ++ ; for ( I = 0 ; I < N ; I ++ ) AA [ J ] [ I ] = AA [ J ] [ I ] - CS [ J ] * AA [ K ] [ I ] ; if ( M != 0 ) for ( I = 0 ; I < M ; I ++ ) BB [ J ] [ I ] = BB [ J ] [ I ] - CS [ J ] * BB [ K ] [ I ] ; } } for ( I = N - 1 ; I >= 0 ; I -- ) { IK = ( int ) PC [ I ] ; if ( IK == I ) goto fin1 ; for ( J = 0 ; J < N ; J ++ ) { TT = AA [ I ] [ J ] ; AA [ I ] [ J ] = AA [ IK ] [ J ] ; AA [ IK ] [ J ] = TT ; } if ( M != 0 ) for ( J = 0 ; J < M ; J ++ ) { TT = BB [ I ] [ J ] ; BB [ I ] [ J ] = BB [ IK ] [ J ] ; BB [ IK ] [ J ] = TT ; } fin1 : ; } for ( J = N - 1 ; J >= 0 ; J -- ) { JK = ( int ) PL [ J ] ; if ( JK == J ) goto fin2 ; for ( I = 0 ; I < N ; I ++ ) { TT = AA [ I ] [ J ] ; AA [ I ] [ J ] = AA [ I ] [ JK ] ; AA [ I ] [ JK ] = TT ; } fin2 : ; } return ; }
SOLVING A LINEAR MATRIX SYSTEM AX = B
with Gauss-Jordan method using full
pivoting at each step.
SOLVING A LINEAR MATRIX SYSTEM AX = B with Gauss - Jordan method using full pivoting at each step .
// Initializations :                                                              // Main loop                                                                          // Searching greatest pivot :                                                // Search terminated, the pivot is in location I=IK, J=JK.                      // Memorizing pivot location:  		   // DETERMINANT DET is actualised                                               // If DET=0, ERROR MESSAGE and STOP                                            // Machine dependant EPSMACH equals here 1e-20 //          fprintf(fp2,"\n  The determinant equals ZERO !!!\n");                                                               // POSITIONNING PIVOT IN K,K:		 // EXCHANGE LINES IK and K of matrix AA:			 // PIVOT is at correct line		   // EXCHANGE COLUMNS JK and K of matrix AA: // The PIVOT is at correct column.                                               // and is located in K,K.                                                    // Column K of matrix AA is stored in CS vector                              // then column K is set to zero. 		   // Line K of matrix AA is modified:		 //          fprintf(fp2,"\n  PIVOT TOO SMALL - STOP\n");                                // Other lines of matrix AA are modified:		   // Line J of matrix AA is modified: // Line K is ready //End of K loop // MATRIX AA INVERSION IS DONE - REARRANGEMENT OF MATRIX AA // EXCHANGE LINES                                                                         // EXCHANGE LINES I and PC(I) of matrix AA: // NO MORE EXCHANGE is NECESSARY                                                       // Go to next line                                                  // for i                                                                      // EXCHANGE COLUMNS:  	   // EXCHANGE COLUMNS J ET PL(J) of matrix AA: // NO MORE EXCHANGE is NECESSARY                                                       // Go to next column. // REARRANGEMENT TERMINATED.                                                        