# Contents of /trunk/esys2/finley/src/finleyC/Util.c

Revision 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (15 years, 2 months ago) by jgs
File MIME type: text/plain
File size: 15987 byte(s)
```*** empty log message ***

```
 1 /* \$Id\$ */ 2 3 /**************************************************************/ 4 5 /* Some utility routines: */ 6 7 /**************************************************************/ 8 9 /* Copyrights by ACcESS Australia, 2003 */ 10 /* author: gross@access.edu.au */ 11 /* Version: \$Id\$ */ 12 13 /**************************************************************/ 14 15 #include "Common.h" 16 #include "Finley.h" 17 #include "Util.h" 18 19 /**************************************************************/ 20 21 /* gathers double values out from in by index: */ 22 23 /* out(1:numData,1:len)=in(1:numData,index(1:len)) */ 24 25 void Finley_Util_Gather_double(int len,maybelong* index,int numData,double* in, double * out){ 26 int s,i; 27 for (s=0;s0. ){ 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 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; 213 sprintf(Finley_ErrorMsg,"Non-regular matrix"); 214 return; 215 } 216 } 217 break; 218 219 case 2: 220 for (q=0;q 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; 236 sprintf(Finley_ErrorMsg,"Non-regular matrix"); 237 return; 238 } 239 } 240 break; 241 242 case 3: 243 for (q=0;q 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; 269 sprintf(Finley_ErrorMsg,"Non-regular matrix"); 270 return; 271 } 272 } 273 break; 274 275 } 276 return; 277 } 278 279 /* sets the derterminate of a set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */ 280 281 void Finley_Util_DetOfSmallMat(int len,int dim,double* A, double* det){ 282 int q; 283 double A11,A12,A13,A21,A22,A23,A31,A32,A33; 284 285 switch(dim) { 286 case 1: 287 for (q=0;q0) { 339 Finley_ErrorCode=ZERO_DIVISION_ERROR; 340 sprintf(Finley_ErrorMsg,"area equals zero."); 341 return; 342 } else { 343 invlength=1./length; 344 Normal[INDEX2(0,q,dim)]=A21*invlength; 345 Normal[INDEX2(1,q,dim)]=-A11*invlength; 346 } 347 } 348 break; 349 case 3: 350 for (q=0;q0) { 362 Finley_ErrorCode=ZERO_DIVISION_ERROR; 363 sprintf(Finley_ErrorMsg,"area equals zero."); 364 return; 365 } else { 366 invlength=1./length; 367 Normal[INDEX2(0,q,dim)]=CO_A13*invlength; 368 Normal[INDEX2(1,q,dim)]=CO_A23*invlength; 369 Normal[INDEX2(2,q,dim)]=CO_A33*invlength; 370 } 371 372 } 373 break; 374 375 } 376 return; 377 } 378 379 /* return the length of the vector which is orthogonal to the vectors A(:,0,q) and A(:,1,q) in the case of dim=3 */ 380 /* or the vector A(:,0,q) in the case of dim=2 */ 381 382 void Finley_LengthOfNormalVector(int len, int dim, int dim1, double* A,double* length) { 383 int q; 384 double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33; 385 386 switch(dim) { 387 case 1: 388 for (q=0;q=0) invMap[Map[i]]=i; 425 } 426 } 427 428 /* orders a Finley_Util_ValueAndIndex array by value */ 429 /* it is assumed that n is large */ 430 431 int Finley_Util_ValueAndIndex_compar(const void *arg1 , const void *arg2 ) { 432 Finley_Util_ValueAndIndex *e1,*e2; 433 e1=(Finley_Util_ValueAndIndex*) arg1; 434 e2=(Finley_Util_ValueAndIndex*) arg2; 435 if (e1->value < e2->value) return -1; 436 if (e1->value > e2->value) return 1; 437 return 0; 438 } 439 void Finley_Util_sortValueAndIndex(int n,Finley_Util_ValueAndIndex* array) { 440 /* OMP : needs parallelization !*/ 441 qsort(array,n,sizeof(Finley_Util_ValueAndIndex),Finley_Util_ValueAndIndex_compar); 442 } 443 444 445 /**************************************************************/ 446 447 /* calculates the minimum value from a dim X N integer array */ 448 449 maybelong Finley_Util_getMinInt(int dim,int N,maybelong* values) { 450 maybelong i,j,out; 451 out=MAYBELONG_MAX; 452 if (values!=NULL && dim*N>0 ) { 453 /* OMP */ 454 out=values[0]; 455 for (j=0;j0 ) { 468 /* OMP */ 469 out=values[0]; 470 for (j=0;j=0) { 485 index[out]=k; 486 out++; 487 } 488 } 489 return out; 490 } 491 492 /* returns true if array contains value */ 493 int Finley_Util_isAny(maybelong N,maybelong* array,maybelong value) { 494 int out=FALSE; 495 maybelong i; 496 #pragma omp parallel for private(i) schedule(static) reduction(||:out) 497 for (i=0;i

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision