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

revision 757 by woo409, Mon Jun 26 13:12:56 2006 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1
1  // \$Id\$
2  /*  /*******************************************************
3   ************************************************************  *
4   *          Copyright 2006 by ACcESS MNRF                   *  * Copyright (c) 2003-2009 by University of Queensland
5   *                                                          *  * Earth Systems Science Computational Center (ESSCC)
6   *              http://www.access.edu.au                    *  * http://www.uq.edu.au/esscc
7   *       Primary Business: Queensland, Australia            *  *
11   ************************************************************  *
12  */  *******************************************************/
13
14
15  #if !defined escript_LocalOps_H  #if !defined escript_LocalOps_H
16  #define escript_LocalOps_H  #define escript_LocalOps_H
17  #ifdef _WIN32 && __INTEL_COMPILER  #if defined(_WIN32) && defined(__INTEL_COMPILER)
18  #include <mathimf.h>  #   include <mathimf.h>
#define M_PI 3.141592653589
19  #else  #else
20  #include <math.h>  #   include <math.h>
21  #endif  #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
# Line 27  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 40  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 60  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 108  void eigenvalues3(const double A00, cons Line 120  void eigenvalues3(const double A00, cons
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 155  void  vectorInKernel2(const double A00,c Line 167  void  vectorInKernel2(const double A00,c
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 187  void  vectorInKernel3__nonZeroA00(const Line 199  void  vectorInKernel3__nonZeroA00(const
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 209  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 228  void  eigenvalues_and_eigenvectors2(cons Line 240  void  eigenvalues_and_eigenvectors2(cons
240              *V00=-TEMP0*scale;              *V00=-TEMP0*scale;
241              *V10=-TEMP1*scale;              *V10=-TEMP1*scale;
242              if (TEMP1<0.) {              if (TEMP1<0.) {
243                 *V01=  *V10;                 *V01=  *V10;
244                 *V11=-(*V00);                 *V11=-(*V00);
245              } else {              } else {
246                 *V01=-(*V10);                 *V01=-(*V10);
# Line 238  void  eigenvalues_and_eigenvectors2(cons Line 250  void  eigenvalues_and_eigenvectors2(cons
250              *V00=TEMP0*scale;              *V00=TEMP0*scale;
251              *V10=TEMP1*scale;              *V10=TEMP1*scale;
252              if (TEMP1<0.) {              if (TEMP1<0.) {
253                 *V01=-(*V10);                 *V01=-(*V10);
254                 *V11= (*V00);                 *V11= (*V00);
255              } else {              } else {
256                 *V01= (*V10);                 *V01= (*V10);
257                 *V11=-(*V00);                 *V11=-(*V00);
258              }              }
259          } else {          } else {
260             *V00=0.;             *V00=0.;
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 290  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 385  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 395  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 414  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.757 changed lines Added in v.2548