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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 97 - (hide annotations)
Tue Dec 14 05:39:33 2004 UTC (14 years, 10 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 jgs 82 /* $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 jgs 97 /* 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 jgs 82 /* 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 jgs 97 if (ABS(D) > 0 ){
208     det[q]=D;
209     D=1./D;
210     invA[INDEX3(0,0,q,dim,dim)]=D;
211     } else {
212 jgs 82 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 jgs 97 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 jgs 82 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 jgs 97 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 jgs 82 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 jgs 97 * Revision 1.2 2004/12/14 05:39:31 jgs
509     * *** empty log message ***
510 jgs 82 *
511 jgs 97 * 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 jgs 82 * 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