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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 123 - (hide annotations)
Fri Jul 8 04:08:13 2005 UTC (14 years, 6 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Util.c
File MIME type: text/plain
File size: 14982 byte(s)
Merge of development branch back to main trunk on 2005-07-08

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 jgs 123 void Finley_Util_Gather_double(dim_t len,index_t* index,dim_t numData,double* in, double * out){
29     dim_t s,i;
30 jgs 82 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 jgs 123 void Finley_Util_Gather_int(dim_t len,index_t* index,dim_t numData, index_t* in, index_t * out){
45     dim_t s,i;
46 jgs 82 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 jgs 123 void Finley_Util_AddScatter(dim_t len,index_t* index,dim_t numData,double* in,double * out){
60     dim_t i,s;
61 jgs 82 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 jgs 123 void Finley_Util_SmallMatMult(dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
74     dim_t i,j,s;
75 jgs 82 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 jgs 123 void Finley_Util_SmallMatSetMult(dim_t len,dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
90     dim_t q,i,j,s;
91 jgs 82 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 jgs 123 void Finley_Util_InvertSmallMat(dim_t len,dim_t dim,double* A,double *invA, double* det){
106     dim_t q;
107     register double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;
108 jgs 82
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 jgs 123 void Finley_Util_DetOfSmallMat(dim_t len,dim_t dim,double* A, double* det){
188     dim_t q;
189     register double A11,A12,A13,A21,A22,A23,A31,A32,A33;
190 jgs 82
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 jgs 123 void Finley_NormalVector(dim_t len, dim_t dim, dim_t dim1, double* A,double* Normal) {
232     dim_t q;
233     register double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33,length,invlength;
234 jgs 82
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 jgs 123 void Finley_LengthOfNormalVector(dim_t len, dim_t dim, dim_t dim1, double* A,double* length) {
289     dim_t q;
290 jgs 82 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 jgs 123 void Finley_Util_InvertMap(dim_t lenInvMap, index_t* invMap,dim_t lenMap, index_t* Map) {
327     dim_t i;
328 jgs 82 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 jgs 123 void Finley_Util_sortValueAndIndex(dim_t n,Finley_Util_ValueAndIndex* array) {
346 jgs 82 /* 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 jgs 123 index_t Finley_Util_getMinInt(dim_t dim,dim_t N,index_t* values) {
356     dim_t i,j;
357     index_t out,out_local;
358     out=INDEX_T_MAX;
359 jgs 82 if (values!=NULL && dim*N>0 ) {
360     out=values[0];
361 jgs 115 #pragma omp parallel private(out_local)
362     {
363     out_local=out;
364     #pragma omp for private(i,j) schedule(static)
365     for (j=0;j<N;j++) {
366     for (i=0;i<dim;i++) out_local=MIN(out_local,values[INDEX2(i,j,dim)]);
367     }
368     #pragma omp critical
369     out=MIN(out_local,out);
370 jgs 82 }
371     }
372     return out;
373     }
374    
375     /* calculates the maximum value from a dim X N integer array */
376    
377 jgs 123 index_t Finley_Util_getMaxInt(dim_t dim,dim_t N,index_t* values) {
378     dim_t i,j;
379     index_t out,out_local;
380     out=-INDEX_T_MAX;
381 jgs 82 if (values!=NULL && dim*N>0 ) {
382     out=values[0];
383 jgs 115 #pragma omp parallel private(out_local)
384     {
385     out_local=out;
386     #pragma omp for private(i,j) schedule(static)
387     for (j=0;j<N;j++) {
388     for (i=0;i<dim;i++) out_local=MAX(out_local,values[INDEX2(i,j,dim)]);
389     }
390     #pragma omp critical
391     out=MAX(out_local,out);
392     }
393 jgs 82 }
394     return out;
395     }
396    
397     /* set the index of the positive entries in mask. The length of index is returned. */
398    
399 jgs 123 dim_t Finley_Util_packMask(dim_t N,index_t* mask,index_t* index) {
400     dim_t out,k;
401 jgs 82 out=0;
402     /*OMP */
403     for (k=0;k<N;k++) {
404     if (mask[k]>=0) {
405     index[out]=k;
406     out++;
407     }
408     }
409     return out;
410     }
411    
412     /* returns true if array contains value */
413 jgs 123 bool_t Finley_Util_isAny(dim_t N,index_t* array,index_t value) {
414     bool_t out=FALSE;
415     dim_t i;
416 jgs 82 #pragma omp parallel for private(i) schedule(static) reduction(||:out)
417 jgs 115 for (i=0;i<N;i++) out = out || (array[i]==value);
418 jgs 82 return out;
419     }
420 jgs 113 /* calculates the cummultative sum in array and returns the total sum */
421 jgs 123 index_t Finley_Util_cumsum(dim_t N,index_t* array) {
422     index_t out=0,tmp;
423     dim_t i;
424 jgs 113 #ifdef _OPENMP
425 jgs 123 index_t partial_sums[omp_get_max_threads()],sum;
426 jgs 113 #pragma omp parallel private(sum,i,tmp)
427     {
428     sum=0;
429 jgs 115 #pragma omp for schedule(static)
430     for (i=0;i<N;++i) sum+=array[i];
431 jgs 113 partial_sums[omp_get_thread_num()]=sum;
432 jgs 115 #pragma omp barrier
433 jgs 113 #pragma omp master
434     {
435     out=0;
436     for (i=0;i<omp_get_max_threads();++i) {
437     tmp=out;
438     out+=partial_sums[i];
439     partial_sums[i]=tmp;
440     }
441     }
442 jgs 115 #pragma omp barrier
443 jgs 113 sum=partial_sums[omp_get_thread_num()];
444 jgs 115 #pragma omp for schedule(static)
445     for (i=0;i<N;++i) {
446     tmp=sum;
447     sum+=array[i];
448     array[i]=tmp;
449     }
450 jgs 113 }
451     #else
452     for (i=0;i<N;++i) {
453     tmp=out;
454     out+=array[i];
455     array[i]=tmp;
456     }
457     #endif
458     return out;
459     }
460 jgs 82
461 jgs 123 void Finley_copyDouble(dim_t n,double* source, double* target) {
462     dim_t i;
463 jgs 82 for (i=0;i<n;i++) target[i]=source[i];
464     }
465 jgs 123
466     /*
467     * $Log$
468     * Revision 1.7 2005/07/08 04:07:59 jgs
469     * Merge of development branch back to main trunk on 2005-07-08
470     *
471     * Revision 1.1.1.1.2.4 2005/06/29 02:34:57 gross
472     * some changes towards 64 integers in finley
473     *
474     * Revision 1.1.1.1.2.3 2005/03/02 23:35:06 gross
475     * reimplementation of the ILU in Finley. block size>1 still needs some testing
476     *
477     * Revision 1.1.1.1.2.2 2005/02/18 02:27:31 gross
478     * two function that will be used for a reimplementation of the ILU preconditioner
479     *
480     * Revision 1.1.1.1.2.1 2004/11/12 06:58:19 gross
481     * 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
482     *
483     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
484     * initial import of project esys2
485     *
486     * Revision 1.3 2004/08/26 12:03:52 gross
487     * Some other bug in Finley_Assemble_gradient fixed.
488     *
489     * Revision 1.2 2004/07/02 04:21:13 gross
490     * Finley C code has been included
491     *
492     * Revision 1.1.1.1 2004/06/24 04:00:40 johng
493     * Initial version of eys using boost-python.
494     *
495     *
496     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26