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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1020 - (hide annotations)
Mon Mar 12 10:12:36 2007 UTC (12 years, 8 months ago) by phornby
File MIME type: text/plain
File size: 14873 byte(s)
Added explicit destructors to all Exception classes.

Fixed an ifdef in TestCase.cpp

Made the conditional definition of M_PI in LocalOps.h
depend only on M_PI being undefined.

Replace dynamically dimensioned arrays in DataFactory & DataTagged with malloc.

sort() method of list does not take a named argument
(despite the manual's claims to the contary).


1 gross 576 // $Id$
2     /*
3 elspeth 615 ************************************************************
4     * Copyright 2006 by ACcESS MNRF *
5     * *
6     * http://www.access.edu.au *
7     * Primary Business: Queensland, Australia *
8     * Licensed under the Open Software License version 3.0 *
9     * http://www.opensource.org/licenses/osl-3.0.php *
10     * *
11     ************************************************************
12 gross 576 */
13    
14     #if !defined escript_LocalOps_H
15     #define escript_LocalOps_H
16 woo409 779 #ifdef __INTEL_COMPILER
17 phornby 1020 # include <mathimf.h>
18 woo409 757 #else
19 phornby 1020 # include <math.h>
20 woo409 757 #endif
21 phornby 1020 #ifndef M_PI
22     # define M_PI 3.14159265358979323846 /* pi */
23     #endif
24 woo409 757
25 gross 576 namespace escript {
26    
27    
28     /**
29     \brief
30     solves a 1x1 eigenvalue A*V=ev*V problem
31    
32     \param A00 Input - A_00
33     \param ev0 Output - eigenvalue
34     */
35     inline
36     void eigenvalues1(const double A00,double* ev0) {
37    
38 gross 580 *ev0=A00;
39 gross 576
40     }
41     /**
42     \brief
43     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A
44    
45     \param A00 Input - A_00
46     \param A01 Input - A_01
47     \param A11 Input - A_11
48     \param ev0 Output - smallest eigenvalue
49     \param ev1 Output - largest eigenvalue
50     */
51     inline
52 gross 583 void eigenvalues2(const double A00,const double A01,const double A11,
53 gross 576 double* ev0, double* ev1) {
54 gross 580 const register double trA=(A00+A11)/2.;
55     const register double A_00=A00-trA;
56     const register double A_11=A11-trA;
57     const register double s=sqrt(A01*A01-A_00*A_11);
58     *ev0=trA-s;
59     *ev1=trA+s;
60 gross 576 }
61     /**
62     \brief
63     solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A
64    
65     \param A00 Input - A_00
66     \param A01 Input - A_01
67     \param A02 Input - A_02
68     \param A11 Input - A_11
69     \param A12 Input - A_12
70     \param A22 Input - A_22
71     \param ev0 Output - smallest eigenvalue
72     \param ev1 Output - eigenvalue
73     \param ev2 Output - largest eigenvalue
74     */
75     inline
76     void eigenvalues3(const double A00, const double A01, const double A02,
77     const double A11, const double A12,
78     const double A22,
79     double* ev0, double* ev1,double* ev2) {
80    
81 gross 580 const register double trA=(A00+A11+A22)/3.;
82     const register double A_00=A00-trA;
83     const register double A_11=A11-trA;
84     const register double A_22=A22-trA;
85     const register double A01_2=A01*A01;
86     const register double A02_2=A02*A02;
87     const register double A12_2=A12*A12;
88     const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.;
89 gross 585 if (p<=0.) {
90     *ev2=trA;
91     *ev1=trA;
92     *ev0=trA;
93    
94     } else {
95     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);
96     const register double sq_p=sqrt(p/3.);
97     register double z=-q/(2*pow(sq_p,3));
98     if (z<-1.) {
99     z=-1.;
100     } else if (z>1.) {
101     z=1.;
102     }
103     const register double alpha_3=acos(z)/3.;
104     *ev2=trA+2.*sq_p*cos(alpha_3);
105     *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.);
106     *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.);
107 gross 580 }
108 gross 576 }
109 gross 583 /**
110     \brief
111     solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A
112    
113     \param A00 Input - A_00
114     \param ev0 Output - eigenvalue
115     \param V00 Output - eigenvector
116     \param tol Input - tolerance to identify to eigenvalues
117     */
118     inline
119     void eigenvalues_and_eigenvectors1(const double A00,double* ev0,double* V00,const double tol)
120     {
121     eigenvalues1(A00,ev0);
122     *V00=1.;
123     return;
124     }
125     /**
126     \brief
127     returns a non-zero vector in the kernel of [[A00,A01],[A01,A11]] assuming that the kernel dimension is at least 1.
128    
129     \param A00 Input - matrix component
130     \param A10 Input - matrix component
131     \param A01 Input - matrix component
132     \param A11 Input - matrix component
133     \param V0 Output - vector component
134     \param V1 Output - vector component
135     */
136     inline
137     void vectorInKernel2(const double A00,const double A10,const double A01,const double A11,
138     double* V0, double*V1)
139     {
140     register double absA00=fabs(A00);
141 gross 587 register double absA10=fabs(A10);
142 gross 583 register double absA01=fabs(A01);
143     register double absA11=fabs(A11);
144 gross 587 register double m=absA11>absA10 ? absA11 : absA10;
145     if (absA00>m || absA01>m) {
146     *V0=-A01;
147 gross 583 *V1=A00;
148     } else {
149     if (m<=0) {
150     *V0=1.;
151     *V1=0.;
152     } else {
153     *V0=A11;
154 gross 587 *V1=-A10;
155 gross 583 }
156     }
157     }
158     /**
159     \brief
160     returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]]
161     assuming that the kernel dimension is at least 1 and A00 is non zero.
162    
163     \param A00 Input - matrix component
164     \param A10 Input - matrix component
165     \param A20 Input - matrix component
166     \param A01 Input - matrix component
167     \param A11 Input - matrix component
168     \param A21 Input - matrix component
169     \param A02 Input - matrix component
170     \param A12 Input - matrix component
171     \param A22 Input - matrix component
172     \param V0 Output - vector component
173     \param V1 Output - vector component
174     \param V2 Output - vector component
175     */
176     inline
177     void vectorInKernel3__nonZeroA00(const double A00,const double A10,const double A20,
178     const double A01,const double A11,const double A21,
179     const double A02,const double A12,const double A22,
180     double* V0,double* V1,double* V2)
181     {
182     double TEMP0,TEMP1;
183     register const double I00=1./A00;
184     register const double IA10=I00*A10;
185     register const double IA20=I00*A20;
186 gross 588 vectorInKernel2(A11-IA10*A01,A12-IA10*A02,
187     A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1);
188 gross 583 *V0=-(A10*TEMP0+A20*TEMP1);
189     *V1=A00*TEMP0;
190     *V2=A00*TEMP1;
191     }
192    
193     /**
194     \brief
195     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
196     ordered by increasing value and eigen vectors are normalizeVector3d such that
197     length is zero and first non-zero component is positive.
198    
199     \param A00 Input - A_00
200     \param A01 Input - A_01
201     \param A11 Input - A_11
202     \param ev0 Output - smallest eigenvalue
203     \param ev1 Output - eigenvalue
204     \param V00 Output - eigenvector componenent coresponding to ev0
205     \param V10 Output - eigenvector componenent coresponding to ev0
206     \param V01 Output - eigenvector componenent coresponding to ev1
207     \param V11 Output - eigenvector componenent coresponding to ev1
208     \param tol Input - tolerance to identify to eigenvalues
209     */
210     inline
211     void eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11,
212     double* ev0, double* ev1,
213     double* V00, double* V10, double* V01, double* V11,
214     const double tol)
215     {
216     double TEMP0,TEMP1;
217     eigenvalues2(A00,A01,A11,ev0,ev1);
218     const register double absev0=fabs(*ev0);
219     const register double absev1=fabs(*ev1);
220     register double max_ev=absev0>absev1 ? absev0 : absev1;
221     if (fabs((*ev0)-(*ev1))<tol*max_ev) {
222     *V00=1.;
223     *V10=0.;
224     *V01=0.;
225     *V11=1.;
226     } else {
227     vectorInKernel2(A00-(*ev0),A01,A01,A11-(*ev0),&TEMP0,&TEMP1);
228     const register double scale=1./sqrt(TEMP0*TEMP0+TEMP1*TEMP1);
229     if (TEMP0<0.) {
230     *V00=-TEMP0*scale;
231     *V10=-TEMP1*scale;
232     if (TEMP1<0.) {
233 gross 587 *V01= *V10;
234 gross 583 *V11=-(*V00);
235     } else {
236     *V01=-(*V10);
237 gross 587 *V11= (*V10);
238 gross 583 }
239     } else if (TEMP0>0.) {
240     *V00=TEMP0*scale;
241     *V10=TEMP1*scale;
242     if (TEMP1<0.) {
243 gross 587 *V01=-(*V10);
244 gross 583 *V11= (*V00);
245     } else {
246 gross 587 *V01= (*V10);
247     *V11=-(*V00);
248 gross 583 }
249     } else {
250     *V00=0.;
251     *V10=1;
252     *V11=0.;
253     *V01=1.;
254     }
255     }
256     }
257     /**
258     \brief
259     nomalizes a 3-d vector such that length is one and first non-zero component is positive.
260    
261     \param V0 - vector componenent
262     \param V1 - vector componenent
263     \param V2 - vector componenent
264     */
265     inline
266     void normalizeVector3(double* V0,double* V1,double* V2)
267     {
268     register double s;
269     if (*V0>0) {
270     s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
271     *V0*=s;
272     *V1*=s;
273     *V2*=s;
274     } else if (*V0<0) {
275     s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2));
276     *V0*=s;
277     *V1*=s;
278     *V2*=s;
279     } else {
280     if (*V1>0) {
281     s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
282     *V1*=s;
283     *V2*=s;
284     } else if (*V1<0) {
285     s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2));
286     *V1*=s;
287     *V2*=s;
288     } else {
289     *V2=1.;
290     }
291     }
292     }
293     /**
294     \brief
295     solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are
296     ordered by increasing value and eigen vectors are normalizeVector3d such that
297     length is zero and first non-zero component is positive.
298    
299     \param A00 Input - A_00
300     \param A01 Input - A_01
301     \param A11 Input - A_11
302     \param ev0 Output - smallest eigenvalue
303     \param ev1 Output - eigenvalue
304     \param V00 Output - eigenvector componenent coresponding to ev0
305     \param V10 Output - eigenvector componenent coresponding to ev0
306     \param V01 Output - eigenvector componenent coresponding to ev1
307     \param V11 Output - eigenvector componenent coresponding to ev1
308     \param tol Input - tolerance to identify to eigenvalues
309     */
310     inline
311     void eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02,
312     const double A11, const double A12, const double A22,
313     double* ev0, double* ev1, double* ev2,
314     double* V00, double* V10, double* V20,
315     double* V01, double* V11, double* V21,
316     double* V02, double* V12, double* V22,
317     const double tol)
318     {
319     register const double absA01=fabs(A01);
320     register const double absA02=fabs(A02);
321     register const double m=absA01>absA02 ? absA01 : absA02;
322     if (m<=0) {
323     double TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1;
324     eigenvalues_and_eigenvectors2(A11,A12,A22,
325     &TEMP_EV0,&TEMP_EV1,
326     &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol);
327     if (A00<=TEMP_EV0) {
328     *V00=1.;
329     *V10=0.;
330     *V20=0.;
331     *V01=0.;
332     *V11=TEMP_V00;
333     *V21=TEMP_V10;
334     *V02=0.;
335     *V12=TEMP_V01;
336     *V22=TEMP_V11;
337     *ev0=A00;
338     *ev1=TEMP_EV0;
339     *ev2=TEMP_EV1;
340     } else if (A00>TEMP_EV1) {
341 gross 588 *V02=1.;
342 gross 583 *V12=0.;
343 gross 588 *V22=0.;
344     *V00=0.;
345     *V10=TEMP_V00;
346     *V20=TEMP_V10;
347     *V01=0.;
348     *V11=TEMP_V01;
349     *V21=TEMP_V11;
350 gross 583 *ev0=TEMP_EV0;
351     *ev1=TEMP_EV1;
352 gross 588 *ev2=A00;
353 gross 583 } else {
354 gross 588 *V01=1.;
355     *V11=0.;
356     *V21=0.;
357     *V00=0.;
358     *V10=TEMP_V00;
359 gross 583 *V20=TEMP_V10;
360 gross 588 *V02=0.;
361     *V12=TEMP_V01;
362 gross 583 *V22=TEMP_V11;
363     *ev0=TEMP_EV0;
364     *ev1=A00;
365     *ev2=TEMP_EV1;
366     }
367     } else {
368     eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2);
369     const register double absev0=fabs(*ev0);
370     const register double absev1=fabs(*ev1);
371     const register double absev2=fabs(*ev2);
372     register double max_ev=absev0>absev1 ? absev0 : absev1;
373     max_ev=max_ev>absev2 ? max_ev : absev2;
374     const register double d_01=fabs((*ev0)-(*ev1));
375     const register double d_12=fabs((*ev1)-(*ev2));
376     const register double max_d=d_01>d_12 ? d_01 : d_12;
377     if (max_d<=tol*max_ev) {
378     *V00=1.;
379     *V10=0;
380     *V20=0;
381     *V01=0;
382     *V11=1.;
383     *V21=0;
384     *V02=0;
385     *V12=0;
386     *V22=1.;
387     } else {
388     const register double S00=A00-(*ev0);
389     const register double absS00=fabs(S00);
390     if (fabs(S00)>m) {
391     vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20);
392     } else if (absA02<m) {
393     vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20);
394     } else {
395     vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20);
396     }
397     normalizeVector3(V00,V10,V20);;
398     const register double T00=A00-(*ev2);
399     const register double absT00=fabs(T00);
400     if (fabs(T00)>m) {
401     vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22);
402     } else if (absA02<m) {
403     vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22);
404     } else {
405     vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22);
406     }
407     const register double dot=(*V02)*(*V00)+(*V12)*(*V10)+(*V22)*(*V20);
408     *V02-=dot*(*V00);
409     *V12-=dot*(*V10);
410     *V22-=dot*(*V20);
411     normalizeVector3(V02,V12,V22);
412     *V01=(*V10)*(*V22)-(*V12)*(*V20);
413     *V11=(*V20)*(*V02)-(*V00)*(*V22);
414     *V21=(*V00)*(*V12)-(*V02)*(*V10);
415     normalizeVector3(V01,V11,V21);
416     }
417     }
418     }
419 ksteube 813
420     // General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR)
421     // SM is the product of the last axis_offset entries in arg_0.getShape().
422     inline
423     void matrix_matrix_product(const int SL, const int SM, const int SR, const double* A, const double* B, double* C, int transpose)
424     {
425     if (transpose == 0) {
426     for (int i=0; i<SL; i++) {
427     for (int j=0; j<SR; j++) {
428     double sum = 0.0;
429     for (int l=0; l<SM; l++) {
430     sum += A[i+SL*l] * B[l+SM*j];
431     }
432     C[i+SL*j] = sum;
433     }
434     }
435     }
436     else if (transpose == 1) {
437     for (int i=0; i<SL; i++) {
438     for (int j=0; j<SR; j++) {
439     double sum = 0.0;
440     for (int l=0; l<SM; l++) {
441     sum += A[i*SM+l] * B[l+SM*j];
442     }
443     C[i+SL*j] = sum;
444     }
445     }
446     }
447     else if (transpose == 2) {
448     for (int i=0; i<SL; i++) {
449     for (int j=0; j<SR; j++) {
450     double sum = 0.0;
451     for (int l=0; l<SM; l++) {
452     sum += A[i+SL*l] * B[l*SR+j];
453     }
454     C[i+SL*j] = sum;
455     }
456     }
457     }
458     }
459    
460 gross 576 } // end of namespace
461     #endif

  ViewVC Help
Powered by ViewVC 1.1.26