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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 97 - (show annotations)
Tue Dec 14 05:39:33 2004 UTC (14 years, 8 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Util.c
File MIME type: text/plain
File size: 16615 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;s<len;s++) {
28 for (i=0;i<numData;i++) {
29 out[INDEX2(i,s,numData)]=in[INDEX2(i,index[s],numData)];
30 }
31 }
32 }
33
34 /**************************************************************/
35
36
37 /* gathers maybelong values out from in by index: */
38
39 /* out(1:numData,1:len)=in(1:numData,index(1:len)) */
40
41 void Finley_Util_Gather_int(int len,maybelong* index,int numData, maybelong* in, maybelong * out){
42 int s,i;
43 for (s=0;s<len;s++) {
44 for (i=0;i<numData;i++) {
45 out[INDEX2(i,s,numData)]=in[INDEX2(i,index[s],numData)];
46 }
47 }
48 }
49
50 /**************************************************************/
51
52 /* adds a vector in into out using and index. */
53
54 /* out(1:numData,index(1:len))+=in(1:numData,1:len) */
55
56 void Finley_Util_AddScatter(int len,maybelong* index,int numData,double* in,double * out){
57 int i,s;
58 for (s=0;s<len;s++) {
59 for(i=0;i<numData;i++) {
60 #pragma omp atomic
61 out[INDEX2(i,index[s],numData)]+=in[INDEX2(i,s,numData)];
62 }
63 }
64 }
65
66 /* multiplies two matrices */
67
68 /* A(1:A1,1:A2)=B(1:A1,1:B2)*C(1:B2,1:A2) */
69
70 void Finley_Util_SmallMatMult(int A1,int A2, double* A, int B2, double*B, double* C) {
71 int i,j,s;
72 for (i=0;i<A1*A2;i++) A[i]=0;
73 for (i=0;i<A1;i++) {
74 for (j=0;j<A2;j++) {
75 for (s=0;s<B2;s++) {
76 A[INDEX2(i,j,A1)]+=B[INDEX2(i,s,A1)]*C[INDEX2(s,j,B2)];
77 }
78 }
79 }
80 }
81
82 /* multiplies a two sets of matries: */
83
84 /* A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2,i) i=1,len */
85
86 void Finley_Util_SmallMatSetMult(int len,int A1,int A2, double* A, int B2, double*B, double* C) {
87 int q,i,j,s;
88 for (i=0;i<A1*A2*len;i++) A[i]=0;
89 for (q=0;q<len;q++) {
90 for (i=0;i<A1;i++) {
91 for (j=0;j<A2;j++) {
92 for (s=0;s<B2;s++) {
93 A[INDEX3(i,j,q,A1,A2)]+=B[INDEX3(i,s,q,A1,B2)]*C[INDEX3(s,j,q,B2,A2)];
94 }
95 }
96 }
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 */
197 /* the determinante is returned. */
198
199 void Finley_Util_InvertSmallMat(int len,int dim,double* A,double *invA, double* det){
200 int q;
201 double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;
202
203 switch(dim) {
204 case 1:
205 for (q=0;q<len;q++) {
206 D=A[INDEX3(0,0,q,dim,dim)];
207 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;
213 sprintf(Finley_ErrorMsg,"Non-regular matrix");
214 return;
215 }
216 }
217 break;
218
219 case 2:
220 for (q=0;q<len;q++) {
221 A11=A[INDEX3(0,0,q,dim,dim)];
222 A12=A[INDEX3(0,1,q,dim,dim)];
223 A21=A[INDEX3(1,0,q,dim,dim)];
224 A22=A[INDEX3(1,1,q,dim,dim)];
225
226 D = A11*A22-A12*A21;
227 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;
236 sprintf(Finley_ErrorMsg,"Non-regular matrix");
237 return;
238 }
239 }
240 break;
241
242 case 3:
243 for (q=0;q<len;q++) {
244 A11=A[INDEX3(0,0,q,dim,dim)];
245 A21=A[INDEX3(1,0,q,dim,dim)];
246 A31=A[INDEX3(2,0,q,dim,dim)];
247 A12=A[INDEX3(0,1,q,dim,dim)];
248 A22=A[INDEX3(1,1,q,dim,dim)];
249 A32=A[INDEX3(2,1,q,dim,dim)];
250 A13=A[INDEX3(0,2,q,dim,dim)];
251 A23=A[INDEX3(1,2,q,dim,dim)];
252 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);
255 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;
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;q<len;q++) {
288 det[q]=A[INDEX3(0,0,q,dim,dim)];
289 }
290 break;
291
292 case 2:
293 for (q=0;q<len;q++) {
294 A11=A[INDEX3(0,0,q,dim,dim)];
295 A12=A[INDEX3(0,1,q,dim,dim)];
296 A21=A[INDEX3(1,0,q,dim,dim)];
297 A22=A[INDEX3(1,1,q,dim,dim)];
298
299 det[q] = A11*A22-A12*A21;
300 }
301 break;
302
303 case 3:
304 for (q=0;q<len;q++) {
305 A11=A[INDEX3(0,0,q,dim,dim)];
306 A21=A[INDEX3(1,0,q,dim,dim)];
307 A31=A[INDEX3(2,0,q,dim,dim)];
308 A12=A[INDEX3(0,1,q,dim,dim)];
309 A22=A[INDEX3(1,1,q,dim,dim)];
310 A32=A[INDEX3(2,1,q,dim,dim)];
311 A13=A[INDEX3(0,2,q,dim,dim)];
312 A23=A[INDEX3(1,2,q,dim,dim)];
313 A33=A[INDEX3(2,2,q,dim,dim)];
314
315 det[q] = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
316 }
317 break;
318
319 }
320 return;
321 }
322 /* returns the normalized vector Normal[dim,len] orthogonal to A(:,0,q) and A(:,1,q) in the case of dim=3 */
323 /* or the vector A(:,0,q) in the case of dim=2 */
324
325 void Finley_NormalVector(int len, int dim, int dim1, double* A,double* Normal) {
326 int q;
327 double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33,length,invlength;
328
329 switch(dim) {
330 case 1:
331 for (q=0;q<len;q++) Normal[INDEX1(q)] =1;
332 break;
333 case 2:
334 for (q=0;q<len;q++) {
335 A11=A[INDEX3(0,0,q,dim,dim1)];
336 A21=A[INDEX3(1,0,q,dim,dim1)];
337 length = sqrt(A11*A11+A21*A21);
338 if (! length>0) {
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;q<len;q++) {
351 A11=A[INDEX3(0,0,q,dim,dim1)];
352 A21=A[INDEX3(1,0,q,dim,dim1)];
353 A31=A[INDEX3(2,0,q,dim,dim1)];
354 A12=A[INDEX3(0,1,q,dim,dim1)];
355 A22=A[INDEX3(1,1,q,dim,dim1)];
356 A32=A[INDEX3(2,1,q,dim,dim1)];
357 CO_A13=A21*A32-A31*A22;
358 CO_A23=A31*A12-A11*A32;
359 CO_A33=A11*A22-A21*A12;
360 length=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
361 if (! length>0) {
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<len;q++) length[INDEX1(q)] =1;
389 break;
390 case 2:
391 for (q=0;q<len;q++) {
392 A11=A[INDEX3(0,0,q,dim,dim1)];
393 A21=A[INDEX3(1,0,q,dim,dim1)];
394 length[q] = sqrt(A11*A11+A21*A21);
395 }
396 break;
397 case 3:
398 for (q=0;q<len;q++) {
399 A11=A[INDEX3(0,0,q,dim,dim1)];
400 A21=A[INDEX3(1,0,q,dim,dim1)];
401 A31=A[INDEX3(2,0,q,dim,dim1)];
402 A12=A[INDEX3(0,1,q,dim,dim1)];
403 A22=A[INDEX3(1,1,q,dim,dim1)];
404 A32=A[INDEX3(2,1,q,dim,dim1)];
405 CO_A13=A21*A32-A31*A22;
406 CO_A23=A31*A12-A11*A32;
407 CO_A33=A11*A22-A21*A12;
408 length[q]=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
409 }
410 break;
411
412 }
413 return;
414 }
415
416 /* inverts the map map of length len */
417 /* there is no range checking! */
418 /* at output Map[invMap[i]]=i for i=0:lenInvMap */
419
420 void Finley_Util_InvertMap(int lenInvMap, maybelong* invMap,int lenMap, maybelong* Map) {
421 int i;
422 for (i=0;i<lenInvMap;i++) invMap[i]=0;
423 for (i=0;i<lenMap;i++) {
424 if (Map[i]>=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;j<N;j++) {
456 for (i=0;i<dim;i++) out=MIN(out,values[INDEX2(i,j,dim)]);
457 }
458 }
459 return out;
460 }
461
462 /* calculates the maximum value from a dim X N integer array */
463
464 maybelong Finley_Util_getMaxInt(int dim,int N,maybelong* values) {
465 maybelong i,j,out;
466 out=-MAYBELONG_MAX;
467 if (values!=NULL && dim*N>0 ) {
468 /* OMP */
469 out=values[0];
470 for (j=0;j<N;j++) {
471 for (i=0;i<dim;i++) out=MAX(out,values[INDEX2(i,j,dim)]);
472 }
473 }
474 return out;
475 }
476
477 /* set the index of the positive entries in mask. The length of index is returned. */
478
479 maybelong Finley_Util_packMask(maybelong N,maybelong* mask,maybelong* index) {
480 maybelong out,k;
481 out=0;
482 /*OMP */
483 for (k=0;k<N;k++) {
484 if (mask[k]>=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<N;i++) out=out || (array[i]==value);
498 return out;
499 }
500
501 void Finley_copyDouble(int n,double* source, double* target) {
502 int i;
503 for (i=0;i<n;i++) target[i]=source[i];
504 }
505
506 /*
507 * $Log$
508 * Revision 1.2 2004/12/14 05:39:31 jgs
509 * *** empty log message ***
510 *
511 * Revision 1.1.1.1.2.1 2004/11/12 06:58:19 gross
512 * 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
513 *
514 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
515 * initial import of project esys2
516 *
517 * Revision 1.3 2004/08/26 12:03:52 gross
518 * Some other bug in Finley_Assemble_gradient fixed.
519 *
520 * Revision 1.2 2004/07/02 04:21:13 gross
521 * Finley C code has been included
522 *
523 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
524 * Initial version of eys using boost-python.
525 *
526 *
527 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26