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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26