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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 576 by gross, Fri Mar 3 08:28:42 2006 UTC revision 587 by gross, Fri Mar 10 02:26:50 2006 UTC
# Line 15  Line 15 
15    
16  #if !defined escript_LocalOps_H  #if !defined escript_LocalOps_H
17  #define escript_LocalOps_H  #define escript_LocalOps_H
18    #include <math.h>
19  namespace escript {  namespace escript {
20    
21    
# Line 29  namespace escript { Line 29  namespace escript {
29  inline  inline
30  void eigenvalues1(const double A00,double* ev0) {  void eigenvalues1(const double A00,double* ev0) {
31    
32     *ev0=1.;     *ev0=A00;
33    
34  }  }
35  /**  /**
# Line 43  void eigenvalues1(const double A00,doubl Line 43  void eigenvalues1(const double A00,doubl
43     \param ev1 Output - largest eigenvalue     \param ev1 Output - largest eigenvalue
44  */  */
45  inline  inline
46  void eigenvalues2(const double A00,const double A01  void eigenvalues2(const double A00,const double A01,const double A11,
                                  ,const double A11,  
47                   double* ev0, double* ev1) {                   double* ev0, double* ev1) {
48          const register double trA=(A00+A11)/2.;
49     *ev0=1.;        const register double A_00=A00-trA;
50     *ev1=2.;        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     \brief
# Line 70  void eigenvalues3(const double A00, cons Line 72  void eigenvalues3(const double A00, cons
72                                                       const double A22,                                                       const double A22,
73                   double* ev0, double* ev1,double* ev2) {                   double* ev0, double* ev1,double* ev2) {
74    
75     *ev0=1.;        const register double trA=(A00+A11+A22)/3.;
76     *ev1=2.;        const register double A_00=A00-trA;
77     *ev2=3.;        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  } // end of namespace
413  #endif  #endif

Legend:
Removed from v.576  
changed lines
  Added in v.587

  ViewVC Help
Powered by ViewVC 1.1.26