/[escript]/trunk/finley/src/Util.c
ViewVC logotype

Diff of /trunk/finley/src/Util.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 97 by jgs, Tue Dec 14 05:39:33 2004 UTC revision 100 by jgs, Wed Dec 15 03:48:48 2004 UTC
# Line 97  void Finley_Util_SmallMatSetMult(int len Line 97  void Finley_Util_SmallMatSetMult(int len
97      }      }
98  }  }
99    
 /* calcultes the LU factorization for a small matrix dimxdim matrix A */  
 /* TODO: use LAPACK */  
   
 int Finley_Util_SmallMatLU(int dim,double* A,double *LU,int* pivot){  
      double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;  
      int info=0;  
      /* LAPACK version */  
      /* memcpy(LU,A,sizeof(douple)); */  
      /* dgetf2_(&dim,&dim,A,&dim,LU,pivot,&info); */  
      switch(dim) {  
       case 1:  
             D=A[INDEX2(0,0,dim)];  
             if (ABS(D) >0. ){  
                LU[INDEX2(0,0,dim)]=1./D;  
             } else {  
                info=2;  
             }  
          break;  
   
       case 2:  
             A11=A[INDEX2(0,0,dim)];  
             A12=A[INDEX2(0,1,dim)];  
             A21=A[INDEX2(1,0,dim)];  
             A22=A[INDEX2(1,1,dim)];  
   
             D = A11*A22-A12*A21;  
             if (ABS(D) > 0 ){  
                D=1./D;  
                LU[INDEX2(0,0,dim)]= A22*D;  
                LU[INDEX2(1,0,dim)]=-A21*D;  
                LU[INDEX2(0,1,dim)]=-A12*D;  
                LU[INDEX2(1,1,dim)]= A11*D;  
             } else {  
                info=2;  
             }  
          break;  
   
       case 3:  
             A11=A[INDEX2(0,0,dim)];  
             A21=A[INDEX2(1,0,dim)];  
             A31=A[INDEX2(2,0,dim)];  
             A12=A[INDEX2(0,1,dim)];  
             A22=A[INDEX2(1,1,dim)];  
             A32=A[INDEX2(2,1,dim)];  
             A13=A[INDEX2(0,2,dim)];  
             A23=A[INDEX2(1,2,dim)];  
             A33=A[INDEX2(2,2,dim)];  
   
             D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);  
             if (ABS(D) > 0 ){  
                D=1./D;  
                LU[INDEX2(0,0,dim)]=(A22*A33-A23*A32)*D;  
                LU[INDEX2(1,0,dim)]=(A31*A23-A21*A33)*D;  
                LU[INDEX2(2,0,dim)]=(A21*A32-A31*A22)*D;  
                LU[INDEX2(0,1,dim)]=(A13*A32-A12*A33)*D;  
                LU[INDEX2(1,1,dim)]=(A11*A33-A31*A13)*D;  
                LU[INDEX2(2,1,dim)]=(A12*A31-A11*A32)*D;  
                LU[INDEX2(0,2,dim)]=(A12*A23-A13*A22)*D;  
                LU[INDEX2(1,2,dim)]=(A13*A21-A11*A23)*D;  
                LU[INDEX2(2,2,dim)]=(A11*A22-A12*A21)*D;  
             } else {  
                info=2;  
             }  
          break;  
        default:  
             info=1;  
    }  
    return info;  
 }  
   
 /* solves LUx=b where LU is a LU factorization calculated by an Finley_Util_SmallMatLU call */  
 void Finley_Util_SmallMatForwardBackwardSolve(int dim ,int nrhs,double* LU,int* pivot,double* x,double* b) {  
      int i;  
      switch(dim) {  
       case 1:  
          for (i=0;i<nrhs;i++) {  
             x[INDEX2(0,i,dim)]=LU[0]*b[INDEX2(0,i,dim)];  
          }  
          break;  
       case 2:  
          for (i=0;i<nrhs;i++) {  
             x[INDEX2(0,i,dim)]=LU[INDEX2(0,0,dim)]*b[INDEX2(0,i,dim)]+LU[INDEX2(0,1,dim)]*b[INDEX2(1,i,dim)];  
             x[INDEX2(1,i,dim)]=LU[INDEX2(1,0,dim)]*b[INDEX2(0,i,dim)]+LU[INDEX2(1,1,dim)]*b[INDEX2(1,i,dim)];  
          }  
          break;  
   
       case 3:  
          for (i=0;i<nrhs;i++) {  
             x[INDEX2(0,i,dim)]=LU[INDEX2(0,0,dim)]*b[INDEX2(0,i,dim)]+LU[INDEX2(0,1,dim)]*b[INDEX2(1,i,dim)]+LU[INDEX2(0,2,dim)]*b[INDEX2(2,i,dim)];  
             x[INDEX2(1,i,dim)]=LU[INDEX2(1,0,dim)]*b[INDEX2(0,i,dim)]+LU[INDEX2(1,1,dim)]*b[INDEX2(1,i,dim)]+LU[INDEX2(1,2,dim)]*b[INDEX2(2,i,dim)];  
             x[INDEX2(2,i,dim)]=LU[INDEX2(2,0,dim)]*b[INDEX2(0,i,dim)]+LU[INDEX2(2,1,dim)]*b[INDEX2(1,i,dim)]+LU[INDEX2(2,2,dim)]*b[INDEX2(2,i,dim)];  
          }  
          break;  
    }  
    return;  
 }  
100  /*    inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */  /*    inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
101  /*    the determinante is returned. */  /*    the determinante is returned. */
102    
# Line 204  void Finley_Util_InvertSmallMat(int len, Line 108  void Finley_Util_InvertSmallMat(int len,
108        case 1:        case 1:
109           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
110              D=A[INDEX3(0,0,q,dim,dim)];              D=A[INDEX3(0,0,q,dim,dim)];
111              if (ABS(D) > 0 ){              if (D == 0 ){
                det[q]=D;  
                D=1./D;  
                invA[INDEX3(0,0,q,dim,dim)]=D;  
             } else {  
112                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
113                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
114                 return;                 return;
115              }              }
116                det[q]=D;
117                D=1./D;
118                invA[INDEX3(0,0,q,dim,dim)]=D;
119           }           }
120           break;           break;
121    
# Line 224  void Finley_Util_InvertSmallMat(int len, Line 127  void Finley_Util_InvertSmallMat(int len,
127              A22=A[INDEX3(1,1,q,dim,dim)];              A22=A[INDEX3(1,1,q,dim,dim)];
128    
129              D = A11*A22-A12*A21;              D = A11*A22-A12*A21;
130              if (ABS(D) > 0 ){              if (D == 0 ){
                det[q]=D;  
                D=1./D;  
                invA[INDEX3(0,0,q,dim,dim)]= A22*D;  
                invA[INDEX3(1,0,q,dim,dim)]=-A21*D;  
                invA[INDEX3(0,1,q,dim,dim)]=-A12*D;  
                invA[INDEX3(1,1,q,dim,dim)]= A11*D;  
             } else {  
131                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
132                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
133                 return;                 return;
134              }              }
135                det[q]=D;
136                D=1./D;
137                invA[INDEX3(0,0,q,dim,dim)]= A22*D;
138                invA[INDEX3(1,0,q,dim,dim)]=-A21*D;
139                invA[INDEX3(0,1,q,dim,dim)]=-A12*D;
140                invA[INDEX3(1,1,q,dim,dim)]= A11*D;
141           }           }
142           break;           break;
143    
# Line 252  void Finley_Util_InvertSmallMat(int len, Line 154  void Finley_Util_InvertSmallMat(int len,
154              A33=A[INDEX3(2,2,q,dim,dim)];              A33=A[INDEX3(2,2,q,dim,dim)];
155    
156              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);              D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
157              if (ABS(D) > 0 ){              if (D == 0 ){
                det[q]    =D;  
                D=1./D;  
                invA[INDEX3(0,0,q,dim,dim)]=(A22*A33-A23*A32)*D;  
                invA[INDEX3(1,0,q,dim,dim)]=(A31*A23-A21*A33)*D;  
                invA[INDEX3(2,0,q,dim,dim)]=(A21*A32-A31*A22)*D;  
                invA[INDEX3(0,1,q,dim,dim)]=(A13*A32-A12*A33)*D;  
                invA[INDEX3(1,1,q,dim,dim)]=(A11*A33-A31*A13)*D;  
                invA[INDEX3(2,1,q,dim,dim)]=(A12*A31-A11*A32)*D;  
                invA[INDEX3(0,2,q,dim,dim)]=(A12*A23-A13*A22)*D;  
                invA[INDEX3(1,2,q,dim,dim)]=(A13*A21-A11*A23)*D;  
                invA[INDEX3(2,2,q,dim,dim)]=(A11*A22-A12*A21)*D;  
             } else {  
158                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
159                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
160                 return;                 return;
161              }              }
162                det[q]    =D;
163                D=1./D;
164                invA[INDEX3(0,0,q,dim,dim)]=(A22*A33-A23*A32)*D;
165                invA[INDEX3(1,0,q,dim,dim)]=(A31*A23-A21*A33)*D;
166                invA[INDEX3(2,0,q,dim,dim)]=(A21*A32-A31*A22)*D;
167                invA[INDEX3(0,1,q,dim,dim)]=(A13*A32-A12*A33)*D;
168                invA[INDEX3(1,1,q,dim,dim)]=(A11*A33-A31*A13)*D;
169                invA[INDEX3(2,1,q,dim,dim)]=(A12*A31-A11*A32)*D;
170                invA[INDEX3(0,2,q,dim,dim)]=(A12*A23-A13*A22)*D;
171                invA[INDEX3(1,2,q,dim,dim)]=(A13*A21-A11*A23)*D;
172                invA[INDEX3(2,2,q,dim,dim)]=(A11*A22-A12*A21)*D;
173           }           }
174           break;           break;
175    
# Line 505  void Finley_copyDouble(int n,double* sou Line 406  void Finley_copyDouble(int n,double* sou
406    
407  /*  /*
408   * $Log$   * $Log$
409   * Revision 1.2  2004/12/14 05:39:31  jgs   * Revision 1.3  2004/12/15 03:48:47  jgs
410   * *** empty log message ***   * *** empty log message ***
411   *   *
  * Revision 1.1.1.1.2.1  2004/11/12 06:58:19  gross  
  * a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry  
  *  
412   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs   * Revision 1.1.1.1  2004/10/26 06:53:57  jgs
413   * initial import of project esys2   * initial import of project esys2
414   *   *

Legend:
Removed from v.97  
changed lines
  Added in v.100

  ViewVC Help
Powered by ViewVC 1.1.26