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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 615 - (show annotations)
Wed Mar 22 02:12:00 2006 UTC (13 years, 6 months ago) by elspeth
File MIME type: text/plain
File size: 13722 byte(s)
More copyright information.

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

  ViewVC Help
Powered by ViewVC 1.1.26