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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 115 - (hide annotations)
Fri Mar 4 07:12:47 2005 UTC (14 years, 7 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Util.c
File MIME type: text/plain
File size: 13771 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 jgs 113 #ifdef _OPENMP
19     #include <omp.h>
20     #endif
21 jgs 82
22     /**************************************************************/
23    
24     /* gathers double values out from in by index: */
25    
26     /* out(1:numData,1:len)=in(1:numData,index(1:len)) */
27    
28     void Finley_Util_Gather_double(int len,maybelong* index,int numData,double* in, double * out){
29     int s,i;
30     for (s=0;s<len;s++) {
31     for (i=0;i<numData;i++) {
32     out[INDEX2(i,s,numData)]=in[INDEX2(i,index[s],numData)];
33     }
34     }
35     }
36    
37     /**************************************************************/
38    
39    
40     /* gathers maybelong values out from in by index: */
41    
42     /* out(1:numData,1:len)=in(1:numData,index(1:len)) */
43    
44     void Finley_Util_Gather_int(int len,maybelong* index,int numData, maybelong* in, maybelong * out){
45     int s,i;
46     for (s=0;s<len;s++) {
47     for (i=0;i<numData;i++) {
48     out[INDEX2(i,s,numData)]=in[INDEX2(i,index[s],numData)];
49     }
50     }
51     }
52    
53     /**************************************************************/
54    
55     /* adds a vector in into out using and index. */
56    
57     /* out(1:numData,index(1:len))+=in(1:numData,1:len) */
58    
59     void Finley_Util_AddScatter(int len,maybelong* index,int numData,double* in,double * out){
60     int i,s;
61     for (s=0;s<len;s++) {
62     for(i=0;i<numData;i++) {
63     #pragma omp atomic
64     out[INDEX2(i,index[s],numData)]+=in[INDEX2(i,s,numData)];
65     }
66     }
67     }
68    
69     /* multiplies two matrices */
70    
71     /* A(1:A1,1:A2)=B(1:A1,1:B2)*C(1:B2,1:A2) */
72    
73     void Finley_Util_SmallMatMult(int A1,int A2, double* A, int B2, double*B, double* C) {
74     int i,j,s;
75     for (i=0;i<A1*A2;i++) A[i]=0;
76     for (i=0;i<A1;i++) {
77     for (j=0;j<A2;j++) {
78     for (s=0;s<B2;s++) {
79     A[INDEX2(i,j,A1)]+=B[INDEX2(i,s,A1)]*C[INDEX2(s,j,B2)];
80     }
81     }
82     }
83     }
84    
85     /* multiplies a two sets of matries: */
86    
87     /* A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2,i) i=1,len */
88    
89     void Finley_Util_SmallMatSetMult(int len,int A1,int A2, double* A, int B2, double*B, double* C) {
90     int q,i,j,s;
91     for (i=0;i<A1*A2*len;i++) A[i]=0;
92     for (q=0;q<len;q++) {
93     for (i=0;i<A1;i++) {
94     for (j=0;j<A2;j++) {
95     for (s=0;s<B2;s++) {
96     A[INDEX3(i,j,q,A1,A2)]+=B[INDEX3(i,s,q,A1,B2)]*C[INDEX3(s,j,q,B2,A2)];
97     }
98     }
99     }
100     }
101     }
102     /* inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
103     /* the determinante is returned. */
104    
105     void Finley_Util_InvertSmallMat(int len,int dim,double* A,double *invA, double* det){
106     int q;
107     double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;
108    
109     switch(dim) {
110     case 1:
111     for (q=0;q<len;q++) {
112 jgs 115 D=A[q];
113 jgs 102 if (ABS(D) > 0 ){
114     det[q]=D;
115     D=1./D;
116 jgs 115 invA[q]=D;
117 jgs 102 } else {
118 jgs 82 Finley_ErrorCode=ZERO_DIVISION_ERROR;
119     sprintf(Finley_ErrorMsg,"Non-regular matrix");
120     return;
121     }
122     }
123     break;
124    
125     case 2:
126     for (q=0;q<len;q++) {
127 jgs 115 A11=A[INDEX3(0,0,q,2,2)];
128     A12=A[INDEX3(0,1,q,2,2)];
129     A21=A[INDEX3(1,0,q,2,2)];
130     A22=A[INDEX3(1,1,q,2,2)];
131 jgs 82
132     D = A11*A22-A12*A21;
133 jgs 102 if (ABS(D) > 0 ){
134     det[q]=D;
135     D=1./D;
136 jgs 115 invA[INDEX3(0,0,q,2,2)]= A22*D;
137     invA[INDEX3(1,0,q,2,2)]=-A21*D;
138     invA[INDEX3(0,1,q,2,2)]=-A12*D;
139     invA[INDEX3(1,1,q,2,2)]= A11*D;
140 jgs 102 } else {
141 jgs 82 Finley_ErrorCode=ZERO_DIVISION_ERROR;
142     sprintf(Finley_ErrorMsg,"Non-regular matrix");
143     return;
144     }
145     }
146     break;
147    
148     case 3:
149     for (q=0;q<len;q++) {
150 jgs 115 A11=A[INDEX3(0,0,q,3,3)];
151     A21=A[INDEX3(1,0,q,3,3)];
152     A31=A[INDEX3(2,0,q,3,3)];
153     A12=A[INDEX3(0,1,q,3,3)];
154     A22=A[INDEX3(1,1,q,3,3)];
155     A32=A[INDEX3(2,1,q,3,3)];
156     A13=A[INDEX3(0,2,q,3,3)];
157     A23=A[INDEX3(1,2,q,3,3)];
158     A33=A[INDEX3(2,2,q,3,3)];
159 jgs 82
160     D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
161 jgs 102 if (ABS(D) > 0 ){
162     det[q] =D;
163     D=1./D;
164 jgs 115 invA[INDEX3(0,0,q,3,3)]=(A22*A33-A23*A32)*D;
165     invA[INDEX3(1,0,q,3,3)]=(A31*A23-A21*A33)*D;
166     invA[INDEX3(2,0,q,3,3)]=(A21*A32-A31*A22)*D;
167     invA[INDEX3(0,1,q,3,3)]=(A13*A32-A12*A33)*D;
168     invA[INDEX3(1,1,q,3,3)]=(A11*A33-A31*A13)*D;
169     invA[INDEX3(2,1,q,3,3)]=(A12*A31-A11*A32)*D;
170     invA[INDEX3(0,2,q,3,3)]=(A12*A23-A13*A22)*D;
171     invA[INDEX3(1,2,q,3,3)]=(A13*A21-A11*A23)*D;
172     invA[INDEX3(2,2,q,3,3)]=(A11*A22-A12*A21)*D;
173 jgs 102 } else {
174 jgs 82 Finley_ErrorCode=ZERO_DIVISION_ERROR;
175     sprintf(Finley_ErrorMsg,"Non-regular matrix");
176     return;
177     }
178     }
179     break;
180    
181     }
182     return;
183     }
184    
185     /* sets the derterminate of a set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
186    
187     void Finley_Util_DetOfSmallMat(int len,int dim,double* A, double* det){
188     int q;
189     double A11,A12,A13,A21,A22,A23,A31,A32,A33;
190    
191     switch(dim) {
192     case 1:
193     for (q=0;q<len;q++) {
194 jgs 115 det[q]=A[q];
195 jgs 82 }
196     break;
197    
198     case 2:
199     for (q=0;q<len;q++) {
200 jgs 115 A11=A[INDEX3(0,0,q,2,2)];
201     A12=A[INDEX3(0,1,q,2,2)];
202     A21=A[INDEX3(1,0,q,2,2)];
203     A22=A[INDEX3(1,1,q,2,2)];
204 jgs 82
205     det[q] = A11*A22-A12*A21;
206     }
207     break;
208    
209     case 3:
210     for (q=0;q<len;q++) {
211 jgs 115 A11=A[INDEX3(0,0,q,3,3)];
212     A21=A[INDEX3(1,0,q,3,3)];
213     A31=A[INDEX3(2,0,q,3,3)];
214     A12=A[INDEX3(0,1,q,3,3)];
215     A22=A[INDEX3(1,1,q,3,3)];
216     A32=A[INDEX3(2,1,q,3,3)];
217     A13=A[INDEX3(0,2,q,3,3)];
218     A23=A[INDEX3(1,2,q,3,3)];
219     A33=A[INDEX3(2,2,q,3,3)];
220 jgs 82
221     det[q] = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
222     }
223     break;
224    
225     }
226     return;
227     }
228     /* returns the normalized vector Normal[dim,len] orthogonal to A(:,0,q) and A(:,1,q) in the case of dim=3 */
229     /* or the vector A(:,0,q) in the case of dim=2 */
230    
231     void Finley_NormalVector(int len, int dim, int dim1, double* A,double* Normal) {
232     int q;
233     double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33,length,invlength;
234    
235     switch(dim) {
236     case 1:
237 jgs 115 for (q=0;q<len;q++) Normal[q] =1;
238 jgs 82 break;
239     case 2:
240     for (q=0;q<len;q++) {
241 jgs 115 A11=A[INDEX3(0,0,q,2,dim1)];
242     A21=A[INDEX3(1,0,q,2,dim1)];
243 jgs 82 length = sqrt(A11*A11+A21*A21);
244     if (! length>0) {
245     Finley_ErrorCode=ZERO_DIVISION_ERROR;
246     sprintf(Finley_ErrorMsg,"area equals zero.");
247     return;
248     } else {
249     invlength=1./length;
250 jgs 115 Normal[INDEX2(0,q,2)]=A21*invlength;
251     Normal[INDEX2(1,q,2)]=-A11*invlength;
252 jgs 82 }
253     }
254     break;
255     case 3:
256     for (q=0;q<len;q++) {
257 jgs 115 A11=A[INDEX3(0,0,q,3,dim1)];
258     A21=A[INDEX3(1,0,q,3,dim1)];
259     A31=A[INDEX3(2,0,q,3,dim1)];
260     A12=A[INDEX3(0,1,q,3,dim1)];
261     A22=A[INDEX3(1,1,q,3,dim1)];
262     A32=A[INDEX3(2,1,q,3,dim1)];
263 jgs 82 CO_A13=A21*A32-A31*A22;
264     CO_A23=A31*A12-A11*A32;
265     CO_A33=A11*A22-A21*A12;
266     length=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
267     if (! length>0) {
268     Finley_ErrorCode=ZERO_DIVISION_ERROR;
269     sprintf(Finley_ErrorMsg,"area equals zero.");
270     return;
271     } else {
272     invlength=1./length;
273 jgs 115 Normal[INDEX2(0,q,3)]=CO_A13*invlength;
274     Normal[INDEX2(1,q,3)]=CO_A23*invlength;
275     Normal[INDEX2(2,q,3)]=CO_A33*invlength;
276 jgs 82 }
277    
278     }
279     break;
280    
281     }
282     return;
283     }
284    
285     /* 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 */
286     /* or the vector A(:,0,q) in the case of dim=2 */
287    
288     void Finley_LengthOfNormalVector(int len, int dim, int dim1, double* A,double* length) {
289     int q;
290     double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33;
291    
292     switch(dim) {
293     case 1:
294 jgs 115 for (q=0;q<len;q++) length[q] =1;
295 jgs 82 break;
296     case 2:
297     for (q=0;q<len;q++) {
298 jgs 115 A11=A[INDEX3(0,0,q,2,dim1)];
299     A21=A[INDEX3(1,0,q,2,dim1)];
300 jgs 82 length[q] = sqrt(A11*A11+A21*A21);
301     }
302     break;
303     case 3:
304     for (q=0;q<len;q++) {
305 jgs 115 A11=A[INDEX3(0,0,q,3,dim1)];
306     A21=A[INDEX3(1,0,q,3,dim1)];
307     A31=A[INDEX3(2,0,q,3,dim1)];
308     A12=A[INDEX3(0,1,q,3,dim1)];
309     A22=A[INDEX3(1,1,q,3,dim1)];
310     A32=A[INDEX3(2,1,q,3,dim1)];
311 jgs 82 CO_A13=A21*A32-A31*A22;
312     CO_A23=A31*A12-A11*A32;
313     CO_A33=A11*A22-A21*A12;
314     length[q]=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
315     }
316     break;
317    
318     }
319     return;
320     }
321    
322     /* inverts the map map of length len */
323     /* there is no range checking! */
324     /* at output Map[invMap[i]]=i for i=0:lenInvMap */
325    
326     void Finley_Util_InvertMap(int lenInvMap, maybelong* invMap,int lenMap, maybelong* Map) {
327     int i;
328     for (i=0;i<lenInvMap;i++) invMap[i]=0;
329     for (i=0;i<lenMap;i++) {
330     if (Map[i]>=0) invMap[Map[i]]=i;
331     }
332     }
333    
334     /* orders a Finley_Util_ValueAndIndex array by value */
335     /* it is assumed that n is large */
336    
337     int Finley_Util_ValueAndIndex_compar(const void *arg1 , const void *arg2 ) {
338     Finley_Util_ValueAndIndex *e1,*e2;
339     e1=(Finley_Util_ValueAndIndex*) arg1;
340     e2=(Finley_Util_ValueAndIndex*) arg2;
341     if (e1->value < e2->value) return -1;
342     if (e1->value > e2->value) return 1;
343     return 0;
344     }
345     void Finley_Util_sortValueAndIndex(int n,Finley_Util_ValueAndIndex* array) {
346     /* OMP : needs parallelization !*/
347     qsort(array,n,sizeof(Finley_Util_ValueAndIndex),Finley_Util_ValueAndIndex_compar);
348     }
349    
350    
351     /**************************************************************/
352    
353     /* calculates the minimum value from a dim X N integer array */
354    
355     maybelong Finley_Util_getMinInt(int dim,int N,maybelong* values) {
356 jgs 115 maybelong i,j,out,out_local;
357 jgs 82 out=MAYBELONG_MAX;
358     if (values!=NULL && dim*N>0 ) {
359     out=values[0];
360 jgs 115 #pragma omp parallel private(out_local)
361     {
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 jgs 82 }
370     }
371     return out;
372     }
373    
374     /* calculates the maximum value from a dim X N integer array */
375    
376     maybelong Finley_Util_getMaxInt(int dim,int N,maybelong* values) {
377 jgs 115 maybelong i,j,out,out_local;
378 jgs 82 out=-MAYBELONG_MAX;
379     if (values!=NULL && dim*N>0 ) {
380     out=values[0];
381 jgs 115 #pragma omp parallel private(out_local)
382     {
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 jgs 82 }
392     return out;
393     }
394    
395     /* set the index of the positive entries in mask. The length of index is returned. */
396    
397     maybelong Finley_Util_packMask(maybelong N,maybelong* mask,maybelong* index) {
398     maybelong out,k;
399     out=0;
400     /*OMP */
401     for (k=0;k<N;k++) {
402     if (mask[k]>=0) {
403     index[out]=k;
404     out++;
405     }
406     }
407     return out;
408     }
409    
410     /* returns true if array contains value */
411     int Finley_Util_isAny(maybelong N,maybelong* array,maybelong value) {
412     int out=FALSE;
413     maybelong i;
414     #pragma omp parallel for private(i) schedule(static) reduction(||:out)
415 jgs 115 for (i=0;i<N;i++) out = out || (array[i]==value);
416 jgs 82 return out;
417     }
418 jgs 113 /* calculates the cummultative sum in array and returns the total sum */
419     maybelong Finley_Util_cumsum(maybelong N,maybelong* array) {
420     maybelong out=0,tmp,i;
421     #ifdef _OPENMP
422     maybelong partial_sums[omp_get_max_threads()],sum;
423     #pragma omp parallel private(sum,i,tmp)
424     {
425     sum=0;
426 jgs 115 #pragma omp for schedule(static)
427     for (i=0;i<N;++i) sum+=array[i];
428 jgs 113 partial_sums[omp_get_thread_num()]=sum;
429 jgs 115 #pragma omp barrier
430 jgs 113 #pragma omp master
431     {
432     out=0;
433     for (i=0;i<omp_get_max_threads();++i) {
434     tmp=out;
435     out+=partial_sums[i];
436     partial_sums[i]=tmp;
437     }
438     }
439 jgs 115 #pragma omp barrier
440 jgs 113 sum=partial_sums[omp_get_thread_num()];
441 jgs 115 #pragma omp for schedule(static)
442     for (i=0;i<N;++i) {
443     tmp=sum;
444     sum+=array[i];
445     array[i]=tmp;
446     }
447 jgs 113 }
448     #else
449     for (i=0;i<N;++i) {
450     tmp=out;
451     out+=array[i];
452     array[i]=tmp;
453     }
454     #endif
455     return out;
456     }
457 jgs 82
458     void Finley_copyDouble(int n,double* source, double* target) {
459     int i;
460     for (i=0;i<n;i++) target[i]=source[i];
461     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26