/[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 82 by jgs, Tue Oct 26 06:53:54 2004 UTC revision 113 by jgs, Mon Feb 28 07:06:33 2005 UTC
# Line 15  Line 15 
15  #include "Common.h"  #include "Common.h"
16  #include "Finley.h"  #include "Finley.h"
17  #include "Util.h"  #include "Util.h"
18    #ifdef _OPENMP
19    #include <omp.h>
20    #endif
21    
22  /**************************************************************/  /**************************************************************/
23    
# Line 97  void Finley_Util_SmallMatSetMult(int len Line 100  void Finley_Util_SmallMatSetMult(int len
100      }      }
101  }  }
102    
103    /* calcultes the LU factorization for a small matrix dimxdim matrix A */
104    /* TODO: use LAPACK */
105    
106    int Finley_Util_SmallMatLU(int dim,double* A,double *LU,int* pivot){
107         double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;
108         int info=0;
109         /* LAPACK version */
110         /* memcpy(LU,A,sizeof(douple)); */
111         /* dgetf2_(&dim,&dim,A,&dim,LU,pivot,&info); */
112         switch(dim) {
113          case 1:
114                D=A[INDEX2(0,0,dim)];
115                if (ABS(D) >0. ){
116                   LU[INDEX2(0,0,dim)]=1./D;
117                } else {
118                   info=2;
119                }
120             break;
121    
122          case 2:
123                A11=A[INDEX2(0,0,dim)];
124                A12=A[INDEX2(0,1,dim)];
125                A21=A[INDEX2(1,0,dim)];
126                A22=A[INDEX2(1,1,dim)];
127    
128                D = A11*A22-A12*A21;
129                if (ABS(D) > 0 ){
130                   D=1./D;
131                   LU[INDEX2(0,0,dim)]= A22*D;
132                   LU[INDEX2(1,0,dim)]=-A21*D;
133                   LU[INDEX2(0,1,dim)]=-A12*D;
134                   LU[INDEX2(1,1,dim)]= A11*D;
135                } else {
136                   info=2;
137                }
138             break;
139    
140          case 3:
141                A11=A[INDEX2(0,0,dim)];
142                A21=A[INDEX2(1,0,dim)];
143                A31=A[INDEX2(2,0,dim)];
144                A12=A[INDEX2(0,1,dim)];
145                A22=A[INDEX2(1,1,dim)];
146                A32=A[INDEX2(2,1,dim)];
147                A13=A[INDEX2(0,2,dim)];
148                A23=A[INDEX2(1,2,dim)];
149                A33=A[INDEX2(2,2,dim)];
150    
151                D  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
152                if (ABS(D) > 0 ){
153                   D=1./D;
154                   LU[INDEX2(0,0,dim)]=(A22*A33-A23*A32)*D;
155                   LU[INDEX2(1,0,dim)]=(A31*A23-A21*A33)*D;
156                   LU[INDEX2(2,0,dim)]=(A21*A32-A31*A22)*D;
157                   LU[INDEX2(0,1,dim)]=(A13*A32-A12*A33)*D;
158                   LU[INDEX2(1,1,dim)]=(A11*A33-A31*A13)*D;
159                   LU[INDEX2(2,1,dim)]=(A12*A31-A11*A32)*D;
160                   LU[INDEX2(0,2,dim)]=(A12*A23-A13*A22)*D;
161                   LU[INDEX2(1,2,dim)]=(A13*A21-A11*A23)*D;
162                   LU[INDEX2(2,2,dim)]=(A11*A22-A12*A21)*D;
163                } else {
164                   info=2;
165                }
166             break;
167           default:
168                info=1;
169       }
170       return info;
171    }
172    
173    /* solves LUx=b where LU is a LU factorization calculated by an Finley_Util_SmallMatLU call */
174    void Finley_Util_SmallMatForwardBackwardSolve(int dim ,int nrhs,double* LU,int* pivot,double* x,double* b) {
175         int i;
176         switch(dim) {
177          case 1:
178             for (i=0;i<nrhs;i++) {
179                x[INDEX2(0,i,dim)]=LU[0]*b[INDEX2(0,i,dim)];
180             }
181             break;
182          case 2:
183             for (i=0;i<nrhs;i++) {
184                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)];
185                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)];
186             }
187             break;
188    
189          case 3:
190             for (i=0;i<nrhs;i++) {
191                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)];
192                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)];
193                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)];
194             }
195             break;
196       }
197       return;
198    }
199  /*    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 */
200  /*    the determinante is returned. */  /*    the determinante is returned. */
201    
# Line 108  void Finley_Util_InvertSmallMat(int len, Line 207  void Finley_Util_InvertSmallMat(int len,
207        case 1:        case 1:
208           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
209              D=A[INDEX3(0,0,q,dim,dim)];              D=A[INDEX3(0,0,q,dim,dim)];
210              if (D == 0 ){              if (ABS(D) > 0 ){
211                   det[q]=D;
212                   D=1./D;
213                   invA[INDEX3(0,0,q,dim,dim)]=D;
214                } else {
215                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
216                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
217                 return;                 return;
218              }              }
             det[q]=D;  
             D=1./D;  
             invA[INDEX3(0,0,q,dim,dim)]=D;  
219           }           }
220           break;           break;
221    
# Line 127  void Finley_Util_InvertSmallMat(int len, Line 227  void Finley_Util_InvertSmallMat(int len,
227              A22=A[INDEX3(1,1,q,dim,dim)];              A22=A[INDEX3(1,1,q,dim,dim)];
228    
229              D = A11*A22-A12*A21;              D = A11*A22-A12*A21;
230              if (D == 0 ){              if (ABS(D) > 0 ){
231                   det[q]=D;
232                   D=1./D;
233                   invA[INDEX3(0,0,q,dim,dim)]= A22*D;
234                   invA[INDEX3(1,0,q,dim,dim)]=-A21*D;
235                   invA[INDEX3(0,1,q,dim,dim)]=-A12*D;
236                   invA[INDEX3(1,1,q,dim,dim)]= A11*D;
237                } else {
238                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
239                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
240                 return;                 return;
241              }              }
             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;  
242           }           }
243           break;           break;
244    
# Line 154  void Finley_Util_InvertSmallMat(int len, Line 255  void Finley_Util_InvertSmallMat(int len,
255              A33=A[INDEX3(2,2,q,dim,dim)];              A33=A[INDEX3(2,2,q,dim,dim)];
256    
257              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);
258              if (D == 0 ){              if (ABS(D) > 0 ){
259                   det[q]    =D;
260                   D=1./D;
261                   invA[INDEX3(0,0,q,dim,dim)]=(A22*A33-A23*A32)*D;
262                   invA[INDEX3(1,0,q,dim,dim)]=(A31*A23-A21*A33)*D;
263                   invA[INDEX3(2,0,q,dim,dim)]=(A21*A32-A31*A22)*D;
264                   invA[INDEX3(0,1,q,dim,dim)]=(A13*A32-A12*A33)*D;
265                   invA[INDEX3(1,1,q,dim,dim)]=(A11*A33-A31*A13)*D;
266                   invA[INDEX3(2,1,q,dim,dim)]=(A12*A31-A11*A32)*D;
267                   invA[INDEX3(0,2,q,dim,dim)]=(A12*A23-A13*A22)*D;
268                   invA[INDEX3(1,2,q,dim,dim)]=(A13*A21-A11*A23)*D;
269                   invA[INDEX3(2,2,q,dim,dim)]=(A11*A22-A12*A21)*D;
270                } else {
271                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
272                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
273                 return;                 return;
274              }              }
             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;  
275           }           }
276           break;           break;
277    
# Line 398  int Finley_Util_isAny(maybelong N,maybel Line 500  int Finley_Util_isAny(maybelong N,maybel
500     for (i=0;i<N;i++) out=out || (array[i]==value);     for (i=0;i<N;i++) out=out || (array[i]==value);
501     return out;     return out;
502  }  }
503    /* calculates the cummultative sum in array and returns the total sum */
504    maybelong Finley_Util_cumsum(maybelong N,maybelong* array) {
505       maybelong out=0,tmp,i;
506       #ifdef _OPENMP
507          maybelong partial_sums[omp_get_max_threads()],sum;
508          #pragma omp parallel private(sum,i,tmp)
509          {
510            sum=0;
511            #pragma omp for
512            for (i=0;i<N;++i) {
513              tmp=sum;
514              sum+=array[i];
515              array[i]=tmp;
516            }
517            #pragma omp critical
518            partial_sums[omp_get_thread_num()]=sum;
519            #pragma omp master
520            {
521              out=0;
522              for (i=0;i<omp_get_max_threads();++i) {
523                 tmp=out;
524                 out+=partial_sums[i];
525                 partial_sums[i]=tmp;
526               }
527            }
528            sum=partial_sums[omp_get_thread_num()];
529            #pragma omp for
530            for (i=0;i<N;++i) array[i]+=sum;
531          }
532       #else
533          for (i=0;i<N;++i) {
534             tmp=out;
535             out+=array[i];
536             array[i]=tmp;
537          }
538       #endif
539       return out;
540    }
541    
542  void Finley_copyDouble(int n,double* source, double* target) {  void Finley_copyDouble(int n,double* source, double* target) {
543    int i;    int i;
544    for (i=0;i<n;i++) target[i]=source[i];    for (i=0;i<n;i++) target[i]=source[i];
545  }  }
   
 /*  
  * $Log$  
  * Revision 1.1  2004/10/26 06:53:57  jgs  
  * Initial revision  
  *  
  * 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.  
  *  
  *  
  */  

Legend:
Removed from v.82  
changed lines
  Added in v.113

  ViewVC Help
Powered by ViewVC 1.1.26