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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26