SALSA Analysis Modules
|
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