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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4286 - (hide annotations)
Thu Mar 7 04:28:11 2013 UTC (6 years, 6 months ago) by caltinay
File MIME type: text/plain
File size: 18730 byte(s)
Assorted spelling fixes.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26