/[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 583 by gross, Wed Mar 8 08:15:34 2006 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
1  // $Id$  
2  /*  /*******************************************************
3   ******************************************************************************  *
4   *                                                                            *  * Copyright (c) 2003-2009 by University of Queensland
5   *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *  * Earth Systems Science Computational Center (ESSCC)
6   *                                                                            *  * http://www.uq.edu.au/esscc
7   * This software is the property of ACcESS. No part of this code              *  *
8   * may be copied in any form or by any means without the expressed written    *  * Primary Business: Queensland, Australia
9   * consent of ACcESS.  Copying, use or modification of this software          *  * Licensed under the Open Software License version 3.0
10   * by any unauthorised person is illegal unless that person has a software    *  * http://www.opensource.org/licenses/osl-3.0.php
11   * license agreement with ACcESS.                                             *  *
12   *                                                                            *  *******************************************************/
13   ******************************************************************************  
 */  
14    
15  #if !defined escript_LocalOps_H  #if !defined escript_LocalOps_H
16  #define escript_LocalOps_H  #define escript_LocalOps_H
17  #include <math.h>  #if defined(_WIN32) && defined(__INTEL_COMPILER)
18    #   include <mathimf.h>
19    #else
20    #   include <math.h>
21    #endif
22    #ifndef M_PI
23    #   define M_PI           3.14159265358979323846  /* pi */
24    #endif
25    
26    
27    /**
28    \file LocalOps.h
29    \brief Describes binary operations performed on double*.
30    
31    For operations on DataAbstract see BinaryOp.h.
32    For operations on DataVector see DataMaths.h.
33    */
34    
35  namespace escript {  namespace escript {
36    
37    
# Line 23  namespace escript { Line 39  namespace escript {
39     \brief     \brief
40     solves a 1x1 eigenvalue A*V=ev*V problem     solves a 1x1 eigenvalue A*V=ev*V problem
41    
42     \param A00 Input - A_00     \param A00 Input - A_00
43     \param ev0 Output - eigenvalue     \param ev0 Output - eigenvalue
44  */  */
45  inline  inline
# Line 36  void eigenvalues1(const double A00,doubl Line 52  void eigenvalues1(const double A00,doubl
52     \brief     \brief
53     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
54    
55     \param A00 Input - A_00     \param A00 Input - A_00
56     \param A01 Input - A_01     \param A01 Input - A_01
57     \param A11 Input - A_11     \param A11 Input - A_11
58     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
59     \param ev1 Output - largest eigenvalue     \param ev1 Output - largest eigenvalue
# Line 56  void eigenvalues2(const double A00,const Line 72  void eigenvalues2(const double A00,const
72     \brief     \brief
73     solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A     solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
74    
75     \param A00 Input - A_00     \param A00 Input - A_00
76     \param A01 Input - A_01     \param A01 Input - A_01
77     \param A02 Input - A_02     \param A02 Input - A_02
78     \param A11 Input - A_11     \param A11 Input - A_11
79     \param A12 Input - A_12     \param A12 Input - A_12
80     \param A22 Input - A_22     \param A22 Input - A_22
81     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
82     \param ev1 Output - eigenvalue     \param ev1 Output - eigenvalue
83     \param ev2 Output - largest eigenvalue     \param ev2 Output - largest eigenvalue
# Line 80  void eigenvalues3(const double A00, cons Line 96  void eigenvalues3(const double A00, cons
96        const register double A02_2=A02*A02;        const register double A02_2=A02*A02;
97        const register double A12_2=A12*A12;        const register double A12_2=A12*A12;
98        const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;        const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
99        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);        if (p<=0.) {
100        const register double sq_p=sqrt(p/3.);           *ev2=trA;
101        register double z=-q/(2*pow(sq_p,3));           *ev1=trA;
102        if (z<-1.) {           *ev0=trA;
103           z=-1.;  
104        } else if (z>1.) {        } else {
105           z=1.;           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);
106             const register double sq_p=sqrt(p/3.);
107             register double z=-q/(2*pow(sq_p,3));
108             if (z<-1.) {
109                z=-1.;
110             } else if (z>1.) {
111                z=1.;
112             }
113             const register double alpha_3=acos(z)/3.;
114             *ev2=trA+2.*sq_p*cos(alpha_3);
115             *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.);
116             *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.);
117        }        }
       const register double alpha_3=acos(z)/3.;  
       *ev2=trA+2.*sq_p*cos(alpha_3);  
       *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.);  
       *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.);  
118  }  }
119  /**  /**
120     \brief     \brief
121     solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A     solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
122    
123     \param A00 Input - A_00     \param A00 Input - A_00
124     \param ev0 Output - eigenvalue     \param ev0 Output - eigenvalue
125     \param V00 Output - eigenvector     \param V00 Output - eigenvector
126     \param tol Input - tolerance to identify to eigenvalues     \param tol Input - tolerance to identify to eigenvalues
# Line 125  void  vectorInKernel2(const double A00,c Line 148  void  vectorInKernel2(const double A00,c
148                        double* V0, double*V1)                        double* V0, double*V1)
149  {  {
150        register double absA00=fabs(A00);        register double absA00=fabs(A00);
       register double absA01=fabs(A01);  
151        register double absA10=fabs(A10);        register double absA10=fabs(A10);
152          register double absA01=fabs(A01);
153        register double absA11=fabs(A11);        register double absA11=fabs(A11);
154        register double m=absA11>absA01 ? absA11 : absA01;        register double m=absA11>absA10 ? absA11 : absA10;
155        if (absA00>m || absA10>m) {        if (absA00>m || absA01>m) {
156           *V0=-A10;           *V0=-A01;
157           *V1=A00;           *V1=A00;
158        } else {        } else {
159           if (m<=0) {           if (m<=0) {
# Line 138  void  vectorInKernel2(const double A00,c Line 161  void  vectorInKernel2(const double A00,c
161             *V1=0.;             *V1=0.;
162           } else {           } else {
163             *V0=A11;             *V0=A11;
164             *V1=-A01;             *V1=-A10;
165           }           }
166       }       }
167  }  }
168  /**  /**
169     \brief     \brief
170     returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]]     returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]]
171     assuming that the kernel dimension is at least 1 and A00 is non zero.     assuming that the kernel dimension is at least 1 and A00 is non zero.
172    
173     \param A00 Input - matrix component     \param A00 Input - matrix component
# Line 170  void  vectorInKernel3__nonZeroA00(const Line 193  void  vectorInKernel3__nonZeroA00(const
193      register const double I00=1./A00;      register const double I00=1./A00;
194      register const double IA10=I00*A10;      register const double IA10=I00*A10;
195      register const double IA20=I00*A20;      register const double IA20=I00*A20;
196      vectorInKernel2(A11-IA10*A01,A21-IA20*A01,A12-IA10*A02,A22-IA20*A02,&TEMP0,&TEMP1);      vectorInKernel2(A11-IA10*A01,A12-IA10*A02,
197                        A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
198      *V0=-(A10*TEMP0+A20*TEMP1);      *V0=-(A10*TEMP0+A20*TEMP1);
199      *V1=A00*TEMP0;      *V1=A00*TEMP0;
200      *V2=A00*TEMP1;      *V2=A00*TEMP1;
201  }  }
202          
203  /**  /**
204     \brief     \brief
205     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
206     ordered by increasing value and eigen vectors are normalizeVector3d such that     ordered by increasing value and eigen vectors are normalizeVector3d such that
207     length is zero and first non-zero component is positive.     length is zero and first non-zero component is positive.
208    
209     \param A00 Input - A_00     \param A00 Input - A_00
210     \param A01 Input - A_01     \param A01 Input - A_01
211     \param A11 Input - A_11     \param A11 Input - A_11
212     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
213     \param ev1 Output - eigenvalue     \param ev1 Output - eigenvalue
214     \param V00 Output - eigenvector componenent coresponding to ev0     \param V00 Output - eigenvector componenent coresponding to ev0
# Line 197  inline Line 221  inline
221  void  eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11,  void  eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11,
222                                      double* ev0, double* ev1,                                      double* ev0, double* ev1,
223                                      double* V00, double* V10, double* V01, double* V11,                                      double* V00, double* V10, double* V01, double* V11,
224                                      const double tol)                                      const double tol)
225  {  {
226       double TEMP0,TEMP1;       double TEMP0,TEMP1;
227       eigenvalues2(A00,A01,A11,ev0,ev1);       eigenvalues2(A00,A01,A11,ev0,ev1);
# Line 220  void  eigenvalues_and_eigenvectors2(cons Line 244  void  eigenvalues_and_eigenvectors2(cons
244                 *V11=-(*V00);                 *V11=-(*V00);
245              } else {              } else {
246                 *V01=-(*V10);                 *V01=-(*V10);
247                 *V11= (*V00);                 *V11= (*V10);
248              }              }
249          } else if (TEMP0>0.) {          } else if (TEMP0>0.) {
250              *V00=TEMP0*scale;              *V00=TEMP0*scale;
# Line 237  void  eigenvalues_and_eigenvectors2(cons Line 261  void  eigenvalues_and_eigenvectors2(cons
261             *V10=1;             *V10=1;
262             *V11=0.;             *V11=0.;
263             *V01=1.;             *V01=1.;
264         }         }
265     }     }
266  }  }
267  /**  /**
268     \brief     \brief
269     nomalizes a 3-d vector such that length is one and first non-zero component is positive.     nomalizes a 3-d vector such that length is one and first non-zero component is positive.
270    
271     \param V0 - vector componenent     \param V0 - vector componenent
272     \param V1 - vector componenent     \param V1 - vector componenent
273     \param V2 - vector componenent     \param V2 - vector componenent
274  */  */
# Line 278  void  normalizeVector3(double* V0,double Line 302  void  normalizeVector3(double* V0,double
302  }  }
303  /**  /**
304     \brief     \brief
305     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
306     ordered by increasing value and eigen vectors are normalizeVector3d such that     ordered by increasing value and eigen vectors are normalizeVector3d such that
307     length is zero and first non-zero component is positive.     length is zero and first non-zero component is positive.
308    
309     \param A00 Input - A_00     \param A00 Input - A_00
310     \param A01 Input - A_01     \param A01 Input - A_01
311     \param A11 Input - A_11     \param A02 Input - A_02
312       \param A11 Input - A_11
313       \param A12 Input - A_12
314       \param A22 Input - A_22
315     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
316     \param ev1 Output - eigenvalue     \param ev1 Output - eigenvalue
317       \param ev2 Output -
318     \param V00 Output - eigenvector componenent coresponding to ev0     \param V00 Output - eigenvector componenent coresponding to ev0
319     \param V10 Output - eigenvector componenent coresponding to ev0     \param V10 Output - eigenvector componenent coresponding to ev0
320       \param V20 Output -
321     \param V01 Output - eigenvector componenent coresponding to ev1     \param V01 Output - eigenvector componenent coresponding to ev1
322     \param V11 Output - eigenvector componenent coresponding to ev1     \param V11 Output - eigenvector componenent coresponding to ev1
323       \param V21 Output -
324       \param V02 Output -
325       \param V12 Output -
326       \param V22 Output -
327     \param tol Input - tolerance to identify to eigenvalues     \param tol Input - tolerance to identify to eigenvalues
328  */  */
329  inline  inline
330  void  eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02,  void  eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02,
331                                      const double A11, const double A12, const double A22,                                      const double A11, const double A12, const double A22,
332                                      double* ev0, double* ev1, double* ev2,                                      double* ev0, double* ev1, double* ev2,
333                                      double* V00, double* V10, double* V20,                                      double* V00, double* V10, double* V20,
334                                      double* V01, double* V11, double* V21,                                      double* V01, double* V11, double* V21,
335                                      double* V02, double* V12, double* V22,                                      double* V02, double* V12, double* V22,
336                                      const double tol)                                      const double tol)
337  {  {
338        register const double absA01=fabs(A01);        register const double absA01=fabs(A01);
# Line 324  void  eigenvalues_and_eigenvectors3(cons Line 357  void  eigenvalues_and_eigenvectors3(cons
357              *ev1=TEMP_EV0;              *ev1=TEMP_EV0;
358              *ev2=TEMP_EV1;              *ev2=TEMP_EV1;
359          } else if (A00>TEMP_EV1) {          } else if (A00>TEMP_EV1) {
360              *V00=TEMP_V00;              *V02=1.;
             *V10=TEMP_V10;  
             *V20=0.;  
             *V01=TEMP_V01;  
             *V11=TEMP_V11;  
             *V21=0.;  
             *V02=0.;  
361              *V12=0.;              *V12=0.;
362              *V22=1.;              *V22=0.;
363                *V00=0.;
364                *V10=TEMP_V00;
365                *V20=TEMP_V10;
366                *V01=0.;
367                *V11=TEMP_V01;
368                *V21=TEMP_V11;
369              *ev0=TEMP_EV0;              *ev0=TEMP_EV0;
370              *ev1=TEMP_EV1;              *ev1=TEMP_EV1;
371              *ev0=A00;              *ev2=A00;
372          } else {          } else {
373              *V00=TEMP_V00;              *V01=1.;
374              *V10=0;              *V11=0.;
             *V20=TEMP_V10;  
             *V01=0.;  
             *V11=1.;  
375              *V21=0.;              *V21=0.;
376              *V02=TEMP_V01;              *V00=0.;
377              *V12=0.;              *V10=TEMP_V00;
378                *V20=TEMP_V10;
379                *V02=0.;
380                *V12=TEMP_V01;
381              *V22=TEMP_V11;              *V22=TEMP_V11;
382              *ev0=TEMP_EV0;              *ev0=TEMP_EV0;
383              *ev1=A00;              *ev1=A00;
# Line 373  void  eigenvalues_and_eigenvectors3(cons Line 406  void  eigenvalues_and_eigenvectors3(cons
406           } else {           } else {
407              const register double S00=A00-(*ev0);              const register double S00=A00-(*ev0);
408              const register double absS00=fabs(S00);              const register double absS00=fabs(S00);
409              if (fabs(S00)>m) {              if (absS00>m) {
410                  vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);                  vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);
411              } else if (absA02<m) {              } else if (absA02<m) {
412                  vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);                  vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);
# Line 383  void  eigenvalues_and_eigenvectors3(cons Line 416  void  eigenvalues_and_eigenvectors3(cons
416              normalizeVector3(V00,V10,V20);;              normalizeVector3(V00,V10,V20);;
417              const register double T00=A00-(*ev2);              const register double T00=A00-(*ev2);
418              const register double absT00=fabs(T00);              const register double absT00=fabs(T00);
419              if (fabs(T00)>m) {              if (absT00>m) {
420                   vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);                   vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);
421              } else if (absA02<m) {              } else if (absA02<m) {
422                   vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);                   vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);
# Line 402  void  eigenvalues_and_eigenvectors3(cons Line 435  void  eigenvalues_and_eigenvectors3(cons
435           }           }
436     }     }
437  }  }
438    
439    // General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
440    // SM is the product of the last axis_offset entries in arg_0.getShape().
441    inline
442    void matrix_matrix_product(const int SL, const int SM, const int SR, const double* A, const double* B, double* C, int transpose)
443    {
444      if (transpose == 0) {
445        for (int i=0; i<SL; i++) {
446          for (int j=0; j<SR; j++) {
447            double sum = 0.0;
448            for (int l=0; l<SM; l++) {
449          sum += A[i+SL*l] * B[l+SM*j];
450            }
451            C[i+SL*j] = sum;
452          }
453        }
454      }
455      else if (transpose == 1) {
456        for (int i=0; i<SL; i++) {
457          for (int j=0; j<SR; j++) {
458            double sum = 0.0;
459            for (int l=0; l<SM; l++) {
460          sum += A[i*SM+l] * B[l+SM*j];
461            }
462            C[i+SL*j] = sum;
463          }
464        }
465      }
466      else if (transpose == 2) {
467        for (int i=0; i<SL; i++) {
468          for (int j=0; j<SR; j++) {
469            double sum = 0.0;
470            for (int l=0; l<SM; l++) {
471          sum += A[i+SL*l] * B[l*SR+j];
472            }
473            C[i+SL*j] = sum;
474          }
475        }
476      }
477    }
478    
479    template <typename UnaryFunction>
480    inline void tensor_unary_operation(const int size,
481                     const double *arg1,
482                     double * argRes,
483                     UnaryFunction operation)
484    {
485      for (int i = 0; i < size; ++i) {
486        argRes[i] = operation(arg1[i]);
487      }
488      return;
489    }
490    
491    template <typename BinaryFunction>
492    inline void tensor_binary_operation(const int size,
493                     const double *arg1,
494                     const double *arg2,
495                     double * argRes,
496                     BinaryFunction operation)
497    {
498      for (int i = 0; i < size; ++i) {
499        argRes[i] = operation(arg1[i], arg2[i]);
500      }
501      return;
502    }
503    
504    template <typename BinaryFunction>
505    inline void tensor_binary_operation(const int size,
506                     double arg1,
507                     const double *arg2,
508                     double *argRes,
509                     BinaryFunction operation)
510    {
511      for (int i = 0; i < size; ++i) {
512        argRes[i] = operation(arg1, arg2[i]);
513      }
514      return;
515    }
516    
517    template <typename BinaryFunction>
518    inline void tensor_binary_operation(const int size,
519                     const double *arg1,
520                     double arg2,
521                     double *argRes,
522                     BinaryFunction operation)
523    {
524      for (int i = 0; i < size; ++i) {
525        argRes[i] = operation(arg1[i], arg2);
526      }
527      return;
528    }
529    
530  } // end of namespace  } // end of namespace
531  #endif  #endif

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

  ViewVC Help
Powered by ViewVC 1.1.26