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

Revision 1811 - (show annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 4 months ago) by ksteube
File MIME type: text/plain
File size: 15808 byte(s)
```Copyright updated in all files

```
 1 2 /******************************************************* 3 * 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 14 15 #if !defined escript_LocalOps_H 16 #define escript_LocalOps_H 17 #ifdef __INTEL_COMPILER 18 # include 19 #else 20 # include 21 #endif 22 #ifndef M_PI 23 # define M_PI 3.14159265358979323846 /* pi */ 24 #endif 25 26 namespace escript { 27 28 29 /** 30 \brief 31 solves a 1x1 eigenvalue A*V=ev*V problem 32 33 \param A00 Input - A_00 34 \param ev0 Output - eigenvalue 35 */ 36 inline 37 void eigenvalues1(const double A00,double* ev0) { 38 39 *ev0=A00; 40 41 } 42 /** 43 \brief 44 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A 45 46 \param A00 Input - A_00 47 \param A01 Input - A_01 48 \param A11 Input - A_11 49 \param ev0 Output - smallest eigenvalue 50 \param ev1 Output - largest eigenvalue 51 */ 52 inline 53 void eigenvalues2(const double A00,const double A01,const double A11, 54 double* ev0, double* ev1) { 55 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 } 62 /** 63 \brief 64 solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A 65 66 \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 \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 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 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 } 109 } 110 /** 111 \brief 112 solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A 113 114 \param A00 Input - A_00 115 \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 register double absA10=fabs(A10); 143 register double absA01=fabs(A01); 144 register double absA11=fabs(A11); 145 register double m=absA11>absA10 ? absA11 : absA10; 146 if (absA00>m || absA01>m) { 147 *V0=-A01; 148 *V1=A00; 149 } else { 150 if (m<=0) { 151 *V0=1.; 152 *V1=0.; 153 } else { 154 *V0=A11; 155 *V1=-A10; 156 } 157 } 158 } 159 /** 160 \brief 161 returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]] 162 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 vectorInKernel2(A11-IA10*A01,A12-IA10*A02, 188 A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1); 189 *V0=-(A10*TEMP0+A20*TEMP1); 190 *V1=A00*TEMP0; 191 *V2=A00*TEMP1; 192 } 193 194 /** 195 \brief 196 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 length is zero and first non-zero component is positive. 199 200 \param A00 Input - A_00 201 \param A01 Input - A_01 202 \param A11 Input - A_11 203 \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 const double tol) 216 { 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))0.) { 241 *V00=TEMP0*scale; 242 *V10=TEMP1*scale; 243 if (TEMP1<0.) { 244 *V01=-(*V10); 245 *V11= (*V00); 246 } else { 247 *V01= (*V10); 248 *V11=-(*V00); 249 } 250 } else { 251 *V00=0.; 252 *V10=1; 253 *V11=0.; 254 *V01=1.; 255 } 256 } 257 } 258 /** 259 \brief 260 nomalizes a 3-d vector such that length is one and first non-zero component is positive. 261 262 \param V0 - vector componenent 263 \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 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 length is zero and first non-zero component is positive. 299 300 \param A00 Input - A_00 301 \param A01 Input - A_01 302 \param A11 Input - A_11 303 \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 double* V00, double* V10, double* V20, 316 double* V01, double* V11, double* V21, 317 double* V02, double* V12, double* V22, 318 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 *V02=1.; 343 *V12=0.; 344 *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 *ev0=TEMP_EV0; 352 *ev1=TEMP_EV1; 353 *ev2=A00; 354 } else { 355 *V01=1.; 356 *V11=0.; 357 *V21=0.; 358 *V00=0.; 359 *V10=TEMP_V00; 360 *V20=TEMP_V10; 361 *V02=0.; 362 *V12=TEMP_V01; 363 *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 if (fabs(S00)>m) { 392 vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20); 393 } else if (absA02m) { 402 vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22); 403 } else if (absA02 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 template 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 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 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 } // end of namespace 513 #endif