ESYS13
Revision_
|
00001 00002 /******************************************************* 00003 * 00004 * Copyright (c) 2003-2012 by University of Queensland 00005 * Earth Systems Science Computational Center (ESSCC) 00006 * http://www.uq.edu.au/esscc 00007 * 00008 * Primary Business: Queensland, Australia 00009 * Licensed under the Open Software License version 3.0 00010 * http://www.opensource.org/licenses/osl-3.0.php 00011 * 00012 *******************************************************/ 00013 00014 00015 #if !defined escript_LocalOps_H 00016 #define escript_LocalOps_H 00017 #if defined(_WIN32) && defined(__INTEL_COMPILER) 00018 # include <mathimf.h> 00019 #else 00020 # include <math.h> 00021 #endif 00022 #ifndef M_PI 00023 # define M_PI 3.14159265358979323846 /* pi */ 00024 #endif 00025 00026 00035 namespace escript { 00036 00041 inline 00042 bool nancheck(double d) 00043 { 00044 // Q: so why not just test d!=d? 00045 // A: Coz it doesn't always work [I've checked]. 00046 // One theory is that the optimizer skips the test. 00047 #ifdef isnan 00048 return isnan(d); 00049 #elif defined _isnan 00050 return _isnan(d); 00051 #else 00052 return false; 00053 #endif 00054 } 00055 00060 inline 00061 double makeNaN() 00062 { 00063 #ifdef nan 00064 return nan(""); 00065 #else 00066 return sqrt(-1.); 00067 #endif 00068 00069 } 00070 00071 00079 inline 00080 void eigenvalues1(const double A00,double* ev0) { 00081 00082 *ev0=A00; 00083 00084 } 00095 inline 00096 void eigenvalues2(const double A00,const double A01,const double A11, 00097 double* ev0, double* ev1) { 00098 const register double trA=(A00+A11)/2.; 00099 const register double A_00=A00-trA; 00100 const register double A_11=A11-trA; 00101 const register double s=sqrt(A01*A01-A_00*A_11); 00102 *ev0=trA-s; 00103 *ev1=trA+s; 00104 } 00119 inline 00120 void eigenvalues3(const double A00, const double A01, const double A02, 00121 const double A11, const double A12, 00122 const double A22, 00123 double* ev0, double* ev1,double* ev2) { 00124 00125 const register double trA=(A00+A11+A22)/3.; 00126 const register double A_00=A00-trA; 00127 const register double A_11=A11-trA; 00128 const register double A_22=A22-trA; 00129 const register double A01_2=A01*A01; 00130 const register double A02_2=A02*A02; 00131 const register double A12_2=A12*A12; 00132 const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.; 00133 if (p<=0.) { 00134 *ev2=trA; 00135 *ev1=trA; 00136 *ev0=trA; 00137 00138 } else { 00139 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); 00140 const register double sq_p=sqrt(p/3.); 00141 register double z=-q/(2*pow(sq_p,3)); 00142 if (z<-1.) { 00143 z=-1.; 00144 } else if (z>1.) { 00145 z=1.; 00146 } 00147 const register double alpha_3=acos(z)/3.; 00148 *ev2=trA+2.*sq_p*cos(alpha_3); 00149 *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.); 00150 *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.); 00151 } 00152 } 00162 inline 00163 void eigenvalues_and_eigenvectors1(const double A00,double* ev0,double* V00,const double tol) 00164 { 00165 eigenvalues1(A00,ev0); 00166 *V00=1.; 00167 return; 00168 } 00180 inline 00181 void vectorInKernel2(const double A00,const double A10,const double A01,const double A11, 00182 double* V0, double*V1) 00183 { 00184 register double absA00=fabs(A00); 00185 register double absA10=fabs(A10); 00186 register double absA01=fabs(A01); 00187 register double absA11=fabs(A11); 00188 register double m=absA11>absA10 ? absA11 : absA10; 00189 if (absA00>m || absA01>m) { 00190 *V0=-A01; 00191 *V1=A00; 00192 } else { 00193 if (m<=0) { 00194 *V0=1.; 00195 *V1=0.; 00196 } else { 00197 *V0=A11; 00198 *V1=-A10; 00199 } 00200 } 00201 } 00220 inline 00221 void vectorInKernel3__nonZeroA00(const double A00,const double A10,const double A20, 00222 const double A01,const double A11,const double A21, 00223 const double A02,const double A12,const double A22, 00224 double* V0,double* V1,double* V2) 00225 { 00226 double TEMP0,TEMP1; 00227 register const double I00=1./A00; 00228 register const double IA10=I00*A10; 00229 register const double IA20=I00*A20; 00230 vectorInKernel2(A11-IA10*A01,A12-IA10*A02, 00231 A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1); 00232 *V0=-(A10*TEMP0+A20*TEMP1); 00233 *V1=A00*TEMP0; 00234 *V2=A00*TEMP1; 00235 } 00236 00254 inline 00255 void eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11, 00256 double* ev0, double* ev1, 00257 double* V00, double* V10, double* V01, double* V11, 00258 const double tol) 00259 { 00260 double TEMP0,TEMP1; 00261 eigenvalues2(A00,A01,A11,ev0,ev1); 00262 const register double absev0=fabs(*ev0); 00263 const register double absev1=fabs(*ev1); 00264 register double max_ev=absev0>absev1 ? absev0 : absev1; 00265 if (fabs((*ev0)-(*ev1))<tol*max_ev) { 00266 *V00=1.; 00267 *V10=0.; 00268 *V01=0.; 00269 *V11=1.; 00270 } else { 00271 vectorInKernel2(A00-(*ev0),A01,A01,A11-(*ev0),&TEMP0,&TEMP1); 00272 const register double scale=1./sqrt(TEMP0*TEMP0+TEMP1*TEMP1); 00273 if (TEMP0<0.) { 00274 *V00=-TEMP0*scale; 00275 *V10=-TEMP1*scale; 00276 if (TEMP1<0.) { 00277 *V01= *V10; 00278 *V11=-(*V00); 00279 } else { 00280 *V01=-(*V10); 00281 *V11= (*V00); 00282 } 00283 } else if (TEMP0>0.) { 00284 *V00=TEMP0*scale; 00285 *V10=TEMP1*scale; 00286 if (TEMP1<0.) { 00287 *V01=-(*V10); 00288 *V11= (*V00); 00289 } else { 00290 *V01= (*V10); 00291 *V11=-(*V00); 00292 } 00293 } else { 00294 *V00=0.; 00295 *V10=1; 00296 *V11=0.; 00297 *V01=1.; 00298 } 00299 } 00300 } 00309 inline 00310 void normalizeVector3(double* V0,double* V1,double* V2) 00311 { 00312 register double s; 00313 if (*V0>0) { 00314 s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2)); 00315 *V0*=s; 00316 *V1*=s; 00317 *V2*=s; 00318 } else if (*V0<0) { 00319 s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2)); 00320 *V0*=s; 00321 *V1*=s; 00322 *V2*=s; 00323 } else { 00324 if (*V1>0) { 00325 s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2)); 00326 *V1*=s; 00327 *V2*=s; 00328 } else if (*V1<0) { 00329 s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2)); 00330 *V1*=s; 00331 *V2*=s; 00332 } else { 00333 *V2=1.; 00334 } 00335 } 00336 } 00363 inline 00364 void eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02, 00365 const double A11, const double A12, const double A22, 00366 double* ev0, double* ev1, double* ev2, 00367 double* V00, double* V10, double* V20, 00368 double* V01, double* V11, double* V21, 00369 double* V02, double* V12, double* V22, 00370 const double tol) 00371 { 00372 register const double absA01=fabs(A01); 00373 register const double absA02=fabs(A02); 00374 register const double m=absA01>absA02 ? absA01 : absA02; 00375 if (m<=0) { 00376 double TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1; 00377 eigenvalues_and_eigenvectors2(A11,A12,A22, 00378 &TEMP_EV0,&TEMP_EV1, 00379 &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol); 00380 if (A00<=TEMP_EV0) { 00381 *V00=1.; 00382 *V10=0.; 00383 *V20=0.; 00384 *V01=0.; 00385 *V11=TEMP_V00; 00386 *V21=TEMP_V10; 00387 *V02=0.; 00388 *V12=TEMP_V01; 00389 *V22=TEMP_V11; 00390 *ev0=A00; 00391 *ev1=TEMP_EV0; 00392 *ev2=TEMP_EV1; 00393 } else if (A00>TEMP_EV1) { 00394 *V02=1.; 00395 *V12=0.; 00396 *V22=0.; 00397 *V00=0.; 00398 *V10=TEMP_V00; 00399 *V20=TEMP_V10; 00400 *V01=0.; 00401 *V11=TEMP_V01; 00402 *V21=TEMP_V11; 00403 *ev0=TEMP_EV0; 00404 *ev1=TEMP_EV1; 00405 *ev2=A00; 00406 } else { 00407 *V01=1.; 00408 *V11=0.; 00409 *V21=0.; 00410 *V00=0.; 00411 *V10=TEMP_V00; 00412 *V20=TEMP_V10; 00413 *V02=0.; 00414 *V12=TEMP_V01; 00415 *V22=TEMP_V11; 00416 *ev0=TEMP_EV0; 00417 *ev1=A00; 00418 *ev2=TEMP_EV1; 00419 } 00420 } else { 00421 eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2); 00422 const register double absev0=fabs(*ev0); 00423 const register double absev1=fabs(*ev1); 00424 const register double absev2=fabs(*ev2); 00425 register double max_ev=absev0>absev1 ? absev0 : absev1; 00426 max_ev=max_ev>absev2 ? max_ev : absev2; 00427 const register double d_01=fabs((*ev0)-(*ev1)); 00428 const register double d_12=fabs((*ev1)-(*ev2)); 00429 const register double max_d=d_01>d_12 ? d_01 : d_12; 00430 if (max_d<=tol*max_ev) { 00431 *V00=1.; 00432 *V10=0; 00433 *V20=0; 00434 *V01=0; 00435 *V11=1.; 00436 *V21=0; 00437 *V02=0; 00438 *V12=0; 00439 *V22=1.; 00440 } else { 00441 const register double S00=A00-(*ev0); 00442 const register double absS00=fabs(S00); 00443 if (absS00>m) { 00444 vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20); 00445 } else if (absA02<m) { 00446 vectorInKernel3__nonZeroA00(A01,A11-(*ev0),A12,S00,A01,A02,A02,A12,A22-(*ev0),V00,V10,V20); 00447 } else { 00448 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev0),S00,A01,A02,A01,A11-(*ev0),A12,V00,V10,V20); 00449 } 00450 normalizeVector3(V00,V10,V20);; 00451 const register double T00=A00-(*ev2); 00452 const register double absT00=fabs(T00); 00453 if (absT00>m) { 00454 vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22); 00455 } else if (absA02<m) { 00456 vectorInKernel3__nonZeroA00(A01,A11-(*ev2),A12,T00,A01,A02,A02,A12,A22-(*ev2),V02,V12,V22); 00457 } else { 00458 vectorInKernel3__nonZeroA00(A02,A12,A22-(*ev2),T00,A01,A02,A01,A11-(*ev2),A12,V02,V12,V22); 00459 } 00460 const register double dot=(*V02)*(*V00)+(*V12)*(*V10)+(*V22)*(*V20); 00461 *V02-=dot*(*V00); 00462 *V12-=dot*(*V10); 00463 *V22-=dot*(*V20); 00464 normalizeVector3(V02,V12,V22); 00465 *V01=(*V10)*(*V22)-(*V12)*(*V20); 00466 *V11=(*V20)*(*V02)-(*V00)*(*V22); 00467 *V21=(*V00)*(*V12)-(*V02)*(*V10); 00468 normalizeVector3(V01,V11,V21); 00469 } 00470 } 00471 } 00472 00473 // General tensor product: arg_2(SL x SR) = arg_0(SL x SM) * arg_1(SM x SR) 00474 // SM is the product of the last axis_offset entries in arg_0.getShape(). 00475 inline 00476 void matrix_matrix_product(const int SL, const int SM, const int SR, const double* A, const double* B, double* C, int transpose) 00477 { 00478 if (transpose == 0) { 00479 for (int i=0; i<SL; i++) { 00480 for (int j=0; j<SR; j++) { 00481 double sum = 0.0; 00482 for (int l=0; l<SM; l++) { 00483 sum += A[i+SL*l] * B[l+SM*j]; 00484 } 00485 C[i+SL*j] = sum; 00486 } 00487 } 00488 } 00489 else if (transpose == 1) { 00490 for (int i=0; i<SL; i++) { 00491 for (int j=0; j<SR; j++) { 00492 double sum = 0.0; 00493 for (int l=0; l<SM; l++) { 00494 sum += A[i*SM+l] * B[l+SM*j]; 00495 } 00496 C[i+SL*j] = sum; 00497 } 00498 } 00499 } 00500 else if (transpose == 2) { 00501 for (int i=0; i<SL; i++) { 00502 for (int j=0; j<SR; j++) { 00503 double sum = 0.0; 00504 for (int l=0; l<SM; l++) { 00505 sum += A[i+SL*l] * B[l*SR+j]; 00506 } 00507 C[i+SL*j] = sum; 00508 } 00509 } 00510 } 00511 } 00512 00513 template <typename UnaryFunction> 00514 inline void tensor_unary_operation(const int size, 00515 const double *arg1, 00516 double * argRes, 00517 UnaryFunction operation) 00518 { 00519 for (int i = 0; i < size; ++i) { 00520 argRes[i] = operation(arg1[i]); 00521 } 00522 return; 00523 } 00524 00525 template <typename BinaryFunction> 00526 inline void tensor_binary_operation(const int size, 00527 const double *arg1, 00528 const double *arg2, 00529 double * argRes, 00530 BinaryFunction operation) 00531 { 00532 for (int i = 0; i < size; ++i) { 00533 argRes[i] = operation(arg1[i], arg2[i]); 00534 } 00535 return; 00536 } 00537 00538 template <typename BinaryFunction> 00539 inline void tensor_binary_operation(const int size, 00540 double arg1, 00541 const double *arg2, 00542 double *argRes, 00543 BinaryFunction operation) 00544 { 00545 for (int i = 0; i < size; ++i) { 00546 argRes[i] = operation(arg1, arg2[i]); 00547 } 00548 return; 00549 } 00550 00551 template <typename BinaryFunction> 00552 inline void tensor_binary_operation(const int size, 00553 const double *arg1, 00554 double arg2, 00555 double *argRes, 00556 BinaryFunction operation) 00557 { 00558 for (int i = 0; i < size; ++i) { 00559 argRes[i] = operation(arg1[i], arg2); 00560 } 00561 return; 00562 } 00563 00564 } // end of namespace 00565 #endif