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.;
      if (p<=0.) {
         *ev2=trA;
         *ev1=trA;
         *ev0=trA;

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