/[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 1020 by phornby, Mon Mar 12 10:12:36 2007 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 __INTEL_COMPILER
17    #   include <mathimf.h>
18    #else
19    #   include <math.h>
20    #endif
21    #ifndef M_PI
22    #   define M_PI           3.14159265358979323846  /* pi */
23    #endif
24    
25  namespace escript {  namespace escript {
26    
# Line 29  namespace escript { Line 35  namespace escript {
35  inline  inline
36  void eigenvalues1(const double A00,double* ev0) {  void eigenvalues1(const double A00,double* ev0) {
37    
38     *ev0=1.;     *ev0=A00;
39    
40  }  }
41  /**  /**
# Line 43  void eigenvalues1(const double A00,doubl Line 49  void eigenvalues1(const double A00,doubl
49     \param ev1 Output - largest eigenvalue     \param ev1 Output - largest eigenvalue
50  */  */
51  inline  inline
52  void eigenvalues2(const double A00,const double A01  void eigenvalues2(const double A00,const double A01,const double A11,
                                  ,const double A11,  
53                   double* ev0, double* ev1) {                   double* ev0, double* ev1) {
54          const register double trA=(A00+A11)/2.;
55     *ev0=1.;        const register double A_00=A00-trA;
56     *ev1=2.;        const register double A_11=A11-trA;
57          const register double s=sqrt(A01*A01-A_00*A_11);
58          *ev0=trA-s;
59          *ev1=trA+s;
60  }  }
61  /**  /**
62     \brief     \brief
# Line 70  void eigenvalues3(const double A00, cons Line 78  void eigenvalues3(const double A00, cons
78                                                       const double A22,                                                       const double A22,
79                   double* ev0, double* ev1,double* ev2) {                   double* ev0, double* ev1,double* ev2) {
80    
81     *ev0=1.;        const register double trA=(A00+A11+A22)/3.;
82     *ev1=2.;        const register double A_00=A00-trA;
83     *ev2=3.;        const register double A_11=A11-trA;
84          const register double A_22=A22-trA;
85          const register double A01_2=A01*A01;
86          const register double A02_2=A02*A02;
87          const register double A12_2=A12*A12;
88          const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
89          if (p<=0.) {
90             *ev2=trA;
91             *ev1=trA;
92             *ev0=trA;
93    
94          } else {
95             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);
96             const register double sq_p=sqrt(p/3.);
97             register double z=-q/(2*pow(sq_p,3));
98             if (z<-1.) {
99                z=-1.;
100             } else if (z>1.) {
101                z=1.;
102             }
103             const register double alpha_3=acos(z)/3.;
104             *ev2=trA+2.*sq_p*cos(alpha_3);
105             *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.);
106             *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.);
107          }
108    }
109    /**
110       \brief
111       solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
112    
113       \param A00 Input - A_00
114       \param ev0 Output - eigenvalue
115       \param V00 Output - eigenvector
116       \param tol Input - tolerance to identify to eigenvalues
117    */
118    inline
119    void  eigenvalues_and_eigenvectors1(const double A00,double* ev0,double* V00,const double tol)
120    {
121          eigenvalues1(A00,ev0);
122          *V00=1.;
123          return;
124    }
125    /**
126       \brief
127       returns a non-zero vector in the kernel of [[A00,A01],[A01,A11]] assuming that the kernel dimension is at least 1.
128    
129       \param A00 Input - matrix component
130       \param A10 Input - matrix component
131       \param A01 Input - matrix component
132       \param A11 Input - matrix component
133       \param V0 Output - vector component
134       \param V1 Output - vector component
135    */
136    inline
137    void  vectorInKernel2(const double A00,const double A10,const double A01,const double A11,
138                          double* V0, double*V1)
139    {
140          register double absA00=fabs(A00);
141          register double absA10=fabs(A10);
142          register double absA01=fabs(A01);
143          register double absA11=fabs(A11);
144          register double m=absA11>absA10 ? absA11 : absA10;
145          if (absA00>m || absA01>m) {
146             *V0=-A01;
147             *V1=A00;
148          } else {
149             if (m<=0) {
150               *V0=1.;
151               *V1=0.;
152             } else {
153               *V0=A11;
154               *V1=-A10;
155             }
156         }
157    }
158    /**
159       \brief
160       returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]]
161       assuming that the kernel dimension is at least 1 and A00 is non zero.
162    
163       \param A00 Input - matrix component
164       \param A10 Input - matrix component
165       \param A20 Input - matrix component
166       \param A01 Input - matrix component
167       \param A11 Input - matrix component
168       \param A21 Input - matrix component
169       \param A02 Input - matrix component
170       \param A12 Input - matrix component
171       \param A22 Input - matrix component
172       \param V0 Output - vector component
173       \param V1 Output - vector component
174       \param V2 Output - vector component
175    */
176    inline
177    void  vectorInKernel3__nonZeroA00(const double A00,const double A10,const double A20,
178                                    const double A01,const double A11,const double A21,
179                                    const double A02,const double A12,const double A22,
180                                    double* V0,double* V1,double* V2)
181    {
182        double TEMP0,TEMP1;
183        register const double I00=1./A00;
184        register const double IA10=I00*A10;
185        register const double IA20=I00*A20;
186        vectorInKernel2(A11-IA10*A01,A12-IA10*A02,
187                        A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
188        *V0=-(A10*TEMP0+A20*TEMP1);
189        *V1=A00*TEMP0;
190        *V2=A00*TEMP1;
191    }
192          
193    /**
194       \brief
195       solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
196       ordered by increasing value and eigen vectors are normalizeVector3d such that
197       length is zero and first non-zero component is positive.
198    
199       \param A00 Input - A_00
200       \param A01 Input - A_01
201       \param A11 Input - A_11
202       \param ev0 Output - smallest eigenvalue
203       \param ev1 Output - eigenvalue
204       \param V00 Output - eigenvector componenent coresponding to ev0
205       \param V10 Output - eigenvector componenent coresponding to ev0
206       \param V01 Output - eigenvector componenent coresponding to ev1
207       \param V11 Output - eigenvector componenent coresponding to ev1
208       \param tol Input - tolerance to identify to eigenvalues
209    */
210    inline
211    void  eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11,
212                                        double* ev0, double* ev1,
213                                        double* V00, double* V10, double* V01, double* V11,
214                                        const double tol)
215    {
216         double TEMP0,TEMP1;
217         eigenvalues2(A00,A01,A11,ev0,ev1);
218         const register double absev0=fabs(*ev0);
219         const register double absev1=fabs(*ev1);
220         register double max_ev=absev0>absev1 ? absev0 : absev1;
221         if (fabs((*ev0)-(*ev1))<tol*max_ev) {
222            *V00=1.;
223            *V10=0.;
224            *V01=0.;
225            *V11=1.;
226         } else {
227            vectorInKernel2(A00-(*ev0),A01,A01,A11-(*ev0),&TEMP0,&TEMP1);
228            const register double scale=1./sqrt(TEMP0*TEMP0+TEMP1*TEMP1);
229            if (TEMP0<0.) {
230                *V00=-TEMP0*scale;
231                *V10=-TEMP1*scale;
232                if (TEMP1<0.) {
233                   *V01=  *V10;
234                   *V11=-(*V00);
235                } else {
236                   *V01=-(*V10);
237                   *V11= (*V10);
238                }
239            } else if (TEMP0>0.) {
240                *V00=TEMP0*scale;
241                *V10=TEMP1*scale;
242                if (TEMP1<0.) {
243                   *V01=-(*V10);
244                   *V11= (*V00);
245                } else {
246                   *V01= (*V10);
247                   *V11=-(*V00);
248                }
249            } else {
250               *V00=0.;
251               *V10=1;
252               *V11=0.;
253               *V01=1.;
254           }
255       }
256    }
257    /**
258       \brief
259       nomalizes a 3-d vector such that length is one and first non-zero component is positive.
260    
261       \param V0 - vector componenent
262       \param V1 - vector componenent
263       \param V2 - vector componenent
264    */
265    inline
266    void  normalizeVector3(double* V0,double* V1,double* V2)
267    {
268        register double s;
269        if (*V0>0) {
270            s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
271            *V0*=s;
272            *V1*=s;
273            *V2*=s;
274        } else if (*V0<0)  {
275            s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
276            *V0*=s;
277            *V1*=s;
278            *V2*=s;
279        } else {
280            if (*V1>0) {
281                s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
282                *V1*=s;
283                *V2*=s;
284            } else if (*V1<0)  {
285                s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
286                *V1*=s;
287                *V2*=s;
288            } else {
289                *V2=1.;
290            }
291        }
292    }
293    /**
294       \brief
295       solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
296       ordered by increasing value and eigen vectors are normalizeVector3d such that
297       length is zero and first non-zero component is positive.
298    
299       \param A00 Input - A_00
300       \param A01 Input - A_01
301       \param A11 Input - A_11
302       \param ev0 Output - smallest eigenvalue
303       \param ev1 Output - eigenvalue
304       \param V00 Output - eigenvector componenent coresponding to ev0
305       \param V10 Output - eigenvector componenent coresponding to ev0
306       \param V01 Output - eigenvector componenent coresponding to ev1
307       \param V11 Output - eigenvector componenent coresponding to ev1
308       \param tol Input - tolerance to identify to eigenvalues
309    */
310    inline
311    void  eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02,
312                                        const double A11, const double A12, const double A22,
313                                        double* ev0, double* ev1, double* ev2,
314                                        double* V00, double* V10, double* V20,
315                                        double* V01, double* V11, double* V21,
316                                        double* V02, double* V12, double* V22,
317                                        const double tol)
318    {
319          register const double absA01=fabs(A01);
320          register const double absA02=fabs(A02);
321          register const double m=absA01>absA02 ? absA01 : absA02;
322          if (m<=0) {
323            double TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1;
324            eigenvalues_and_eigenvectors2(A11,A12,A22,
325                                          &TEMP_EV0,&TEMP_EV1,
326                                          &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
327            if (A00<=TEMP_EV0) {
328                *V00=1.;
329                *V10=0.;
330                *V20=0.;
331                *V01=0.;
332                *V11=TEMP_V00;
333                *V21=TEMP_V10;
334                *V02=0.;
335                *V12=TEMP_V01;
336                *V22=TEMP_V11;
337                *ev0=A00;
338                *ev1=TEMP_EV0;
339                *ev2=TEMP_EV1;
340            } else if (A00>TEMP_EV1) {
341                *V02=1.;
342                *V12=0.;
343                *V22=0.;
344                *V00=0.;
345                *V10=TEMP_V00;
346                *V20=TEMP_V10;
347                *V01=0.;
348                *V11=TEMP_V01;
349                *V21=TEMP_V11;
350                *ev0=TEMP_EV0;
351                *ev1=TEMP_EV1;
352                *ev2=A00;
353            } else {
354                *V01=1.;
355                *V11=0.;
356                *V21=0.;
357                *V00=0.;
358                *V10=TEMP_V00;
359                *V20=TEMP_V10;
360                *V02=0.;
361                *V12=TEMP_V01;
362                *V22=TEMP_V11;
363                *ev0=TEMP_EV0;
364                *ev1=A00;
365                *ev2=TEMP_EV1;
366            }
367          } else {
368             eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2);
369             const register double absev0=fabs(*ev0);
370             const register double absev1=fabs(*ev1);
371             const register double absev2=fabs(*ev2);
372             register double max_ev=absev0>absev1 ? absev0 : absev1;
373             max_ev=max_ev>absev2 ? max_ev : absev2;
374             const register double d_01=fabs((*ev0)-(*ev1));
375             const register double d_12=fabs((*ev1)-(*ev2));
376             const register double max_d=d_01>d_12 ? d_01 : d_12;
377             if (max_d<=tol*max_ev) {
378                 *V00=1.;
379                 *V10=0;
380                 *V20=0;
381                 *V01=0;
382                 *V11=1.;
383                 *V21=0;
384                 *V02=0;
385                 *V12=0;
386                 *V22=1.;
387             } else {
388                const register double S00=A00-(*ev0);
389                const register double absS00=fabs(S00);
390                if (fabs(S00)>m) {
391                    vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);
392                } else if (absA02<m) {
393                    vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);
394                } else {
395                    vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
396                }
397                normalizeVector3(V00,V10,V20);;
398                const register double T00=A00-(*ev2);
399                const register double absT00=fabs(T00);
400                if (fabs(T00)>m) {
401                     vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);
402                } else if (absA02<m) {
403                     vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);
404                } else {
405                     vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
406                }
407                const register double dot=(*V02)*(*V00)+(*V12)*(*V10)+(*V22)*(*V20);
408                *V02-=dot*(*V00);
409                *V12-=dot*(*V10);
410                *V22-=dot*(*V20);
411                normalizeVector3(V02,V12,V22);
412                *V01=(*V10)*(*V22)-(*V12)*(*V20);
413                *V11=(*V20)*(*V02)-(*V00)*(*V22);
414                *V21=(*V00)*(*V12)-(*V02)*(*V10);
415                normalizeVector3(V01,V11,V21);
416             }
417       }
418    }
419    
420    // General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
421    // SM is the product of the last axis_offset entries in arg_0.getShape().
422    inline
423    void matrix_matrix_product(const int SL, const int SM, const int SR, const double* A, const double* B, double* C, int transpose)
424    {
425      if (transpose == 0) {
426        for (int i=0; i<SL; i++) {
427          for (int j=0; j<SR; j++) {
428            double sum = 0.0;
429            for (int l=0; l<SM; l++) {
430          sum += A[i+SL*l] * B[l+SM*j];
431            }
432            C[i+SL*j] = sum;
433          }
434        }
435      }
436      else if (transpose == 1) {
437        for (int i=0; i<SL; i++) {
438          for (int j=0; j<SR; j++) {
439            double sum = 0.0;
440            for (int l=0; l<SM; l++) {
441          sum += A[i*SM+l] * B[l+SM*j];
442            }
443            C[i+SL*j] = sum;
444          }
445        }
446      }
447      else if (transpose == 2) {
448        for (int i=0; i<SL; i++) {
449          for (int j=0; j<SR; j++) {
450            double sum = 0.0;
451            for (int l=0; l<SM; l++) {
452          sum += A[i+SL*l] * B[l*SR+j];
453            }
454            C[i+SL*j] = sum;
455          }
456        }
457      }
458  }  }
459    
460  } // end of namespace  } // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26