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

Revision 583 - (hide annotations)
Wed Mar 8 08:15:34 2006 UTC (16 years, 11 months ago) by gross
File MIME type: text/plain
File size: 13882 byte(s)
```_eigenvalues_and_eigenvector method added of data object. the algorithm has been tested on floats in python but not on data objects.
```
 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 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); 84 const register double sq_p=sqrt(p/3.); 85 register double z=-q/(2*pow(sq_p,3)); 86 if (z<-1.) { 87 z=-1.; 88 } else if (z>1.) { 89 z=1.; 90 } 91 const register double alpha_3=acos(z)/3.; 92 *ev2=trA+2.*sq_p*cos(alpha_3); 93 *ev1=trA-2.*sq_p*cos(alpha_3+M_PI/3.); 94 *ev0=trA-2.*sq_p*cos(alpha_3-M_PI/3.); 95 gross 576 } 96 gross 583 /** 97 \brief 98 solves a 1x1 eigenvalue A*V=ev*V problem for symmetric A 99 100 \param A00 Input - A_00 101 \param ev0 Output - eigenvalue 102 \param V00 Output - eigenvector 103 \param tol Input - tolerance to identify to eigenvalues 104 */ 105 inline 106 void eigenvalues_and_eigenvectors1(const double A00,double* ev0,double* V00,const double tol) 107 { 108 eigenvalues1(A00,ev0); 109 *V00=1.; 110 return; 111 } 112 /** 113 \brief 114 returns a non-zero vector in the kernel of [[A00,A01],[A01,A11]] assuming that the kernel dimension is at least 1. 115 116 \param A00 Input - matrix component 117 \param A10 Input - matrix component 118 \param A01 Input - matrix component 119 \param A11 Input - matrix component 120 \param V0 Output - vector component 121 \param V1 Output - vector component 122 */ 123 inline 124 void vectorInKernel2(const double A00,const double A10,const double A01,const double A11, 125 double* V0, double*V1) 126 { 127 register double absA00=fabs(A00); 128 register double absA01=fabs(A01); 129 register double absA10=fabs(A10); 130 register double absA11=fabs(A11); 131 register double m=absA11>absA01 ? absA11 : absA01; 132 if (absA00>m || absA10>m) { 133 *V0=-A10; 134 *V1=A00; 135 } else { 136 if (m<=0) { 137 *V0=1.; 138 *V1=0.; 139 } else { 140 *V0=A11; 141 *V1=-A01; 142 } 143 } 144 } 145 /** 146 \brief 147 returns a non-zero vector in the kernel of [[A00,A01,A02],[A10,A11,A12],[A20,A21,A22]] 148 assuming that the kernel dimension is at least 1 and A00 is non zero. 149 150 \param A00 Input - matrix component 151 \param A10 Input - matrix component 152 \param A20 Input - matrix component 153 \param A01 Input - matrix component 154 \param A11 Input - matrix component 155 \param A21 Input - matrix component 156 \param A02 Input - matrix component 157 \param A12 Input - matrix component 158 \param A22 Input - matrix component 159 \param V0 Output - vector component 160 \param V1 Output - vector component 161 \param V2 Output - vector component 162 */ 163 inline 164 void vectorInKernel3__nonZeroA00(const double A00,const double A10,const double A20, 165 const double A01,const double A11,const double A21, 166 const double A02,const double A12,const double A22, 167 double* V0,double* V1,double* V2) 168 { 169 double TEMP0,TEMP1; 170 register const double I00=1./A00; 171 register const double IA10=I00*A10; 172 register const double IA20=I00*A20; 173 vectorInKernel2(A11-IA10*A01,A21-IA20*A01,A12-IA10*A02,A22-IA20*A02,&TEMP0,&TEMP1); 174 *V0=-(A10*TEMP0+A20*TEMP1); 175 *V1=A00*TEMP0; 176 *V2=A00*TEMP1; 177 } 178 179 /** 180 \brief 181 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are 182 ordered by increasing value and eigen vectors are normalizeVector3d such that 183 length is zero and first non-zero component is positive. 184 185 \param A00 Input - A_00 186 \param A01 Input - A_01 187 \param A11 Input - A_11 188 \param ev0 Output - smallest eigenvalue 189 \param ev1 Output - eigenvalue 190 \param V00 Output - eigenvector componenent coresponding to ev0 191 \param V10 Output - eigenvector componenent coresponding to ev0 192 \param V01 Output - eigenvector componenent coresponding to ev1 193 \param V11 Output - eigenvector componenent coresponding to ev1 194 \param tol Input - tolerance to identify to eigenvalues 195 */ 196 inline 197 void eigenvalues_and_eigenvectors2(const double A00,const double A01,const double A11, 198 double* ev0, double* ev1, 199 double* V00, double* V10, double* V01, double* V11, 200 const double tol) 201 { 202 double TEMP0,TEMP1; 203 eigenvalues2(A00,A01,A11,ev0,ev1); 204 const register double absev0=fabs(*ev0); 205 const register double absev1=fabs(*ev1); 206 register double max_ev=absev0>absev1 ? absev0 : absev1; 207 if (fabs((*ev0)-(*ev1))0.) { 226 *V00=TEMP0*scale; 227 *V10=TEMP1*scale; 228 if (TEMP1<0.) { 229 *V01=-(*V10); 230 *V11= (*V00); 231 } else { 232 *V01= (*V10); 233 *V11=-(*V00); 234 } 235 } else { 236 *V00=0.; 237 *V10=1; 238 *V11=0.; 239 *V01=1.; 240 } 241 } 242 } 243 /** 244 \brief 245 nomalizes a 3-d vector such that length is one and first non-zero component is positive. 246 247 \param V0 - vector componenent 248 \param V1 - vector componenent 249 \param V2 - vector componenent 250 */ 251 inline 252 void normalizeVector3(double* V0,double* V1,double* V2) 253 { 254 register double s; 255 if (*V0>0) { 256 s=1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2)); 257 *V0*=s; 258 *V1*=s; 259 *V2*=s; 260 } else if (*V0<0) { 261 s=-1./sqrt((*V0)*(*V0)+(*V1)*(*V1)+(*V2)*(*V2)); 262 *V0*=s; 263 *V1*=s; 264 *V2*=s; 265 } else { 266 if (*V1>0) { 267 s=1./sqrt((*V1)*(*V1)+(*V2)*(*V2)); 268 *V1*=s; 269 *V2*=s; 270 } else if (*V1<0) { 271 s=-1./sqrt((*V1)*(*V1)+(*V2)*(*V2)); 272 *V1*=s; 273 *V2*=s; 274 } else { 275 *V2=1.; 276 } 277 } 278 } 279 /** 280 \brief 281 solves a 2x2 eigenvalue A*V=ev*V problem for symmetric A. Eigenvectors are 282 ordered by increasing value and eigen vectors are normalizeVector3d such that 283 length is zero and first non-zero component is positive. 284 285 \param A00 Input - A_00 286 \param A01 Input - A_01 287 \param A11 Input - A_11 288 \param ev0 Output - smallest eigenvalue 289 \param ev1 Output - eigenvalue 290 \param V00 Output - eigenvector componenent coresponding to ev0 291 \param V10 Output - eigenvector componenent coresponding to ev0 292 \param V01 Output - eigenvector componenent coresponding to ev1 293 \param V11 Output - eigenvector componenent coresponding to ev1 294 \param tol Input - tolerance to identify to eigenvalues 295 */ 296 inline 297 void eigenvalues_and_eigenvectors3(const double A00, const double A01, const double A02, 298 const double A11, const double A12, const double A22, 299 double* ev0, double* ev1, double* ev2, 300 double* V00, double* V10, double* V20, 301 double* V01, double* V11, double* V21, 302 double* V02, double* V12, double* V22, 303 const double tol) 304 { 305 register const double absA01=fabs(A01); 306 register const double absA02=fabs(A02); 307 register const double m=absA01>absA02 ? absA01 : absA02; 308 if (m<=0) { 309 double TEMP_V00,TEMP_V10,TEMP_V01,TEMP_V11,TEMP_EV0,TEMP_EV1; 310 eigenvalues_and_eigenvectors2(A11,A12,A22, 311 &TEMP_EV0,&TEMP_EV1, 312 &TEMP_V00,&TEMP_V10,&TEMP_V01,&TEMP_V11,tol); 313 if (A00<=TEMP_EV0) { 314 *V00=1.; 315 *V10=0.; 316 *V20=0.; 317 *V01=0.; 318 *V11=TEMP_V00; 319 *V21=TEMP_V10; 320 *V02=0.; 321 *V12=TEMP_V01; 322 *V22=TEMP_V11; 323 *ev0=A00; 324 *ev1=TEMP_EV0; 325 *ev2=TEMP_EV1; 326 } else if (A00>TEMP_EV1) { 327 *V00=TEMP_V00; 328 *V10=TEMP_V10; 329 *V20=0.; 330 *V01=TEMP_V01; 331 *V11=TEMP_V11; 332 *V21=0.; 333 *V02=0.; 334 *V12=0.; 335 *V22=1.; 336 *ev0=TEMP_EV0; 337 *ev1=TEMP_EV1; 338 *ev0=A00; 339 } else { 340 *V00=TEMP_V00; 341 *V10=0; 342 *V20=TEMP_V10; 343 *V01=0.; 344 *V11=1.; 345 *V21=0.; 346 *V02=TEMP_V01; 347 *V12=0.; 348 *V22=TEMP_V11; 349 *ev0=TEMP_EV0; 350 *ev1=A00; 351 *ev2=TEMP_EV1; 352 } 353 } else { 354 eigenvalues3(A00,A01,A02,A11,A12,A22,ev0,ev1,ev2); 355 const register double absev0=fabs(*ev0); 356 const register double absev1=fabs(*ev1); 357 const register double absev2=fabs(*ev2); 358 register double max_ev=absev0>absev1 ? absev0 : absev1; 359 max_ev=max_ev>absev2 ? max_ev : absev2; 360 const register double d_01=fabs((*ev0)-(*ev1)); 361 const register double d_12=fabs((*ev1)-(*ev2)); 362 const register double max_d=d_01>d_12 ? d_01 : d_12; 363 if (max_d<=tol*max_ev) { 364 *V00=1.; 365 *V10=0; 366 *V20=0; 367 *V01=0; 368 *V11=1.; 369 *V21=0; 370 *V02=0; 371 *V12=0; 372 *V22=1.; 373 } else { 374 const register double S00=A00-(*ev0); 375 const register double absS00=fabs(S00); 376 if (fabs(S00)>m) { 377 vectorInKernel3__nonZeroA00(S00,A01,A02,A01,A11-(*ev0),A12,A02,A12,A22-(*ev0),V00,V10,V20); 378 } else if (absA02m) { 387 vectorInKernel3__nonZeroA00(T00,A01,A02,A01,A11-(*ev2),A12,A02,A12,A22-(*ev2),V02,V12,V22); 388 } else if (absA02