/[escript]/branches/domexper/dudley/src/Util.c
ViewVC logotype

Annotation of /branches/domexper/dudley/src/Util.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3114 - (hide annotations)
Fri Aug 27 05:26:25 2010 UTC (9 years, 5 months ago) by jfenwick
File MIME type: text/plain
File size: 19272 byte(s)
It doesn't pass all tests but this is major progress

1 jgs 82
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13 jgs 82
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17 ksteube 1312 /* Some utility routines: */
18 jgs 82
19     /**************************************************************/
20    
21     #include "Util.h"
22 woo409 757
23 jgs 113 #ifdef _OPENMP
24     #include <omp.h>
25     #endif
26 jgs 82
27     /**************************************************************/
28    
29 jgs 147 /* returns true if any of the values in the short array values is not equalt to Zero */
30    
31 jfenwick 3086 bool_t Dudley_Util_anyNonZeroDouble(dim_t N, double* values) {
32 jgs 147 dim_t q;
33     for (q=0;q<N;++q) if (ABS(values[q])>0) return TRUE;
34     return FALSE;
35     }
36     /**************************************************************/
37    
38 jgs 82 /* gathers double values out from in by index: */
39    
40     /* out(1:numData,1:len)=in(1:numData,index(1:len)) */
41    
42 jfenwick 3086 void Dudley_Util_Gather_double(dim_t len,index_t* index,dim_t numData,double* in, double * out){
43 jgs 123 dim_t s,i;
44 jgs 82 for (s=0;s<len;s++) {
45     for (i=0;i<numData;i++) {
46     out[INDEX2(i,s,numData)]=in[INDEX2(i,index[s],numData)];
47     }
48     }
49     }
50    
51     /**************************************************************/
52    
53    
54     /* gathers maybelong values out from in by index: */
55    
56     /* out(1:numData,1:len)=in(1:numData,index(1:len)) */
57    
58 jfenwick 3086 void Dudley_Util_Gather_int(dim_t len,index_t* index,dim_t numData, index_t* in, index_t * out){
59 jgs 123 dim_t s,i;
60 jgs 82 for (s=0;s<len;s++) {
61     for (i=0;i<numData;i++) {
62     out[INDEX2(i,s,numData)]=in[INDEX2(i,index[s],numData)];
63     }
64     }
65     }
66    
67     /**************************************************************/
68    
69     /* adds a vector in into out using and index. */
70    
71 gross 798 /* out(1:numData,index[p])+=in(1:numData,p) where p = {k=1...len , index[k]<upperBound}*/
72 jgs 82
73 gross 798
74 jfenwick 3086 void Dudley_Util_AddScatter(dim_t len,index_t* index,dim_t numData,double* in,double * out, index_t upperBound){
75 jgs 123 dim_t i,s;
76 jgs 82 for (s=0;s<len;s++) {
77     for(i=0;i<numData;i++) {
78 ksteube 799 if( index[s]<upperBound ) {
79 bcumming 751 out[INDEX2(i,index[s],numData)]+=in[INDEX2(i,s,numData)];
80 ksteube 799 }
81 bcumming 751 }
82     }
83 gross 798 }
84 bcumming 751
85 jgs 82 /* multiplies two matrices */
86    
87     /* A(1:A1,1:A2)=B(1:A1,1:B2)*C(1:B2,1:A2) */
88    
89 jfenwick 3086 void Dudley_Util_SmallMatMult(dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
90 jgs 123 dim_t i,j,s;
91 gross 2748 register double rtmp;
92 jgs 82 for (i=0;i<A1;i++) {
93     for (j=0;j<A2;j++) {
94 gross 2748 rtmp=0;
95     for (s=0;s<B2;s++) rtmp+=B[INDEX2(i,s,A1)]*C[INDEX2(s,j,B2)];
96     A[INDEX2(i,j,A1)]=rtmp;
97 jgs 82 }
98     }
99     }
100    
101     /* multiplies a two sets of matries: */
102    
103     /* A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2,i) i=1,len */
104    
105 jfenwick 3086 void Dudley_Util_SmallMatSetMult(dim_t len,dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
106 jgs 123 dim_t q,i,j,s;
107 gross 2748 register double rtmp;
108 jgs 82 for (q=0;q<len;q++) {
109     for (i=0;i<A1;i++) {
110     for (j=0;j<A2;j++) {
111 gross 2748 rtmp=0;
112     for (s=0;s<B2;s++) rtmp+=B[INDEX3(i,s,q,A1,B2)]*C[INDEX3(s,j,q,B2,A2)];
113     A[INDEX3(i,j,q, A1,A2)]=rtmp;
114 jgs 82 }
115     }
116     }
117     }
118 gross 2748 /* multiplies a set of matries with a single matrix: */
119    
120     /* A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2) i=1,len */
121    
122 jfenwick 3086 void Dudley_Util_SmallMatSetMult1(dim_t len,dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
123 gross 2748 dim_t q,i,j,s;
124     register double rtmp;
125     for (q=0;q<len;q++) {
126     for (i=0;i<A1;i++) {
127     for (j=0;j<A2;j++) {
128     rtmp=0;
129     for (s=0;s<B2;s++) rtmp+=B[INDEX3(i,s,q, A1,B2)]*C[INDEX2(s,j,B2)];
130     A[INDEX3(i,j,q,A1,A2)]=rtmp;
131     }
132     }
133     }
134     }
135 jgs 82 /* inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
136     /* the determinante is returned. */
137    
138 jfenwick 3086 void Dudley_Util_InvertSmallMat(dim_t len,dim_t dim,double* A,double *invA, double* det){
139 jgs 123 dim_t q;
140     register double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;
141 jgs 82
142     switch(dim) {
143     case 1:
144     for (q=0;q<len;q++) {
145 jgs 115 D=A[q];
146 jgs 102 if (ABS(D) > 0 ){
147     det[q]=D;
148     D=1./D;
149 jgs 115 invA[q]=D;
150 jgs 102 } else {
151 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
152 jgs 82 return;
153     }
154     }
155     break;
156    
157     case 2:
158     for (q=0;q<len;q++) {
159 jgs 115 A11=A[INDEX3(0,0,q,2,2)];
160     A12=A[INDEX3(0,1,q,2,2)];
161     A21=A[INDEX3(1,0,q,2,2)];
162     A22=A[INDEX3(1,1,q,2,2)];
163 jgs 82
164     D = A11*A22-A12*A21;
165 jgs 102 if (ABS(D) > 0 ){
166     det[q]=D;
167     D=1./D;
168 jgs 115 invA[INDEX3(0,0,q,2,2)]= A22*D;
169     invA[INDEX3(1,0,q,2,2)]=-A21*D;
170     invA[INDEX3(0,1,q,2,2)]=-A12*D;
171     invA[INDEX3(1,1,q,2,2)]= A11*D;
172 jgs 102 } else {
173 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
174 jgs 82 return;
175     }
176     }
177     break;
178    
179     case 3:
180     for (q=0;q<len;q++) {
181 jgs 115 A11=A[INDEX3(0,0,q,3,3)];
182     A21=A[INDEX3(1,0,q,3,3)];
183     A31=A[INDEX3(2,0,q,3,3)];
184     A12=A[INDEX3(0,1,q,3,3)];
185     A22=A[INDEX3(1,1,q,3,3)];
186     A32=A[INDEX3(2,1,q,3,3)];
187     A13=A[INDEX3(0,2,q,3,3)];
188     A23=A[INDEX3(1,2,q,3,3)];
189     A33=A[INDEX3(2,2,q,3,3)];
190 jgs 82
191     D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
192 jgs 102 if (ABS(D) > 0 ){
193     det[q] =D;
194     D=1./D;
195 jgs 115 invA[INDEX3(0,0,q,3,3)]=(A22*A33-A23*A32)*D;
196     invA[INDEX3(1,0,q,3,3)]=(A31*A23-A21*A33)*D;
197     invA[INDEX3(2,0,q,3,3)]=(A21*A32-A31*A22)*D;
198     invA[INDEX3(0,1,q,3,3)]=(A13*A32-A12*A33)*D;
199     invA[INDEX3(1,1,q,3,3)]=(A11*A33-A31*A13)*D;
200     invA[INDEX3(2,1,q,3,3)]=(A12*A31-A11*A32)*D;
201     invA[INDEX3(0,2,q,3,3)]=(A12*A23-A13*A22)*D;
202     invA[INDEX3(1,2,q,3,3)]=(A13*A21-A11*A23)*D;
203     invA[INDEX3(2,2,q,3,3)]=(A11*A22-A12*A21)*D;
204 jgs 102 } else {
205 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
206 jgs 82 return;
207     }
208     }
209     break;
210    
211     }
212     return;
213     }
214    
215     /* sets the derterminate of a set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
216    
217 jfenwick 3086 void Dudley_Util_DetOfSmallMat(dim_t len,dim_t dim,double* A, double* det){
218 jgs 123 dim_t q;
219     register double A11,A12,A13,A21,A22,A23,A31,A32,A33;
220 jgs 82
221     switch(dim) {
222     case 1:
223     for (q=0;q<len;q++) {
224 jgs 115 det[q]=A[q];
225 jgs 82 }
226     break;
227    
228     case 2:
229     for (q=0;q<len;q++) {
230 jgs 115 A11=A[INDEX3(0,0,q,2,2)];
231     A12=A[INDEX3(0,1,q,2,2)];
232     A21=A[INDEX3(1,0,q,2,2)];
233     A22=A[INDEX3(1,1,q,2,2)];
234 jgs 82
235     det[q] = A11*A22-A12*A21;
236     }
237     break;
238    
239     case 3:
240     for (q=0;q<len;q++) {
241 jgs 115 A11=A[INDEX3(0,0,q,3,3)];
242     A21=A[INDEX3(1,0,q,3,3)];
243     A31=A[INDEX3(2,0,q,3,3)];
244     A12=A[INDEX3(0,1,q,3,3)];
245     A22=A[INDEX3(1,1,q,3,3)];
246     A32=A[INDEX3(2,1,q,3,3)];
247     A13=A[INDEX3(0,2,q,3,3)];
248     A23=A[INDEX3(1,2,q,3,3)];
249     A33=A[INDEX3(2,2,q,3,3)];
250 jgs 82
251     det[q] = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
252     }
253     break;
254    
255     }
256     return;
257     }
258     /* returns the normalized vector Normal[dim,len] orthogonal to A(:,0,q) and A(:,1,q) in the case of dim=3 */
259     /* or the vector A(:,0,q) in the case of dim=2 */
260    
261 jfenwick 3086 void Dudley_NormalVector(dim_t len, dim_t dim, dim_t dim1, double* A,double* Normal) {
262 jgs 123 dim_t q;
263     register double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33,length,invlength;
264 jgs 82
265     switch(dim) {
266     case 1:
267 jgs 115 for (q=0;q<len;q++) Normal[q] =1;
268 jgs 82 break;
269     case 2:
270     for (q=0;q<len;q++) {
271 jgs 115 A11=A[INDEX3(0,0,q,2,dim1)];
272     A21=A[INDEX3(1,0,q,2,dim1)];
273 jgs 82 length = sqrt(A11*A11+A21*A21);
274     if (! length>0) {
275 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,"__FILE__: area equals zero.");
276 jgs 82 return;
277     } else {
278     invlength=1./length;
279 jgs 115 Normal[INDEX2(0,q,2)]=A21*invlength;
280     Normal[INDEX2(1,q,2)]=-A11*invlength;
281 jgs 82 }
282     }
283     break;
284     case 3:
285     for (q=0;q<len;q++) {
286 jgs 115 A11=A[INDEX3(0,0,q,3,dim1)];
287     A21=A[INDEX3(1,0,q,3,dim1)];
288     A31=A[INDEX3(2,0,q,3,dim1)];
289     A12=A[INDEX3(0,1,q,3,dim1)];
290     A22=A[INDEX3(1,1,q,3,dim1)];
291     A32=A[INDEX3(2,1,q,3,dim1)];
292 jgs 82 CO_A13=A21*A32-A31*A22;
293     CO_A23=A31*A12-A11*A32;
294     CO_A33=A11*A22-A21*A12;
295     length=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
296     if (! length>0) {
297 jfenwick 3086 Dudley_setError(ZERO_DIVISION_ERROR,"__FILE__: area equals zero.");
298 jgs 82 return;
299     } else {
300     invlength=1./length;
301 jgs 115 Normal[INDEX2(0,q,3)]=CO_A13*invlength;
302     Normal[INDEX2(1,q,3)]=CO_A23*invlength;
303     Normal[INDEX2(2,q,3)]=CO_A33*invlength;
304 jgs 82 }
305    
306     }
307     break;
308    
309     }
310     return;
311     }
312    
313     /* 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 */
314     /* or the vector A(:,0,q) in the case of dim=2 */
315    
316 jfenwick 3086 void Dudley_LengthOfNormalVector(dim_t len, dim_t dim, dim_t dim1, double* A,double* length) {
317 jgs 123 dim_t q;
318 jgs 82 double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33;
319    
320     switch(dim) {
321     case 1:
322 jgs 115 for (q=0;q<len;q++) length[q] =1;
323 jgs 82 break;
324     case 2:
325     for (q=0;q<len;q++) {
326 jgs 115 A11=A[INDEX3(0,0,q,2,dim1)];
327     A21=A[INDEX3(1,0,q,2,dim1)];
328 jgs 82 length[q] = sqrt(A11*A11+A21*A21);
329     }
330     break;
331     case 3:
332     for (q=0;q<len;q++) {
333 jgs 115 A11=A[INDEX3(0,0,q,3,dim1)];
334     A21=A[INDEX3(1,0,q,3,dim1)];
335     A31=A[INDEX3(2,0,q,3,dim1)];
336     A12=A[INDEX3(0,1,q,3,dim1)];
337     A22=A[INDEX3(1,1,q,3,dim1)];
338     A32=A[INDEX3(2,1,q,3,dim1)];
339 jgs 82 CO_A13=A21*A32-A31*A22;
340     CO_A23=A31*A12-A11*A32;
341     CO_A33=A11*A22-A21*A12;
342     length[q]=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
343     }
344     break;
345    
346     }
347     return;
348     }
349    
350     /* inverts the map map of length len */
351     /* there is no range checking! */
352     /* at output Map[invMap[i]]=i for i=0:lenInvMap */
353    
354 jfenwick 3086 void Dudley_Util_InvertMap(dim_t lenInvMap, index_t* invMap,dim_t lenMap, index_t* Map) {
355 jgs 123 dim_t i;
356 jgs 82 for (i=0;i<lenInvMap;i++) invMap[i]=0;
357     for (i=0;i<lenMap;i++) {
358     if (Map[i]>=0) invMap[Map[i]]=i;
359     }
360     }
361    
362 jfenwick 3086 /* orders a Dudley_Util_ValueAndIndex array by value */
363 jgs 82 /* it is assumed that n is large */
364    
365 jfenwick 3086 int Dudley_Util_ValueAndIndex_compar(const void *arg1 , const void *arg2 ) {
366     Dudley_Util_ValueAndIndex *e1,*e2;
367     e1=(Dudley_Util_ValueAndIndex*) arg1;
368     e2=(Dudley_Util_ValueAndIndex*) arg2;
369 jgs 82 if (e1->value < e2->value) return -1;
370     if (e1->value > e2->value) return 1;
371 gross 763 if (e1->index < e2->index) return -1;
372     if (e1->index > e2->index) return 1;
373 jgs 82 return 0;
374     }
375 woo409 757
376 jfenwick 3086 void Dudley_Util_sortValueAndIndex(dim_t n,Dudley_Util_ValueAndIndex* array) {
377 jgs 82 /* OMP : needs parallelization !*/
378 jfenwick 3086 qsort(array,n,sizeof(Dudley_Util_ValueAndIndex),Dudley_Util_ValueAndIndex_compar);
379 jgs 82 }
380    
381    
382     /**************************************************************/
383    
384     /* calculates the minimum value from a dim X N integer array */
385    
386 jfenwick 3086 index_t Dudley_Util_getMinInt(dim_t dim,dim_t N,index_t* values) {
387 jgs 123 dim_t i,j;
388     index_t out,out_local;
389     out=INDEX_T_MAX;
390 jgs 82 if (values!=NULL && dim*N>0 ) {
391     out=values[0];
392 jgs 115 #pragma omp parallel private(out_local)
393     {
394     out_local=out;
395     #pragma omp for private(i,j) schedule(static)
396     for (j=0;j<N;j++) {
397     for (i=0;i<dim;i++) out_local=MIN(out_local,values[INDEX2(i,j,dim)]);
398     }
399     #pragma omp critical
400     out=MIN(out_local,out);
401 jgs 82 }
402     }
403     return out;
404     }
405    
406     /* calculates the maximum value from a dim X N integer array */
407    
408 jfenwick 3086 index_t Dudley_Util_getMaxInt(dim_t dim,dim_t N,index_t* values) {
409 jgs 123 dim_t i,j;
410     index_t out,out_local;
411     out=-INDEX_T_MAX;
412 jgs 82 if (values!=NULL && dim*N>0 ) {
413     out=values[0];
414 jgs 115 #pragma omp parallel private(out_local)
415     {
416     out_local=out;
417     #pragma omp for private(i,j) schedule(static)
418     for (j=0;j<N;j++) {
419 jfenwick 3114 for (i=0;i<dim;i++)
420     {
421     //printf("%d,%d,%d[%d] %d\n",i,j,dim,INDEX2(i,j,dim), values[INDEX2(i,j,dim)]);
422     out_local=MAX(out_local,values[INDEX2(i,j,dim)]);
423    
424     }
425 jgs 115 }
426     #pragma omp critical
427     out=MAX(out_local,out);
428     }
429 jgs 82 }
430     return out;
431     }
432 ksteube 1312 /**************************************************************/
433 jgs 82
434 ksteube 1312 /* calculates the minimum value from a dim X N integer array */
435    
436 jfenwick 3086 index_t Dudley_Util_getFlaggedMinInt(dim_t dim,dim_t N,index_t* values, index_t ignore) {
437 ksteube 1312 dim_t i,j;
438     index_t out,out_local;
439     out=INDEX_T_MAX;
440     if (values!=NULL && dim*N>0 ) {
441     out=values[0];
442     #pragma omp parallel private(out_local)
443     {
444     out_local=out;
445     #pragma omp for private(i,j) schedule(static)
446     for (j=0;j<N;j++) {
447     for (i=0;i<dim;i++) if (values[INDEX2(i,j,dim)]!=ignore) out_local=MIN(out_local,values[INDEX2(i,j,dim)]);
448     }
449     #pragma omp critical
450     out=MIN(out_local,out);
451     }
452     }
453     return out;
454     }
455    
456     /* calculates the maximum value from a dim X N integer array */
457    
458 jfenwick 3086 index_t Dudley_Util_getFlaggedMaxInt(dim_t dim,dim_t N,index_t* values, index_t ignore) {
459 ksteube 1312 dim_t i,j;
460     index_t out,out_local;
461     out=-INDEX_T_MAX;
462     if (values!=NULL && dim*N>0 ) {
463     out=values[0];
464     #pragma omp parallel private(out_local)
465     {
466     out_local=out;
467     #pragma omp for private(i,j) schedule(static)
468     for (j=0;j<N;j++) {
469     for (i=0;i<dim;i++) if (values[INDEX2(i,j,dim)]!=ignore) out_local=MAX(out_local,values[INDEX2(i,j,dim)]);
470     }
471     #pragma omp critical
472     out=MAX(out_local,out);
473     }
474     }
475     return out;
476     }
477    
478 jgs 82 /* set the index of the positive entries in mask. The length of index is returned. */
479    
480 jfenwick 3086 dim_t Dudley_Util_packMask(dim_t N,index_t* mask,index_t* index) {
481 jgs 123 dim_t out,k;
482 jgs 82 out=0;
483     /*OMP */
484     for (k=0;k<N;k++) {
485     if (mask[k]>=0) {
486     index[out]=k;
487     out++;
488     }
489     }
490     return out;
491     }
492    
493     /* returns true if array contains value */
494 jfenwick 3086 bool_t Dudley_Util_isAny(dim_t N,index_t* array,index_t value) {
495 jgs 123 bool_t out=FALSE;
496     dim_t i;
497 jgs 82 #pragma omp parallel for private(i) schedule(static) reduction(||:out)
498 jgs 115 for (i=0;i<N;i++) out = out || (array[i]==value);
499 jgs 82 return out;
500     }
501 jgs 113 /* calculates the cummultative sum in array and returns the total sum */
502 jfenwick 3086 index_t Dudley_Util_cumsum(dim_t N,index_t* array) {
503 jgs 123 index_t out=0,tmp;
504     dim_t i;
505 jgs 113 #ifdef _OPENMP
506 gross 1564 index_t *partial_sums=NULL, sum;
507     partial_sums=TMPMEMALLOC(omp_get_max_threads(),index_t);
508 jgs 113 #pragma omp parallel private(sum,i,tmp)
509     {
510     sum=0;
511 jgs 115 #pragma omp for schedule(static)
512     for (i=0;i<N;++i) sum+=array[i];
513 jgs 113 partial_sums[omp_get_thread_num()]=sum;
514 jgs 115 #pragma omp barrier
515 jgs 113 #pragma omp master
516     {
517     out=0;
518     for (i=0;i<omp_get_max_threads();++i) {
519     tmp=out;
520     out+=partial_sums[i];
521     partial_sums[i]=tmp;
522     }
523     }
524 jgs 115 #pragma omp barrier
525 jgs 113 sum=partial_sums[omp_get_thread_num()];
526 jgs 115 #pragma omp for schedule(static)
527     for (i=0;i<N;++i) {
528     tmp=sum;
529     sum+=array[i];
530     array[i]=tmp;
531     }
532 jgs 113 }
533 gross 1564 TMPMEMFREE(partial_sums);
534 jgs 113 #else
535     for (i=0;i<N;++i) {
536     tmp=out;
537     out+=array[i];
538     array[i]=tmp;
539     }
540     #endif
541     return out;
542     }
543 jfenwick 3086 void Dudley_Util_setValuesInUse(const index_t *values, const dim_t numValues, dim_t *numValuesInUse, index_t **valuesInUse, Paso_MPIInfo* mpiinfo)
544 gross 1716 {
545     dim_t i;
546     index_t lastFoundValue=INDEX_T_MIN, minFoundValue, local_minFoundValue, *newValuesInUse=NULL;
547     register index_t itmp;
548     bool_t allFound=FALSE;
549     dim_t nv=0;
550 jgs 82
551 gross 1716 while (! allFound) {
552     /*
553     * find smallest value bigger than lastFoundValue
554     */
555     minFoundValue=INDEX_T_MAX;
556     #pragma omp parallel private(local_minFoundValue)
557     {
558     local_minFoundValue=minFoundValue;
559     #pragma omp for private(i,itmp) schedule(static)
560     for (i=0;i< numValues;i++) {
561     itmp=values[i];
562     if ((itmp>lastFoundValue) && (itmp<local_minFoundValue)) local_minFoundValue=itmp;
563     }
564     #pragma omp critical
565 gross 2425 {
566     if (local_minFoundValue<minFoundValue) minFoundValue=local_minFoundValue;
567     }
568    
569 gross 1716 }
570     #ifdef PASO_MPI
571     local_minFoundValue=minFoundValue;
572 gross 2425 MPI_Allreduce(&local_minFoundValue,&minFoundValue, 1, MPI_INT, MPI_MIN, mpiinfo->comm );
573 gross 1716 #endif
574     /* if we found a new tag we need to add this too the valuesInUseList */
575    
576     if (minFoundValue < INDEX_T_MAX) {
577     newValuesInUse=MEMALLOC(nv+1,index_t);
578     if (*valuesInUse!=NULL) {
579     memcpy(newValuesInUse,*valuesInUse,sizeof(index_t)*nv);
580     MEMFREE(*valuesInUse);
581     }
582     newValuesInUse[nv]=minFoundValue;
583     *valuesInUse=newValuesInUse;
584     newValuesInUse=NULL;
585     nv++;
586     lastFoundValue=minFoundValue;
587     } else {
588     allFound=TRUE;
589     }
590     }
591     *numValuesInUse=nv;
592     }
593    
594    
595 bcumming 751 #ifdef PASO_MPI
596 jfenwick 3086 void Dudley_printDoubleArray( FILE *fid, dim_t n, double *array, char *name )
597 bcumming 751 {
598     index_t i;
599    
600     if( name )
601     fprintf( fid, "%s [ ", name );
602     else
603     fprintf( fid, "[ " );
604 bcumming 782 for( i=0; i<(n<60 ? n : 60); i++ )
605 bcumming 751 fprintf( fid, "%g ", array[i] );
606     if( n>=30 )
607     fprintf( fid, "... " );
608     fprintf( fid, "]\n" );
609     }
610 jfenwick 3086 void Dudley_printIntArray( FILE *fid, dim_t n, int *array, char *name )
611 bcumming 751 {
612     index_t i;
613    
614     if( name )
615     fprintf( fid, "%s [ ", name );
616     else
617     fprintf( fid, "[ " );
618 bcumming 782 for( i=0; i<(n<60 ? n : 60); i++ )
619 bcumming 751 fprintf( fid, "%d ", array[i] );
620     if( n>=30 )
621     fprintf( fid, "... " );
622     fprintf( fid, "]\n" );
623     }
624 jfenwick 3086 void Dudley_printMaskArray( FILE *fid, dim_t n, int *array, char *name )
625 bcumming 751 {
626     index_t i;
627    
628     if( name )
629     fprintf( fid, "%s [ ", name );
630     else
631     fprintf( fid, "[ " );
632 bcumming 782 for( i=0; i<(n<60 ? n : 60); i++ )
633 bcumming 751 if( array[i]!=-1 )
634 bcumming 782 fprintf( fid, "%3d ", array[i] );
635 bcumming 751 else
636 bcumming 782 fprintf( fid, " * " );
637 bcumming 751 if( n>=30 )
638     fprintf( fid, "... " );
639     fprintf( fid, "]\n" );
640     }
641     #endif

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26