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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.576  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26