/[escript]/trunk/escript/src/LocalOps.h
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 580 - (show annotations)
Wed Mar 8 05:45:51 2006 UTC (13 years, 6 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 <math.h>
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

  ViewVC Help
Powered by ViewVC 1.1.26