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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 580 - (hide annotations)
Wed Mar 8 05:45:51 2006 UTC (13 years, 8 months ago) by gross
File MIME type: text/plain
File size: 3326 byte(s)
faster version of the local eigenvalue calculation
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 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 gross 580 *ev0=A00;
33 gross 576
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 gross 580 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 gross 576 }
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 gross 580 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 gross 576 }
97     } // end of namespace
98     #endif

  ViewVC Help
Powered by ViewVC 1.1.26