00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
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
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(1,"DGEESX Problem with parameter %d",INFO);
00078 } else if (INFO>0) {
00079 SETERRQ1(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
00102 if (db<dt) {
00103
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
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(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; PetscTruth 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 *dummy, PetscTruth *flg)
00197 {
00198 Mat A = (Mat)prob;
00199 PetscReal v=0; PetscTruth 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(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 *dummy, PetscTruth *flg)
00223 {
00224 Mat A = (Mat)prob;
00225 PetscTruth 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(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 *dummy,PetscTruth *flg)
00248 {
00249 Mat A = (Mat)prob;
00250 PetscTruth 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(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 *dummy,PetscTruth *flg)
00273 {
00274 Mat A = (Mat)prob;
00275 PetscTruth 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(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 *dummy,PetscTruth *flg)
00298 {
00299 Mat A = (Mat)prob;
00300 PetscTruth 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(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 *dummy,PetscTruth *flg)
00323 {
00324 Mat A = (Mat)prob;
00325 PetscTruth 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(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 *dummy,PetscTruth *flg)
00348 {
00349 Mat A = (Mat)prob;
00350 PetscTruth 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(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 *dummy,PetscTruth *flg)
00373 {
00374 Mat A = (Mat)prob;
00375 PetscTruth 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(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 *dummy,PetscTruth *flg)
00398 {
00399 Mat A = (Mat)prob;
00400 PetscTruth 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(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
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 }