SALSA Analysis Modules
|
00001 /*! \file lapack.c \ingroup categories 00002 \brief Dense routines from Lapack: eigenvalue and schur stuff 00003 00004 \section lapack Dense calculations from Lapack 00005 00006 \section Usage 00007 00008 Activate this module with 00009 00010 PetscErrorCode RegisterLapackModules(); 00011 00012 Compute these elements with 00013 00014 ComputeQuantity("lapack",<element>,A,(AnalysisItem*)&res,&flg); 00015 00016 Available elements are: 00017 - something 00018 */ 00019 00020 #include <stdlib.h> 00021 #include "anamod.h" 00022 #include "anamodsalsamodules.h" 00023 #include "petscerror.h" 00024 #include "petscmat.h" 00025 00026 extern "C" { 00027 extern void LAPACK_DGEESX(const char* jobvs, const char* sort, 00028 char* sel, 00029 const char* sense, const int* n, double* a, 00030 const int* lda, int* sdim, double* wr, double* wi, 00031 double* vs, const int* ldvs, double* rconde, 00032 double* rcondv, double* work, const int* lwork, 00033 int* iwork, const int* liwork, int* bwork, int* info); 00034 } 00035 00036 #undef __FUNCT__ 00037 #define __FUNCT__ "eigenvaluecomp" 00038 /*! 00039 This is the computational routine for lapack eigenvalue calculation. 00040 */ 00041 static PetscErrorCode eigenvaluecomp(Mat A) 00042 { 00043 MPI_Comm comm; Mat Adense; 00044 int maxbyreset=0,maxbyimset=0,maxbymagset=0,minbymagset=0; 00045 PetscReal dep=0, 00046 maxbymagre,maxbymagim, minbymagre,minbymagim, 00047 maxbyrere,maxbyreim, maxbyimre,maxbyimim, 00048 maxmag,minmag; 00049 PetscErrorCode ierr; 00050 00051 PetscFunctionBegin; 00052 00053 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00054 ierr = MatConvert(A,MATSEQDENSE,MAT_INITIAL_MATRIX,&Adense); CHKERRQ(ierr); 00055 00056 { 00057 char JOBVS = 'N', SORT = 'N', SENSE = 'N', SELECT='X'; 00058 int N,LDA,SDIM=0,LDVS=1,LWORK,LIWORK,*IWORK,INFO; 00059 PetscScalar *AA,*VS=NULL,*RCONDE=NULL,*RCONDV=NULL; 00060 PetscReal *WR,*WI,*WORK,*BWORK=NULL; 00061 00062 ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr); 00063 LDA = N; 00064 ierr = MatGetArray(Adense,&AA); CHKERRQ(ierr); 00065 ierr = PetscMalloc(N*sizeof(PetscReal),&WR); CHKERRQ(ierr); 00066 ierr = PetscMalloc(N*sizeof(PetscReal),&WI); CHKERRQ(ierr); 00067 LWORK = 3*N+1; 00068 ierr = PetscMalloc(LWORK*sizeof(PetscReal),&WORK); CHKERRQ(ierr); 00069 LIWORK = 10; 00070 ierr = PetscMalloc(LIWORK*sizeof(int),&IWORK); CHKERRQ(ierr); 00071 CHKMEMQ; 00072 LAPACK_DGEESX(&JOBVS, &SORT, &SELECT, &SENSE, &N, AA, &LDA, &SDIM, 00073 WR, WI, VS, &LDVS, RCONDE, RCONDV, WORK, &LWORK, 00074 IWORK, &LIWORK, (int*)BWORK, &INFO); 00075 if (INFO<0) { 00076 INFO = -INFO; 00077 SETERRQ1(MPI_COMM_WORLD,1,"DGEESX Problem with parameter %d",INFO); 00078 } else if (INFO>0) { 00079 SETERRQ1(MPI_COMM_WORLD,1,"DGEESX error %d",INFO); 00080 } 00081 CHKMEMQ; 00082 { 00083 int i,j,rind=0; PetscReal f=0; 00084 for (i=0; i<N; i++) { 00085 #define ABS(x) (x>0?x:-x) 00086 #define EIGENVALUE(re,im) { \ 00087 PetscReal mag = sqrt(re*re+im*im); \ 00088 if (!maxbyreset) {maxbyreset = 1; maxbyrere = re; maxbyreim = im; } \ 00089 else if (re>maxbyrere) { maxbyrere = re; maxbyreim = im; } \ 00090 if (!maxbyimset) {maxbyimset = 1; maxbyimre = re; maxbyimim = im; } \ 00091 else if (im>maxbyimre) { maxbyimre = re; maxbyimim = im; } \ 00092 if (!maxbymagset) {maxbymagset = 1; maxmag = mag; \ 00093 maxbymagre = re; maxbymagim = im; } \ 00094 else if (mag>maxmag) { maxbymagre = re; maxbymagim = im; } \ 00095 if (!minbymagset) {minbymagset = 1; minmag = mag; \ 00096 minbymagre = re; minbymagim = im; } \ 00097 else if (mag<minmag) { minbymagre = re; minbymagim = im; } \ 00098 } 00099 PetscReal dd = AA[rind],db = ABS(AA[rind+1]),dt; 00100 dt = 1.e-14 * ABS(dd); 00101 /*printf("row %d, diagonal %e, below %e, test %e\n",i,dd,db,dt);*/ 00102 if (db<dt) { 00103 /* case of 1x1 eigenvalue block */ 00104 int ind=rind+N; 00105 for (j=i+1; j<N; j++) { 00106 PetscReal aa = AA[ind]; 00107 f += aa*aa; ind += N; 00108 } 00109 EIGENVALUE(dd,0.); 00110 } else { 00111 /* case of 2x2 eigenvalue block */ 00112 int ind; 00113 PetscReal dd = AA[rind],do1=AA[rind+1],do2=AA[rind+N],dn=AA[rind+N+1]; 00114 if (ABS((dd-dn)/dd)>1.e-13) 00115 SETERRQ4(PETSC_COMM_SELF,1,"Strange 2x2 block %e,%e,%e,%e\n",dd,do1,do2,dn); 00116 EIGENVALUE(dd,sqrt(do1*do2)); 00117 ind = rind+2*N; 00118 for (j=i+2; j<N; j++) { 00119 PetscReal aa = AA[ind]; 00120 f += aa*aa; ind += N; 00121 } 00122 i++; rind += N+1; ind = rind+N; 00123 for (j=i+1; j<N; j++) { 00124 PetscReal aa = AA[ind]; 00125 f += aa*aa; ind += N; 00126 } 00127 } 00128 rind += N+1; 00129 } 00130 dep = sqrt(f); 00131 } 00132 ierr = PetscFree(WR); CHKERRQ(ierr); 00133 ierr = PetscFree(WI); CHKERRQ(ierr); 00134 } 00135 00136 ierr = MatDestroy(&Adense); CHKERRQ(ierr); 00137 00138 { 00139 int id; PetscBool has; 00140 ierr = GetDataID("lapack","departure",&id,&has); CHKERRQ(ierr); 00141 HASTOEXIST(has); 00142 ierr = PetscObjectComposedDataSetScalar 00143 ((PetscObject)A,id,dep); CHKERRQ(ierr); 00144 00145 ierr = GetDataID 00146 ("lapack","lambda-max-by-magnitude-re",&id,&has); CHKERRQ(ierr); 00147 HASTOEXIST(has); 00148 ierr = PetscObjectComposedDataSetReal 00149 ((PetscObject)A,id,maxbymagre); CHKERRQ(ierr); 00150 ierr = GetDataID 00151 ("lapack","lambda-max-by-magnitude-im",&id,&has); CHKERRQ(ierr); 00152 HASTOEXIST(has); 00153 ierr = PetscObjectComposedDataSetReal 00154 ((PetscObject)A,id,maxbymagim); CHKERRQ(ierr); 00155 00156 ierr = GetDataID 00157 ("lapack","lambda-min-by-magnitude-re",&id,&has); CHKERRQ(ierr); 00158 HASTOEXIST(has); 00159 ierr = PetscObjectComposedDataSetReal 00160 ((PetscObject)A,id,minbymagre); CHKERRQ(ierr); 00161 ierr = GetDataID 00162 ("lapack","lambda-min-by-magnitude-im",&id,&has); CHKERRQ(ierr); 00163 HASTOEXIST(has); 00164 ierr = PetscObjectComposedDataSetReal 00165 ((PetscObject)A,id,minbymagim); CHKERRQ(ierr); 00166 00167 ierr = GetDataID 00168 ("lapack","lambda-max-by-real-part-re",&id,&has); CHKERRQ(ierr); 00169 HASTOEXIST(has); 00170 ierr = PetscObjectComposedDataSetReal 00171 ((PetscObject)A,id,maxbyrere); CHKERRQ(ierr); 00172 ierr = GetDataID 00173 ("lapack","lambda-max-by-real-part-im",&id,&has); CHKERRQ(ierr); 00174 HASTOEXIST(has); 00175 ierr = PetscObjectComposedDataSetReal 00176 ((PetscObject)A,id,maxbyreim); CHKERRQ(ierr); 00177 00178 ierr = GetDataID 00179 ("lapack","lambda-max-by-im-part-re",&id,&has); CHKERRQ(ierr); 00180 HASTOEXIST(has); 00181 ierr = PetscObjectComposedDataSetReal 00182 ((PetscObject)A,id,maxbyimre); CHKERRQ(ierr); 00183 ierr = GetDataID 00184 ("lapack","lambda-max-by-im-part-im",&id,&has); CHKERRQ(ierr); 00185 HASTOEXIST(has); 00186 ierr = PetscObjectComposedDataSetReal 00187 ((PetscObject)A,id,maxbyimim); CHKERRQ(ierr); 00188 } 00189 00190 PetscFunctionReturn(0); 00191 } 00192 00193 #undef __FUNCT__ 00194 #define __FUNCT__ "Departure" 00195 static PetscErrorCode Departure 00196 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00197 { 00198 Mat A = (Mat)prob; 00199 PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; 00200 PetscFunctionBegin; 00201 ierr = GetDataID("lapack","departure",&id,&has); CHKERRQ(ierr); 00202 HASTOEXIST(has); 00203 ierr = PetscObjectComposedDataGetScalar 00204 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00205 if (*flg) { 00206 rv->r = v; 00207 } else { 00208 ierr = eigenvaluecomp(A); CHKERRQ(ierr); 00209 ierr = PetscObjectComposedDataGetScalar 00210 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00211 if (*flg) { 00212 rv->r = v; 00213 } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed departure by now"); 00214 } 00215 00216 PetscFunctionReturn(0); 00217 } 00218 00219 #undef __FUNCT__ 00220 #define __FUNCT__ "MaxEVbyMagRe" 00221 static PetscErrorCode MaxEVbyMagRe 00222 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00223 { 00224 Mat A = (Mat)prob; 00225 PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; 00226 PetscFunctionBegin; 00227 ierr = GetDataID 00228 ("lapack","lambda-max-by-magnitude-re",&id,&has); CHKERRQ(ierr); 00229 HASTOEXIST(has); 00230 ierr = PetscObjectComposedDataGetReal 00231 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00232 if (*flg) rv->r = v; 00233 else { 00234 ierr = eigenvaluecomp(A); CHKERRQ(ierr); 00235 ierr = PetscObjectComposedDataGetScalar 00236 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00237 if (*flg) { 00238 rv->r = v; 00239 } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-magnitude-re by now"); 00240 } 00241 PetscFunctionReturn(0); 00242 } 00243 00244 #undef __FUNCT__ 00245 #define __FUNCT__ "MaxEVbyMagIm" 00246 static PetscErrorCode MaxEVbyMagIm 00247 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00248 { 00249 Mat A = (Mat)prob; 00250 PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; 00251 PetscFunctionBegin; 00252 ierr = GetDataID 00253 ("lapack","lambda-max-by-magnitude-im",&id,&has); CHKERRQ(ierr); 00254 HASTOEXIST(has); 00255 ierr = PetscObjectComposedDataGetReal 00256 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00257 if (*flg) rv->r = v; 00258 else { 00259 ierr = eigenvaluecomp(A); CHKERRQ(ierr); 00260 ierr = PetscObjectComposedDataGetScalar 00261 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00262 if (*flg) { 00263 rv->r = v; 00264 } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-magnitude-im by now"); 00265 } 00266 PetscFunctionReturn(0); 00267 } 00268 00269 #undef __FUNCT__ 00270 #define __FUNCT__ "MinEVbyMagRe" 00271 static PetscErrorCode MinEVbyMagRe 00272 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00273 { 00274 Mat A = (Mat)prob; 00275 PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; 00276 PetscFunctionBegin; 00277 ierr = GetDataID 00278 ("lapack","lambda-min-by-magnitude-re",&id,&has); CHKERRQ(ierr); 00279 HASTOEXIST(has); 00280 ierr = PetscObjectComposedDataGetReal 00281 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00282 if (*flg) rv->r = v; 00283 else { 00284 ierr = eigenvaluecomp(A); CHKERRQ(ierr); 00285 ierr = PetscObjectComposedDataGetScalar 00286 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00287 if (*flg) { 00288 rv->r = v; 00289 } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-min-by-magnitude-re by now"); 00290 } 00291 PetscFunctionReturn(0); 00292 } 00293 00294 #undef __FUNCT__ 00295 #define __FUNCT__ "MinEVbyMagIm" 00296 static PetscErrorCode MinEVbyMagIm 00297 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00298 { 00299 Mat A = (Mat)prob; 00300 PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; 00301 PetscFunctionBegin; 00302 ierr = GetDataID 00303 ("lapack","lambda-min-by-magnitude-im",&id,&has); CHKERRQ(ierr); 00304 HASTOEXIST(has); 00305 ierr = PetscObjectComposedDataGetReal 00306 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00307 if (*flg) rv->r = v; 00308 else { 00309 ierr = eigenvaluecomp(A); CHKERRQ(ierr); 00310 ierr = PetscObjectComposedDataGetScalar 00311 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00312 if (*flg) { 00313 rv->r = v; 00314 } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-min-by-magnitude-im by now"); 00315 } 00316 PetscFunctionReturn(0); 00317 } 00318 00319 #undef __FUNCT__ 00320 #define __FUNCT__ "MaxEVbyRealRe" 00321 static PetscErrorCode MaxEVbyRealRe 00322 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00323 { 00324 Mat A = (Mat)prob; 00325 PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; 00326 PetscFunctionBegin; 00327 ierr = GetDataID 00328 ("lapack","lambda-max-by-real-part-re",&id,&has); CHKERRQ(ierr); 00329 HASTOEXIST(has); 00330 ierr = PetscObjectComposedDataGetReal 00331 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00332 if (*flg) rv->r = v; 00333 else { 00334 ierr = eigenvaluecomp(A); CHKERRQ(ierr); 00335 ierr = PetscObjectComposedDataGetScalar 00336 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00337 if (*flg) { 00338 rv->r = v; 00339 } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-real-part-re by now"); 00340 } 00341 PetscFunctionReturn(0); 00342 } 00343 00344 #undef __FUNCT__ 00345 #define __FUNCT__ "MaxEVbyRealIm" 00346 static PetscErrorCode MaxEVbyRealIm 00347 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00348 { 00349 Mat A = (Mat)prob; 00350 PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; 00351 PetscFunctionBegin; 00352 ierr = GetDataID 00353 ("lapack","lambda-max-by-real-part-im",&id,&has); CHKERRQ(ierr); 00354 HASTOEXIST(has); 00355 ierr = PetscObjectComposedDataGetReal 00356 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00357 if (*flg) rv->r = v; 00358 else { 00359 ierr = eigenvaluecomp(A); CHKERRQ(ierr); 00360 ierr = PetscObjectComposedDataGetScalar 00361 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00362 if (*flg) { 00363 rv->r = v; 00364 } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-real-part-im by now"); 00365 } 00366 PetscFunctionReturn(0); 00367 } 00368 00369 #undef __FUNCT__ 00370 #define __FUNCT__ "MaxEVbyImRe" 00371 static PetscErrorCode MaxEVbyImRe 00372 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00373 { 00374 Mat A = (Mat)prob; 00375 PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; 00376 PetscFunctionBegin; 00377 ierr = GetDataID 00378 ("lapack","lambda-max-by-im-part-re",&id,&has); CHKERRQ(ierr); 00379 HASTOEXIST(has); 00380 ierr = PetscObjectComposedDataGetReal 00381 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00382 if (*flg) rv->r = v; 00383 else { 00384 ierr = eigenvaluecomp(A); CHKERRQ(ierr); 00385 ierr = PetscObjectComposedDataGetScalar 00386 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00387 if (*flg) { 00388 rv->r = v; 00389 } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-im-part-re by now"); 00390 } 00391 PetscFunctionReturn(0); 00392 } 00393 00394 #undef __FUNCT__ 00395 #define __FUNCT__ "MaxEVbyImIm" 00396 static PetscErrorCode MaxEVbyImIm 00397 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00398 { 00399 Mat A = (Mat)prob; 00400 PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0; 00401 PetscFunctionBegin; 00402 ierr = GetDataID 00403 ("lapack","lambda-max-by-im-part-im",&id,&has); CHKERRQ(ierr); 00404 HASTOEXIST(has); 00405 ierr = PetscObjectComposedDataGetReal 00406 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00407 if (*flg) rv->r = v; 00408 else { 00409 ierr = eigenvaluecomp(A); CHKERRQ(ierr); 00410 ierr = PetscObjectComposedDataGetScalar 00411 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00412 if (*flg) { 00413 rv->r = v; 00414 } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed lambda-max-by-im-part-im by now"); 00415 } 00416 PetscFunctionReturn(0); 00417 } 00418 00419 #undef __FUNCT__ 00420 #define __FUNCT__ "RegisterLapackModules" 00421 /*! 00422 Declare norm-like modules that can be performed with lapack calculations. 00423 */ 00424 PetscErrorCode RegisterLapackModules() 00425 { 00426 PetscErrorCode ierr; 00427 PetscFunctionBegin; 00428 00429 ierr = RegisterModule 00430 ("lapack","departure",ANALYSISDOUBLE,&Departure); CHKERRQ(ierr); 00431 00432 ierr = RegisterModule 00433 ("lapack","lambda-max-by-magnitude-re", 00434 ANALYSISDOUBLE,&MaxEVbyMagRe); CHKERRQ(ierr); 00435 ierr = RegisterModule 00436 ("lapack","lambda-max-by-magnitude-im", 00437 ANALYSISDOUBLE,&MaxEVbyMagIm); CHKERRQ(ierr); 00438 ierr = RegisterModule 00439 ("lapack","lambda-min-by-magnitude-re", 00440 ANALYSISDOUBLE,&MinEVbyMagRe); CHKERRQ(ierr); 00441 ierr = RegisterModule 00442 ("lapack","lambda-min-by-magnitude-im", 00443 ANALYSISDOUBLE,&MinEVbyMagIm); CHKERRQ(ierr); 00444 ierr = RegisterModule 00445 ("lapack","lambda-max-by-real-part-re", 00446 ANALYSISDOUBLE,&MaxEVbyRealRe); CHKERRQ(ierr); 00447 ierr = RegisterModule 00448 ("lapack","lambda-max-by-real-part-im", 00449 ANALYSISDOUBLE,&MaxEVbyRealIm); CHKERRQ(ierr); 00450 ierr = RegisterModule 00451 ("lapack","lambda-max-by-im-part-re", 00452 ANALYSISDOUBLE,&MaxEVbyImRe); CHKERRQ(ierr); 00453 ierr = RegisterModule 00454 ("lapack","lambda-max-by-im-part-im", 00455 ANALYSISDOUBLE,&MaxEVbyImIm); CHKERRQ(ierr); 00456 00457 PetscFunctionReturn(0); 00458 }