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