/[escript]/branches/intelc_win32/finley/src/Util.c
ViewVC logotype

Contents of /branches/intelc_win32/finley/src/Util.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 752 - (show annotations)
Mon Jun 26 02:25:41 2006 UTC (12 years, 11 months ago) by woo409
File MIME type: text/plain
File size: 16030 byte(s)
+ Added a qsort.c file which contains a drop in replacement for qsort (call it as qsortG). This one appears to be a stable implementation and the test .msh files on windows have been set up to be the same as unix again except for the exponent digits (3 instead of 2).
With ALL the qsorts replaced with qsortG only two tests fail now on win32:
test_normal_onFunctionOnContactOne
test_normal_onFunctionOnContactZero

Both give wrong result errors.

I will check this same code on the altix (including the use of qsortG and see if Altix has the same problem.
1 /*
2 ************************************************************
3 * Copyright 2006 by ACcESS MNRF *
4 * *
5 * http://www.access.edu.au *
6 * Primary Business: Queensland, Australia *
7 * Licensed under the Open Software License version 3.0 *
8 * http://www.opensource.org/licenses/osl-3.0.php *
9 * *
10 ************************************************************
11 */
12
13 /**************************************************************/
14
15 /* Some utility routines: */
16
17 /**************************************************************/
18
19 /* author: gross@access.edu.au */
20 /* Version: $Id$ */
21
22 /**************************************************************/
23
24 #include "Finley.h"
25 #include "Util.h"
26 #include "qsortG.h"
27 #ifdef _OPENMP
28 #include <omp.h>
29 #endif
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 Finley_Util_anyNonZeroDouble(dim_t N, double* values) {
36 dim_t q;
37 for (q=0;q<N;++q) if (ABS(values[q])>0) return TRUE;
38 return FALSE;
39 }
40 /**************************************************************/
41
42 /* gathers double values out from in by index: */
43
44 /* out(1:numData,1:len)=in(1:numData,index(1:len)) */
45
46 void Finley_Util_Gather_double(dim_t len,index_t* index,dim_t numData,double* in, double * out){
47 dim_t s,i;
48 for (s=0;s<len;s++) {
49 for (i=0;i<numData;i++) {
50 out[INDEX2(i,s,numData)]=in[INDEX2(i,index[s],numData)];
51 }
52 }
53 }
54
55 /**************************************************************/
56
57
58 /* gathers maybelong values out from in by index: */
59
60 /* out(1:numData,1:len)=in(1:numData,index(1:len)) */
61
62 void Finley_Util_Gather_int(dim_t len,index_t* index,dim_t numData, index_t* in, index_t * out){
63 dim_t s,i;
64 for (s=0;s<len;s++) {
65 for (i=0;i<numData;i++) {
66 out[INDEX2(i,s,numData)]=in[INDEX2(i,index[s],numData)];
67 }
68 }
69 }
70
71 /**************************************************************/
72
73 /* adds a vector in into out using and index. */
74
75 /* out(1:numData,index(1:len))+=in(1:numData,1:len) */
76
77 void Finley_Util_AddScatter(dim_t len,index_t* index,dim_t numData,double* in,double * out){
78 dim_t i,s;
79 for (s=0;s<len;s++) {
80 for(i=0;i<numData;i++) {
81 #pragma omp atomic
82 out[INDEX2(i,index[s],numData)]+=in[INDEX2(i,s,numData)];
83 }
84 }
85 }
86
87 /* multiplies two matrices */
88
89 /* A(1:A1,1:A2)=B(1:A1,1:B2)*C(1:B2,1:A2) */
90
91 void Finley_Util_SmallMatMult(dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
92 dim_t i,j,s;
93 for (i=0;i<A1*A2;i++) A[i]=0;
94 for (i=0;i<A1;i++) {
95 for (j=0;j<A2;j++) {
96 for (s=0;s<B2;s++) {
97 A[INDEX2(i,j,A1)]+=B[INDEX2(i,s,A1)]*C[INDEX2(s,j,B2)];
98 }
99 }
100 }
101 }
102
103 /* multiplies a two sets of matries: */
104
105 /* A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2,i) i=1,len */
106
107 void Finley_Util_SmallMatSetMult(dim_t len,dim_t A1,dim_t A2, double* A, dim_t B2, double*B, double* C) {
108 dim_t q,i,j,s;
109 for (i=0;i<A1*A2*len;i++) A[i]=0;
110 for (q=0;q<len;q++) {
111 for (i=0;i<A1;i++) {
112 for (j=0;j<A2;j++) {
113 for (s=0;s<B2;s++) {
114 A[INDEX3(i,j,q,A1,A2)]+=B[INDEX3(i,s,q,A1,B2)]*C[INDEX3(s,j,q,B2,A2)];
115 }
116 }
117 }
118 }
119 }
120 /* inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
121 /* the determinante is returned. */
122
123 void Finley_Util_InvertSmallMat(dim_t len,dim_t dim,double* A,double *invA, double* det){
124 dim_t q;
125 register double D,A11,A12,A13,A21,A22,A23,A31,A32,A33;
126
127 switch(dim) {
128 case 1:
129 for (q=0;q<len;q++) {
130 D=A[q];
131 if (ABS(D) > 0 ){
132 det[q]=D;
133 D=1./D;
134 invA[q]=D;
135 } else {
136 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
137 return;
138 }
139 }
140 break;
141
142 case 2:
143 for (q=0;q<len;q++) {
144 A11=A[INDEX3(0,0,q,2,2)];
145 A12=A[INDEX3(0,1,q,2,2)];
146 A21=A[INDEX3(1,0,q,2,2)];
147 A22=A[INDEX3(1,1,q,2,2)];
148
149 D = A11*A22-A12*A21;
150 if (ABS(D) > 0 ){
151 det[q]=D;
152 D=1./D;
153 invA[INDEX3(0,0,q,2,2)]= A22*D;
154 invA[INDEX3(1,0,q,2,2)]=-A21*D;
155 invA[INDEX3(0,1,q,2,2)]=-A12*D;
156 invA[INDEX3(1,1,q,2,2)]= A11*D;
157 } else {
158 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
159 return;
160 }
161 }
162 break;
163
164 case 3:
165 for (q=0;q<len;q++) {
166 A11=A[INDEX3(0,0,q,3,3)];
167 A21=A[INDEX3(1,0,q,3,3)];
168 A31=A[INDEX3(2,0,q,3,3)];
169 A12=A[INDEX3(0,1,q,3,3)];
170 A22=A[INDEX3(1,1,q,3,3)];
171 A32=A[INDEX3(2,1,q,3,3)];
172 A13=A[INDEX3(0,2,q,3,3)];
173 A23=A[INDEX3(1,2,q,3,3)];
174 A33=A[INDEX3(2,2,q,3,3)];
175
176 D = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
177 if (ABS(D) > 0 ){
178 det[q] =D;
179 D=1./D;
180 invA[INDEX3(0,0,q,3,3)]=(A22*A33-A23*A32)*D;
181 invA[INDEX3(1,0,q,3,3)]=(A31*A23-A21*A33)*D;
182 invA[INDEX3(2,0,q,3,3)]=(A21*A32-A31*A22)*D;
183 invA[INDEX3(0,1,q,3,3)]=(A13*A32-A12*A33)*D;
184 invA[INDEX3(1,1,q,3,3)]=(A11*A33-A31*A13)*D;
185 invA[INDEX3(2,1,q,3,3)]=(A12*A31-A11*A32)*D;
186 invA[INDEX3(0,2,q,3,3)]=(A12*A23-A13*A22)*D;
187 invA[INDEX3(1,2,q,3,3)]=(A13*A21-A11*A23)*D;
188 invA[INDEX3(2,2,q,3,3)]=(A11*A22-A12*A21)*D;
189 } else {
190 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: Non-regular matrix");
191 return;
192 }
193 }
194 break;
195
196 }
197 return;
198 }
199
200 /* sets the derterminate of a set of dim x dim matrices A(:,:,1:len) with dim=1,2,3 */
201
202 void Finley_Util_DetOfSmallMat(dim_t len,dim_t dim,double* A, double* det){
203 dim_t q;
204 register double A11,A12,A13,A21,A22,A23,A31,A32,A33;
205
206 switch(dim) {
207 case 1:
208 for (q=0;q<len;q++) {
209 det[q]=A[q];
210 }
211 break;
212
213 case 2:
214 for (q=0;q<len;q++) {
215 A11=A[INDEX3(0,0,q,2,2)];
216 A12=A[INDEX3(0,1,q,2,2)];
217 A21=A[INDEX3(1,0,q,2,2)];
218 A22=A[INDEX3(1,1,q,2,2)];
219
220 det[q] = A11*A22-A12*A21;
221 }
222 break;
223
224 case 3:
225 for (q=0;q<len;q++) {
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 det[q] = A11*(A22*A33-A23*A32)+ A12*(A31*A23-A21*A33)+A13*(A21*A32-A31*A22);
237 }
238 break;
239
240 }
241 return;
242 }
243 /* returns the normalized vector Normal[dim,len] orthogonal to A(:,0,q) and A(:,1,q) in the case of dim=3 */
244 /* or the vector A(:,0,q) in the case of dim=2 */
245
246 void Finley_NormalVector(dim_t len, dim_t dim, dim_t dim1, double* A,double* Normal) {
247 dim_t q;
248 register double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33,length,invlength;
249
250 switch(dim) {
251 case 1:
252 for (q=0;q<len;q++) Normal[q] =1;
253 break;
254 case 2:
255 for (q=0;q<len;q++) {
256 A11=A[INDEX3(0,0,q,2,dim1)];
257 A21=A[INDEX3(1,0,q,2,dim1)];
258 length = sqrt(A11*A11+A21*A21);
259 if (! length>0) {
260 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: area equals zero.");
261 return;
262 } else {
263 invlength=1./length;
264 Normal[INDEX2(0,q,2)]=A21*invlength;
265 Normal[INDEX2(1,q,2)]=-A11*invlength;
266 }
267 }
268 break;
269 case 3:
270 for (q=0;q<len;q++) {
271 A11=A[INDEX3(0,0,q,3,dim1)];
272 A21=A[INDEX3(1,0,q,3,dim1)];
273 A31=A[INDEX3(2,0,q,3,dim1)];
274 A12=A[INDEX3(0,1,q,3,dim1)];
275 A22=A[INDEX3(1,1,q,3,dim1)];
276 A32=A[INDEX3(2,1,q,3,dim1)];
277 CO_A13=A21*A32-A31*A22;
278 CO_A23=A31*A12-A11*A32;
279 CO_A33=A11*A22-A21*A12;
280 length=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
281 if (! length>0) {
282 Finley_setError(ZERO_DIVISION_ERROR,"__FILE__: area equals zero.");
283 return;
284 } else {
285 invlength=1./length;
286 Normal[INDEX2(0,q,3)]=CO_A13*invlength;
287 Normal[INDEX2(1,q,3)]=CO_A23*invlength;
288 Normal[INDEX2(2,q,3)]=CO_A33*invlength;
289 }
290
291 }
292 break;
293
294 }
295 return;
296 }
297
298 /* 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 */
299 /* or the vector A(:,0,q) in the case of dim=2 */
300
301 void Finley_LengthOfNormalVector(dim_t len, dim_t dim, dim_t dim1, double* A,double* length) {
302 dim_t q;
303 double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33;
304
305 switch(dim) {
306 case 1:
307 for (q=0;q<len;q++) length[q] =1;
308 break;
309 case 2:
310 for (q=0;q<len;q++) {
311 A11=A[INDEX3(0,0,q,2,dim1)];
312 A21=A[INDEX3(1,0,q,2,dim1)];
313 length[q] = sqrt(A11*A11+A21*A21);
314 }
315 break;
316 case 3:
317 for (q=0;q<len;q++) {
318 A11=A[INDEX3(0,0,q,3,dim1)];
319 A21=A[INDEX3(1,0,q,3,dim1)];
320 A31=A[INDEX3(2,0,q,3,dim1)];
321 A12=A[INDEX3(0,1,q,3,dim1)];
322 A22=A[INDEX3(1,1,q,3,dim1)];
323 A32=A[INDEX3(2,1,q,3,dim1)];
324 CO_A13=A21*A32-A31*A22;
325 CO_A23=A31*A12-A11*A32;
326 CO_A33=A11*A22-A21*A12;
327 length[q]=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
328 }
329 break;
330
331 }
332 return;
333 }
334
335 /* inverts the map map of length len */
336 /* there is no range checking! */
337 /* at output Map[invMap[i]]=i for i=0:lenInvMap */
338
339 void Finley_Util_InvertMap(dim_t lenInvMap, index_t* invMap,dim_t lenMap, index_t* Map) {
340 dim_t i;
341 for (i=0;i<lenInvMap;i++) invMap[i]=0;
342 for (i=0;i<lenMap;i++) {
343 if (Map[i]>=0) invMap[Map[i]]=i;
344 }
345 }
346
347 /* orders a Finley_Util_ValueAndIndex array by value */
348 /* it is assumed that n is large */
349
350 int Finley_Util_ValueAndIndex_compar(const void *arg1 , const void *arg2 ) {
351 Finley_Util_ValueAndIndex *e1,*e2;
352 e1=(Finley_Util_ValueAndIndex*) arg1;
353 e2=(Finley_Util_ValueAndIndex*) arg2;
354 if (e1->value < e2->value) return -1;
355 if (e1->value > e2->value) return 1;
356 return 0;
357 }
358
359 void Finley_Util_sortValueAndIndex(dim_t n,Finley_Util_ValueAndIndex* array) {
360 /* OMP : needs parallelization !*/
361 qsortG(array,n,sizeof(Finley_Util_ValueAndIndex),Finley_Util_ValueAndIndex_compar);
362 }
363
364
365 /**************************************************************/
366
367 /* calculates the minimum value from a dim X N integer array */
368
369 index_t Finley_Util_getMinInt(dim_t dim,dim_t N,index_t* values) {
370 dim_t i,j;
371 index_t out,out_local;
372 out=INDEX_T_MAX;
373 if (values!=NULL && dim*N>0 ) {
374 out=values[0];
375 #pragma omp parallel private(out_local)
376 {
377 out_local=out;
378 #pragma omp for private(i,j) schedule(static)
379 for (j=0;j<N;j++) {
380 for (i=0;i<dim;i++) out_local=MIN(out_local,values[INDEX2(i,j,dim)]);
381 }
382 #pragma omp critical
383 out=MIN(out_local,out);
384 }
385 }
386 return out;
387 }
388
389 /* calculates the maximum value from a dim X N integer array */
390
391 index_t Finley_Util_getMaxInt(dim_t dim,dim_t N,index_t* values) {
392 dim_t i,j;
393 index_t out,out_local;
394 out=-INDEX_T_MAX;
395 if (values!=NULL && dim*N>0 ) {
396 out=values[0];
397 #pragma omp parallel private(out_local)
398 {
399 out_local=out;
400 #pragma omp for private(i,j) schedule(static)
401 for (j=0;j<N;j++) {
402 for (i=0;i<dim;i++) out_local=MAX(out_local,values[INDEX2(i,j,dim)]);
403 }
404 #pragma omp critical
405 out=MAX(out_local,out);
406 }
407 }
408 return out;
409 }
410
411 /* set the index of the positive entries in mask. The length of index is returned. */
412
413 dim_t Finley_Util_packMask(dim_t N,index_t* mask,index_t* index) {
414 dim_t out,k;
415 out=0;
416 /*OMP */
417 for (k=0;k<N;k++) {
418 if (mask[k]>=0) {
419 index[out]=k;
420 out++;
421 }
422 }
423 return out;
424 }
425
426 /* returns true if array contains value */
427 bool_t Finley_Util_isAny(dim_t N,index_t* array,index_t value) {
428 bool_t out=FALSE;
429 dim_t i;
430 #pragma omp parallel for private(i) schedule(static) reduction(||:out)
431 for (i=0;i<N;i++) out = out || (array[i]==value);
432 return out;
433 }
434 /* calculates the cummultative sum in array and returns the total sum */
435 index_t Finley_Util_cumsum(dim_t N,index_t* array) {
436 index_t out=0,tmp;
437 dim_t i;
438 #ifdef _OPENMP
439 index_t partial_sums[omp_get_max_threads()],sum;
440 #pragma omp parallel private(sum,i,tmp)
441 {
442 sum=0;
443 #pragma omp for schedule(static)
444 for (i=0;i<N;++i) sum+=array[i];
445 partial_sums[omp_get_thread_num()]=sum;
446 #pragma omp barrier
447 #pragma omp master
448 {
449 out=0;
450 for (i=0;i<omp_get_max_threads();++i) {
451 tmp=out;
452 out+=partial_sums[i];
453 partial_sums[i]=tmp;
454 }
455 }
456 #pragma omp barrier
457 sum=partial_sums[omp_get_thread_num()];
458 #pragma omp for schedule(static)
459 for (i=0;i<N;++i) {
460 tmp=sum;
461 sum+=array[i];
462 array[i]=tmp;
463 }
464 }
465 #else
466 for (i=0;i<N;++i) {
467 tmp=out;
468 out+=array[i];
469 array[i]=tmp;
470 }
471 #endif
472 return out;
473 }
474
475 void Finley_copyDouble(dim_t n,double* source, double* target) {
476 dim_t i;
477 for (i=0;i<n;i++) target[i]=source[i];
478 }
479
480 /*
481 * Revision 1.8 2005/08/12 01:45:43 jgs
482 * erge of development branch dev-02 back to main trunk on 2005-08-12
483 *
484 * Revision 1.7.2.2 2005/09/07 06:26:22 gross
485 * the solver from finley are put into the standalone package paso now
486 *
487 * Revision 1.7.2.1 2005/08/04 22:41:11 gross
488 * some extra routines for finley that might speed-up RHS assembling in some cases (not actived right now)
489 *
490 * Revision 1.7 2005/07/08 04:07:59 jgs
491 * Merge of development branch back to main trunk on 2005-07-08
492 *
493 * Revision 1.1.1.1.2.4 2005/06/29 02:34:57 gross
494 * some changes towards 64 integers in finley
495 *
496 * Revision 1.1.1.1.2.3 2005/03/02 23:35:06 gross
497 * reimplementation of the ILU in Finley. block size>1 still needs some testing
498 *
499 * Revision 1.1.1.1.2.2 2005/02/18 02:27:31 gross
500 * two function that will be used for a reimplementation of the ILU preconditioner
501 *
502 * Revision 1.1.1.1.2.1 2004/11/12 06:58:19 gross
503 * 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
504 *
505 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
506 * initial import of project esys2
507 *
508 * Revision 1.3 2004/08/26 12:03:52 gross
509 * Some other bug in Finley_Assemble_gradient fixed.
510 *
511 * Revision 1.2 2004/07/02 04:21:13 gross
512 * Finley C code has been included
513 *
514 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
515 * Initial version of eys using boost-python.
516 *
517 *
518 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26