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

revision 114 by jgs, Mon Feb 28 07:06:33 2005 UTC revision 115 by jgs, Fri Mar 4 07:12:47 2005 UTC
# Line 99  void Finley_Util_SmallMatSetMult(int len Line 99  void Finley_Util_SmallMatSetMult(int len
99         }         }
100      }      }
101  }  }

/* 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;
}
102  /*    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 */
103  /*    the determinante is returned. */  /*    the determinante is returned. */
104
# Line 206  void Finley_Util_InvertSmallMat(int len, Line 109  void Finley_Util_InvertSmallMat(int len,
109     switch(dim) {     switch(dim) {
110        case 1:        case 1:
111           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
112              D=A[INDEX3(0,0,q,dim,dim)];              D=A[q];
113              if (ABS(D) > 0 ){              if (ABS(D) > 0 ){
114                 det[q]=D;                 det[q]=D;
115                 D=1./D;                 D=1./D;
116                 invA[INDEX3(0,0,q,dim,dim)]=D;                 invA[q]=D;
117              } else {              } else {
118                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
119                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
# Line 221  void Finley_Util_InvertSmallMat(int len, Line 124  void Finley_Util_InvertSmallMat(int len,
124
125        case 2:        case 2:
126           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
127              A11=A[INDEX3(0,0,q,dim,dim)];              A11=A[INDEX3(0,0,q,2,2)];
128              A12=A[INDEX3(0,1,q,dim,dim)];              A12=A[INDEX3(0,1,q,2,2)];
129              A21=A[INDEX3(1,0,q,dim,dim)];              A21=A[INDEX3(1,0,q,2,2)];
130              A22=A[INDEX3(1,1,q,dim,dim)];              A22=A[INDEX3(1,1,q,2,2)];
131
132              D = A11*A22-A12*A21;              D = A11*A22-A12*A21;
133              if (ABS(D) > 0 ){              if (ABS(D) > 0 ){
134                 det[q]=D;                 det[q]=D;
135                 D=1./D;                 D=1./D;
136                 invA[INDEX3(0,0,q,dim,dim)]= A22*D;                 invA[INDEX3(0,0,q,2,2)]= A22*D;
137                 invA[INDEX3(1,0,q,dim,dim)]=-A21*D;                 invA[INDEX3(1,0,q,2,2)]=-A21*D;
138                 invA[INDEX3(0,1,q,dim,dim)]=-A12*D;                 invA[INDEX3(0,1,q,2,2)]=-A12*D;
139                 invA[INDEX3(1,1,q,dim,dim)]= A11*D;                 invA[INDEX3(1,1,q,2,2)]= A11*D;
140              } else {              } else {
141                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
142                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
# Line 244  void Finley_Util_InvertSmallMat(int len, Line 147  void Finley_Util_InvertSmallMat(int len,
147
148        case 3:        case 3:
149       for (q=0;q<len;q++) {       for (q=0;q<len;q++) {
150              A11=A[INDEX3(0,0,q,dim,dim)];              A11=A[INDEX3(0,0,q,3,3)];
151              A21=A[INDEX3(1,0,q,dim,dim)];              A21=A[INDEX3(1,0,q,3,3)];
152              A31=A[INDEX3(2,0,q,dim,dim)];              A31=A[INDEX3(2,0,q,3,3)];
153              A12=A[INDEX3(0,1,q,dim,dim)];              A12=A[INDEX3(0,1,q,3,3)];
154              A22=A[INDEX3(1,1,q,dim,dim)];              A22=A[INDEX3(1,1,q,3,3)];
155              A32=A[INDEX3(2,1,q,dim,dim)];              A32=A[INDEX3(2,1,q,3,3)];
156              A13=A[INDEX3(0,2,q,dim,dim)];              A13=A[INDEX3(0,2,q,3,3)];
157              A23=A[INDEX3(1,2,q,dim,dim)];              A23=A[INDEX3(1,2,q,3,3)];
158              A33=A[INDEX3(2,2,q,dim,dim)];              A33=A[INDEX3(2,2,q,3,3)];
159
160              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);
161              if (ABS(D) > 0 ){              if (ABS(D) > 0 ){
162                 det[q]    =D;                 det[q]    =D;
163                 D=1./D;                 D=1./D;
164                 invA[INDEX3(0,0,q,dim,dim)]=(A22*A33-A23*A32)*D;                 invA[INDEX3(0,0,q,3,3)]=(A22*A33-A23*A32)*D;
165                 invA[INDEX3(1,0,q,dim,dim)]=(A31*A23-A21*A33)*D;                 invA[INDEX3(1,0,q,3,3)]=(A31*A23-A21*A33)*D;
166                 invA[INDEX3(2,0,q,dim,dim)]=(A21*A32-A31*A22)*D;                 invA[INDEX3(2,0,q,3,3)]=(A21*A32-A31*A22)*D;
167                 invA[INDEX3(0,1,q,dim,dim)]=(A13*A32-A12*A33)*D;                 invA[INDEX3(0,1,q,3,3)]=(A13*A32-A12*A33)*D;
168                 invA[INDEX3(1,1,q,dim,dim)]=(A11*A33-A31*A13)*D;                 invA[INDEX3(1,1,q,3,3)]=(A11*A33-A31*A13)*D;
169                 invA[INDEX3(2,1,q,dim,dim)]=(A12*A31-A11*A32)*D;                 invA[INDEX3(2,1,q,3,3)]=(A12*A31-A11*A32)*D;
170                 invA[INDEX3(0,2,q,dim,dim)]=(A12*A23-A13*A22)*D;                 invA[INDEX3(0,2,q,3,3)]=(A12*A23-A13*A22)*D;
171                 invA[INDEX3(1,2,q,dim,dim)]=(A13*A21-A11*A23)*D;                 invA[INDEX3(1,2,q,3,3)]=(A13*A21-A11*A23)*D;
172                 invA[INDEX3(2,2,q,dim,dim)]=(A11*A22-A12*A21)*D;                 invA[INDEX3(2,2,q,3,3)]=(A11*A22-A12*A21)*D;
173              } else {              } else {
174                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
175                 sprintf(Finley_ErrorMsg,"Non-regular matrix");                 sprintf(Finley_ErrorMsg,"Non-regular matrix");
# Line 288  void Finley_Util_DetOfSmallMat(int len,i Line 191  void Finley_Util_DetOfSmallMat(int len,i
191     switch(dim) {     switch(dim) {
192        case 1:        case 1:
193           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
194              det[q]=A[INDEX3(0,0,q,dim,dim)];              det[q]=A[q];
195           }           }
196           break;           break;
197
198        case 2:        case 2:
199           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
200              A11=A[INDEX3(0,0,q,dim,dim)];              A11=A[INDEX3(0,0,q,2,2)];
201              A12=A[INDEX3(0,1,q,dim,dim)];              A12=A[INDEX3(0,1,q,2,2)];
202              A21=A[INDEX3(1,0,q,dim,dim)];              A21=A[INDEX3(1,0,q,2,2)];
203              A22=A[INDEX3(1,1,q,dim,dim)];              A22=A[INDEX3(1,1,q,2,2)];
204
205              det[q] = A11*A22-A12*A21;              det[q] = A11*A22-A12*A21;
206           }           }
# Line 305  void Finley_Util_DetOfSmallMat(int len,i Line 208  void Finley_Util_DetOfSmallMat(int len,i
208
209        case 3:        case 3:
210       for (q=0;q<len;q++) {       for (q=0;q<len;q++) {
211              A11=A[INDEX3(0,0,q,dim,dim)];              A11=A[INDEX3(0,0,q,3,3)];
212              A21=A[INDEX3(1,0,q,dim,dim)];              A21=A[INDEX3(1,0,q,3,3)];
213              A31=A[INDEX3(2,0,q,dim,dim)];              A31=A[INDEX3(2,0,q,3,3)];
214              A12=A[INDEX3(0,1,q,dim,dim)];              A12=A[INDEX3(0,1,q,3,3)];
215              A22=A[INDEX3(1,1,q,dim,dim)];              A22=A[INDEX3(1,1,q,3,3)];
216              A32=A[INDEX3(2,1,q,dim,dim)];              A32=A[INDEX3(2,1,q,3,3)];
217              A13=A[INDEX3(0,2,q,dim,dim)];              A13=A[INDEX3(0,2,q,3,3)];
218              A23=A[INDEX3(1,2,q,dim,dim)];              A23=A[INDEX3(1,2,q,3,3)];
219              A33=A[INDEX3(2,2,q,dim,dim)];              A33=A[INDEX3(2,2,q,3,3)];
220
221              det[q]  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);              det[q]  =  A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
222           }           }
# Line 331  void  Finley_NormalVector(int len, int d Line 234  void  Finley_NormalVector(int len, int d
234
235     switch(dim) {     switch(dim) {
236        case 1:        case 1:
237           for (q=0;q<len;q++) Normal[INDEX1(q)]    =1;           for (q=0;q<len;q++) Normal[q]    =1;
238           break;           break;
239        case 2:        case 2:
240           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
241              A11=A[INDEX3(0,0,q,dim,dim1)];              A11=A[INDEX3(0,0,q,2,dim1)];
242              A21=A[INDEX3(1,0,q,dim,dim1)];              A21=A[INDEX3(1,0,q,2,dim1)];
243              length = sqrt(A11*A11+A21*A21);              length = sqrt(A11*A11+A21*A21);
244              if (! length>0) {              if (! length>0) {
245                 Finley_ErrorCode=ZERO_DIVISION_ERROR;                 Finley_ErrorCode=ZERO_DIVISION_ERROR;
# Line 344  void  Finley_NormalVector(int len, int d Line 247  void  Finley_NormalVector(int len, int d
247                 return;                 return;
248              } else {              } else {
249                 invlength=1./length;                 invlength=1./length;
250                 Normal[INDEX2(0,q,dim)]=A21*invlength;                 Normal[INDEX2(0,q,2)]=A21*invlength;
251                 Normal[INDEX2(1,q,dim)]=-A11*invlength;                 Normal[INDEX2(1,q,2)]=-A11*invlength;
252              }              }
253           }           }
254           break;           break;
255        case 3:        case 3:
256           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
257              A11=A[INDEX3(0,0,q,dim,dim1)];              A11=A[INDEX3(0,0,q,3,dim1)];
258              A21=A[INDEX3(1,0,q,dim,dim1)];              A21=A[INDEX3(1,0,q,3,dim1)];
259              A31=A[INDEX3(2,0,q,dim,dim1)];              A31=A[INDEX3(2,0,q,3,dim1)];
260              A12=A[INDEX3(0,1,q,dim,dim1)];              A12=A[INDEX3(0,1,q,3,dim1)];
261              A22=A[INDEX3(1,1,q,dim,dim1)];              A22=A[INDEX3(1,1,q,3,dim1)];
262              A32=A[INDEX3(2,1,q,dim,dim1)];              A32=A[INDEX3(2,1,q,3,dim1)];
263              CO_A13=A21*A32-A31*A22;              CO_A13=A21*A32-A31*A22;
264              CO_A23=A31*A12-A11*A32;              CO_A23=A31*A12-A11*A32;
265              CO_A33=A11*A22-A21*A12;              CO_A33=A11*A22-A21*A12;
# Line 367  void  Finley_NormalVector(int len, int d Line 270  void  Finley_NormalVector(int len, int d
270                 return;                 return;
271              } else {              } else {
272                 invlength=1./length;                 invlength=1./length;
273                 Normal[INDEX2(0,q,dim)]=CO_A13*invlength;                 Normal[INDEX2(0,q,3)]=CO_A13*invlength;
274                 Normal[INDEX2(1,q,dim)]=CO_A23*invlength;                 Normal[INDEX2(1,q,3)]=CO_A23*invlength;
275                 Normal[INDEX2(2,q,dim)]=CO_A33*invlength;                 Normal[INDEX2(2,q,3)]=CO_A33*invlength;
276             }             }
277
278        }        }
# Line 388  void  Finley_LengthOfNormalVector(int le Line 291  void  Finley_LengthOfNormalVector(int le
291
292     switch(dim) {     switch(dim) {
293        case 1:        case 1:
294           for (q=0;q<len;q++) length[INDEX1(q)]    =1;           for (q=0;q<len;q++) length[q]    =1;
295           break;           break;
296        case 2:        case 2:
297           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
298              A11=A[INDEX3(0,0,q,dim,dim1)];              A11=A[INDEX3(0,0,q,2,dim1)];
299              A21=A[INDEX3(1,0,q,dim,dim1)];              A21=A[INDEX3(1,0,q,2,dim1)];
300              length[q] = sqrt(A11*A11+A21*A21);              length[q] = sqrt(A11*A11+A21*A21);
301           }           }
302           break;           break;
303        case 3:        case 3:
304           for (q=0;q<len;q++) {           for (q=0;q<len;q++) {
305              A11=A[INDEX3(0,0,q,dim,dim1)];              A11=A[INDEX3(0,0,q,3,dim1)];
306              A21=A[INDEX3(1,0,q,dim,dim1)];              A21=A[INDEX3(1,0,q,3,dim1)];
307              A31=A[INDEX3(2,0,q,dim,dim1)];              A31=A[INDEX3(2,0,q,3,dim1)];
308              A12=A[INDEX3(0,1,q,dim,dim1)];              A12=A[INDEX3(0,1,q,3,dim1)];
309              A22=A[INDEX3(1,1,q,dim,dim1)];              A22=A[INDEX3(1,1,q,3,dim1)];
310              A32=A[INDEX3(2,1,q,dim,dim1)];              A32=A[INDEX3(2,1,q,3,dim1)];
311              CO_A13=A21*A32-A31*A22;              CO_A13=A21*A32-A31*A22;
312              CO_A23=A31*A12-A11*A32;              CO_A23=A31*A12-A11*A32;
313              CO_A33=A11*A22-A21*A12;              CO_A33=A11*A22-A21*A12;
# Line 450  void Finley_Util_sortValueAndIndex(int n Line 353  void Finley_Util_sortValueAndIndex(int n
353  /* calculates the minimum value from a dim X N integer array */  /* calculates the minimum value from a dim X N integer array */
354
355  maybelong Finley_Util_getMinInt(int dim,int N,maybelong* values) {  maybelong Finley_Util_getMinInt(int dim,int N,maybelong* values) {
356     maybelong i,j,out;     maybelong i,j,out,out_local;
357     out=MAYBELONG_MAX;     out=MAYBELONG_MAX;
358     if (values!=NULL && dim*N>0 ) {     if (values!=NULL && dim*N>0 ) {
/* OMP */
359       out=values[0];       out=values[0];
360       for (j=0;j<N;j++) {       #pragma omp parallel private(out_local)
361         for (i=0;i<dim;i++) out=MIN(out,values[INDEX2(i,j,dim)]);       {
362             out_local=out;
363             #pragma omp for private(i,j) schedule(static)
364             for (j=0;j<N;j++) {
365               for (i=0;i<dim;i++) out_local=MIN(out_local,values[INDEX2(i,j,dim)]);
366             }
367             #pragma omp critical
368             out=MIN(out_local,out);
369       }       }
370     }     }
371     return out;     return out;
# Line 465  maybelong Finley_Util_getMinInt(int dim, Line 374  maybelong Finley_Util_getMinInt(int dim,
374  /* calculates the maximum value from a dim X N integer array */  /* calculates the maximum value from a dim X N integer array */
375
376  maybelong Finley_Util_getMaxInt(int dim,int N,maybelong* values) {  maybelong Finley_Util_getMaxInt(int dim,int N,maybelong* values) {
377     maybelong i,j,out;     maybelong i,j,out,out_local;
378     out=-MAYBELONG_MAX;     out=-MAYBELONG_MAX;
379     if (values!=NULL && dim*N>0 ) {     if (values!=NULL && dim*N>0 ) {
/* OMP */
380       out=values[0];       out=values[0];
381       for (j=0;j<N;j++) {       #pragma omp parallel private(out_local)
382         for (i=0;i<dim;i++) out=MAX(out,values[INDEX2(i,j,dim)]);       {
383       }           out_local=out;
384             #pragma omp for private(i,j) schedule(static)
385             for (j=0;j<N;j++) {
386                 for (i=0;i<dim;i++) out_local=MAX(out_local,values[INDEX2(i,j,dim)]);
387             }
388             #pragma omp critical
389             out=MAX(out_local,out);
390          }
391     }     }
392     return out;     return out;
393  }  }
# Line 497  int Finley_Util_isAny(maybelong N,maybel Line 412  int Finley_Util_isAny(maybelong N,maybel
412     int out=FALSE;     int out=FALSE;
413     maybelong i;     maybelong i;
414     #pragma omp parallel for private(i) schedule(static) reduction(||:out)     #pragma omp parallel for private(i) schedule(static) reduction(||:out)
415     for (i=0;i<N;i++) out=out || (array[i]==value);     for (i=0;i<N;i++) out = out || (array[i]==value);
416     return out;     return out;
417  }  }
418  /* calculates the cummultative sum in array and returns the total sum */  /* calculates the cummultative sum in array and returns the total sum */
# Line 508  maybelong Finley_Util_cumsum(maybelong N Line 423  maybelong Finley_Util_cumsum(maybelong N
423        #pragma omp parallel private(sum,i,tmp)        #pragma omp parallel private(sum,i,tmp)
424        {        {
425          sum=0;          sum=0;
426          #pragma omp for          #pragma omp for schedule(static)
427          for (i=0;i<N;++i) {          for (i=0;i<N;++i) sum+=array[i];
tmp=sum;
sum+=array[i];
array[i]=tmp;
}
#pragma omp critical
429            #pragma omp barrier
430          #pragma omp master          #pragma omp master
431          {          {
432            out=0;            out=0;
# Line 525  maybelong Finley_Util_cumsum(maybelong N Line 436  maybelong Finley_Util_cumsum(maybelong N
436               partial_sums[i]=tmp;               partial_sums[i]=tmp;
437             }             }
438          }          }
439            #pragma omp barrier
441          #pragma omp for          #pragma omp for schedule(static)
442          for (i=0;i<N;++i) array[i]+=sum;          for (i=0;i<N;++i) {
443              tmp=sum;
444              sum+=array[i];
445              array[i]=tmp;
446            }
447        }        }
448     #else     #else
449        for (i=0;i<N;++i) {        for (i=0;i<N;++i) {

Legend:
 Removed from v.114 changed lines Added in v.115