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 #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