escript/src/LocalOps.h

// $Id$
/* 
 ******************************************************************************
 *                                                                            *
 *       COPYRIGHT  ACcESS 2004 -  All Rights Reserved                        *
 *                                                                            *
 * This software is the property of ACcESS. No part of this code              *
 * may be copied in any form or by any means without the expressed written    *
 * consent of ACcESS.  Copying, use or modification of this software          *
 * by any unauthorised person is illegal unless that person has a software    *
 * license agreement with ACcESS.                                             *
 *                                                                            *
 ******************************************************************************
*/

#if !defined escript_LocalOps_H
#define escript_LocalOps_H
#include <math.h>
namespace escript {


/**
   \brief
   solves a 1x1 eigenvalue A*V=ev*V problem

   \param A00 Input - A_00 
   \param ev0 Output - eigenvalue
*/
inline
void eigenvalues1(const double A00,double* ev0) {

   *ev0=A00;

}
/**
   \brief
   solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A

   \param A00 Input - A_00 
   \param A01 Input - A_01 
   \param A11 Input - A_11
   \param ev0 Output - smallest eigenvalue
   \param ev1 Output - largest eigenvalue
*/
inline
void eigenvalues2(const double A00,const double A01
                                 ,const double A11,
                 double* ev0, double* ev1) {
      const register double trA=(A00+A11)/2.;
      const register double A_00=A00-trA;
      const register double A_11=A11-trA;
      const register double s=sqrt(A01*A01-A_00*A_11);
      *ev0=trA-s;
      *ev1=trA+s;
}
/**
   \brief
   solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A

   \param A00 Input - A_00 
   \param A01 Input - A_01 
   \param A02 Input - A_02 
   \param A11 Input - A_11 
   \param A12 Input - A_12 
   \param A22 Input - A_22 
   \param ev0 Output - smallest eigenvalue
   \param ev1 Output - eigenvalue
   \param ev2 Output - largest eigenvalue
*/
inline
void eigenvalues3(const double A00, const double A01, const double A02,
                                   const double A11, const double A12,
                                                     const double A22,
                 double* ev0, double* ev1,double* ev2) {

      const register double trA=(A00+A11+A22)/3.;
      const register double A_00=A00-trA;
      const register double A_11=A11-trA;
      const register double A_22=A22-trA;
      const register double A01_2=A01*A01;
      const register double A02_2=A02*A02;
      const register double A12_2=A12*A12;
      const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
      const register double q=(A02_2*A_11+A12_2*A_00+A01_2*A_22)-(A_00*A_11*A_22+2*A01*A12*A02);
      const register double sq_p=sqrt(p/3.);
      register double z=-q/(2*pow(sq_p,3));
      if (z<-1.) {
         z=-1.;
      } else if (z>1.) {
         z=1.;
      }
      const register double alpha_3=acos(z)/3.;
      *ev2=trA+2.*sq_p*cos(alpha_3);
      *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.);
      *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.);
}
} // end of namespace
#endif
1	gross	576	// $Id$
2			/*
3			******************************************************************************
4			* *
5			* COPYRIGHT ACcESS 2004 - All Rights Reserved *
6			* *
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 *
9			* consent of ACcESS. Copying, use or modification of this software *
10			* by any unauthorised person is illegal unless that person has a software *
11			* license agreement with ACcESS. *
12			* *
13			******************************************************************************
14			*/
15
16			#if !defined escript_LocalOps_H
17			#define escript_LocalOps_H
18	gross	580	#include <math.h>
19	gross	576	namespace escript {
20
21
22			/**
23			\brief
24			solves a 1x1 eigenvalue AV=evV problem
25
26			\param A00 Input - A_00
27			\param ev0 Output - eigenvalue
28			*/
29			inline
30			void eigenvalues1(const double A00,double* ev0) {
31
32	gross	580	*ev0=A00;
33	gross	576
34			}
35			/**
36			\brief
37			solves a 2x2 eigenvalue AV=evV problem for symmetric A
38
39			\param A00 Input - A_00
40			\param A01 Input - A_01
41			\param A11 Input - A_11
42			\param ev0 Output - smallest eigenvalue
43			\param ev1 Output - largest eigenvalue
44			*/
45			inline
46			void eigenvalues2(const double A00,const double A01
47			,const double A11,
48			double* ev0, double* ev1) {
49	gross	580	const register double trA=(A00+A11)/2.;
50			const register double A_00=A00-trA;
51			const register double A_11=A11-trA;
52			const register double s=sqrt(A01A01-A_00A_11);
53			*ev0=trA-s;
54			*ev1=trA+s;
55	gross	576	}
56			/**
57			\brief
58			solves a 3x3 eigenvalue AV=evV problem for symmetric A
59
60			\param A00 Input - A_00
61			\param A01 Input - A_01
62			\param A02 Input - A_02
63			\param A11 Input - A_11
64			\param A12 Input - A_12
65			\param A22 Input - A_22
66			\param ev0 Output - smallest eigenvalue
67			\param ev1 Output - eigenvalue
68			\param ev2 Output - largest eigenvalue
69			*/
70			inline
71			void eigenvalues3(const double A00, const double A01, const double A02,
72			const double A11, const double A12,
73			const double A22,
74			double* ev0, double* ev1,double* ev2) {
75
76	gross	580	const register double trA=(A00+A11+A22)/3.;
77			const register double A_00=A00-trA;
78			const register double A_11=A11-trA;
79			const register double A_22=A22-trA;
80			const register double A01_2=A01*A01;
81			const register double A02_2=A02*A02;
82			const register double A12_2=A12*A12;
83			const register double p=A02_2+A12_2+A01_2+(A_00A_00+A_11A_11+A_22*A_22)/2.;
84			const register double q=(A02_2A_11+A12_2A_00+A01_2A_22)-(A_00A_11A_22+2A01A12A02);
85			const register double sq_p=sqrt(p/3.);
86			register double z=-q/(2*pow(sq_p,3));
87			if (z<-1.) {
88			z=-1.;
89			} else if (z>1.) {
90			z=1.;
91			}
92			const register double alpha_3=acos(z)/3.;
93			ev2=trA+2.sq_p*cos(alpha_3);
94			ev1=trA-2.sq_p*cos(alpha_3+M_PI/3.);
95			ev0=trA-2.sq_p*cos(alpha_3-M_PI/3.);
96	gross	576	}
97			} // end of namespace
98			#endif