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 |
|
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 |
qsort(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 |
*/ |