/[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 100 by jgs, Wed Dec 15 03:48:48 2004 UTC revision 102 by jgs, Wed Dec 15 07:08:39 2004 UTC
# Line 97  void Finley_Util_SmallMatSetMult(int len Line 97  void Finley_Util_SmallMatSetMult(int len
97      }      }
98  }  }
99    
100    /* calcultes the LU factorization for a small matrix dimxdim matrix A */
101    /* TODO: use LAPACK */
102    
103    int Finley_Util_SmallMatLU(int dim,double* A,double *LU,int* pivot){
104         double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;
105         int info=0;
106         /* LAPACK version */
107         /* memcpy(LU,A,sizeof(douple)); */
108         /* dgetf2_(&dim,&dim,A,&dim,LU,pivot,&info); */
109         switch(dim) {
110          case 1:
111                D=A[INDEX2(0,0,dim)];
112                if (ABS(D) >0. ){
113                   LU[INDEX2(0,0,dim)]=1./D;
114                } else {
115                   info=2;
116                }
117             break;
118    
119          case 2:
120                A11=A[INDEX2(0,0,dim)];
121                A12=A[INDEX2(0,1,dim)];
122                A21=A[INDEX2(1,0,dim)];
123                A22=A[INDEX2(1,1,dim)];
124    
125                D = A11*A22-A12*A21;
126                if (ABS(D) > 0 ){
127                   D=1./D;
128                   LU[INDEX2(0,0,dim)]= A22*D;
129                   LU[INDEX2(1,0,dim)]=-A21*D;
130                   LU[INDEX2(0,1,dim)]=-A12*D;
131                   LU[INDEX2(1,1,dim)]= A11*D;
132                } else {
133                   info=2;
134                }
135             break;
136    
137          case 3:
138                A11=A[INDEX2(0,0,dim)];
139                A21=A[INDEX2(1,0,dim)];
140                A31=A[INDEX2(2,0,dim)];
141                A12=A[INDEX2(0,1,dim)];
142                A22=A[INDEX2(1,1,dim)];
143                A32=A[INDEX2(2,1,dim)];
144                A13=A[INDEX2(0,2,dim)];
145                A23=A[INDEX2(1,2,dim)];
146                A33=A[INDEX2(2,2,dim)];
147    
148                D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
149                if (ABS(D) > 0 ){
150                   D=1./D;
151                   LU[INDEX2(0,0,dim)]=(A22*A33-A23*A32)*D;
152                   LU[INDEX2(1,0,dim)]=(A31*A23-A21*A33)*D;
153                   LU[INDEX2(2,0,dim)]=(A21*A32-A31*A22)*D;
154                   LU[INDEX2(0,1,dim)]=(A13*A32-A12*A33)*D;
155                   LU[INDEX2(1,1,dim)]=(A11*A33-A31*A13)*D;
156                   LU[INDEX2(2,1,dim)]=(A12*A31-A11*A32)*D;
157                   LU[INDEX2(0,2,dim)]=(A12*A23-A13*A22)*D;
158                   LU[INDEX2(1,2,dim)]=(A13*A21-A11*A23)*D;
159                   LU[INDEX2(2,2,dim)]=(A11*A22-A12*A21)*D;
160                } else {
161                   info=2;
162                }
163             break;
164           default:
165                info=1;
166       }
167       return info;
168    }
169    
170    /* solves LUx=b where LU is a LU factorization calculated by an Finley_Util_SmallMatLU call */
171    void Finley_Util_SmallMatForwardBackwardSolve(int dim ,int nrhs,double* LU,int* pivot,double* x,double* b) {
172         int i;
173         switch(dim) {
174          case 1:
175             for (i=0;i<nrhs;i++) {
176                x[INDEX2(0,i,dim)]=LU[0]*b[INDEX2(0,i,dim)];
177             }
178             break;
179          case 2:
180             for (i=0;i<nrhs;i++) {
181                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)];
182                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)];
183             }
184             break;
185    
186          case 3:
187             for (i=0;i<nrhs;i++) {
188                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)];
189                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)];
190                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)];
191             }
192             break;
193       }
194       return;
195    }
196  /*    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 */
197  /*    the determinante is returned. */  /*    the determinante is returned. */
198    
# Line 108  void Finley_Util_InvertSmallMat(int len, Line 204  void Finley_Util_InvertSmallMat(int len,
204        case 1:        case 1:
205           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
206              D=A[INDEX3(0,0,q,dim,dim)];              D=A[INDEX3(0,0,q,dim,dim)];
207              if (D == 0 ){              if (ABS(D) > 0 ){
208                   det[q]=D;
209                   D=1./D;
210                   invA[INDEX3(0,0,q,dim,dim)]=D;
211                } else {
212                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
213                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
214                 return;                 return;
215              }              }
             det[q]=D;  
             D=1./D;  
             invA[INDEX3(0,0,q,dim,dim)]=D;  
216           }           }
217           break;           break;
218    
# Line 127  void Finley_Util_InvertSmallMat(int len, Line 224  void Finley_Util_InvertSmallMat(int len,
224              A22=A[INDEX3(1,1,q,dim,dim)];              A22=A[INDEX3(1,1,q,dim,dim)];
225    
226              D = A11*A22-A12*A21;              D = A11*A22-A12*A21;
227              if (D == 0 ){              if (ABS(D) > 0 ){
228                   det[q]=D;
229                   D=1./D;
230                   invA[INDEX3(0,0,q,dim,dim)]= A22*D;
231                   invA[INDEX3(1,0,q,dim,dim)]=-A21*D;
232                   invA[INDEX3(0,1,q,dim,dim)]=-A12*D;
233                   invA[INDEX3(1,1,q,dim,dim)]= A11*D;
234                } else {
235                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
236                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
237                 return;                 return;
238              }              }
             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;  
239           }           }
240           break;           break;
241    
# Line 154  void Finley_Util_InvertSmallMat(int len, Line 252  void Finley_Util_InvertSmallMat(int len,
252              A33=A[INDEX3(2,2,q,dim,dim)];              A33=A[INDEX3(2,2,q,dim,dim)];
253    
254              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);
255              if (D == 0 ){              if (ABS(D) > 0 ){
256                   det[q]    =D;
257                   D=1./D;
258                   invA[INDEX3(0,0,q,dim,dim)]=(A22*A33-A23*A32)*D;
259                   invA[INDEX3(1,0,q,dim,dim)]=(A31*A23-A21*A33)*D;
260                   invA[INDEX3(2,0,q,dim,dim)]=(A21*A32-A31*A22)*D;
261                   invA[INDEX3(0,1,q,dim,dim)]=(A13*A32-A12*A33)*D;
262                   invA[INDEX3(1,1,q,dim,dim)]=(A11*A33-A31*A13)*D;
263                   invA[INDEX3(2,1,q,dim,dim)]=(A12*A31-A11*A32)*D;
264                   invA[INDEX3(0,2,q,dim,dim)]=(A12*A23-A13*A22)*D;
265                   invA[INDEX3(1,2,q,dim,dim)]=(A13*A21-A11*A23)*D;
266                   invA[INDEX3(2,2,q,dim,dim)]=(A11*A22-A12*A21)*D;
267                } else {
268                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
269                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
270                 return;                 return;
271              }              }
             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;  
272           }           }
273           break;           break;
274    
# Line 406  void Finley_copyDouble(int n,double* sou Line 505  void Finley_copyDouble(int n,double* sou
505    
506  /*  /*
507   * $Log$   * $Log$
508   * Revision 1.3  2004/12/15 03:48:47  jgs   * Revision 1.4  2004/12/15 07:08:34  jgs
509   * *** empty log message ***   * *** empty log message ***
510   *   *
  * Revision 1.1.1.1  2004/10/26 06:53:57  jgs  
  * initial import of project esys2  
  *  
  * Revision 1.3  2004/08/26 12:03:52  gross  
  * Some other bug in Finley_Assemble_gradient fixed.  
  *  
  * Revision 1.2  2004/07/02 04:21:13  gross  
  * Finley C code has been included  
  *  
  * Revision 1.1.1.1  2004/06/24 04:00:40  johng  
  * Initial version of eys using boost-python.  
511   *   *
512   *   *
513   */   */

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

  ViewVC Help
Powered by ViewVC 1.1.26