/[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 588 by gross, Fri Mar 10 04:45:04 2006 UTC revision 2777 by jfenwick, Thu Nov 26 01:06:00 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    /**
38    \brief acts as a wrapper to isnan.
39    \warning if compiler does not support FP_NAN this function will always return false.
40    */
41    inline
42    bool nancheck(double d)
43    {
44    #ifndef isnan       // Q: so why not just test d!=d?
45        return false;   // A: Coz it doesn't always work [I've checked].
46    #else           // One theory is that the optimizer skips the test.
47        return isnan(d);
48    #endif
49    }
50    
51    /**
52    \brief returns a NaN.
53    \warning Should probably only used where you know you can test for NaNs
54    */
55    inline
56    double makeNaN()
57    {
58    #ifdef isnan
59        return nan("");
60    #else
61        return sqrt(-1);
62    #endif
63    
64    }
65    
66    
67  /**  /**
68     \brief     \brief
69     solves a 1x1 eigenvalue A*V=ev*V problem     solves a 1x1 eigenvalue A*V=ev*V problem
70    
71     \param A00 Input - A_00     \param A00 Input - A_00
72     \param ev0 Output - eigenvalue     \param ev0 Output - eigenvalue
73  */  */
74  inline  inline
# Line 36  void eigenvalues1(const double A00,doubl Line 81  void eigenvalues1(const double A00,doubl
81     \brief     \brief
82     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
83    
84     \param A00 Input - A_00     \param A00 Input - A_00
85     \param A01 Input - A_01     \param A01 Input - A_01
86     \param A11 Input - A_11     \param A11 Input - A_11
87     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
88     \param ev1 Output - largest eigenvalue     \param ev1 Output - largest eigenvalue
# Line 56  void eigenvalues2(const double A00,const Line 101  void eigenvalues2(const double A00,const
101     \brief     \brief
102     solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A     solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
103    
104     \param A00 Input - A_00     \param A00 Input - A_00
105     \param A01 Input - A_01     \param A01 Input - A_01
106     \param A02 Input - A_02     \param A02 Input - A_02
107     \param A11 Input - A_11     \param A11 Input - A_11
108     \param A12 Input - A_12     \param A12 Input - A_12
109     \param A22 Input - A_22     \param A22 Input - A_22
110     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
111     \param ev1 Output - eigenvalue     \param ev1 Output - eigenvalue
112     \param ev2 Output - largest eigenvalue     \param ev2 Output - largest eigenvalue
# Line 104  void eigenvalues3(const double A00, cons Line 149  void eigenvalues3(const double A00, cons
149     \brief     \brief
150     solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A     solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
151    
152     \param A00 Input - A_00     \param A00 Input - A_00
153     \param ev0 Output - eigenvalue     \param ev0 Output - eigenvalue
154     \param V00 Output - eigenvector     \param V00 Output - eigenvector
155     \param tol Input - tolerance to identify to eigenvalues     \param tol Input - tolerance to identify to eigenvalues
# Line 151  void  vectorInKernel2(const double A00,c Line 196  void  vectorInKernel2(const double A00,c
196  }  }
197  /**  /**
198     \brief     \brief
199     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]]
200     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.
201    
202     \param A00 Input - matrix component     \param A00 Input - matrix component
# Line 183  void  vectorInKernel3__nonZeroA00(const Line 228  void  vectorInKernel3__nonZeroA00(const
228      *V1=A00*TEMP0;      *V1=A00*TEMP0;
229      *V2=A00*TEMP1;      *V2=A00*TEMP1;
230  }  }
231          
232  /**  /**
233     \brief     \brief
234     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
235     ordered by increasing value and eigen vectors are normalizeVector3d such that     ordered by increasing value and eigen vectors are normalizeVector3d such that
236     length is zero and first non-zero component is positive.     length is zero and first non-zero component is positive.
237    
238     \param A00 Input - A_00     \param A00 Input - A_00
239     \param A01 Input - A_01     \param A01 Input - A_01
240     \param A11 Input - A_11     \param A11 Input - A_11
241     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
242     \param ev1 Output - eigenvalue     \param ev1 Output - eigenvalue
243     \param V00 Output - eigenvector componenent coresponding to ev0     \param V00 Output - eigenvector componenent coresponding to ev0
# Line 205  inline Line 250  inline
250  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,
251                                      double* ev0, double* ev1,                                      double* ev0, double* ev1,
252                                      double* V00, double* V10, double* V01, double* V11,                                      double* V00, double* V10, double* V01, double* V11,
253                                      const double tol)                                      const double tol)
254  {  {
255       double TEMP0,TEMP1;       double TEMP0,TEMP1;
256       eigenvalues2(A00,A01,A11,ev0,ev1);       eigenvalues2(A00,A01,A11,ev0,ev1);
# Line 224  void  eigenvalues_and_eigenvectors2(cons Line 269  void  eigenvalues_and_eigenvectors2(cons
269              *V00=-TEMP0*scale;              *V00=-TEMP0*scale;
270              *V10=-TEMP1*scale;              *V10=-TEMP1*scale;
271              if (TEMP1<0.) {              if (TEMP1<0.) {
272                 *V01=  *V10;                 *V01=  *V10;
273                 *V11=-(*V00);                 *V11=-(*V00);
274              } else {              } else {
275                 *V01=-(*V10);                 *V01=-(*V10);
# Line 234  void  eigenvalues_and_eigenvectors2(cons Line 279  void  eigenvalues_and_eigenvectors2(cons
279              *V00=TEMP0*scale;              *V00=TEMP0*scale;
280              *V10=TEMP1*scale;              *V10=TEMP1*scale;
281              if (TEMP1<0.) {              if (TEMP1<0.) {
282                 *V01=-(*V10);                 *V01=-(*V10);
283                 *V11= (*V00);                 *V11= (*V00);
284              } else {              } else {
285                 *V01= (*V10);                 *V01= (*V10);
286                 *V11=-(*V00);                 *V11=-(*V00);
287              }              }
288          } else {          } else {
289             *V00=0.;             *V00=0.;
290             *V10=1;             *V10=1;
291             *V11=0.;             *V11=0.;
292             *V01=1.;             *V01=1.;
293         }         }
294     }     }
295  }  }
296  /**  /**
297     \brief     \brief
298     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.
299    
300     \param V0 - vector componenent     \param V0 - vector componenent
301     \param V1 - vector componenent     \param V1 - vector componenent
302     \param V2 - vector componenent     \param V2 - vector componenent
303  */  */
# Line 286  void  normalizeVector3(double* V0,double Line 331  void  normalizeVector3(double* V0,double
331  }  }
332  /**  /**
333     \brief     \brief
334     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
335     ordered by increasing value and eigen vectors are normalizeVector3d such that     ordered by increasing value and eigen vectors are normalizeVector3d such that
336     length is zero and first non-zero component is positive.     length is zero and first non-zero component is positive.
337    
338     \param A00 Input - A_00     \param A00 Input - A_00
339     \param A01 Input - A_01     \param A01 Input - A_01
340     \param A11 Input - A_11     \param A02 Input - A_02
341       \param A11 Input - A_11
342       \param A12 Input - A_12
343       \param A22 Input - A_22
344     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
345     \param ev1 Output - eigenvalue     \param ev1 Output - eigenvalue
346       \param ev2 Output -
347     \param V00 Output - eigenvector componenent coresponding to ev0     \param V00 Output - eigenvector componenent coresponding to ev0
348     \param V10 Output - eigenvector componenent coresponding to ev0     \param V10 Output - eigenvector componenent coresponding to ev0
349       \param V20 Output -
350     \param V01 Output - eigenvector componenent coresponding to ev1     \param V01 Output - eigenvector componenent coresponding to ev1
351     \param V11 Output - eigenvector componenent coresponding to ev1     \param V11 Output - eigenvector componenent coresponding to ev1
352       \param V21 Output -
353       \param V02 Output -
354       \param V12 Output -
355       \param V22 Output -
356     \param tol Input - tolerance to identify to eigenvalues     \param tol Input - tolerance to identify to eigenvalues
357  */  */
358  inline  inline
359  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,
360                                      const double A11, const double A12, const double A22,                                      const double A11, const double A12, const double A22,
361                                      double* ev0, double* ev1, double* ev2,                                      double* ev0, double* ev1, double* ev2,
362                                      double* V00, double* V10, double* V20,                                      double* V00, double* V10, double* V20,
363                                      double* V01, double* V11, double* V21,                                      double* V01, double* V11, double* V21,
364                                      double* V02, double* V12, double* V22,                                      double* V02, double* V12, double* V22,
365                                      const double tol)                                      const double tol)
366  {  {
367        register const double absA01=fabs(A01);        register const double absA01=fabs(A01);
# Line 381  void  eigenvalues_and_eigenvectors3(cons Line 435  void  eigenvalues_and_eigenvectors3(cons
435           } else {           } else {
436              const register double S00=A00-(*ev0);              const register double S00=A00-(*ev0);
437              const register double absS00=fabs(S00);              const register double absS00=fabs(S00);
438              if (fabs(S00)>m) {              if (absS00>m) {
439                  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);
440              } else if (absA02<m) {              } else if (absA02<m) {
441                  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 391  void  eigenvalues_and_eigenvectors3(cons Line 445  void  eigenvalues_and_eigenvectors3(cons
445              normalizeVector3(V00,V10,V20);;              normalizeVector3(V00,V10,V20);;
446              const register double T00=A00-(*ev2);              const register double T00=A00-(*ev2);
447              const register double absT00=fabs(T00);              const register double absT00=fabs(T00);
448              if (fabs(T00)>m) {              if (absT00>m) {
449                   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);
450              } else if (absA02<m) {              } else if (absA02<m) {
451                   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 410  void  eigenvalues_and_eigenvectors3(cons Line 464  void  eigenvalues_and_eigenvectors3(cons
464           }           }
465     }     }
466  }  }
467    
468    // General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
469    // SM is the product of the last axis_offset entries in arg_0.getShape().
470    inline
471    void matrix_matrix_product(const int SL, const int SM, const int SR, const double* A, const double* B, double* C, int transpose)
472    {
473      if (transpose == 0) {
474        for (int i=0; i<SL; i++) {
475          for (int j=0; j<SR; j++) {
476            double sum = 0.0;
477            for (int l=0; l<SM; l++) {
478          sum += A[i+SL*l] * B[l+SM*j];
479            }
480            C[i+SL*j] = sum;
481          }
482        }
483      }
484      else if (transpose == 1) {
485        for (int i=0; i<SL; i++) {
486          for (int j=0; j<SR; j++) {
487            double sum = 0.0;
488            for (int l=0; l<SM; l++) {
489          sum += A[i*SM+l] * B[l+SM*j];
490            }
491            C[i+SL*j] = sum;
492          }
493        }
494      }
495      else if (transpose == 2) {
496        for (int i=0; i<SL; i++) {
497          for (int j=0; j<SR; j++) {
498            double sum = 0.0;
499            for (int l=0; l<SM; l++) {
500          sum += A[i+SL*l] * B[l*SR+j];
501            }
502            C[i+SL*j] = sum;
503          }
504        }
505      }
506    }
507    
508    template <typename UnaryFunction>
509    inline void tensor_unary_operation(const int size,
510                     const double *arg1,
511                     double * argRes,
512                     UnaryFunction operation)
513    {
514      for (int i = 0; i < size; ++i) {
515        argRes[i] = operation(arg1[i]);
516      }
517      return;
518    }
519    
520    template <typename BinaryFunction>
521    inline void tensor_binary_operation(const int size,
522                     const double *arg1,
523                     const double *arg2,
524                     double * argRes,
525                     BinaryFunction operation)
526    {
527      for (int i = 0; i < size; ++i) {
528        argRes[i] = operation(arg1[i], arg2[i]);
529      }
530      return;
531    }
532    
533    template <typename BinaryFunction>
534    inline void tensor_binary_operation(const int size,
535                     double arg1,
536                     const double *arg2,
537                     double *argRes,
538                     BinaryFunction operation)
539    {
540      for (int i = 0; i < size; ++i) {
541        argRes[i] = operation(arg1, arg2[i]);
542      }
543      return;
544    }
545    
546    template <typename BinaryFunction>
547    inline void tensor_binary_operation(const int size,
548                     const double *arg1,
549                     double arg2,
550                     double *argRes,
551                     BinaryFunction operation)
552    {
553      for (int i = 0; i < size; ++i) {
554        argRes[i] = operation(arg1[i], arg2);
555      }
556      return;
557    }
558    
559  } // end of namespace  } // end of namespace
560  #endif  #endif

Legend:
Removed from v.588  
changed lines
  Added in v.2777

  ViewVC Help
Powered by ViewVC 1.1.26