/[escript]/trunk/escript/src/LocalOps.h
ViewVC logotype

Contents of /trunk/escript/src/LocalOps.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (10 years, 3 months ago) by ksteube
File MIME type: text/plain
File size: 15808 byte(s)
Copyright updated in all files

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2008 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 #if !defined escript_LocalOps_H
16 #define escript_LocalOps_H
17 #ifdef __INTEL_COMPILER
18 # include <mathimf.h>
19 #else
20 # include <math.h>
21 #endif
22 #ifndef M_PI
23 # define M_PI 3.14159265358979323846 /* pi */
24 #endif
25
26 namespace escript {
27
28
29 /**
30 \brief
31 solves a 1x1 eigenvalue A*V=ev*V problem
32
33 \param A00 Input - A_00
34 \param ev0 Output - eigenvalue
35 */
36 inline
37 void eigenvalues1(const double A00,double* ev0) {
38
39 *ev0=A00;
40
41 }
42 /**
43 \brief
44 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
45
46 \param A00 Input - A_00
47 \param A01 Input - A_01
48 \param A11 Input - A_11
49 \param ev0 Output - smallest eigenvalue
50 \param ev1 Output - largest eigenvalue
51 */
52 inline
53 void eigenvalues2(const double A00,const double A01,const double A11,
54 double* ev0, double* ev1) {
55 const register double trA=(A00+A11)/2.;
56 const register double A_00=A00-trA;
57 const register double A_11=A11-trA;
58 const register double s=sqrt(A01*A01-A_00*A_11);
59 *ev0=trA-s;
60 *ev1=trA+s;
61 }
62 /**
63 \brief
64 solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
65
66 \param A00 Input - A_00
67 \param A01 Input - A_01
68 \param A02 Input - A_02
69 \param A11 Input - A_11
70 \param A12 Input - A_12
71 \param A22 Input - A_22
72 \param ev0 Output - smallest eigenvalue
73 \param ev1 Output - eigenvalue
74 \param ev2 Output - largest eigenvalue
75 */
76 inline
77 void eigenvalues3(const double A00, const double A01, const double A02,
78 const double A11, const double A12,
79 const double A22,
80 double* ev0, double* ev1,double* ev2) {
81
82 const register double trA=(A00+A11+A22)/3.;
83 const register double A_00=A00-trA;
84 const register double A_11=A11-trA;
85 const register double A_22=A22-trA;
86 const register double A01_2=A01*A01;
87 const register double A02_2=A02*A02;
88 const register double A12_2=A12*A12;
89 const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
90 if (p<=0.) {
91 *ev2=trA;
92 *ev1=trA;
93 *ev0=trA;
94
95 } else {
96 const register double q=(A02_2*A_11+A12_2*A_00+A01_2*A_22)-(A_00*A_11*A_22+2*A01*A12*A02);
97 const register double sq_p=sqrt(p/3.);
98 register double z=-q/(2*pow(sq_p,3));
99 if (z<-1.) {
100 z=-1.;
101 } else if (z>1.) {
102 z=1.;
103 }
104 const register double alpha_3=acos(z)/3.;
105 *ev2=trA+2.*sq_p*cos(alpha_3);
106 *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.);
107 *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.);
108 }
109 }
110 /**
111 \brief
112 solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
113
114 \param A00 Input - A_00
115 \param ev0 Output - eigenvalue
116 \param V00 Output - eigenvector
117 \param tol Input - tolerance to identify to eigenvalues
118 */
119 inline
120 void eigenvalues_and_eigenvectors1(const double A00,double* ev0,double* V00,const double tol)
121 {
122 eigenvalues1(A00,ev0);
123 *V00=1.;
124 return;
125 }
126 /**
127 \brief
128 returns a non-zero vector in the kernel of [[A00,A01],[A01,A11]] assuming that the kernel dimension is at least 1.
129
130 \param A00 Input - matrix component
131 \param A10 Input - matrix component
132 \param A01 Input - matrix component
133 \param A11 Input - matrix component
134 \param V0 Output - vector component
135 \param V1 Output - vector component
136 */
137 inline
138 void vectorInKernel2(const double A00,const double A10,const double A01,const double A11,
139 double* V0, double*V1)
140 {
141 register double absA00=fabs(A00);
142 register double absA10=fabs(A10);
143 register double absA01=fabs(A01);
144 register double absA11=fabs(A11);
145 register double m=absA11>absA10 ? absA11 : absA10;
146 if (absA00>m || absA01>m) {
147 *V0=-A01;
148 *V1=A00;
149 } else {
150 if (m<=0) {
151 *V0=1.;
152 *V1=0.;
153 } else {
154 *V0=A11;
155 *V1=-A10;
156 }
157 }
158 }
159 /**
160 \brief
161 returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]]
162 assuming that the kernel dimension is at least 1 and A00 is non zero.
163
164 \param A00 Input - matrix component
165 \param A10 Input - matrix component
166 \param A20 Input - matrix component
167 \param A01 Input - matrix component
168 \param A11 Input - matrix component
169 \param A21 Input - matrix component
170 \param A02 Input - matrix component
171 \param A12 Input - matrix component
172 \param A22 Input - matrix component
173 \param V0 Output - vector component
174 \param V1 Output - vector component
175 \param V2 Output - vector component
176 */
177 inline
178 void vectorInKernel3__nonZeroA00(const double A00,const double A10,const double A20,
179 const double A01,const double A11,const double A21,
180 const double A02,const double A12,const double A22,
181 double* V0,double* V1,double* V2)
182 {
183 double TEMP0,TEMP1;
184 register const double I00=1./A00;
185 register const double IA10=I00*A10;
186 register const double IA20=I00*A20;
187 vectorInKernel2(A11-IA10*A01,A12-IA10*A02,
188 A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
189 *V0=-(A10*TEMP0+A20*TEMP1);
190 *V1=A00*TEMP0;
191 *V2=A00*TEMP1;
192 }
193
194 /**
195 \brief
196 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
197 ordered by increasing value and eigen vectors are normalizeVector3d such that
198 length is zero and first non-zero component is positive.
199
200 \param A00 Input - A_00
201 \param A01 Input - A_01
202 \param A11 Input - A_11
203 \param ev0 Output - smallest eigenvalue
204 \param ev1 Output - eigenvalue
205 \param V00 Output - eigenvector componenent coresponding to ev0
206 \param V10 Output - eigenvector componenent coresponding to ev0
207 \param V01 Output - eigenvector componenent coresponding to ev1
208 \param V11 Output - eigenvector componenent coresponding to ev1
209 \param tol Input - tolerance to identify to eigenvalues
210 */
211 inline
212 void eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11,
213 double* ev0, double* ev1,
214 double* V00, double* V10, double* V01, double* V11,
215 const double tol)
216 {
217 double TEMP0,TEMP1;
218 eigenvalues2(A00,A01,A11,ev0,ev1);
219 const register double absev0=fabs(*ev0);
220 const register double absev1=fabs(*ev1);
221 register double max_ev=absev0>absev1 ? absev0 : absev1;
222 if (fabs((*ev0)-(*ev1))<tol*max_ev) {
223 *V00=1.;
224 *V10=0.;
225 *V01=0.;
226 *V11=1.;
227 } else {
228 vectorInKernel2(A00-(*ev0),A01,A01,A11-(*ev0),&TEMP0,&TEMP1);
229 const register double scale=1./sqrt(TEMP0*TEMP0+TEMP1*TEMP1);
230 if (TEMP0<0.) {
231 *V00=-TEMP0*scale;
232 *V10=-TEMP1*scale;
233 if (TEMP1<0.) {
234 *V01= *V10;
235 *V11=-(*V00);
236 } else {
237 *V01=-(*V10);
238 *V11= (*V10);
239 }
240 } else if (TEMP0>0.) {
241 *V00=TEMP0*scale;
242 *V10=TEMP1*scale;
243 if (TEMP1<0.) {
244 *V01=-(*V10);
245 *V11= (*V00);
246 } else {
247 *V01= (*V10);
248 *V11=-(*V00);
249 }
250 } else {
251 *V00=0.;
252 *V10=1;
253 *V11=0.;
254 *V01=1.;
255 }
256 }
257 }
258 /**
259 \brief
260 nomalizes a 3-d vector such that length is one and first non-zero component is positive.
261
262 \param V0 - vector componenent
263 \param V1 - vector componenent
264 \param V2 - vector componenent
265 */
266 inline
267 void normalizeVector3(double* V0,double* V1,double* V2)
268 {
269 register double s;
270 if (*V0>0) {
271 s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
272 *V0*=s;
273 *V1*=s;
274 *V2*=s;
275 } else if (*V0<0) {
276 s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
277 *V0*=s;
278 *V1*=s;
279 *V2*=s;
280 } else {
281 if (*V1>0) {
282 s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
283 *V1*=s;
284 *V2*=s;
285 } else if (*V1<0) {
286 s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
287 *V1*=s;
288 *V2*=s;
289 } else {
290 *V2=1.;
291 }
292 }
293 }
294 /**
295 \brief
296 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
297 ordered by increasing value and eigen vectors are normalizeVector3d such that
298 length is zero and first non-zero component is positive.
299
300 \param A00 Input - A_00
301 \param A01 Input - A_01
302 \param A11 Input - A_11
303 \param ev0 Output - smallest eigenvalue
304 \param ev1 Output - eigenvalue
305 \param V00 Output - eigenvector componenent coresponding to ev0
306 \param V10 Output - eigenvector componenent coresponding to ev0
307 \param V01 Output - eigenvector componenent coresponding to ev1
308 \param V11 Output - eigenvector componenent coresponding to ev1
309 \param tol Input - tolerance to identify to eigenvalues
310 */
311 inline
312 void eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02,
313 const double A11, const double A12, const double A22,
314 double* ev0, double* ev1, double* ev2,
315 double* V00, double* V10, double* V20,
316 double* V01, double* V11, double* V21,
317 double* V02, double* V12, double* V22,
318 const double tol)
319 {
320 register const double absA01=fabs(A01);
321 register const double absA02=fabs(A02);
322 register const double m=absA01>absA02 ? absA01 : absA02;
323 if (m<=0) {
324 double TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1;
325 eigenvalues_and_eigenvectors2(A11,A12,A22,
326 &TEMP_EV0,&TEMP_EV1,
327 &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
328 if (A00<=TEMP_EV0) {
329 *V00=1.;
330 *V10=0.;
331 *V20=0.;
332 *V01=0.;
333 *V11=TEMP_V00;
334 *V21=TEMP_V10;
335 *V02=0.;
336 *V12=TEMP_V01;
337 *V22=TEMP_V11;
338 *ev0=A00;
339 *ev1=TEMP_EV0;
340 *ev2=TEMP_EV1;
341 } else if (A00>TEMP_EV1) {
342 *V02=1.;
343 *V12=0.;
344 *V22=0.;
345 *V00=0.;
346 *V10=TEMP_V00;
347 *V20=TEMP_V10;
348 *V01=0.;
349 *V11=TEMP_V01;
350 *V21=TEMP_V11;
351 *ev0=TEMP_EV0;
352 *ev1=TEMP_EV1;
353 *ev2=A00;
354 } else {
355 *V01=1.;
356 *V11=0.;
357 *V21=0.;
358 *V00=0.;
359 *V10=TEMP_V00;
360 *V20=TEMP_V10;
361 *V02=0.;
362 *V12=TEMP_V01;
363 *V22=TEMP_V11;
364 *ev0=TEMP_EV0;
365 *ev1=A00;
366 *ev2=TEMP_EV1;
367 }
368 } else {
369 eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2);
370 const register double absev0=fabs(*ev0);
371 const register double absev1=fabs(*ev1);
372 const register double absev2=fabs(*ev2);
373 register double max_ev=absev0>absev1 ? absev0 : absev1;
374 max_ev=max_ev>absev2 ? max_ev : absev2;
375 const register double d_01=fabs((*ev0)-(*ev1));
376 const register double d_12=fabs((*ev1)-(*ev2));
377 const register double max_d=d_01>d_12 ? d_01 : d_12;
378 if (max_d<=tol*max_ev) {
379 *V00=1.;
380 *V10=0;
381 *V20=0;
382 *V01=0;
383 *V11=1.;
384 *V21=0;
385 *V02=0;
386 *V12=0;
387 *V22=1.;
388 } else {
389 const register double S00=A00-(*ev0);
390 const register double absS00=fabs(S00);
391 if (fabs(S00)>m) {
392 vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);
393 } else if (absA02<m) {
394 vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);
395 } else {
396 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
397 }
398 normalizeVector3(V00,V10,V20);;
399 const register double T00=A00-(*ev2);
400 const register double absT00=fabs(T00);
401 if (fabs(T00)>m) {
402 vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);
403 } else if (absA02<m) {
404 vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);
405 } else {
406 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
407 }
408 const register double dot=(*V02)*(*V00)+(*V12)*(*V10)+(*V22)*(*V20);
409 *V02-=dot*(*V00);
410 *V12-=dot*(*V10);
411 *V22-=dot*(*V20);
412 normalizeVector3(V02,V12,V22);
413 *V01=(*V10)*(*V22)-(*V12)*(*V20);
414 *V11=(*V20)*(*V02)-(*V00)*(*V22);
415 *V21=(*V00)*(*V12)-(*V02)*(*V10);
416 normalizeVector3(V01,V11,V21);
417 }
418 }
419 }
420
421 // General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
422 // SM is the product of the last axis_offset entries in arg_0.getShape().
423 inline
424 void matrix_matrix_product(const int SL, const int SM, const int SR, const double* A, const double* B, double* C, int transpose)
425 {
426 if (transpose == 0) {
427 for (int i=0; i<SL; i++) {
428 for (int j=0; j<SR; j++) {
429 double sum = 0.0;
430 for (int l=0; l<SM; l++) {
431 sum += A[i+SL*l] * B[l+SM*j];
432 }
433 C[i+SL*j] = sum;
434 }
435 }
436 }
437 else if (transpose == 1) {
438 for (int i=0; i<SL; i++) {
439 for (int j=0; j<SR; j++) {
440 double sum = 0.0;
441 for (int l=0; l<SM; l++) {
442 sum += A[i*SM+l] * B[l+SM*j];
443 }
444 C[i+SL*j] = sum;
445 }
446 }
447 }
448 else if (transpose == 2) {
449 for (int i=0; i<SL; i++) {
450 for (int j=0; j<SR; j++) {
451 double sum = 0.0;
452 for (int l=0; l<SM; l++) {
453 sum += A[i+SL*l] * B[l*SR+j];
454 }
455 C[i+SL*j] = sum;
456 }
457 }
458 }
459 }
460
461 template <typename UnaryFunction>
462 inline void tensor_unary_operation(const int size,
463 const double *arg1,
464 double * argRes,
465 UnaryFunction operation)
466 {
467 for (int i = 0; i < size; ++i) {
468 argRes[i] = operation(arg1[i]);
469 }
470 return;
471 }
472
473 template <typename BinaryFunction>
474 inline void tensor_binary_operation(const int size,
475 const double *arg1,
476 const double *arg2,
477 double * argRes,
478 BinaryFunction operation)
479 {
480 for (int i = 0; i < size; ++i) {
481 argRes[i] = operation(arg1[i], arg2[i]);
482 }
483 return;
484 }
485
486 template <typename BinaryFunction>
487 inline void tensor_binary_operation(const int size,
488 double arg1,
489 const double *arg2,
490 double *argRes,
491 BinaryFunction operation)
492 {
493 for (int i = 0; i < size; ++i) {
494 argRes[i] = operation(arg1, arg2[i]);
495 }
496 return;
497 }
498
499 template <typename BinaryFunction>
500 inline void tensor_binary_operation(const int size,
501 const double *arg1,
502 double arg2,
503 double *argRes,
504 BinaryFunction operation)
505 {
506 for (int i = 0; i < size; ++i) {
507 argRes[i] = operation(arg1[i], arg2);
508 }
509 return;
510 }
511
512 } // end of namespace
513 #endif

  ViewVC Help
Powered by ViewVC 1.1.26