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	// $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	#include <math.h>
19	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	*ev0=A00;
33
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	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	}
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	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	}
97	} // end of namespace
98	#endif