/[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 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC revision 2777 by jfenwick, Thu Nov 26 01:06:00 2009 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2009 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 14  Line 14 
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 23  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    /**
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
# Line 299  void  normalizeVector3(double* V0,double Line 337  void  normalizeVector3(double* V0,double
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 A02 Input - A_02
341     \param A11 Input - A_11     \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
# Line 388  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 398  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);

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

  ViewVC Help
Powered by ViewVC 1.1.26