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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 588 - (hide annotations)
Fri Mar 10 04:45:04 2006 UTC (13 years, 6 months ago) by gross
File MIME type: text/plain
File size: 14044 byte(s)
1D and 3D tests for eigenvalues_and_eigenvector added.
1 gross 576 // $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 gross 580 #include <math.h>
19 gross 576 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 gross 580 *ev0=A00;
33 gross 576
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 gross 583 void eigenvalues2(const double A00,const double A01,const double A11,
47 gross 576 double* ev0, double* ev1) {
48 gross 580 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 gross 576 }
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 gross 580 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 gross 585 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 gross 580 }
102 gross 576 }
103 gross 583 /**
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 gross 587 register double absA10=fabs(A10);
136 gross 583 register double absA01=fabs(A01);
137     register double absA11=fabs(A11);
138 gross 587 register double m=absA11>absA10 ? absA11 : absA10;
139     if (absA00>m || absA01>m) {
140     *V0=-A01;
141 gross 583 *V1=A00;
142     } else {
143     if (m<=0) {
144     *V0=1.;
145     *V1=0.;
146     } else {
147     *V0=A11;
148 gross 587 *V1=-A10;
149 gross 583 }
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 gross 588 vectorInKernel2(A11-IA10*A01,A12-IA10*A02,
181     A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
182 gross 583 *V0=-(A10*TEMP0+A20*TEMP1);
183     *V1=A00*TEMP0;
184     *V2=A00*TEMP1;
185     }
186    
187     /**
188     \brief
189     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
190     ordered by increasing value and eigen vectors are normalizeVector3d such that
191     length is zero and first non-zero component is positive.
192    
193     \param A00 Input - A_00
194     \param A01 Input - A_01
195     \param A11 Input - A_11
196     \param ev0 Output - smallest eigenvalue
197     \param ev1 Output - eigenvalue
198     \param V00 Output - eigenvector componenent coresponding to ev0
199     \param V10 Output - eigenvector componenent coresponding to ev0
200     \param V01 Output - eigenvector componenent coresponding to ev1
201     \param V11 Output - eigenvector componenent coresponding to ev1
202     \param tol Input - tolerance to identify to eigenvalues
203     */
204     inline
205     void eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11,
206     double* ev0, double* ev1,
207     double* V00, double* V10, double* V01, double* V11,
208     const double tol)
209     {
210     double TEMP0,TEMP1;
211     eigenvalues2(A00,A01,A11,ev0,ev1);
212     const register double absev0=fabs(*ev0);
213     const register double absev1=fabs(*ev1);
214     register double max_ev=absev0>absev1 ? absev0 : absev1;
215     if (fabs((*ev0)-(*ev1))<tol*max_ev) {
216     *V00=1.;
217     *V10=0.;
218     *V01=0.;
219     *V11=1.;
220     } else {
221     vectorInKernel2(A00-(*ev0),A01,A01,A11-(*ev0),&TEMP0,&TEMP1);
222     const register double scale=1./sqrt(TEMP0*TEMP0+TEMP1*TEMP1);
223     if (TEMP0<0.) {
224     *V00=-TEMP0*scale;
225     *V10=-TEMP1*scale;
226     if (TEMP1<0.) {
227 gross 587 *V01= *V10;
228 gross 583 *V11=-(*V00);
229     } else {
230     *V01=-(*V10);
231 gross 587 *V11= (*V10);
232 gross 583 }
233     } else if (TEMP0>0.) {
234     *V00=TEMP0*scale;
235     *V10=TEMP1*scale;
236     if (TEMP1<0.) {
237 gross 587 *V01=-(*V10);
238 gross 583 *V11= (*V00);
239     } else {
240 gross 587 *V01= (*V10);
241     *V11=-(*V00);
242 gross 583 }
243     } else {
244     *V00=0.;
245     *V10=1;
246     *V11=0.;
247     *V01=1.;
248     }
249     }
250     }
251     /**
252     \brief
253     nomalizes a 3-d vector such that length is one and first non-zero component is positive.
254    
255     \param V0 - vector componenent
256     \param V1 - vector componenent
257     \param V2 - vector componenent
258     */
259     inline
260     void normalizeVector3(double* V0,double* V1,double* V2)
261     {
262     register double s;
263     if (*V0>0) {
264     s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
265     *V0*=s;
266     *V1*=s;
267     *V2*=s;
268     } else if (*V0<0) {
269     s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
270     *V0*=s;
271     *V1*=s;
272     *V2*=s;
273     } else {
274     if (*V1>0) {
275     s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
276     *V1*=s;
277     *V2*=s;
278     } else if (*V1<0) {
279     s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
280     *V1*=s;
281     *V2*=s;
282     } else {
283     *V2=1.;
284     }
285     }
286     }
287     /**
288     \brief
289     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
290     ordered by increasing value and eigen vectors are normalizeVector3d such that
291     length is zero and first non-zero component is positive.
292    
293     \param A00 Input - A_00
294     \param A01 Input - A_01
295     \param A11 Input - A_11
296     \param ev0 Output - smallest eigenvalue
297     \param ev1 Output - eigenvalue
298     \param V00 Output - eigenvector componenent coresponding to ev0
299     \param V10 Output - eigenvector componenent coresponding to ev0
300     \param V01 Output - eigenvector componenent coresponding to ev1
301     \param V11 Output - eigenvector componenent coresponding to ev1
302     \param tol Input - tolerance to identify to eigenvalues
303     */
304     inline
305     void eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02,
306     const double A11, const double A12, const double A22,
307     double* ev0, double* ev1, double* ev2,
308     double* V00, double* V10, double* V20,
309     double* V01, double* V11, double* V21,
310     double* V02, double* V12, double* V22,
311     const double tol)
312     {
313     register const double absA01=fabs(A01);
314     register const double absA02=fabs(A02);
315     register const double m=absA01>absA02 ? absA01 : absA02;
316     if (m<=0) {
317     double TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1;
318     eigenvalues_and_eigenvectors2(A11,A12,A22,
319     &TEMP_EV0,&TEMP_EV1,
320     &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
321     if (A00<=TEMP_EV0) {
322     *V00=1.;
323     *V10=0.;
324     *V20=0.;
325     *V01=0.;
326     *V11=TEMP_V00;
327     *V21=TEMP_V10;
328     *V02=0.;
329     *V12=TEMP_V01;
330     *V22=TEMP_V11;
331     *ev0=A00;
332     *ev1=TEMP_EV0;
333     *ev2=TEMP_EV1;
334     } else if (A00>TEMP_EV1) {
335 gross 588 *V02=1.;
336 gross 583 *V12=0.;
337 gross 588 *V22=0.;
338     *V00=0.;
339     *V10=TEMP_V00;
340     *V20=TEMP_V10;
341     *V01=0.;
342     *V11=TEMP_V01;
343     *V21=TEMP_V11;
344 gross 583 *ev0=TEMP_EV0;
345     *ev1=TEMP_EV1;
346 gross 588 *ev2=A00;
347 gross 583 } else {
348 gross 588 *V01=1.;
349     *V11=0.;
350     *V21=0.;
351     *V00=0.;
352     *V10=TEMP_V00;
353 gross 583 *V20=TEMP_V10;
354 gross 588 *V02=0.;
355     *V12=TEMP_V01;
356 gross 583 *V22=TEMP_V11;
357     *ev0=TEMP_EV0;
358     *ev1=A00;
359     *ev2=TEMP_EV1;
360     }
361     } else {
362     eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2);
363     const register double absev0=fabs(*ev0);
364     const register double absev1=fabs(*ev1);
365     const register double absev2=fabs(*ev2);
366     register double max_ev=absev0>absev1 ? absev0 : absev1;
367     max_ev=max_ev>absev2 ? max_ev : absev2;
368     const register double d_01=fabs((*ev0)-(*ev1));
369     const register double d_12=fabs((*ev1)-(*ev2));
370     const register double max_d=d_01>d_12 ? d_01 : d_12;
371     if (max_d<=tol*max_ev) {
372     *V00=1.;
373     *V10=0;
374     *V20=0;
375     *V01=0;
376     *V11=1.;
377     *V21=0;
378     *V02=0;
379     *V12=0;
380     *V22=1.;
381     } else {
382     const register double S00=A00-(*ev0);
383     const register double absS00=fabs(S00);
384     if (fabs(S00)>m) {
385     vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);
386     } else if (absA02<m) {
387     vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);
388     } else {
389     vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
390     }
391     normalizeVector3(V00,V10,V20);;
392     const register double T00=A00-(*ev2);
393     const register double absT00=fabs(T00);
394     if (fabs(T00)>m) {
395     vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);
396     } else if (absA02<m) {
397     vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);
398     } else {
399     vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
400     }
401     const register double dot=(*V02)*(*V00)+(*V12)*(*V10)+(*V22)*(*V20);
402     *V02-=dot*(*V00);
403     *V12-=dot*(*V10);
404     *V22-=dot*(*V20);
405     normalizeVector3(V02,V12,V22);
406     *V01=(*V10)*(*V22)-(*V12)*(*V20);
407     *V11=(*V20)*(*V02)-(*V00)*(*V22);
408     *V21=(*V00)*(*V12)-(*V02)*(*V10);
409     normalizeVector3(V01,V11,V21);
410     }
411     }
412     }
413 gross 576 } // end of namespace
414     #endif

  ViewVC Help
Powered by ViewVC 1.1.26