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

Revision 588 - (hide annotations)
Fri Mar 10 04:45:04 2006 UTC (13 years, 6 months ago) by gross
File MIME type: text/plain
File size: 14044 byte(s)
```1D and 3D tests for eigenvalues_and_eigenvector added.
```
 1 gross 576 // \$Id\$ 2 /* 3 ****************************************************************************** 4 * * 5 * COPYRIGHT ACcESS 2004 - All Rights Reserved * 6 * * 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 * 9 * consent of ACcESS. Copying, use or modification of this software * 10 * by any unauthorised person is illegal unless that person has a software * 11 * license agreement with ACcESS. * 12 * * 13 ****************************************************************************** 14 */ 15 16 #if !defined escript_LocalOps_H 17 #define escript_LocalOps_H 18 gross 580 #include 19 gross 576 namespace escript { 20 21 22 /** 23 \brief 24 solves a 1x1 eigenvalue A*V=ev*V problem 25 26 \param A00 Input - A_00 27 \param ev0 Output - eigenvalue 28 */ 29 inline 30 void eigenvalues1(const double A00,double* ev0) { 31 32 gross 580 *ev0=A00; 33 gross 576 34 } 35 /** 36 \brief 37 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A 38 39 \param A00 Input - A_00 40 \param A01 Input - A_01 41 \param A11 Input - A_11 42 \param ev0 Output - smallest eigenvalue 43 \param ev1 Output - largest eigenvalue 44 */ 45 inline 46 gross 583 void eigenvalues2(const double A00,const double A01,const double A11, 47 gross 576 double* ev0, double* ev1) { 48 gross 580 const register double trA=(A00+A11)/2.; 49 const register double A_00=A00-trA; 50 const register double A_11=A11-trA; 51 const register double s=sqrt(A01*A01-A_00*A_11); 52 *ev0=trA-s; 53 *ev1=trA+s; 54 gross 576 } 55 /** 56 \brief 57 solves a 3x3 eigenvalue A*V=ev*V problem for symmetric A 58 59 \param A00 Input - A_00 60 \param A01 Input - A_01 61 \param A02 Input - A_02 62 \param A11 Input - A_11 63 \param A12 Input - A_12 64 \param A22 Input - A_22 65 \param ev0 Output - smallest eigenvalue 66 \param ev1 Output - eigenvalue 67 \param ev2 Output - largest eigenvalue 68 */ 69 inline 70 void eigenvalues3(const double A00, const double A01, const double A02, 71 const double A11, const double A12, 72 const double A22, 73 double* ev0, double* ev1,double* ev2) { 74 75 gross 580 const register double trA=(A00+A11+A22)/3.; 76 const register double A_00=A00-trA; 77 const register double A_11=A11-trA; 78 const register double A_22=A22-trA; 79 const register double A01_2=A01*A01; 80 const register double A02_2=A02*A02; 81 const register double A12_2=A12*A12; 82 const register double p=A02_2+A12_2+A01_2+(A_00*A_00+A_11*A_11+A_22*A_22)/2.; 83 gross 585 if (p<=0.) { 84 *ev2=trA; 85 *ev1=trA; 86 *ev0=trA; 87 88 } else { 89 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); 90 const register double sq_p=sqrt(p/3.); 91 register double z=-q/(2*pow(sq_p,3)); 92 if (z<-1.) { 93 z=-1.; 94 } else if (z>1.) { 95 z=1.; 96 } 97 const register double alpha_3=acos(z)/3.; 98 *ev2=trA+2.*sq_p*cos(alpha_3); 99 *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.); 100 *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.); 101 gross 580 } 102 gross 576 } 103 gross 583 /** 104 \brief 105 solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A 106 107 \param A00 Input - A_00 108 \param ev0 Output - eigenvalue 109 \param V00 Output - eigenvector 110 \param tol Input - tolerance to identify to eigenvalues 111 */ 112 inline 113 void eigenvalues_and_eigenvectors1(const double A00,double* ev0,double* V00,const double tol) 114 { 115 eigenvalues1(A00,ev0); 116 *V00=1.; 117 return; 118 } 119 /** 120 \brief 121 returns a non-zero vector in the kernel of [[A00,A01],[A01,A11]] assuming that the kernel dimension is at least 1. 122 123 \param A00 Input - matrix component 124 \param A10 Input - matrix component 125 \param A01 Input - matrix component 126 \param A11 Input - matrix component 127 \param V0 Output - vector component 128 \param V1 Output - vector component 129 */ 130 inline 131 void vectorInKernel2(const double A00,const double A10,const double A01,const double A11, 132 double* V0, double*V1) 133 { 134 register double absA00=fabs(A00); 135 gross 587 register double absA10=fabs(A10); 136 gross 583 register double absA01=fabs(A01); 137 register double absA11=fabs(A11); 138 gross 587 register double m=absA11>absA10 ? absA11 : absA10; 139 if (absA00>m || absA01>m) { 140 *V0=-A01; 141 gross 583 *V1=A00; 142 } else { 143 if (m<=0) { 144 *V0=1.; 145 *V1=0.; 146 } else { 147 *V0=A11; 148 gross 587 *V1=-A10; 149 gross 583 } 150 } 151 } 152 /** 153 \brief 154 returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]] 155 assuming that the kernel dimension is at least 1 and A00 is non zero. 156 157 \param A00 Input - matrix component 158 \param A10 Input - matrix component 159 \param A20 Input - matrix component 160 \param A01 Input - matrix component 161 \param A11 Input - matrix component 162 \param A21 Input - matrix component 163 \param A02 Input - matrix component 164 \param A12 Input - matrix component 165 \param A22 Input - matrix component 166 \param V0 Output - vector component 167 \param V1 Output - vector component 168 \param V2 Output - vector component 169 */ 170 inline 171 void vectorInKernel3__nonZeroA00(const double A00,const double A10,const double A20, 172 const double A01,const double A11,const double A21, 173 const double A02,const double A12,const double A22, 174 double* V0,double* V1,double* V2) 175 { 176 double TEMP0,TEMP1; 177 register const double I00=1./A00; 178 register const double IA10=I00*A10; 179 register const double IA20=I00*A20; 180 gross 588 vectorInKernel2(A11-IA10*A01,A12-IA10*A02, 181 A21-IA20*A01,A22-IA20*A02,&TEMP0,&TEMP1); 182 gross 583 *V0=-(A10*TEMP0+A20*TEMP1); 183 *V1=A00*TEMP0; 184 *V2=A00*TEMP1; 185 } 186 187 /** 188 \brief 189 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are 190 ordered by increasing value and eigen vectors are normalizeVector3d such that 191 length is zero and first non-zero component is positive. 192 193 \param A00 Input - A_00 194 \param A01 Input - A_01 195 \param A11 Input - A_11 196 \param ev0 Output - smallest eigenvalue 197 \param ev1 Output - eigenvalue 198 \param V00 Output - eigenvector componenent coresponding to ev0 199 \param V10 Output - eigenvector componenent coresponding to ev0 200 \param V01 Output - eigenvector componenent coresponding to ev1 201 \param V11 Output - eigenvector componenent coresponding to ev1 202 \param tol Input - tolerance to identify to eigenvalues 203 */ 204 inline 205 void eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11, 206 double* ev0, double* ev1, 207 double* V00, double* V10, double* V01, double* V11, 208 const double tol) 209 { 210 double TEMP0,TEMP1; 211 eigenvalues2(A00,A01,A11,ev0,ev1); 212 const register double absev0=fabs(*ev0); 213 const register double absev1=fabs(*ev1); 214 register double max_ev=absev0>absev1 ? absev0 : absev1; 215 if (fabs((*ev0)-(*ev1))0.) { 234 *V00=TEMP0*scale; 235 *V10=TEMP1*scale; 236 if (TEMP1<0.) { 237 gross 587 *V01=-(*V10); 238 gross 583 *V11= (*V00); 239 } else { 240 gross 587 *V01= (*V10); 241 *V11=-(*V00); 242 gross 583 } 243 } else { 244 *V00=0.; 245 *V10=1; 246 *V11=0.; 247 *V01=1.; 248 } 249 } 250 } 251 /** 252 \brief 253 nomalizes a 3-d vector such that length is one and first non-zero component is positive. 254 255 \param V0 - vector componenent 256 \param V1 - vector componenent 257 \param V2 - vector componenent 258 */ 259 inline 260 void normalizeVector3(double* V0,double* V1,double* V2) 261 { 262 register double s; 263 if (*V0>0) { 264 s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2)); 265 *V0*=s; 266 *V1*=s; 267 *V2*=s; 268 } else if (*V0<0) { 269 s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2)); 270 *V0*=s; 271 *V1*=s; 272 *V2*=s; 273 } else { 274 if (*V1>0) { 275 s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2)); 276 *V1*=s; 277 *V2*=s; 278 } else if (*V1<0) { 279 s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2)); 280 *V1*=s; 281 *V2*=s; 282 } else { 283 *V2=1.; 284 } 285 } 286 } 287 /** 288 \brief 289 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are 290 ordered by increasing value and eigen vectors are normalizeVector3d such that 291 length is zero and first non-zero component is positive. 292 293 \param A00 Input - A_00 294 \param A01 Input - A_01 295 \param A11 Input - A_11 296 \param ev0 Output - smallest eigenvalue 297 \param ev1 Output - eigenvalue 298 \param V00 Output - eigenvector componenent coresponding to ev0 299 \param V10 Output - eigenvector componenent coresponding to ev0 300 \param V01 Output - eigenvector componenent coresponding to ev1 301 \param V11 Output - eigenvector componenent coresponding to ev1 302 \param tol Input - tolerance to identify to eigenvalues 303 */ 304 inline 305 void eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02, 306 const double A11, const double A12, const double A22, 307 double* ev0, double* ev1, double* ev2, 308 double* V00, double* V10, double* V20, 309 double* V01, double* V11, double* V21, 310 double* V02, double* V12, double* V22, 311 const double tol) 312 { 313 register const double absA01=fabs(A01); 314 register const double absA02=fabs(A02); 315 register const double m=absA01>absA02 ? absA01 : absA02; 316 if (m<=0) { 317 double TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1; 318 eigenvalues_and_eigenvectors2(A11,A12,A22, 319 &TEMP_EV0,&TEMP_EV1, 320 &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol); 321 if (A00<=TEMP_EV0) { 322 *V00=1.; 323 *V10=0.; 324 *V20=0.; 325 *V01=0.; 326 *V11=TEMP_V00; 327 *V21=TEMP_V10; 328 *V02=0.; 329 *V12=TEMP_V01; 330 *V22=TEMP_V11; 331 *ev0=A00; 332 *ev1=TEMP_EV0; 333 *ev2=TEMP_EV1; 334 } else if (A00>TEMP_EV1) { 335 gross 588 *V02=1.; 336 gross 583 *V12=0.; 337 gross 588 *V22=0.; 338 *V00=0.; 339 *V10=TEMP_V00; 340 *V20=TEMP_V10; 341 *V01=0.; 342 *V11=TEMP_V01; 343 *V21=TEMP_V11; 344 gross 583 *ev0=TEMP_EV0; 345 *ev1=TEMP_EV1; 346 gross 588 *ev2=A00; 347 gross 583 } else { 348 gross 588 *V01=1.; 349 *V11=0.; 350 *V21=0.; 351 *V00=0.; 352 *V10=TEMP_V00; 353 gross 583 *V20=TEMP_V10; 354 gross 588 *V02=0.; 355 *V12=TEMP_V01; 356 gross 583 *V22=TEMP_V11; 357 *ev0=TEMP_EV0; 358 *ev1=A00; 359 *ev2=TEMP_EV1; 360 } 361 } else { 362 eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2); 363 const register double absev0=fabs(*ev0); 364 const register double absev1=fabs(*ev1); 365 const register double absev2=fabs(*ev2); 366 register double max_ev=absev0>absev1 ? absev0 : absev1; 367 max_ev=max_ev>absev2 ? max_ev : absev2; 368 const register double d_01=fabs((*ev0)-(*ev1)); 369 const register double d_12=fabs((*ev1)-(*ev2)); 370 const register double max_d=d_01>d_12 ? d_01 : d_12; 371 if (max_d<=tol*max_ev) { 372 *V00=1.; 373 *V10=0; 374 *V20=0; 375 *V01=0; 376 *V11=1.; 377 *V21=0; 378 *V02=0; 379 *V12=0; 380 *V22=1.; 381 } else { 382 const register double S00=A00-(*ev0); 383 const register double absS00=fabs(S00); 384 if (fabs(S00)>m) { 385 vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20); 386 } else if (absA02m) { 395 vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22); 396 } else if (absA02