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

revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1
1
/* \$Id\$ */

2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
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 __INTEL_COMPILER  #if defined(_WIN32) && defined(__INTEL_COMPILER)
18  #   include <mathimf.h>  #   include <mathimf.h>
19  #else  #else
20  #   include <math.h>  #   include <math.h>
# Line 24  Line 23
23  #   define M_PI           3.14159265358979323846  /* pi */  #   define M_PI           3.14159265358979323846  /* pi */
24  #endif  #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 31  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 44  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 64  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 112  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 159  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 191  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 213  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 232  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 242  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 294  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 389  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 399  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 459  void matrix_matrix_product(const int SL, Line 476  void matrix_matrix_product(const int SL,
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.1312 changed lines Added in v.2548