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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 100 - (hide annotations)
Wed Dec 15 03:48:48 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: 12835 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     /* inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
101     /* the determinante is returned. */
102    
103     void Finley_Util_InvertSmallMat(int len,int dim,double* A,double *invA, double* det){
104     int q;
105     double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;
106    
107     switch(dim) {
108     case 1:
109     for (q=0;q<len;q++) {
110     D=A[INDEX3(0,0,q,dim,dim)];
111 jgs 100 if (D == 0 ){
112 jgs 82 Finley_ErrorCode=ZERO_DIVISION_ERROR;
113     sprintf(Finley_ErrorMsg,"Non-regular matrix");
114     return;
115     }
116 jgs 100 det[q]=D;
117     D=1./D;
118     invA[INDEX3(0,0,q,dim,dim)]=D;
119 jgs 82 }
120     break;
121    
122     case 2:
123     for (q=0;q<len;q++) {
124     A11=A[INDEX3(0,0,q,dim,dim)];
125     A12=A[INDEX3(0,1,q,dim,dim)];
126     A21=A[INDEX3(1,0,q,dim,dim)];
127     A22=A[INDEX3(1,1,q,dim,dim)];
128    
129     D = A11*A22-A12*A21;
130 jgs 100 if (D == 0 ){
131 jgs 82 Finley_ErrorCode=ZERO_DIVISION_ERROR;
132     sprintf(Finley_ErrorMsg,"Non-regular matrix");
133     return;
134     }
135 jgs 100 det[q]=D;
136     D=1./D;
137     invA[INDEX3(0,0,q,dim,dim)]= A22*D;
138     invA[INDEX3(1,0,q,dim,dim)]=-A21*D;
139     invA[INDEX3(0,1,q,dim,dim)]=-A12*D;
140     invA[INDEX3(1,1,q,dim,dim)]= A11*D;
141 jgs 82 }
142     break;
143    
144     case 3:
145     for (q=0;q<len;q++) {
146     A11=A[INDEX3(0,0,q,dim,dim)];
147     A21=A[INDEX3(1,0,q,dim,dim)];
148     A31=A[INDEX3(2,0,q,dim,dim)];
149     A12=A[INDEX3(0,1,q,dim,dim)];
150     A22=A[INDEX3(1,1,q,dim,dim)];
151     A32=A[INDEX3(2,1,q,dim,dim)];
152     A13=A[INDEX3(0,2,q,dim,dim)];
153     A23=A[INDEX3(1,2,q,dim,dim)];
154     A33=A[INDEX3(2,2,q,dim,dim)];
155    
156     D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
157 jgs 100 if (D == 0 ){
158 jgs 82 Finley_ErrorCode=ZERO_DIVISION_ERROR;
159     sprintf(Finley_ErrorMsg,"Non-regular matrix");
160     return;
161     }
162 jgs 100 det[q] =D;
163     D=1./D;
164     invA[INDEX3(0,0,q,dim,dim)]=(A22*A33-A23*A32)*D;
165     invA[INDEX3(1,0,q,dim,dim)]=(A31*A23-A21*A33)*D;
166     invA[INDEX3(2,0,q,dim,dim)]=(A21*A32-A31*A22)*D;
167     invA[INDEX3(0,1,q,dim,dim)]=(A13*A32-A12*A33)*D;
168     invA[INDEX3(1,1,q,dim,dim)]=(A11*A33-A31*A13)*D;
169     invA[INDEX3(2,1,q,dim,dim)]=(A12*A31-A11*A32)*D;
170     invA[INDEX3(0,2,q,dim,dim)]=(A12*A23-A13*A22)*D;
171     invA[INDEX3(1,2,q,dim,dim)]=(A13*A21-A11*A23)*D;
172     invA[INDEX3(2,2,q,dim,dim)]=(A11*A22-A12*A21)*D;
173 jgs 82 }
174     break;
175    
176     }
177     return;
178     }
179    
180     /* sets the derterminate of a set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
181    
182     void Finley_Util_DetOfSmallMat(int len,int dim,double* A, double* det){
183     int q;
184     double A11,A12,A13,A21,A22,A23,A31,A32,A33;
185    
186     switch(dim) {
187     case 1:
188     for (q=0;q<len;q++) {
189     det[q]=A[INDEX3(0,0,q,dim,dim)];
190     }
191     break;
192    
193     case 2:
194     for (q=0;q<len;q++) {
195     A11=A[INDEX3(0,0,q,dim,dim)];
196     A12=A[INDEX3(0,1,q,dim,dim)];
197     A21=A[INDEX3(1,0,q,dim,dim)];
198     A22=A[INDEX3(1,1,q,dim,dim)];
199    
200     det[q] = A11*A22-A12*A21;
201     }
202     break;
203    
204     case 3:
205     for (q=0;q<len;q++) {
206     A11=A[INDEX3(0,0,q,dim,dim)];
207     A21=A[INDEX3(1,0,q,dim,dim)];
208     A31=A[INDEX3(2,0,q,dim,dim)];
209     A12=A[INDEX3(0,1,q,dim,dim)];
210     A22=A[INDEX3(1,1,q,dim,dim)];
211     A32=A[INDEX3(2,1,q,dim,dim)];
212     A13=A[INDEX3(0,2,q,dim,dim)];
213     A23=A[INDEX3(1,2,q,dim,dim)];
214     A33=A[INDEX3(2,2,q,dim,dim)];
215    
216     det[q] = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
217     }
218     break;
219    
220     }
221     return;
222     }
223     /* returns the normalized vector Normal[dim,len] orthogonal to A(:,0,q) and A(:,1,q) in the case of dim=3 */
224     /* or the vector A(:,0,q) in the case of dim=2 */
225    
226     void Finley_NormalVector(int len, int dim, int dim1, double* A,double* Normal) {
227     int q;
228     double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33,length,invlength;
229    
230     switch(dim) {
231     case 1:
232     for (q=0;q<len;q++) Normal[INDEX1(q)] =1;
233     break;
234     case 2:
235     for (q=0;q<len;q++) {
236     A11=A[INDEX3(0,0,q,dim,dim1)];
237     A21=A[INDEX3(1,0,q,dim,dim1)];
238     length = sqrt(A11*A11+A21*A21);
239     if (! length>0) {
240     Finley_ErrorCode=ZERO_DIVISION_ERROR;
241     sprintf(Finley_ErrorMsg,"area equals zero.");
242     return;
243     } else {
244     invlength=1./length;
245     Normal[INDEX2(0,q,dim)]=A21*invlength;
246     Normal[INDEX2(1,q,dim)]=-A11*invlength;
247     }
248     }
249     break;
250     case 3:
251     for (q=0;q<len;q++) {
252     A11=A[INDEX3(0,0,q,dim,dim1)];
253     A21=A[INDEX3(1,0,q,dim,dim1)];
254     A31=A[INDEX3(2,0,q,dim,dim1)];
255     A12=A[INDEX3(0,1,q,dim,dim1)];
256     A22=A[INDEX3(1,1,q,dim,dim1)];
257     A32=A[INDEX3(2,1,q,dim,dim1)];
258     CO_A13=A21*A32-A31*A22;
259     CO_A23=A31*A12-A11*A32;
260     CO_A33=A11*A22-A21*A12;
261     length=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
262     if (! length>0) {
263     Finley_ErrorCode=ZERO_DIVISION_ERROR;
264     sprintf(Finley_ErrorMsg,"area equals zero.");
265     return;
266     } else {
267     invlength=1./length;
268     Normal[INDEX2(0,q,dim)]=CO_A13*invlength;
269     Normal[INDEX2(1,q,dim)]=CO_A23*invlength;
270     Normal[INDEX2(2,q,dim)]=CO_A33*invlength;
271     }
272    
273     }
274     break;
275    
276     }
277     return;
278     }
279    
280     /* 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 */
281     /* or the vector A(:,0,q) in the case of dim=2 */
282    
283     void Finley_LengthOfNormalVector(int len, int dim, int dim1, double* A,double* length) {
284     int q;
285     double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33;
286    
287     switch(dim) {
288     case 1:
289     for (q=0;q<len;q++) length[INDEX1(q)] =1;
290     break;
291     case 2:
292     for (q=0;q<len;q++) {
293     A11=A[INDEX3(0,0,q,dim,dim1)];
294     A21=A[INDEX3(1,0,q,dim,dim1)];
295     length[q] = sqrt(A11*A11+A21*A21);
296     }
297     break;
298     case 3:
299     for (q=0;q<len;q++) {
300     A11=A[INDEX3(0,0,q,dim,dim1)];
301     A21=A[INDEX3(1,0,q,dim,dim1)];
302     A31=A[INDEX3(2,0,q,dim,dim1)];
303     A12=A[INDEX3(0,1,q,dim,dim1)];
304     A22=A[INDEX3(1,1,q,dim,dim1)];
305     A32=A[INDEX3(2,1,q,dim,dim1)];
306     CO_A13=A21*A32-A31*A22;
307     CO_A23=A31*A12-A11*A32;
308     CO_A33=A11*A22-A21*A12;
309     length[q]=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
310     }
311     break;
312    
313     }
314     return;
315     }
316    
317     /* inverts the map map of length len */
318     /* there is no range checking! */
319     /* at output Map[invMap[i]]=i for i=0:lenInvMap */
320    
321     void Finley_Util_InvertMap(int lenInvMap, maybelong* invMap,int lenMap, maybelong* Map) {
322     int i;
323     for (i=0;i<lenInvMap;i++) invMap[i]=0;
324     for (i=0;i<lenMap;i++) {
325     if (Map[i]>=0) invMap[Map[i]]=i;
326     }
327     }
328    
329     /* orders a Finley_Util_ValueAndIndex array by value */
330     /* it is assumed that n is large */
331    
332     int Finley_Util_ValueAndIndex_compar(const void *arg1 , const void *arg2 ) {
333     Finley_Util_ValueAndIndex *e1,*e2;
334     e1=(Finley_Util_ValueAndIndex*) arg1;
335     e2=(Finley_Util_ValueAndIndex*) arg2;
336     if (e1->value < e2->value) return -1;
337     if (e1->value > e2->value) return 1;
338     return 0;
339     }
340     void Finley_Util_sortValueAndIndex(int n,Finley_Util_ValueAndIndex* array) {
341     /* OMP : needs parallelization !*/
342     qsort(array,n,sizeof(Finley_Util_ValueAndIndex),Finley_Util_ValueAndIndex_compar);
343     }
344    
345    
346     /**************************************************************/
347    
348     /* calculates the minimum value from a dim X N integer array */
349    
350     maybelong Finley_Util_getMinInt(int dim,int N,maybelong* values) {
351     maybelong i,j,out;
352     out=MAYBELONG_MAX;
353     if (values!=NULL && dim*N>0 ) {
354     /* OMP */
355     out=values[0];
356     for (j=0;j<N;j++) {
357     for (i=0;i<dim;i++) out=MIN(out,values[INDEX2(i,j,dim)]);
358     }
359     }
360     return out;
361     }
362    
363     /* calculates the maximum value from a dim X N integer array */
364    
365     maybelong Finley_Util_getMaxInt(int dim,int N,maybelong* values) {
366     maybelong i,j,out;
367     out=-MAYBELONG_MAX;
368     if (values!=NULL && dim*N>0 ) {
369     /* OMP */
370     out=values[0];
371     for (j=0;j<N;j++) {
372     for (i=0;i<dim;i++) out=MAX(out,values[INDEX2(i,j,dim)]);
373     }
374     }
375     return out;
376     }
377    
378     /* set the index of the positive entries in mask. The length of index is returned. */
379    
380     maybelong Finley_Util_packMask(maybelong N,maybelong* mask,maybelong* index) {
381     maybelong out,k;
382     out=0;
383     /*OMP */
384     for (k=0;k<N;k++) {
385     if (mask[k]>=0) {
386     index[out]=k;
387     out++;
388     }
389     }
390     return out;
391     }
392    
393     /* returns true if array contains value */
394     int Finley_Util_isAny(maybelong N,maybelong* array,maybelong value) {
395     int out=FALSE;
396     maybelong i;
397     #pragma omp parallel for private(i) schedule(static) reduction(||:out)
398     for (i=0;i<N;i++) out=out || (array[i]==value);
399     return out;
400     }
401    
402     void Finley_copyDouble(int n,double* source, double* target) {
403     int i;
404     for (i=0;i<n;i++) target[i]=source[i];
405     }
406    
407     /*
408     * $Log$
409 jgs 100 * Revision 1.3 2004/12/15 03:48:47 jgs
410 jgs 97 * *** empty log message ***
411 jgs 82 *
412 jgs 97 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
413     * initial import of project esys2
414     *
415 jgs 82 * Revision 1.3 2004/08/26 12:03:52 gross
416     * Some other bug in Finley_Assemble_gradient fixed.
417     *
418     * Revision 1.2 2004/07/02 04:21:13 gross
419     * Finley C code has been included
420     *
421     * Revision 1.1.1.1 2004/06/24 04:00:40 johng
422     * Initial version of eys using boost-python.
423     *
424     *
425     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26