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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1020 - (show annotations)
Mon Mar 12 10:12:36 2007 UTC (12 years, 6 months ago) by phornby
File MIME type: text/plain
File size: 14873 byte(s)
Added explicit destructors to all Exception classes.

Fixed an ifdef in TestCase.cpp

Made the conditional definition of M_PI in LocalOps.h
depend only on M_PI being undefined.

Replace dynamically dimensioned arrays in DataFactory & DataTagged with malloc.

sort() method of list does not take a named argument
(despite the manual's claims to the contary).


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

  ViewVC Help
Powered by ViewVC 1.1.26