SALSA Analysis Modules
spectrum.c
Go to the documentation of this file.
00001 /*! \file spectrum.c \ingroup categories 
00002   \brief Various estimates of spectrum related quantities
00003 
00004   \section spectrum Spectral properties 
00005 
00006   The spectrum (eigenvalues and singular values) are both very
00007   informative about a matrix, and hard to compute. The spectrum module
00008   provides various estimates that are derived from running the system
00009   through a GMRES iteration for a few dozen iterations.
00010 
00011   \subsection Usage
00012 
00013   Activate this module with
00014 
00015   PetscErrorCode RegisterSpectrumModules();
00016 
00017   Compute these elements with
00018 
00019   ComputeQuantity("spectrum",<element>,A,(AnalysisItem*)&res,&reslen,&flg);
00020 
00021   These elements can be computed for both a preconditioned and 
00022   unpreconditioned matrix. Switch between the two with
00023 
00024   PetscErrorCode SpectrumComputePreconditionedSpectrum();
00025   PetscErrorCode SpectrumComputeUnpreconditionedSpectrum();
00026 
00027   In the preconditioned case a preconditioner has to have been defined
00028   in the Petsc options database with prefix "-ana"; for instance
00029   "-ana_pc_type jacobi".
00030 
00031   Available elements are:
00032   - "n-ritz-values" : number of stored Ritz values
00033   - "ritz-values-r" : real parts of stored Ritz values
00034   - "ritz-values-c" : complex parts of stored Ritz values
00035   - "ellipse-ax" : size of \f$x\f$-axis of the enclosing ellipse
00036   - "ellipse-ay" : size of \f$y\f$-axis of the enclosing ellipse
00037   - "ellipse-cx" : \f$x\f$-coordinate of the center of the enclosing ellipse
00038   - "ellipse-cy" : \f$y\f$-coordinate of the center of the enclosing ellipse
00039   - "kappa" : estimated condition number
00040   - "positive-fraction" : fraction of computed eigenvalues that has
00041     positive real part
00042   - "sigma-max" : maximum singular value
00043   - "sigma-min" : minimum singular value
00044   - "lambda-max-by-magnitude-re" : real part of maximum lambda by magnitude
00045   - "lambda-max-by-magnitude-im" : imaginary part of maximum lambda by magnitude
00046   - "lambda-min-by-magnitude-re" : real part of minimum lambda by magnitude
00047   - "lambda-min-by-magnitude-im" : imaginary part of minimum lambda by magnitude
00048   - "lambda-max-by-real-part-re" : real part of maximum lambda by real-part
00049   - "lambda-max-by-real-part-im" : imaginary part of maximum lambda by real-part
00050 
00051   \subsection Theory
00052 
00053   Jacobi methods, such as GMRES, give a reduction $AV=VH$ of $A$ to
00054   upper Hessenberg form $H$. The reduction is defined by the Krylov
00055   basis $V$. While a full reduction gives a Hessenberg matrix that has
00056   the exact same eigenvalues as $A$, using $V$ with a limited number
00057   of columns $n<N$) gives an $H$ with eigenvalues that approximate
00058   those of $A$. In particular, the outer eigenvalues are known to be
00059   fairly accurate even for modest values of $n$.
00060 
00061   This process is used in the Petsc calls KSPSetComputeEigenvalues()
00062   and KSPSetComputeSingularValues(), which we use in the
00063   compute_eigenvalues() function. Alternatively, if the SLEPc library is 
00064   available, that can be used for a similar but probably more accurate 
00065   computation.
00066 
00067   \subsection Tracing
00068 
00069   See SpectrumOptions().
00070 
00071 */
00072 #include <stdlib.h>
00073 #include "string.h"
00074 #include "anamod.h"
00075 #include "anamodsalsamodules.h"
00076 #include "petscerror.h"
00077 #ifdef PETSC_INTERNALS
00078 // #include "src/ksp/ksp/impls/gmres/gmresp.h"
00079 #include "src/ksp/ksp/impls/gmres/gmresimpl.h"
00080 #endif
00081 #include "petscksp.h"
00082 #include "petscmat.h"
00083 #if defined(HAVE_SLEPC)
00084 #include "slepc.h"
00085 #include "slepceps.h"
00086 #include "slepcsvd.h"
00087 #endif
00088 
00089 static PetscBool trace_spectrum = PETSC_FALSE,
00090   trace_hessenberg = PETSC_FALSE;
00091 
00092 static int use_prec = 0;
00093 #undef __FUNCT__
00094 #define __FUNCT__ "SpectrumComputePreconditionedSpectrum"
00095 PetscErrorCode SpectrumComputePreconditionedSpectrum()
00096 {
00097   PetscFunctionBegin;
00098   use_prec = 1;
00099   PetscFunctionReturn(0);
00100 }
00101 
00102 #undef __FUNCT__
00103 #define __FUNCT__ "SpectrumComputeUnpreconditionedSpectrum"
00104 PetscErrorCode SpectrumComputeUnpreconditionedSpectrum()
00105 {
00106   PetscFunctionBegin;
00107   use_prec = 0;
00108   PetscFunctionReturn(0);
00109 }
00110 
00111 #undef __FUNCT__
00112 #define __FUNCT__ "fivestepplan"
00113 static PetscErrorCode fivestepplan
00114 (KSP solver,PetscInt it,PetscReal err,KSPConvergedReason *reason,void *ctx)
00115 {
00116   PetscErrorCode ierr;
00117   PetscFunctionBegin;
00118   if (it<5) {
00119     if (!it) {
00120       ierr = KSPDefaultConverged(solver,it,err,reason,ctx); CHKERRQ(ierr);
00121     }
00122     *reason = KSP_CONVERGED_ITERATING;
00123   } else {
00124     ierr = KSPDefaultConverged(solver,it,err,reason,ctx); CHKERRQ(ierr);
00125   }
00126   PetscFunctionReturn(0);
00127 }
00128 
00129 /*
00130  * Ritz values and ellipse
00131  */
00132 #if defined(HAVE_SLEPC)
00133 #undef __FUNCT__
00134 #define __FUNCT__ "compute_eigenvalues_slepc"
00135 static PetscErrorCode compute_eigenvalues_slepc
00136 (Mat A,PetscReal **rr,PetscReal **rc,int *neig,
00137  PetscBool *success,char **reason)
00138 {
00139   EPS eps; MPI_Comm comm;
00140   EPSWhich whichs[] = {EPS_LARGEST_MAGNITUDE,EPS_SMALLEST_MAGNITUDE,EPS_LARGEST_REAL,EPS_SMALLEST_REAL,EPS_LARGEST_IMAGINARY,EPS_SMALLEST_IMAGINARY};
00141   Vec xi,xr; PetscReal *r,*c;
00142   int M,iwhich,nwhich=6,aimconv,spaceconv,evspace,loc; PetscErrorCode ierr;
00143 
00144   PetscFunctionBegin;
00145   printf("invoking slepc\n");
00146   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00147   ierr = MatGetSize(A,&M,PETSC_NULL); CHKERRQ(ierr);
00148   ierr = EPSCreate(comm,&eps); CHKERRQ(ierr);
00149   if (use_prec) {
00150     Mat B;
00151     ierr = PetscObjectQuery
00152       ((PetscObject)A,"approximation",(PetscObject*)&B); CHKERRQ(ierr);
00153     if (B) {
00154       ierr = EPSSetProblemType(eps,EPS_GNHEP); CHKERRQ(ierr);
00155       ierr = EPSSetOperators(eps,A,B); CHKERRQ(ierr);
00156       goto prec;
00157     }
00158   }
00159   ierr = EPSSetProblemType(eps,EPS_NHEP); CHKERRQ(ierr);
00160   ierr = EPSSetType(eps,EPSKRYLOVSCHUR); CHKERRQ(ierr);
00161   ierr = EPSSetOperators(eps,A,PETSC_NULL); CHKERRQ(ierr);
00162   {
00163     int maxit;
00164     ierr = EPSGetTolerances(eps,PETSC_NULL,&maxit); CHKERRQ(ierr);
00165     maxit *= M/5; if (maxit<100) maxit=100;
00166     printf("maxit=%d\n",maxit);
00167     ierr = EPSSetTolerances(eps,1.e-2,maxit); CHKERRQ(ierr);
00168   }
00169  prec:
00170 
00171   ierr = MatGetVecs(A,PETSC_NULL,&xr); CHKERRQ(ierr);
00172   ierr = MatGetVecs(A,PETSC_NULL,&xi); CHKERRQ(ierr);
00173   /*
00174    * allocate space for the eigenvalues of all problems
00175    * we are going to run
00176    */
00177   aimconv = 6; spaceconv = (int)(5*aimconv);
00178   {
00179     char epstype[20]; PetscBool flg;
00180     ierr = PetscOptionsGetString
00181       (PETSC_NULL,"-eps_type",epstype,20,&flg); CHKERRQ(ierr);
00182     /*
00183       Several space values if lapack is specified through a runtime option
00184     */
00185     if ( flg && strcmp(epstype,EPSLAPACK)==0 ) {
00186       printf("Lapack eps!\n");
00187       evspace = M+1; nwhich = 1;
00188     } else
00189       evspace = nwhich*spaceconv;
00190   }
00191   ierr = EPSSetDimensions(eps,aimconv,spaceconv,spaceconv); CHKERRQ(ierr);
00192   ierr = PetscMalloc((evspace+1)*sizeof(PetscReal),&r); CHKERRQ(ierr);
00193   ierr = PetscMalloc((evspace+1)*sizeof(PetscReal),&c); CHKERRQ(ierr);
00194   loc = 1;
00195   /* allow slepc to fail. maybe */
00196   /*  ierr = PetscPushErrorHandler(PetscReturnErrorHandler,NULL); CHKERRQ(ierr);*/
00197   for (iwhich=0; iwhich<nwhich; iwhich++) {
00198     EPSWhich which = whichs[iwhich]; int nconv,i;
00199     ierr = EPSSetWhichEigenpairs(eps,which); CHKERRQ(ierr);
00200     ierr = EPSSetFromOptions(eps); CHKERRQ(ierr);
00201     ierr = EPSSolve(eps); 
00202     if (ierr) {
00203       printf(">>>> Error %d in EPSSolve\n",ierr);
00204       PetscFunctionReturn(ierr);
00205     }
00206     ierr = EPSGetConverged(eps,&nconv); CHKERRQ(ierr);
00207     if (nconv>0) {
00208       for (i=0; i<nconv; i++) {
00209         if (loc>=evspace+1) SETERRQ(MPI_COMM_WORLD,1,"Eigenvalue space overflow");
00210         ierr = EPSGetEigenpair(eps,i,r+loc,c+loc,xr,xi);CHKERRQ(ierr);
00211         loc++;
00212       }
00213     }
00214   }
00215   /*  ierr = PetscPopErrorHandler(); CHKERRQ(ierr);*/
00216   ierr = VecDestroy(&xr); CHKERRQ(ierr);
00217   ierr = VecDestroy(&xi); CHKERRQ(ierr);
00218   ierr = EPSDestroy(&eps); CHKERRQ(ierr);
00219   loc--; *neig =  loc;
00220   if (loc>0) {
00221     *rr = r; *rc = c;
00222     *success = PETSC_TRUE;
00223     *reason = "normal convergence";
00224   } else {
00225     ierr = PetscFree(r); CHKERRQ(ierr);
00226     ierr = PetscFree(c); CHKERRQ(ierr);
00227     *rr = NULL; *rc = NULL;
00228     *success = PETSC_FALSE;
00229     *reason = "slepc iteration failed to converge";
00230   }
00231   PetscFunctionReturn(0);
00232 }
00233 #undef __FUNCT__
00234 #define __FUNCT__ "compute_sigmas_slepc"
00235 static PetscErrorCode compute_sigmas_slepc
00236 (Mat A,PetscReal *sx,PetscReal *sm,PetscBool *success,char **reason)
00237 {
00238   SVD svd; PetscReal sigmamax,sigmamin;
00239   int nconv; MPI_Comm comm; PetscErrorCode ierr;
00240 
00241   PetscFunctionBegin;
00242 
00243   if (use_prec) {
00244     *success = PETSC_FALSE;
00245     *reason = "can not compute singular values of preconditioned system";
00246     PetscFunctionReturn(0);
00247   }
00248 
00249   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00250   ierr = SVDCreate(PETSC_COMM_WORLD,&svd);CHKERRQ(ierr);
00251   ierr = SVDSetOperator(svd,A);CHKERRQ(ierr);
00252   {
00253     int maxit,M;
00254     ierr = MatGetSize(A,&M,PETSC_NULL); CHKERRQ(ierr);
00255     ierr = SVDGetTolerances(svd,PETSC_NULL,&maxit); CHKERRQ(ierr);
00256     maxit = M/3; if (maxit<60) maxit = 60;
00257     ierr = SVDSetTolerances(svd,.1,maxit); CHKERRQ(ierr);
00258   }
00259   ierr = SVDSetDimensions(svd,1,20,20);CHKERRQ(ierr);
00260   ierr = SVDSetType(svd,SVDTRLANCZOS); CHKERRQ(ierr);
00261 
00262   ierr = SVDSetWhichSingularTriplets(svd,SVD_LARGEST);CHKERRQ(ierr);
00263   ierr = SVDSetFromOptions(svd); CHKERRQ(ierr);
00264   ierr = SVDSolve(svd);CHKERRQ(ierr);
00265   ierr = SVDGetConverged(svd,&nconv);CHKERRQ(ierr);
00266   if (nconv > 0) {
00267     ierr = SVDGetSingularTriplet
00268       (svd,0,&sigmamax,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
00269   } else {
00270     *success = PETSC_FALSE; *reason = "slepc no convergence on sigma max";
00271     goto cleanup;
00272   }
00273 
00274   ierr = SVDSetWhichSingularTriplets(svd,SVD_SMALLEST); CHKERRQ(ierr);
00275   ierr = SVDSetFromOptions(svd); CHKERRQ(ierr);
00276   ierr = SVDSolve(svd); CHKERRQ(ierr);
00277   ierr = SVDGetConverged(svd,&nconv); CHKERRQ(ierr);
00278   if (nconv > 0) {
00279     ierr = SVDGetSingularTriplet
00280       (svd,0,&sigmamin,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
00281   } else {
00282     *success = PETSC_FALSE; *reason = "slepc no convergence on sigma min";
00283     goto cleanup;
00284   }
00285   *sx = sigmamax; *sm = sigmamin;
00286   *success = PETSC_TRUE; *reason = "slepc normal completion";
00287  cleanup:
00288   ierr = SVDDestroy(&svd);CHKERRQ(ierr);
00289 
00290   PetscFunctionReturn(0);
00291 }
00292 #else
00293 #undef __FUNCT__
00294 #define __FUNCT__ "compute_eigenvalues_petsc"
00295 static PetscErrorCode compute_eigenvalues_petsc
00296 (Mat A,PetscReal **r,PetscReal **c,int *rneig,PetscReal *rx,PetscReal *rm,
00297  PetscBool *success,const char **reason)
00298 {
00299   KSP solver; Mat B; PC pc; Vec x=NULL,b=NULL; MPI_Comm comm;
00300   int maxit=60,local_size,neig; 
00301   const PCType pctype; PetscErrorCode ierr;
00302   PetscFunctionBegin;
00303 
00304   *success = PETSC_TRUE;
00305   if (use_prec) {
00306     ierr = AnaModTraceMessage
00307       ("computing preconditioned spectrum\n"); CHKERRQ(ierr);
00308   } else {
00309     ierr = AnaModTraceMessage
00310       ("computing plain spectrum\n"); CHKERRQ(ierr);
00311   }
00312 
00313   /*
00314    * Create solver around A
00315    */
00316   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00317   ierr = KSPCreate(comm,&solver); CHKERRQ(ierr);
00318   ierr = KSPAppendOptionsPrefix(solver,"ana_"); CHKERRQ(ierr);
00319   ierr = KSPGetPC(solver,&pc); CHKERRQ(ierr);
00320   ierr = PetscObjectQuery
00321     ((PetscObject)A,"approximation",(PetscObject*)&B); CHKERRQ(ierr);
00322   if (!B) B = A;
00323   ierr = PCSetOperators(pc,A,B,DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
00324   
00325 
00326   /*
00327    * No pc, or get the current one from the options database
00328    */
00329   if (!use_prec) {
00330     ierr = PetscOptionsClearValue("-ana_pc_type"); CHKERRQ(ierr);
00331     ierr = PCSetType(pc,PCNONE); CHKERRQ(ierr);
00332     pctype = "none";
00333   }
00334   ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
00335   { /* we can not compute the spectrum for direct methods */
00336     PC pc; const PCType type; PetscBool flg;
00337     ierr = KSPGetPC(solver,&pc); CHKERRQ(ierr);
00338     ierr = PCGetType(pc,&type); CHKERRQ(ierr);
00339     ierr = PetscStrcmp(type,PCLU,&flg); CHKERRQ(ierr);
00340     if (flg) {
00341       *success = PETSC_FALSE; *reason = "direct method"; 
00342       goto exit;}
00343   }
00344 
00345   ierr = MatGetLocalSize(A,&local_size,&ierr); CHKERRQ(ierr);
00346   ierr = VecCreateMPI(comm,local_size,PETSC_DECIDE,&x); CHKERRQ(ierr);
00347   ierr = VecDuplicate(x,&b); CHKERRQ(ierr);
00348 
00349   /*
00350    * GMRES method for eigenvalues
00351    */
00352   ierr = KSPSetType(solver,KSPGMRES); CHKERRQ(ierr);
00353   ierr = KSPGMRESSetRestart(solver,maxit); CHKERRQ(ierr);
00354   ierr = KSPSetTolerances(solver,1.e-10,1.e-30,1.e+6,maxit); CHKERRQ(ierr);
00355   ierr = KSPSetComputeEigenvalues(solver,PETSC_TRUE); CHKERRQ(ierr);
00356   ierr = KSPSetComputeSingularValues(solver,PETSC_TRUE); CHKERRQ(ierr);
00357   ierr = PetscOptionsSetValue
00358     ("-ana_ksp_gmres_modifiedgramschmidt",PETSC_NULL); CHKERRQ(ierr);
00359   {
00360     void *ctx;
00361     ierr = KSPDefaultConvergedCreate(&ctx); CHKERRQ(ierr);
00362     ierr = KSPSetConvergenceTest
00363       (solver,&fivestepplan,ctx,PETSC_NULL); CHKERRQ(ierr);
00364   }
00365 
00366   /*
00367    * run one system solution
00368    */
00369   {
00370     PetscScalar one=1.; KSPConvergedReason convergereason; int its;
00371     ierr = VecSet(x,one); CHKERRQ(ierr);
00372     ierr = PetscPushErrorHandler(&PetscReturnErrorHandler,NULL); CHKERRQ(ierr);
00373     ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
00374     ierr = KSPSetUp(solver);
00375     if (!ierr) {ierr = KSPSolve(solver,x,b);}
00376     ierr = PetscPopErrorHandler(); CHKERRQ(ierr);
00377     if (ierr) {
00378       *success = PETSC_FALSE; *reason  = "KSPSetup/Solve breakdown";
00379       goto exit;}
00380     ierr = KSPGetConvergedReason(solver,&convergereason); CHKERRQ(ierr);
00381     ierr = KSPGetIterationNumber(solver,&its); CHKERRQ(ierr);
00382       printf("%d\n",its);
00383     if (convergereason < 0 && convergereason != KSP_DIVERGED_ITS) {
00384       *success = PETSC_FALSE; *reason = "KSPSolve divergence";
00385       goto exit;
00386     }
00387   }
00388 
00389   /*
00390    * analyse the returns
00391    */
00392   {
00393     PetscReal emax,emin, *real_part,*im_part;
00394     emax = 10; emin = 1;
00395     ierr = KSPComputeExtremeSingularValues(solver,&emax,&emin); CHKERRQ(ierr);
00396     *rx = emax; *rm = emin;
00397 
00398     ierr = PetscMalloc((maxit+1)*sizeof(PetscReal),&real_part); CHKERRQ(ierr);
00399     ierr = PetscMalloc((maxit+1)*sizeof(PetscReal),&im_part); CHKERRQ(ierr);
00400     neig = 0;
00401     ierr = KSPComputeEigenvalues
00402       (solver,maxit,real_part,im_part,&neig); CHKERRQ(ierr);
00403     *rneig = neig;
00404 
00405     if (trace_spectrum) {
00406       int i;
00407       printf("Spectrum calculated:\n");
00408       for (i=0; i<neig; i++)
00409         printf("%d: %e+%e i\n",i,real_part[i],im_part[i]);
00410     }
00411     *r = real_part;
00412     *c = im_part;
00413 
00414 #ifdef PETSC_INTERNALS
00415     if (trace_hessenberg) {
00416       Mat H;
00417       void           *solverdata = solver->data;
00418       KSP_GMRES      *gmres = (KSP_GMRES*)solverdata;
00419       PetscInt       /* n = gmres->it + 1,i,*/ N = gmres->max_k + 2;
00420       // PetscScalar    *hess;
00421       ierr = MatCreateSeqDense
00422         (MPI_COMM_SELF,N,N,gmres->hh_origin,&H); CHKERRQ(ierr);
00423       if (trace_hessenberg) {
00424         ierr = MatView(H,0); CHKERRQ(ierr);
00425       }
00426       ierr = MatDestroy(H); CHKERRQ(ierr);
00427       /*
00428       ierr = PetscMalloc(N*N*sizeof(PetscScalar),&hess); CHKERRQ(ierr);
00429       ierr = PetscMemcpy
00430         (hess,gmres->hh_origin,N*N*sizeof(PetscScalar)); CHKERRQ(ierr);
00431       ierr = GetDataID("spectrum","hessenberg",&id,&has); CHKERRQ(ierr);
00432       HASTOEXIST(has);
00433       ierr = PetscObjectComposedDataSetRealstar
00434         ((PetscObject)A,id,hess); CHKERRQ(ierr);
00435       *h = hess;
00436       */
00437     }
00438 #endif
00439   }
00440 
00441  exit:
00442   if (x) {ierr = VecDestroy(&x); CHKERRQ(ierr);}
00443   if (b) {ierr = VecDestroy(&b); CHKERRQ(ierr);}
00444   ierr = KSPDestroy(&solver); CHKERRQ(ierr);
00445   
00446   PetscFunctionReturn(0);
00447 }
00448 #endif
00449 
00450 #undef __FUNCT__
00451 #define __FUNCT__ "compute_eigenvalues"
00452 /*! Compute a number of eigenvalues (returned as separate real and imaginary
00453   parts) of a matrix.
00454 
00455   Parameters:
00456   - success : PETSC_TRUE or PETSC_FALSE
00457   - reason : explanation of the success parameter. This is a static string;
00458   does not need to be freed.
00459   - r,c : arrays of real and imaginary part of the eigenvalues. If success is
00460   false, these are NULL. Otherwise, the calling environment needs to
00461   free them at some point.
00462   - neig : length of the \c r and \c c arrays. 
00463   This is zero if success is false.
00464 */
00465 static PetscErrorCode compute_eigenvalues
00466 (Mat A,PetscReal **r,PetscReal **c,int *neig,
00467  PetscBool *success,const char **reason)
00468 {
00469 #if !defined(HAVE_SLEPC)
00470   PetscReal emax,emin;
00471 #endif
00472   int id; PetscBool has; PetscErrorCode ierr;
00473   PetscFunctionBegin;
00474 #if defined(HAVE_SLEPC)
00475   ierr = compute_eigenvalues_slepc(A,r,c,neig,success,(char**)reason); CHKERRQ(ierr);
00476 #else
00477   ierr = compute_eigenvalues_petsc
00478     (A,r,c,neig,&emax,&emin,success,reason); CHKERRQ(ierr);
00479 #endif
00480 
00481   if (*success) {
00482     ierr = GetDataID("spectrum","n-ritz-values",&id,&has); CHKERRQ(ierr);
00483     HASTOEXIST(has);
00484     ierr = PetscObjectComposedDataSetInt
00485       ((PetscObject)A,id,*neig); CHKERRQ(ierr);
00486     
00487     ierr = GetDataID("spectrum","ritz-values-r",&id,&has); CHKERRQ(ierr);
00488     HASTOEXIST(has);
00489     ierr = PetscObjectComposedDataSetRealstar
00490       ((PetscObject)A,id,*r); CHKERRQ(ierr);
00491     
00492     ierr = GetDataID("spectrum","ritz-values-c",&id,&has); CHKERRQ(ierr);
00493     HASTOEXIST(has);
00494     ierr = PetscObjectComposedDataSetRealstar
00495       ((PetscObject)A,id,*c); CHKERRQ(ierr);
00496     
00497   }
00498 
00499   PetscFunctionReturn(0);
00500 }
00501 
00502 #undef __FUNCT__
00503 #define __FUNCT__ "compute_singularvalues"
00504 /*! Compute a number of singular values of a matrix.
00505 
00506   Parameters:
00507   - success : PETSC_TRUE or PETSC_FALSE
00508   - reason : explanation of the success parameter. This is a static string;
00509   does not need to be freed.
00510 */
00511 static PetscErrorCode compute_singularvalues
00512 (Mat A,PetscBool *success,const char **reason)
00513 {
00514   PetscReal emax,emin; int id; PetscBool has; PetscErrorCode ierr;
00515   PetscFunctionBegin;
00516 #if defined(HAVE_SLEPC)
00517   ierr = compute_sigmas_slepc(A,&emax,&emin,success,(char**)reason); CHKERRQ(ierr);
00518 #else
00519  {
00520    PetscReal *r,*c; int neig;
00521    ierr = compute_eigenvalues_petsc
00522      (A,&r,&c,&neig,&emax,&emin,success,reason); CHKERRQ(ierr);
00523  }
00524 #endif
00525 
00526   if (*success) {
00527     ierr = GetDataID("spectrum","sigma-max",&id,&has); CHKERRQ(ierr);
00528     HASTOEXIST(has);
00529     ierr = PetscObjectComposedDataSetReal
00530       ((PetscObject)A,id,emax); CHKERRQ(ierr);
00531     
00532     ierr = GetDataID("spectrum","sigma-min",&id,&has); CHKERRQ(ierr);
00533     HASTOEXIST(has);
00534     ierr = PetscObjectComposedDataSetReal
00535       ((PetscObject)A,id,emin); CHKERRQ(ierr);
00536   }
00537 
00538   PetscFunctionReturn(0);
00539 }
00540 
00541 #undef __FUNCT__
00542 #define __FUNCT__ "compute_ellipse_from_Ritz_values"
00543 static PetscErrorCode compute_ellipse_from_Ritz_values
00544 (Mat A,PetscReal *r,PetscReal *c,int neig)
00545 {
00546   PetscReal np,ax,ay,cx,cy;
00547   PetscReal mx,Mx,my,My, lambda_min,lambda_max;
00548   int npos,nneg, i,id; PetscBool has; PetscErrorCode ierr;
00549 
00550   PetscFunctionBegin;
00551   /* #positive vs #negative */
00552   nneg = 0;
00553   for (i=0; i<neig; i++) if (r[i]<0) nneg++;
00554   npos = neig-nneg;
00555   if (nneg==0) np = 1; else np = npos/(1.*neig);
00556   ierr = GetDataID("spectrum","positive-fraction",&id,&has); CHKERRQ(ierr);
00557   HASTOEXIST(has);
00558   ierr = PetscObjectComposedDataSetReal
00559     ((PetscObject)A,id,np); CHKERRQ(ierr);
00560   
00561   /* shape of enclosing ellipse */
00562   Mx = mx = My = my = 0.; lambda_max = lambda_min = 0;
00563   for (i=0; i<neig; i++) {
00564     double lambda = sqrt(r[i]*r[i]+c[i]*c[i]);
00565     if (lambda_min==0 || lambda<lambda_min) lambda_min = lambda;
00566     if (lambda>lambda_max) lambda_max = lambda;
00567     if (r[i]>Mx) {
00568       if (mx==Mx && mx==0) mx=Mx=r[i]; else Mx=r[i];
00569     }
00570     if (r[i]<mx) {
00571       if (mx==Mx && mx==0) Mx=mx=r[i]; else mx=r[i];
00572     }
00573     if (c[i]>My) {
00574       if (my==My && my==0) my=My=c[i]; else My=c[i];
00575     }
00576     if (c[i]<my) {
00577       if (my==My && my==0) My=my=c[i]; else my=c[i];
00578     }
00579   }
00580 
00581   ax = (Mx-mx)/2;
00582   ierr = GetDataID("spectrum","ellipse-ax",&id,&has); CHKERRQ(ierr);
00583   HASTOEXIST(has);
00584   ierr = PetscObjectComposedDataSetReal
00585     ((PetscObject)A,id,ax); CHKERRQ(ierr);
00586 
00587   ay = (My-my)/2;
00588   ierr = GetDataID("spectrum","ellipse-ay",&id,&has); CHKERRQ(ierr);
00589   HASTOEXIST(has);
00590   ierr = PetscObjectComposedDataSetReal
00591     ((PetscObject)A,id,ay); CHKERRQ(ierr);
00592 
00593   cx = (Mx+mx)/2;
00594   ierr = GetDataID("spectrum","ellipse-cx",&id,&has); CHKERRQ(ierr);
00595   HASTOEXIST(has);
00596   /*printf("set cx=%e @ %d from %d\n",cx,id,(int)A);*/
00597   ierr = PetscObjectComposedDataSetReal
00598     ((PetscObject)A,id,cx); CHKERRQ(ierr);
00599 
00600   cy = (My+my)/2;
00601   ierr = GetDataID("spectrum","ellipse-cy",&id,&has); CHKERRQ(ierr);
00602   HASTOEXIST(has);
00603   ierr = PetscObjectComposedDataSetReal
00604     ((PetscObject)A,id,cy); CHKERRQ(ierr);
00605 
00606   PetscFunctionReturn(0);
00607 }
00608 
00609 #undef __FUNCT__
00610 #define __FUNCT__ "NRitzValues"
00611 static PetscErrorCode NRitzValues
00612 (AnaModNumericalProblem prob,AnalysisItem *iv,int *lv,PetscBool *flg)
00613 {
00614   Mat A = (Mat)prob;
00615   int v = 0; PetscBool has; int id; PetscErrorCode ierr;
00616   PetscFunctionBegin;
00617   ierr = GetDataID("spectrum","n-ritz-values",&id,&has); CHKERRQ(ierr);
00618   HASTOEXIST(has);
00619   ierr = PetscObjectComposedDataGetInt
00620     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00621   if (*flg) iv->i = v;
00622   else {
00623     PetscReal *r,*c; int neig; PetscBool success; const char *reason;
00624     ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00625     if (!success) {
00626       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00627       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00628     if (neig==0) {
00629       printf("Spectrum computation failed, reason unknown\n");
00630       *flg = PETSC_FALSE;
00631     } else {
00632       ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00633       ierr = PetscObjectComposedDataGetInt
00634         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00635       if (*flg) iv->i = v;
00636     }
00637   }
00638   
00639   PetscFunctionReturn(0);
00640 }
00641 
00642 #undef __FUNCT__
00643 #define __FUNCT__ "RitzValuesR"
00644 static PetscErrorCode RitzValuesR
00645 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00646 {
00647   Mat A = (Mat)prob;
00648   PetscReal *v = NULL; PetscBool has; int id,id2; PetscErrorCode ierr;
00649   PetscFunctionBegin;
00650 
00651   ierr = GetDataID("spectrum","ritz-values-r",&id,&has); CHKERRQ(ierr);
00652   HASTOEXIST(has);
00653   ierr = PetscObjectComposedDataGetRealstar
00654     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00655 
00656   ierr = GetDataID("spectrum","n-ritz-values",&id2,&has); CHKERRQ(ierr);
00657   HASTOEXIST(has);
00658   ierr = PetscObjectComposedDataGetInt
00659     ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00660   if (*flg) rv->rr = v;
00661   else {
00662     PetscReal *r,*c; int neig; PetscBool success; const char *reason;
00663     ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00664     if (!success) {
00665       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00666       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00667     if (neig==0) {
00668       printf("Spectrum computation failed, reason unknown\n");
00669       *flg = PETSC_FALSE;
00670     } else {
00671       ierr = compute_ellipse_from_Ritz_values(A,r+1,c+1,neig); CHKERRQ(ierr);
00672       ierr = PetscObjectComposedDataGetRealstar
00673         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00674       ierr = PetscObjectComposedDataGetInt
00675         ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00676       if (*flg) rv->rr = v;
00677     }
00678   }
00679   PetscFunctionReturn(0);
00680 }
00681 
00682 #undef __FUNCT__
00683 #define __FUNCT__ "RitzValuesC"
00684 static PetscErrorCode RitzValuesC
00685 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00686 {
00687   Mat A = (Mat)prob;
00688   PetscReal *v = NULL; PetscBool has; int id,id2; PetscErrorCode ierr;
00689   PetscFunctionBegin;
00690   ierr = GetDataID("spectrum","ritz-values-c",&id,&has); CHKERRQ(ierr);
00691   HASTOEXIST(has);
00692   ierr = PetscObjectComposedDataGetRealstar
00693     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00694   ierr = GetDataID("spectrum","n-ritz-values",&id2,&has); CHKERRQ(ierr);
00695   HASTOEXIST(has);
00696   ierr = PetscObjectComposedDataGetInt
00697     ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00698   if (*flg) rv->rr = v;
00699   else {
00700     PetscReal *r,*c; int neig; PetscBool success; const char *reason;
00701     ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00702     if (!success) {
00703       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00704       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00705     if (neig==0) {
00706       printf("Spectrum computation failed, reason unknown\n");
00707       *flg = PETSC_FALSE;
00708     } else {
00709       ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00710       ierr = PetscObjectComposedDataGetRealstar
00711         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00712       ierr = PetscObjectComposedDataGetInt
00713         ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00714       if (*flg) rv->rr = v;
00715     }
00716   }
00717   PetscFunctionReturn(0);
00718 }
00719 
00720 #if 0
00721 #undef __FUNCT__
00722 #define __FUNCT__ "Hessenberg"
00723 static PetscErrorCode Hessenberg
00724 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00725 {
00726   Mat A = (Mat)prob;
00727   PetscReal *v = NULL; PetscBool has; int id,id2; PetscErrorCode ierr;
00728   PetscFunctionBegin;
00729   ierr = GetDataID("spectrum","hessenberg",&id,&has); CHKERRQ(ierr);
00730   HASTOEXIST(has);
00731   ierr = PetscObjectComposedDataGetRealstar
00732     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00733   if (*flg) rv->rr = v;
00734   else {
00735     PetscReal *r,*c; int neig; PetscBool success; const char *reason;
00736     ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00737     if (!success) {
00738       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00739       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00740     if (neig==0) {
00741       printf("Spectrum computation failed, reason unknown\n");
00742       *flg = PETSC_FALSE;
00743     } else {
00744       ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00745       ierr = PetscObjectComposedDataGetRealstar
00746         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00747       if (*flg) rv->rr = v;
00748     }
00749   }
00750   PetscFunctionReturn(0);
00751 }
00752 #endif
00753 
00754 #undef __FUNCT__
00755 #define __FUNCT__ "SpectrumAX"
00756 static PetscErrorCode SpectrumAX
00757 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00758 {
00759   Mat A = (Mat)prob;
00760   PetscReal v = 0; PetscBool has; int id; PetscErrorCode ierr;
00761   PetscFunctionBegin;
00762   ierr = GetDataID("spectrum","ellipse-ax",&id,&has); CHKERRQ(ierr);
00763   HASTOEXIST(has);
00764   ierr = PetscObjectComposedDataGetReal
00765     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00766   if (*flg) rv->r = v;
00767   else {
00768     PetscReal *r,*c; int neig; PetscBool success; const char *reason;
00769     ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00770     if (!success) {
00771       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00772       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00773     if (neig==0) {
00774       printf("Spectrum computation failed, reason unknown\n");
00775       *flg = PETSC_FALSE;
00776     } else {
00777       ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00778       ierr = PetscObjectComposedDataGetReal
00779         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00780       if (*flg) rv->r = v;
00781     }
00782   }
00783   PetscFunctionReturn(0);
00784 }
00785 
00786 #undef __FUNCT__
00787 #define __FUNCT__ "SpectrumAY"
00788 static PetscErrorCode SpectrumAY
00789 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00790 {
00791   Mat A = (Mat)prob;
00792   PetscReal v = 0; PetscBool has; int id; PetscErrorCode ierr;
00793   PetscFunctionBegin;
00794   ierr = GetDataID("spectrum","ellipse-ay",&id,&has); CHKERRQ(ierr);
00795   HASTOEXIST(has);
00796   ierr = PetscObjectComposedDataGetReal
00797     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00798   if (*flg) rv->r = v;
00799   else {
00800     PetscReal *r,*c; int neig; PetscBool success; const char *reason;
00801     ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00802     if (!success) {
00803       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00804       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00805     if (neig==0) {
00806       printf("Spectrum computation failed, reason unknown\n");
00807       *flg = PETSC_FALSE;
00808     } else {
00809       ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00810       ierr = PetscObjectComposedDataGetReal
00811         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00812       if (*flg) rv->r = v;
00813     }
00814   }
00815   PetscFunctionReturn(0);
00816 }
00817 
00818 #undef __FUNCT__
00819 #define __FUNCT__ "SpectrumCX"
00820 static PetscErrorCode SpectrumCX
00821 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00822 {
00823   Mat A = (Mat)prob;
00824   PetscReal v = 0; PetscBool has; int id; PetscErrorCode ierr;
00825   PetscFunctionBegin;
00826   ierr = GetDataID("spectrum","ellipse-cx",&id,&has); CHKERRQ(ierr);
00827   HASTOEXIST(has);
00828   ierr = PetscObjectComposedDataGetReal
00829     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00830   /*printf("finding cx @ %d from %d: %d\n",id,(int)A,*flg);*/
00831   if (*flg) rv->r = v;
00832   else {
00833     PetscReal *r,*c; int neig; PetscBool success; const char *reason;
00834     ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00835     if (!success) {
00836       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00837       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00838     if (neig==0) {
00839       printf("Spectrum computation failed, reason unknown\n");
00840       *flg = PETSC_FALSE;
00841     } else {
00842       ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00843       ierr = PetscObjectComposedDataGetReal
00844         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00845       if (*flg) rv->r = v;
00846     }
00847   }
00848   PetscFunctionReturn(0);
00849 }
00850 
00851 #undef __FUNCT__
00852 #define __FUNCT__ "SpectrumCY"
00853 static PetscErrorCode SpectrumCY
00854 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00855 {
00856   Mat A = (Mat)prob;
00857   PetscReal v = 0; PetscBool has; int id; PetscErrorCode ierr;
00858   PetscFunctionBegin;
00859   ierr = GetDataID("spectrum","ellipse-cy",&id,&has); CHKERRQ(ierr);
00860   HASTOEXIST(has);
00861   ierr = PetscObjectComposedDataGetReal
00862     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00863   if (*flg) rv->r = v;
00864   else {
00865     PetscReal *r,*c; int neig; PetscBool success; const char *reason;
00866     ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00867     if (!success) {
00868       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00869       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00870     if (neig==0) {
00871       printf("Spectrum computation failed, reason unknown\n");
00872       *flg = PETSC_FALSE;
00873     } else {
00874       ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00875       ierr = PetscObjectComposedDataGetReal
00876         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00877       if (*flg) rv->r = v;
00878     }
00879   }
00880   PetscFunctionReturn(0);
00881 }
00882 
00883 #undef __FUNCT__
00884 #define __FUNCT__ "PosFraction"
00885 static PetscErrorCode PosFraction
00886 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00887 {
00888   Mat A = (Mat)prob;
00889   PetscReal v = 0; PetscBool has; int id; PetscErrorCode ierr;
00890   PetscFunctionBegin;
00891   ierr = GetDataID("spectrum","positive-fraction",&id,&has); CHKERRQ(ierr);
00892   HASTOEXIST(has);
00893   ierr = PetscObjectComposedDataGetReal
00894     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00895   if (*flg) rv->r = v;
00896   else {
00897     PetscReal *r,*c; int neig; PetscBool success; const char *reason;
00898     ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00899     if (!success) {
00900       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00901       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00902     if (neig==0) {
00903       printf("Spectrum computation failed, reason unknown\n");
00904       *flg = PETSC_FALSE;
00905     } else {
00906       ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00907       ierr = PetscObjectComposedDataGetReal
00908         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00909       if (*flg) rv->r = v;
00910     }
00911   }
00912   PetscFunctionReturn(0);
00913 }
00914 
00915 #undef __FUNCT__
00916 #define __FUNCT__ "SigmaMax"
00917 static PetscErrorCode SigmaMax
00918 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00919 {
00920   Mat A = (Mat)prob;
00921   PetscReal v = 0; PetscBool has; int id; PetscErrorCode ierr;
00922   PetscFunctionBegin;
00923   ierr = GetDataID("spectrum","sigma-max",&id,&has); CHKERRQ(ierr);
00924   HASTOEXIST(has);
00925   ierr = PetscObjectComposedDataGetReal
00926     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00927   if (*flg) rv->r = v;
00928   else {
00929     PetscBool success; const char *reason;
00930     ierr = compute_singularvalues(A,&success,&reason); CHKERRQ(ierr);
00931     if (!success) {
00932       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00933       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00934     ierr = PetscObjectComposedDataGetReal
00935       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00936     if (*flg) rv->r = v;
00937   }
00938   PetscFunctionReturn(0);
00939 }
00940 
00941 #undef __FUNCT__
00942 #define __FUNCT__ "SigmaMin"
00943 static PetscErrorCode SigmaMin
00944 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00945 {
00946   Mat A = (Mat)prob;
00947   PetscReal v = 0; PetscBool has; int id; PetscErrorCode ierr;
00948   PetscFunctionBegin;
00949   ierr = GetDataID("spectrum","sigma-min",&id,&has); CHKERRQ(ierr);
00950   HASTOEXIST(has);
00951   ierr = PetscObjectComposedDataGetReal
00952     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00953   if (*flg) rv->r = v;
00954   else {
00955     PetscBool success; const char *reason;
00956     ierr = compute_singularvalues(A,&success,&reason); CHKERRQ(ierr);
00957     if (!success) {
00958       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00959       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00960     ierr = PetscObjectComposedDataGetReal
00961       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00962     if (*flg) rv->r = v;
00963   }
00964   PetscFunctionReturn(0);
00965 }
00966 
00967 #undef __FUNCT__
00968 #define __FUNCT__ "Kappa"
00969 static PetscErrorCode Kappa
00970 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00971 {
00972   Mat A = (Mat)prob;
00973   PetscReal v = 0; PetscBool has; int id; PetscErrorCode ierr;
00974   PetscFunctionBegin;
00975   ierr = GetDataID("spectrum","kappa",&id,&has); CHKERRQ(ierr);
00976   HASTOEXIST(has);
00977   ierr = PetscObjectComposedDataGetReal
00978     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00979   if (*flg) rv->r = v;
00980   else {
00981     AnalysisItem ret; PetscReal sigmin,sigmax;
00982     ierr = SigmaMin(prob,&ret,lv,flg); CHKERRQ(ierr);
00983     if (!*flg) PetscFunctionReturn(0);
00984     sigmin = ret.r;
00985     ierr = SigmaMax(prob,&ret,lv,flg); CHKERRQ(ierr);
00986     if (!*flg) PetscFunctionReturn(0);
00987     sigmax = ret.r;
00988     rv->r = sigmax/sigmin;
00989     /*
00990     PetscReal *r,*c; int neig; PetscBool success; const char *reason;
00991     ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00992     if (!success) {
00993       printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00994       *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00995     if (neig==0) {
00996       printf("Spectrum computation failed, reason unknown\n");
00997       *flg = PETSC_FALSE;
00998     } else {
00999       ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
01000       ierr = PetscObjectComposedDataGetReal
01001         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01002       if (*flg) rv->r = v;
01003     }
01004     */
01005   }
01006   PetscFunctionReturn(0);
01007 }
01008 
01009 #undef __FUNCT__
01010 #define __FUNCT__ "MaxEVbyMagRe"
01011 static PetscErrorCode MaxEVbyMagRe
01012 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
01013 {
01014   Mat A = (Mat)prob;
01015   PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0;
01016   PetscFunctionBegin;
01017   ierr = GetDataID
01018     ("spectrum","lambda-max-by-magnitude-re",&id,&has); CHKERRQ(ierr);
01019   HASTOEXIST(has);
01020   ierr = PetscObjectComposedDataGetReal
01021     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01022   if (*flg) rv->r = v;
01023   else {
01024     AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01025 
01026     ierr = ComputeQuantity
01027       (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01028     if (!*flg) PetscFunctionReturn(0);
01029     re = q.rr;
01030     ierr = ComputeQuantity
01031       (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01032     if (!*flg) PetscFunctionReturn(0);
01033     im = q.rr;
01034     *lv = n;
01035 
01036     N = 0; mag = 0;
01037     for (i=0; i<n; i++) {
01038       PetscReal m = sqrt(re[i]*re[i]+im[i]*im[i]);
01039       if (m>mag) {mag = m; N = i;}
01040     }
01041     v = re[N];
01042     ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01043     rv->r = v;
01044   }
01045   PetscFunctionReturn(0);
01046 }
01047 
01048 #undef __FUNCT__
01049 #define __FUNCT__ "MaxEVbyMagIm"
01050 static PetscErrorCode MaxEVbyMagIm
01051 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
01052 {
01053   Mat A = (Mat)prob;
01054   PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0;
01055   PetscFunctionBegin;
01056   ierr = GetDataID
01057     ("spectrum","lambda-max-by-magnitude-im",&id,&has); CHKERRQ(ierr);
01058   HASTOEXIST(has);
01059   ierr = PetscObjectComposedDataGetReal
01060     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01061   if (*flg) rv->r = v;
01062   else {
01063     AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01064 
01065     ierr = ComputeQuantity
01066       (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01067     if (!*flg) PetscFunctionReturn(0);
01068     re = q.rr;
01069     ierr = ComputeQuantity
01070       (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01071     if (!*flg) PetscFunctionReturn(0);
01072     im = q.rr;
01073     *lv = n;
01074 
01075     N = 0; mag = 0;
01076     for (i=0; i<n; i++) {
01077       PetscReal m = sqrt(re[i]*re[i]+im[i]*im[i]);
01078       if (m>mag) {mag = m; N = i;}
01079     }
01080     v = im[N];
01081     ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01082     rv->r = v;
01083   }
01084   PetscFunctionReturn(0);
01085 }
01086 
01087 #undef __FUNCT__
01088 #define __FUNCT__ "MinEVbyMagRe"
01089 static PetscErrorCode MinEVbyMagRe
01090 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
01091 {
01092   Mat A = (Mat)prob;
01093   PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0;
01094   PetscFunctionBegin;
01095   ierr = GetDataID
01096     ("spectrum","lambda-min-by-magnitude-re",&id,&has); CHKERRQ(ierr);
01097   HASTOEXIST(has);
01098   ierr = PetscObjectComposedDataGetReal
01099     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01100   if (*flg) rv->r = v;
01101   else {
01102     AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01103 
01104     ierr = ComputeQuantity
01105       (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01106     if (!*flg) PetscFunctionReturn(0);
01107     re = q.rr;
01108     ierr = ComputeQuantity
01109       (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01110     if (!*flg) PetscFunctionReturn(0);
01111     im = q.rr;
01112     *lv = n;
01113 
01114     N = 0; mag = 0;
01115     for (i=0; i<n; i++) {
01116       PetscReal m = sqrt(re[i]*re[i]+im[i]*im[i]);
01117       if (mag==0 || m<mag) {mag = m; N = i;}
01118     }
01119     v = re[N];
01120     ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01121     rv->r = v;
01122   }
01123   PetscFunctionReturn(0);
01124 }
01125 
01126 #undef __FUNCT__
01127 #define __FUNCT__ "MinEVbyMagIm"
01128 static PetscErrorCode MinEVbyMagIm
01129 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
01130 {
01131   Mat A = (Mat)prob;
01132   PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0;
01133   PetscFunctionBegin;
01134   ierr = GetDataID
01135     ("spectrum","lambda-min-by-magnitude-im",&id,&has); CHKERRQ(ierr);
01136   HASTOEXIST(has);
01137   ierr = PetscObjectComposedDataGetReal
01138     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01139   if (*flg) rv->r = v;
01140   else {
01141     AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01142 
01143     ierr = ComputeQuantity
01144       (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01145     if (!*flg) PetscFunctionReturn(0);
01146     re = q.rr;
01147     ierr = ComputeQuantity
01148       (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01149     if (!*flg) PetscFunctionReturn(0);
01150     im = q.rr;
01151     *lv = n;
01152 
01153     N = 0; mag = 0;
01154     for (i=0; i<n; i++) {
01155       PetscReal m = sqrt(re[i]*re[i]+im[i]*im[i]);
01156       if (mag==0 || m<mag) {mag = m; N = i;}
01157     }
01158     v = im[N];
01159     ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01160     rv->r = v;
01161   }
01162   PetscFunctionReturn(0);
01163 }
01164 
01165 #undef __FUNCT__
01166 #define __FUNCT__ "MaxEVbyRealRe"
01167 static PetscErrorCode MaxEVbyRealRe
01168 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
01169 {
01170   Mat A = (Mat)prob;
01171   PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0;
01172   PetscFunctionBegin;
01173   ierr = GetDataID
01174     ("spectrum","lambda-max-by-real-part-re",&id,&has); CHKERRQ(ierr);
01175   HASTOEXIST(has);
01176   ierr = PetscObjectComposedDataGetReal
01177     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01178   if (*flg) rv->r = v;
01179   else {
01180     AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01181 
01182     ierr = ComputeQuantity
01183       (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01184     if (!*flg) PetscFunctionReturn(0);
01185     re = q.rr;
01186     ierr = ComputeQuantity
01187       (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01188     if (!*flg) PetscFunctionReturn(0);
01189     im = q.rr;
01190     *lv = n;
01191 
01192     N = 0; mag = 0;
01193     for (i=0; i<n; i++) {
01194       if (PetscAbsScalar(re[i])>mag) {mag = PetscAbsScalar(re[i]); N = i;}
01195     }
01196     v = re[N];
01197     ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01198     rv->r = v;
01199   }
01200   PetscFunctionReturn(0);
01201 }
01202 
01203 #undef __FUNCT__
01204 #define __FUNCT__ "MaxEVbyRealIm"
01205 static PetscErrorCode MaxEVbyRealIm
01206 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
01207 {
01208   Mat A = (Mat)prob;
01209   PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0;
01210   PetscFunctionBegin;
01211   ierr = GetDataID
01212     ("spectrum","lambda-max-by-real-part-im",&id,&has); CHKERRQ(ierr);
01213   HASTOEXIST(has);
01214   ierr = PetscObjectComposedDataGetReal
01215     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01216   if (*flg) rv->r = v;
01217   else {
01218     AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01219 
01220     ierr = ComputeQuantity
01221       (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01222     if (!*flg) PetscFunctionReturn(0);
01223     re = q.rr;
01224     ierr = ComputeQuantity
01225       (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01226     if (!*flg) PetscFunctionReturn(0);
01227     im = q.rr;
01228     *lv = n;
01229 
01230     N = 0; mag = 0;
01231     for (i=0; i<n; i++) {
01232       if (PetscAbsScalar(re[i])>mag) {mag = PetscAbsScalar(re[i]); N = i;}
01233     }
01234     v = im[N];
01235     ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01236     rv->r = v;
01237   }
01238   PetscFunctionReturn(0);
01239 }
01240 
01241 #undef __FUNCT__
01242 #define __FUNCT__ "MaxEVbyImRe"
01243 static PetscErrorCode MaxEVbyImRe
01244 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
01245 {
01246   Mat A = (Mat)prob;
01247   PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0;
01248   PetscFunctionBegin;
01249   ierr = GetDataID
01250     ("spectrum","lambda-max-by-im-part-re",&id,&has); CHKERRQ(ierr);
01251   HASTOEXIST(has);
01252   ierr = PetscObjectComposedDataGetReal
01253     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01254   if (*flg) rv->r = v;
01255   else {
01256     AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01257 
01258     ierr = ComputeQuantity
01259       (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01260     if (!*flg) PetscFunctionReturn(0);
01261     re = q.rr;
01262     ierr = ComputeQuantity
01263       (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01264     if (!*flg) PetscFunctionReturn(0);
01265     im = q.rr;
01266     *lv = n;
01267 
01268     N = 0; mag = 0;
01269     for (i=0; i<n; i++) {
01270       if (PetscAbsScalar(im[i])>mag) {mag = PetscAbsScalar(im[i]); N = i;}
01271     }
01272     v = re[N];
01273     ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01274     rv->r = v;
01275   }
01276   PetscFunctionReturn(0);
01277 }
01278 
01279 #undef __FUNCT__
01280 #define __FUNCT__ "MaxEVbyImIm"
01281 static PetscErrorCode MaxEVbyImIm
01282 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
01283 {
01284   Mat A = (Mat)prob;
01285   PetscBool has; int id; PetscErrorCode ierr; PetscReal v = 0;
01286   PetscFunctionBegin;
01287   ierr = GetDataID
01288     ("spectrum","lambda-max-by-im-part-im",&id,&has); CHKERRQ(ierr);
01289   HASTOEXIST(has);
01290   ierr = PetscObjectComposedDataGetReal
01291     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01292   if (*flg) rv->r = v;
01293   else {
01294     AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01295 
01296     ierr = ComputeQuantity
01297       (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01298     if (!*flg) PetscFunctionReturn(0);
01299     re = q.rr;
01300     ierr = ComputeQuantity
01301       (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01302     if (!*flg) PetscFunctionReturn(0);
01303     im = q.rr;
01304     *lv = n;
01305 
01306     N = 0; mag = 0;
01307     for (i=0; i<n; i++) {
01308       if (PetscAbsScalar(im[i])>mag) {mag = PetscAbsScalar(im[i]); N = i;}
01309     }
01310     v = im[N];
01311     ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01312     rv->r = v;
01313   }
01314   PetscFunctionReturn(0);
01315 }
01316 
01317 #undef __FUNCT__
01318 #define __FUNCT__ "SpectrumOptions"
01319 /*!
01320   Spectrum options can be set with \c "-anamod_spectrum <someoption>".
01321   Here we process the options. The possibilities are:
01322   - \c "trace_spectrum" : list the computed spectrum values on the screen
01323   - \c "dump_hessenberg" : generate a Matlab file of the Hessenberg matrix
01324 
01325   See \ref optionsfile.
01326 */
01327 static PetscErrorCode SpectrumOptions(char *opt)
01328 {
01329   PetscBool flg; PetscErrorCode ierr;
01330   PetscFunctionBegin;
01331   ierr = PetscStrcmp(opt,"trace_spectrum",&flg); CHKERRQ(ierr);
01332   if (flg)
01333     trace_spectrum = PETSC_TRUE;
01334   ierr = PetscStrcmp(opt,"trace_hessenberg",&flg); CHKERRQ(ierr);
01335   if (flg)
01336     trace_hessenberg = PETSC_TRUE;
01337   PetscFunctionReturn(0);
01338 }
01339 
01340 
01341 #undef __FUNCT__
01342 #define __FUNCT__ "RegisterSpectrumModules"
01343 PetscErrorCode RegisterSpectrumModules()
01344 {
01345   PetscErrorCode ierr;
01346   PetscFunctionBegin;
01347 #if defined(HAVE_SLEPC)
01348   ierr = SlepcInitialize(0,0,0,0); CHKERRQ(ierr);
01349 #endif
01350 
01351   ierr = RegisterModule
01352     ("spectrum","n-ritz-values",ANALYSISINTEGER,&NRitzValues); CHKERRQ(ierr);
01353   ierr = RegisterModule
01354     ("spectrum","ritz-values-r",ANALYSISDBLARRAY,&RitzValuesR); CHKERRQ(ierr);
01355   ierr = RegisterModule
01356     ("spectrum","ritz-values-c",ANALYSISDBLARRAY,&RitzValuesC); CHKERRQ(ierr);
01357   /*
01358   ierr = RegisterModule
01359     ("spectrum","hessenberg",ANALYSISDBLARRAY,&Hessenberg); CHKERRQ(ierr);
01360   */
01361 
01362   ierr = RegisterModule
01363     ("spectrum","ellipse-ax",ANALYSISDOUBLE,&SpectrumAX); CHKERRQ(ierr);
01364   ierr = RegisterModule
01365     ("spectrum","ellipse-ay",ANALYSISDOUBLE,&SpectrumAY); CHKERRQ(ierr);
01366   ierr = RegisterModule
01367     ("spectrum","ellipse-cx",ANALYSISDOUBLE,&SpectrumCX); CHKERRQ(ierr);
01368   ierr = RegisterModule
01369     ("spectrum","ellipse-cy",ANALYSISDOUBLE,&SpectrumCY); CHKERRQ(ierr);
01370 
01371   ierr = RegisterModule
01372     ("spectrum","kappa",ANALYSISDOUBLE,&Kappa); CHKERRQ(ierr);
01373   ierr = RegisterModule
01374     ("spectrum","positive-fraction",ANALYSISDOUBLE,&PosFraction); CHKERRQ(ierr);
01375 
01376   ierr = RegisterModule
01377     ("spectrum","sigma-max",ANALYSISDOUBLE,&SigmaMax); CHKERRQ(ierr);
01378   ierr = RegisterModule
01379     ("spectrum","sigma-min",ANALYSISDOUBLE,&SigmaMin); CHKERRQ(ierr);
01380 
01381   ierr = RegisterModule
01382     ("spectrum","lambda-max-by-magnitude-re",
01383      ANALYSISDOUBLE,&MaxEVbyMagRe); CHKERRQ(ierr);
01384   ierr = RegisterModule
01385     ("spectrum","lambda-max-by-magnitude-im",
01386      ANALYSISDOUBLE,&MaxEVbyMagIm); CHKERRQ(ierr);
01387   ierr = RegisterModule
01388     ("spectrum","lambda-min-by-magnitude-re",
01389      ANALYSISDOUBLE,&MinEVbyMagRe); CHKERRQ(ierr);
01390   ierr = RegisterModule
01391     ("spectrum","lambda-min-by-magnitude-im",
01392      ANALYSISDOUBLE,&MinEVbyMagIm); CHKERRQ(ierr);
01393   ierr = RegisterModule
01394     ("spectrum","lambda-max-by-real-part-re",
01395      ANALYSISDOUBLE,&MaxEVbyRealRe); CHKERRQ(ierr);
01396   ierr = RegisterModule
01397     ("spectrum","lambda-max-by-real-part-im",
01398      ANALYSISDOUBLE,&MaxEVbyRealIm); CHKERRQ(ierr);
01399   ierr = RegisterModule
01400     ("spectrum","lambda-max-by-im-part-re",
01401      ANALYSISDOUBLE,&MaxEVbyImRe); CHKERRQ(ierr);
01402   ierr = RegisterModule
01403     ("spectrum","lambda-max-by-im-part-im",
01404      ANALYSISDOUBLE,&MaxEVbyImIm); CHKERRQ(ierr);
01405 
01406   ierr = DeclareCategoryOptionFunction
01407     ("spectrum",&SpectrumOptions); CHKERRQ(ierr);
01408 
01409   PetscFunctionReturn(0);
01410 }
01411 
01412 #undef __FUNCT__
01413 #define __FUNCT__ "DeregisterSpectrumModules"
01414 PetscErrorCode DeregisterSpectrumModules()
01415 {
01416 #if defined(HAVE_SLEPC)
01417   PetscErrorCode ierr;
01418   PetscFunctionBegin;
01419   ierr = SlepcFinalize(); CHKERRQ(ierr);
01420   PetscFunctionReturn(0);
01421 #else
01422   PetscFunctionBegin;
01423   PetscFunctionReturn(0);
01424 #endif
01425 }
01426 
01427