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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 gross 576 // $Id$
2     /*
3 elspeth 615 ************************************************************
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 gross 576 */
13    
14     #if !defined escript_LocalOps_H
15     #define escript_LocalOps_H
16 gross 580 #include <math.h>
17 gross 576 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 gross 580 *ev0=A00;
31 gross 576
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 gross 583 void eigenvalues2(const double A00,const double A01,const double A11,
45 gross 576 double* ev0, double* ev1) {
46 gross 580 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 gross 576 }
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 gross 580 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 gross 585 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 gross 580 }
100 gross 576 }
101 gross 583 /**
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 gross 587 register double absA10=fabs(A10);
134 gross 583 register double absA01=fabs(A01);
135     register double absA11=fabs(A11);
136 gross 587 register double m=absA11>absA10 ? absA11 : absA10;
137     if (absA00>m || absA01>m) {
138     *V0=-A01;
139 gross 583 *V1=A00;
140     } else {
141     if (m<=0) {
142     *V0=1.;
143     *V1=0.;
144     } else {
145     *V0=A11;
146 gross 587 *V1=-A10;
147 gross 583 }
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 gross 588 vectorInKernel2(A11-IA10*A01,A12-IA10*A02,
179     A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
180 gross 583 *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 gross 587 *V01= *V10;
226 gross 583 *V11=-(*V00);
227     } else {
228     *V01=-(*V10);
229 gross 587 *V11= (*V10);
230 gross 583 }
231     } else if (TEMP0>0.) {
232     *V00=TEMP0*scale;
233     *V10=TEMP1*scale;
234     if (TEMP1<0.) {
235 gross 587 *V01=-(*V10);
236 gross 583 *V11= (*V00);
237     } else {
238 gross 587 *V01= (*V10);
239     *V11=-(*V00);
240 gross 583 }
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 gross 588 *V02=1.;
334 gross 583 *V12=0.;
335 gross 588 *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 gross 583 *ev0=TEMP_EV0;
343     *ev1=TEMP_EV1;
344 gross 588 *ev2=A00;
345 gross 583 } else {
346 gross 588 *V01=1.;
347     *V11=0.;
348     *V21=0.;
349     *V00=0.;
350     *V10=TEMP_V00;
351 gross 583 *V20=TEMP_V10;
352 gross 588 *V02=0.;
353     *V12=TEMP_V01;
354 gross 583 *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 gross 576 } // end of namespace
412     #endif

  ViewVC Help
Powered by ViewVC 1.1.26