SALSA Analysis Modules
lapack.c
Go to the documentation of this file.
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 }