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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1946 - (hide annotations)
Wed Oct 29 05:48:53 2008 UTC (11 years, 4 months ago) by jfenwick
File MIME type: text/plain
File size: 15802 byte(s)
A cleanup of some of the problems I found doing a Wall compile.

Removed some commented out lines.
Swapped some member initialisers.
Removed virtual qualifiers from some methods in FunctionSpace.
Fixed some unused or (possibly) uninitialised variables.


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

  ViewVC Help
Powered by ViewVC 1.1.26