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.);
}
/**
   \brief
   solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A

   \param A00 Input - A_00 
   \param ev0 Output - eigenvalue
   \param V00 Output - eigenvector
   \param tol Input - tolerance to identify to eigenvalues
*/
inline
void  eigenvalues_and_eigenvectors1(const double A00,double* ev0,double* V00,const double tol)
{
      eigenvalues1(A00,ev0);
      *V00=1.;
      return;
}
/**
   \brief
   returns a non-zero vector in the kernel of [[A00,A01],[A01,A11]] assuming that the kernel dimension is at least 1.

   \param A00 Input - matrix component
   \param A10 Input - matrix component
   \param A01 Input - matrix component
   \param A11 Input - matrix component
   \param V0 Output - vector component
   \param V1 Output - vector component
*/
inline
void  vectorInKernel2(const double A00,const double A10,const double A01,const double A11,
                      double* V0, double*V1)
{
      register double absA00=fabs(A00);
      register double absA01=fabs(A01);
      register double absA10=fabs(A10);
      register double absA11=fabs(A11);
      register double m=absA11>absA01 ? absA11 : absA01;
      if (absA00>m || absA10>m) {
         *V0=-A10;
         *V1=A00;
      } else {
         if (m<=0) {
           *V0=1.;
           *V1=0.;
         } else {
           *V0=A11;
           *V1=-A01;
         }
     }
}
/**
   \brief
   returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]] 
   assuming that the kernel dimension is at least 1 and A00 is non zero.

   \param A00 Input - matrix component
   \param A10 Input - matrix component
   \param A20 Input - matrix component
   \param A01 Input - matrix component
   \param A11 Input - matrix component
   \param A21 Input - matrix component
   \param A02 Input - matrix component
   \param A12 Input - matrix component
   \param A22 Input - matrix component
   \param V0 Output - vector component
   \param V1 Output - vector component
   \param V2 Output - vector component
*/
inline
void  vectorInKernel3__nonZeroA00(const double A00,const double A10,const double A20,
                                const double A01,const double A11,const double A21,
                                const double A02,const double A12,const double A22,
                                double* V0,double* V1,double* V2)
{
    double TEMP0,TEMP1;
    register const double I00=1./A00;
    register const double IA10=I00*A10;
    register const double IA20=I00*A20;
    vectorInKernel2(A11-IA10*A01,A21-IA20*A01,A12-IA10*A02,A22-IA20*A02,&TEMP0,&TEMP1);
    *V0=-(A10*TEMP0+A20*TEMP1);
    *V1=A00*TEMP0;
    *V2=A00*TEMP1;
}
       
/**
   \brief
   solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are 
   ordered by increasing value and eigen vectors are normalizeVector3d such that 
   length is zero and first non-zero component is positive.

   \param A00 Input - A_00 
   \param A01 Input - A_01 
   \param A11 Input - A_11 
   \param ev0 Output - smallest eigenvalue
   \param ev1 Output - eigenvalue
   \param V00 Output - eigenvector componenent coresponding to ev0
   \param V10 Output - eigenvector componenent coresponding to ev0
   \param V01 Output - eigenvector componenent coresponding to ev1
   \param V11 Output - eigenvector componenent coresponding to ev1
   \param tol Input - tolerance to identify to eigenvalues
*/
inline
void  eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11,
                                    double* ev0, double* ev1,
                                    double* V00, double* V10, double* V01, double* V11,
                                    const double tol) 
{
     double TEMP0,TEMP1;
     eigenvalues2(A00,A01,A11,ev0,ev1);
     const register double absev0=fabs(*ev0);
     const register double absev1=fabs(*ev1);
     register double max_ev=absev0>absev1 ? absev0 : absev1;
     if (fabs((*ev0)-(*ev1))<tol*max_ev) {
        *V00=1.;
        *V10=0.;
        *V01=0.;
        *V11=1.;
     } else {
        vectorInKernel2(A00-(*ev0),A01,A01,A11-(*ev0),&TEMP0,&TEMP1);
        const register double scale=1./sqrt(TEMP0*TEMP0+TEMP1*TEMP1);
        if (TEMP0<0.) {
            *V00=-TEMP0*scale;
            *V10=-TEMP1*scale;
            if (TEMP1<0.) {
               *V01=  *V10;
               *V11=-(*V00);
            } else {
               *V01=-(*V10);
               *V11= (*V00);
            }
        } else if (TEMP0>0.) {
            *V00=TEMP0*scale;
            *V10=TEMP1*scale;
            if (TEMP1<0.) {
               *V01=-(*V10);
               *V11= (*V00);
            } else {
               *V01= (*V10);
               *V11=-(*V00);
            }
        } else {
           *V00=0.;
           *V10=1;
           *V11=0.;
           *V01=1.;
       } 
   }
}
/**
   \brief
   nomalizes a 3-d vector such that length is one and first non-zero component is positive.

   \param V0 - vector componenent 
   \param V1 - vector componenent
   \param V2 - vector componenent
*/
inline
void  normalizeVector3(double* V0,double* V1,double* V2)
{
    register double s;
    if (*V0>0) {
        s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
        *V0*=s;
        *V1*=s;
        *V2*=s;
    } else if (*V0<0)  {
        s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
        *V0*=s;
        *V1*=s;
        *V2*=s;
    } else {
        if (*V1>0) {
            s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
            *V1*=s;
            *V2*=s;
        } else if (*V1<0)  {
            s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
            *V1*=s;
            *V2*=s;
        } else {
            *V2=1.;
        }
    }
}
/**
   \brief
   solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are 
   ordered by increasing value and eigen vectors are normalizeVector3d such that 
   length is zero and first non-zero component is positive.

   \param A00 Input - A_00 
   \param A01 Input - A_01 
   \param A11 Input - A_11 
   \param ev0 Output - smallest eigenvalue
   \param ev1 Output - eigenvalue
   \param V00 Output - eigenvector componenent coresponding to ev0
   \param V10 Output - eigenvector componenent coresponding to ev0
   \param V01 Output - eigenvector componenent coresponding to ev1
   \param V11 Output - eigenvector componenent coresponding to ev1
   \param tol Input - tolerance to identify to eigenvalues
*/
inline
void  eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02,
                                    const double A11, const double A12, const double A22,
                                    double* ev0, double* ev1, double* ev2,
                                    double* V00, double* V10, double* V20, 
                                    double* V01, double* V11, double* V21, 
                                    double* V02, double* V12, double* V22, 
                                    const double tol)
{
      register const double absA01=fabs(A01);
      register const double absA02=fabs(A02);
      register const double m=absA01>absA02 ? absA01 : absA02;
      if (m<=0) {
        double TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1;
        eigenvalues_and_eigenvectors2(A11,A12,A22,
                                      &TEMP_EV0,&TEMP_EV1,
                                      &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
        if (A00<=TEMP_EV0) {
            *V00=1.;
            *V10=0.;
            *V20=0.;
            *V01=0.;
            *V11=TEMP_V00;
            *V21=TEMP_V10;
            *V02=0.;
            *V12=TEMP_V01;
            *V22=TEMP_V11;
            *ev0=A00;
            *ev1=TEMP_EV0;
            *ev2=TEMP_EV1;
        } else if (A00>TEMP_EV1) {
            *V00=TEMP_V00;
            *V10=TEMP_V10;
            *V20=0.;
            *V01=TEMP_V01;
            *V11=TEMP_V11;
            *V21=0.;
            *V02=0.;
            *V12=0.;
            *V22=1.;
            *ev0=TEMP_EV0;
            *ev1=TEMP_EV1;
            *ev0=A00;
        } else {
            *V00=TEMP_V00;
            *V10=0;
            *V20=TEMP_V10;
            *V01=0.;
            *V11=1.;
            *V21=0.;
            *V02=TEMP_V01;
            *V12=0.;
            *V22=TEMP_V11;
            *ev0=TEMP_EV0;
            *ev1=A00;
            *ev2=TEMP_EV1;
        }
      } else {
         eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2);
         const register double absev0=fabs(*ev0);
         const register double absev1=fabs(*ev1);
         const register double absev2=fabs(*ev2);
         register double max_ev=absev0>absev1 ? absev0 : absev1;
         max_ev=max_ev>absev2 ? max_ev : absev2;
         const register double d_01=fabs((*ev0)-(*ev1));
         const register double d_12=fabs((*ev1)-(*ev2));
         const register double max_d=d_01>d_12 ? d_01 : d_12;
         if (max_d<=tol*max_ev) {
             *V00=1.;
             *V10=0;
             *V20=0;
             *V01=0;
             *V11=1.;
             *V21=0;
             *V02=0;
             *V12=0;
             *V22=1.;
         } else {
            const register double S00=A00-(*ev0);
            const register double absS00=fabs(S00);
            if (fabs(S00)>m) {
                vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);
            } else if (absA02<m) {
                vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);
            } else {
                vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
            }
            normalizeVector3(V00,V10,V20);;
            const register double T00=A00-(*ev2);
            const register double absT00=fabs(T00);
            if (fabs(T00)>m) {
                 vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);
            } else if (absA02<m) {
                 vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);
            } else {
                 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
            }
            const register double dot=(*V02)*(*V00)+(*V12)*(*V10)+(*V22)*(*V20);
            *V02-=dot*(*V00);
            *V12-=dot*(*V10);
            *V22-=dot*(*V20);
            normalizeVector3(V02,V12,V22);
            *V01=(*V10)*(*V22)-(*V12)*(*V20);
            *V11=(*V20)*(*V02)-(*V00)*(*V22);
            *V21=(*V00)*(*V12)-(*V02)*(*V10);
            normalizeVector3(V01,V11,V21);
         }
   }
}
} // 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	gross	583	void eigenvalues2(const double A00,const double A01,const double A11,
47	gross	576	double* ev0, double* ev1) {
48	gross	580	const register double trA=(A00+A11)/2.;
49			const register double A_00=A00-trA;
50			const register double A_11=A11-trA;
51			const register double s=sqrt(A01A01-A_00A_11);
52			*ev0=trA-s;
53			*ev1=trA+s;
54	gross	576	}
55			/**
56			\brief
57			solves a 3x3 eigenvalue AV=evV problem for symmetric A
58
59			\param A00 Input - A_00
60			\param A01 Input - A_01
61			\param A02 Input - A_02
62			\param A11 Input - A_11
63			\param A12 Input - A_12
64			\param A22 Input - A_22
65			\param ev0 Output - smallest eigenvalue
66			\param ev1 Output - eigenvalue
67			\param ev2 Output - largest eigenvalue
68			*/
69			inline
70			void eigenvalues3(const double A00, const double A01, const double A02,
71			const double A11, const double A12,
72			const double A22,
73			double* ev0, double* ev1,double* ev2) {
74
75	gross	580	const register double trA=(A00+A11+A22)/3.;
76			const register double A_00=A00-trA;
77			const register double A_11=A11-trA;
78			const register double A_22=A22-trA;
79			const register double A01_2=A01*A01;
80			const register double A02_2=A02*A02;
81			const register double A12_2=A12*A12;
82			const register double p=A02_2+A12_2+A01_2+(A_00A_00+A_11A_11+A_22*A_22)/2.;
83			const register double q=(A02_2A_11+A12_2A_00+A01_2A_22)-(A_00A_11A_22+2A01A12A02);
84			const register double sq_p=sqrt(p/3.);
85			register double z=-q/(2*pow(sq_p,3));
86			if (z<-1.) {
87			z=-1.;
88			} else if (z>1.) {
89			z=1.;
90			}
91			const register double alpha_3=acos(z)/3.;
92			ev2=trA+2.sq_p*cos(alpha_3);
93			ev1=trA-2.sq_p*cos(alpha_3+M_PI/3.);
94			ev0=trA-2.sq_p*cos(alpha_3-M_PI/3.);
95	gross	576	}
96	gross	583	/**
97			\brief
98			solves a 1x1 eigenvalue AV=evV problem for symmetric A
99
100			\param A00 Input - A_00
101			\param ev0 Output - eigenvalue
102			\param V00 Output - eigenvector
103			\param tol Input - tolerance to identify to eigenvalues
104			*/
105			inline
106			void eigenvalues_and_eigenvectors1(const double A00,double* ev0,double* V00,const double tol)
107			{
108			eigenvalues1(A00,ev0);
109			*V00=1.;
110			return;
111			}
112			/**
113			\brief
114			returns a non-zero vector in the kernel of [[A00,A01],[A01,A11]] assuming that the kernel dimension is at least 1.
115
116			\param A00 Input - matrix component
117			\param A10 Input - matrix component
118			\param A01 Input - matrix component
119			\param A11 Input - matrix component
120			\param V0 Output - vector component
121			\param V1 Output - vector component
122			*/
123			inline
124			void vectorInKernel2(const double A00,const double A10,const double A01,const double A11,
125			double* V0, double*V1)
126			{
127			register double absA00=fabs(A00);
128			register double absA01=fabs(A01);
129			register double absA10=fabs(A10);
130			register double absA11=fabs(A11);
131			register double m=absA11>absA01 ? absA11 : absA01;
132			if (absA00>m \|\| absA10>m) {
133			*V0=-A10;
134			*V1=A00;
135			} else {
136			if (m<=0) {
137			*V0=1.;
138			*V1=0.;
139			} else {
140			*V0=A11;
141			*V1=-A01;
142			}
143			}
144			}
145			/**
146			\brief
147			returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]]
148			assuming that the kernel dimension is at least 1 and A00 is non zero.
149
150			\param A00 Input - matrix component
151			\param A10 Input - matrix component
152			\param A20 Input - matrix component
153			\param A01 Input - matrix component
154			\param A11 Input - matrix component
155			\param A21 Input - matrix component
156			\param A02 Input - matrix component
157			\param A12 Input - matrix component
158			\param A22 Input - matrix component
159			\param V0 Output - vector component
160			\param V1 Output - vector component
161			\param V2 Output - vector component
162			*/
163			inline
164			void vectorInKernel3__nonZeroA00(const double A00,const double A10,const double A20,
165			const double A01,const double A11,const double A21,
166			const double A02,const double A12,const double A22,
167			double* V0,double* V1,double* V2)
168			{
169			double TEMP0,TEMP1;
170			register const double I00=1./A00;
171			register const double IA10=I00*A10;
172			register const double IA20=I00*A20;
173			vectorInKernel2(A11-IA10A01,A21-IA20A01,A12-IA10A02,A22-IA20A02,&TEMP0,&TEMP1);
174			V0=-(A10TEMP0+A20*TEMP1);
175			V1=A00TEMP0;
176			V2=A00TEMP1;
177			}
178
179			/**
180			\brief
181			solves a 2x2 eigenvalue AV=evV problem for symmetric A. Eigenvectors are
182			ordered by increasing value and eigen vectors are normalizeVector3d such that
183			length is zero and first non-zero component is positive.
184
185			\param A00 Input - A_00
186			\param A01 Input - A_01
187			\param A11 Input - A_11
188			\param ev0 Output - smallest eigenvalue
189			\param ev1 Output - eigenvalue
190			\param V00 Output - eigenvector componenent coresponding to ev0
191			\param V10 Output - eigenvector componenent coresponding to ev0
192			\param V01 Output - eigenvector componenent coresponding to ev1
193			\param V11 Output - eigenvector componenent coresponding to ev1
194			\param tol Input - tolerance to identify to eigenvalues
195			*/
196			inline
197			void eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11,
198			double* ev0, double* ev1,
199			double* V00, double* V10, double* V01, double* V11,
200			const double tol)
201			{
202			double TEMP0,TEMP1;
203			eigenvalues2(A00,A01,A11,ev0,ev1);
204			const register double absev0=fabs(*ev0);
205			const register double absev1=fabs(*ev1);
206			register double max_ev=absev0>absev1 ? absev0 : absev1;
207			if (fabs((ev0)-(ev1))<tol*max_ev) {
208			*V00=1.;
209			*V10=0.;
210			*V01=0.;
211			*V11=1.;
212			} else {
213			vectorInKernel2(A00-(ev0),A01,A01,A11-(ev0),&TEMP0,&TEMP1);
214			const register double scale=1./sqrt(TEMP0TEMP0+TEMP1TEMP1);
215			if (TEMP0<0.) {
216			V00=-TEMP0scale;
217			V10=-TEMP1scale;
218			if (TEMP1<0.) {
219			V01= V10;
220			V11=-(V00);
221			} else {
222			V01=-(V10);
223			V11= (V00);
224			}
225			} else if (TEMP0>0.) {
226			V00=TEMP0scale;
227			V10=TEMP1scale;
228			if (TEMP1<0.) {
229			V01=-(V10);
230			V11= (V00);
231			} else {
232			V01= (V10);
233			V11=-(V00);
234			}
235			} else {
236			*V00=0.;
237			*V10=1;
238			*V11=0.;
239			*V01=1.;
240			}
241			}
242			}
243			/**
244			\brief
245			nomalizes a 3-d vector such that length is one and first non-zero component is positive.
246
247			\param V0 - vector componenent
248			\param V1 - vector componenent
249			\param V2 - vector componenent
250			*/
251			inline
252			void normalizeVector3(double* V0,double* V1,double* V2)
253			{
254			register double s;
255			if (*V0>0) {
256			s=1./sqrt((V0)(V0)+(V1)(V1)+(V2)(*V2));
257			V0=s;
258			V1=s;
259			V2=s;
260			} else if (*V0<0) {
261			s=-1./sqrt((V0)(V0)+(V1)(V1)+(V2)(*V2));
262			V0=s;
263			V1=s;
264			V2=s;
265			} else {
266			if (*V1>0) {
267			s=1./sqrt((V1)(V1)+(V2)(V2));
268			V1=s;
269			V2=s;
270			} else if (*V1<0) {
271			s=-1./sqrt((V1)(V1)+(V2)(V2));
272			V1=s;
273			V2=s;
274			} else {
275			*V2=1.;
276			}
277			}
278			}
279			/**
280			\brief
281			solves a 2x2 eigenvalue AV=evV problem for symmetric A. Eigenvectors are
282			ordered by increasing value and eigen vectors are normalizeVector3d such that
283			length is zero and first non-zero component is positive.
284
285			\param A00 Input - A_00
286			\param A01 Input - A_01
287			\param A11 Input - A_11
288			\param ev0 Output - smallest eigenvalue
289			\param ev1 Output - eigenvalue
290			\param V00 Output - eigenvector componenent coresponding to ev0
291			\param V10 Output - eigenvector componenent coresponding to ev0
292			\param V01 Output - eigenvector componenent coresponding to ev1
293			\param V11 Output - eigenvector componenent coresponding to ev1
294			\param tol Input - tolerance to identify to eigenvalues
295			*/
296			inline
297			void eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02,
298			const double A11, const double A12, const double A22,
299			double* ev0, double* ev1, double* ev2,
300			double* V00, double* V10, double* V20,
301			double* V01, double* V11, double* V21,
302			double* V02, double* V12, double* V22,
303			const double tol)
304			{
305			register const double absA01=fabs(A01);
306			register const double absA02=fabs(A02);
307			register const double m=absA01>absA02 ? absA01 : absA02;
308			if (m<=0) {
309			double TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1;
310			eigenvalues_and_eigenvectors2(A11,A12,A22,
311			&TEMP_EV0,&TEMP_EV1,
312			&TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
313			if (A00<=TEMP_EV0) {
314			*V00=1.;
315			*V10=0.;
316			*V20=0.;
317			*V01=0.;
318			*V11=TEMP_V00;
319			*V21=TEMP_V10;
320			*V02=0.;
321			*V12=TEMP_V01;
322			*V22=TEMP_V11;
323			*ev0=A00;
324			*ev1=TEMP_EV0;
325			*ev2=TEMP_EV1;
326			} else if (A00>TEMP_EV1) {
327			*V00=TEMP_V00;
328			*V10=TEMP_V10;
329			*V20=0.;
330			*V01=TEMP_V01;
331			*V11=TEMP_V11;
332			*V21=0.;
333			*V02=0.;
334			*V12=0.;
335			*V22=1.;
336			*ev0=TEMP_EV0;
337			*ev1=TEMP_EV1;
338			*ev0=A00;
339			} else {
340			*V00=TEMP_V00;
341			*V10=0;
342			*V20=TEMP_V10;
343			*V01=0.;
344			*V11=1.;
345			*V21=0.;
346			*V02=TEMP_V01;
347			*V12=0.;
348			*V22=TEMP_V11;
349			*ev0=TEMP_EV0;
350			*ev1=A00;
351			*ev2=TEMP_EV1;
352			}
353			} else {
354			eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2);
355			const register double absev0=fabs(*ev0);
356			const register double absev1=fabs(*ev1);
357			const register double absev2=fabs(*ev2);
358			register double max_ev=absev0>absev1 ? absev0 : absev1;
359			max_ev=max_ev>absev2 ? max_ev : absev2;
360			const register double d_01=fabs((ev0)-(ev1));
361			const register double d_12=fabs((ev1)-(ev2));
362			const register double max_d=d_01>d_12 ? d_01 : d_12;
363			if (max_d<=tol*max_ev) {
364			*V00=1.;
365			*V10=0;
366			*V20=0;
367			*V01=0;
368			*V11=1.;
369			*V21=0;
370			*V02=0;
371			*V12=0;
372			*V22=1.;
373			} else {
374			const register double S00=A00-(*ev0);
375			const register double absS00=fabs(S00);
376			if (fabs(S00)>m) {
377			vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(ev0),A12,A02,A12,A22-(ev0),V00,V10,V20);
378			} else if (absA02<m) {
379			vectorInKernel3__nonZeroA00(A01,A11-(ev0),A12,S00,A01,A02,A02,A12,A22-(ev0),V00,V10,V20);
380			} else {
381			vectorInKernel3__nonZeroA00(A02,A12,A22-(ev0),S00,A01,A02,A01,A11-(ev0),A12,V00,V10,V20);
382			}
383			normalizeVector3(V00,V10,V20);;
384			const register double T00=A00-(*ev2);
385			const register double absT00=fabs(T00);
386			if (fabs(T00)>m) {
387			vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(ev2),A12,A02,A12,A22-(ev2),V02,V12,V22);
388			} else if (absA02<m) {
389			vectorInKernel3__nonZeroA00(A01,A11-(ev2),A12,T00,A01,A02,A02,A12,A22-(ev2),V02,V12,V22);
390			} else {
391			vectorInKernel3__nonZeroA00(A02,A12,A22-(ev2),T00,A01,A02,A01,A11-(ev2),A12,V02,V12,V22);
392			}
393			const register double dot=(V02)(V00)+(V12)(V10)+(V22)(*V20);
394			V02-=dot(*V00);
395			V12-=dot(*V10);
396			V22-=dot(*V20);
397			normalizeVector3(V02,V12,V22);
398			V01=(V10)(V22)-(V12)(*V20);
399			V11=(V20)(V02)-(V00)(*V22);
400			V21=(V00)(V12)-(V02)(*V10);
401			normalizeVector3(V01,V11,V21);
402			}
403			}
404			}
405	gross	576	} // end of namespace
406			#endif