/[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 1326 by ksteube, Mon Sep 24 06:18:44 2007 UTC revision 1327 by matt, Fri Oct 12 07:10:40 2007 UTC
# Line 31  namespace escript { Line 31  namespace escript {
31     \brief     \brief
32     solves a 1x1 eigenvalue A*V=ev*V problem     solves a 1x1 eigenvalue A*V=ev*V problem
33    
34     \param A00 Input - A_00     \param A00 Input - A_00
35     \param ev0 Output - eigenvalue     \param ev0 Output - eigenvalue
36  */  */
37  inline  inline
# Line 44  void eigenvalues1(const double A00,doubl Line 44  void eigenvalues1(const double A00,doubl
44     \brief     \brief
45     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
46    
47     \param A00 Input - A_00     \param A00 Input - A_00
48     \param A01 Input - A_01     \param A01 Input - A_01
49     \param A11 Input - A_11     \param A11 Input - A_11
50     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
51     \param ev1 Output - largest eigenvalue     \param ev1 Output - largest eigenvalue
# Line 64  void eigenvalues2(const double A00,const Line 64  void eigenvalues2(const double A00,const
64     \brief     \brief
65     solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A     solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
66    
67     \param A00 Input - A_00     \param A00 Input - A_00
68     \param A01 Input - A_01     \param A01 Input - A_01
69     \param A02 Input - A_02     \param A02 Input - A_02
70     \param A11 Input - A_11     \param A11 Input - A_11
71     \param A12 Input - A_12     \param A12 Input - A_12
72     \param A22 Input - A_22     \param A22 Input - A_22
73     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
74     \param ev1 Output - eigenvalue     \param ev1 Output - eigenvalue
75     \param ev2 Output - largest eigenvalue     \param ev2 Output - largest eigenvalue
# Line 112  void eigenvalues3(const double A00, cons Line 112  void eigenvalues3(const double A00, cons
112     \brief     \brief
113     solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A     solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
114    
115     \param A00 Input - A_00     \param A00 Input - A_00
116     \param ev0 Output - eigenvalue     \param ev0 Output - eigenvalue
117     \param V00 Output - eigenvector     \param V00 Output - eigenvector
118     \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 159  void  vectorInKernel2(const double A00,c
159  }  }
160  /**  /**
161     \brief     \brief
162     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]]
163     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.
164    
165     \param A00 Input - matrix component     \param A00 Input - matrix component
# Line 191  void  vectorInKernel3__nonZeroA00(const Line 191  void  vectorInKernel3__nonZeroA00(const
191      *V1=A00*TEMP0;      *V1=A00*TEMP0;
192      *V2=A00*TEMP1;      *V2=A00*TEMP1;
193  }  }
194          
195  /**  /**
196     \brief     \brief
197     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
198     ordered by increasing value and eigen vectors are normalizeVector3d such that     ordered by increasing value and eigen vectors are normalizeVector3d such that
199     length is zero and first non-zero component is positive.     length is zero and first non-zero component is positive.
200    
201     \param A00 Input - A_00     \param A00 Input - A_00
202     \param A01 Input - A_01     \param A01 Input - A_01
203     \param A11 Input - A_11     \param A11 Input - A_11
204     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
205     \param ev1 Output - eigenvalue     \param ev1 Output - eigenvalue
206     \param V00 Output - eigenvector componenent coresponding to ev0     \param V00 Output - eigenvector componenent coresponding to ev0
# Line 213  inline Line 213  inline
213  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,
214                                      double* ev0, double* ev1,                                      double* ev0, double* ev1,
215                                      double* V00, double* V10, double* V01, double* V11,                                      double* V00, double* V10, double* V01, double* V11,
216                                      const double tol)                                      const double tol)
217  {  {
218       double TEMP0,TEMP1;       double TEMP0,TEMP1;
219       eigenvalues2(A00,A01,A11,ev0,ev1);       eigenvalues2(A00,A01,A11,ev0,ev1);
# Line 232  void  eigenvalues_and_eigenvectors2(cons Line 232  void  eigenvalues_and_eigenvectors2(cons
232              *V00=-TEMP0*scale;              *V00=-TEMP0*scale;
233              *V10=-TEMP1*scale;              *V10=-TEMP1*scale;
234              if (TEMP1<0.) {              if (TEMP1<0.) {
235                 *V01=  *V10;                 *V01=  *V10;
236                 *V11=-(*V00);                 *V11=-(*V00);
237              } else {              } else {
238                 *V01=-(*V10);                 *V01=-(*V10);
# Line 242  void  eigenvalues_and_eigenvectors2(cons Line 242  void  eigenvalues_and_eigenvectors2(cons
242              *V00=TEMP0*scale;              *V00=TEMP0*scale;
243              *V10=TEMP1*scale;              *V10=TEMP1*scale;
244              if (TEMP1<0.) {              if (TEMP1<0.) {
245                 *V01=-(*V10);                 *V01=-(*V10);
246                 *V11= (*V00);                 *V11= (*V00);
247              } else {              } else {
248                 *V01= (*V10);                 *V01= (*V10);
249                 *V11=-(*V00);                 *V11=-(*V00);
250              }              }
251          } else {          } else {
252             *V00=0.;             *V00=0.;
253             *V10=1;             *V10=1;
254             *V11=0.;             *V11=0.;
255             *V01=1.;             *V01=1.;
256         }         }
257     }     }
258  }  }
259  /**  /**
260     \brief     \brief
261     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.
262    
263     \param V0 - vector componenent     \param V0 - vector componenent
264     \param V1 - vector componenent     \param V1 - vector componenent
265     \param V2 - vector componenent     \param V2 - vector componenent
266  */  */
# Line 294  void  normalizeVector3(double* V0,double Line 294  void  normalizeVector3(double* V0,double
294  }  }
295  /**  /**
296     \brief     \brief
297     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
298     ordered by increasing value and eigen vectors are normalizeVector3d such that     ordered by increasing value and eigen vectors are normalizeVector3d such that
299     length is zero and first non-zero component is positive.     length is zero and first non-zero component is positive.
300    
301     \param A00 Input - A_00     \param A00 Input - A_00
302     \param A01 Input - A_01     \param A01 Input - A_01
303     \param A11 Input - A_11     \param A11 Input - A_11
304     \param ev0 Output - smallest eigenvalue     \param ev0 Output - smallest eigenvalue
305     \param ev1 Output - eigenvalue     \param ev1 Output - eigenvalue
306     \param V00 Output - eigenvector componenent coresponding to ev0     \param V00 Output - eigenvector componenent coresponding to ev0
# Line 313  inline Line 313  inline
313  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,
314                                      const double A11, const double A12, const double A22,                                      const double A11, const double A12, const double A22,
315                                      double* ev0, double* ev1, double* ev2,                                      double* ev0, double* ev1, double* ev2,
316                                      double* V00, double* V10, double* V20,                                      double* V00, double* V10, double* V20,
317                                      double* V01, double* V11, double* V21,                                      double* V01, double* V11, double* V21,
318                                      double* V02, double* V12, double* V22,                                      double* V02, double* V12, double* V22,
319                                      const double tol)                                      const double tol)
320  {  {
321        register const double absA01=fabs(A01);        register const double absA01=fabs(A01);
# Line 459  void matrix_matrix_product(const int SL, Line 459  void matrix_matrix_product(const int SL,
459    }    }
460  }  }
461    
462    template <typename BinaryFunction>
463    inline void tensor_binary_operation(const int size,
464                     const double *arg1,
465                     const double *arg2,
466                     double * argRes,
467                     BinaryFunction operation)
468    {
469      for (int i = 0; i < size; ++i) {
470        argRes[i] = operation(arg1[i], arg2[i]);
471      }
472      return;
473    }
474    
475    template <typename BinaryFunction>
476    inline void tensor_binary_operation(const int size,
477                     double arg1,
478                     const double *arg2,
479                     double *argRes,
480                     BinaryFunction operation)
481    {
482      for (int i = 0; i < size; ++i) {
483        argRes[i] = operation(arg1, arg2[i]);
484      }
485      return;
486    }
487    
488    template <typename BinaryFunction>
489    inline void tensor_binary_operation(const int size,
490                     const double *arg1,
491                     double arg2,
492                     double *argRes,
493                     BinaryFunction operation)
494    {
495      for (int i = 0; i < size; ++i) {
496        argRes[i] = operation(arg1[i], arg2);
497      }
498      return;
499    }
500    
501  } // end of namespace  } // end of namespace
502  #endif  #endif
Legend:
Removed from v.1326  
changed lines
  Added in v.1327
  ViewVC Help
Powered by ViewVC 1.1.26