# Contents of /trunk/escript/src/LocalOps.h

Revision 580 - (show annotations)
Wed Mar 8 05:45:51 2006 UTC (16 years, 8 months ago) by gross
File MIME type: text/plain
File size: 3326 byte(s)
```faster version of the local eigenvalue calculation
```
 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 19 namespace escript { 20 21 22 /** 23 \brief 24 solves a 1x1 eigenvalue A*V=ev*V 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 A*V=ev*V 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(A01*A01-A_00*A_11); 53 *ev0=trA-s; 54 *ev1=trA+s; 55 } 56 /** 57 \brief 58 solves a 3x3 eigenvalue A*V=ev*V 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_00*A_00+A_11*A_11+A_22*A_22)/2.; 84 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); 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