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

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

  ViewVC Help
Powered by ViewVC 1.1.26