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

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

Legend:
 Removed from v.580 changed lines Added in v.583