00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
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
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
00174
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
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,PETSC_DECIDE); 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
00195
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
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,PETSC_DECIDE);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 *neig,PetscReal *rx,PetscReal *rm,
00296 PetscTruth *success,char **reason)
00297 {
00298 KSP solver; Mat B; PC pc; Vec x=NULL,b=NULL; MPI_Comm comm;
00299 int maxit=60,local_size;
00300 PCType pctype; PetscErrorCode ierr;
00301 PetscFunctionBegin;
00302
00303 *success = PETSC_TRUE; *neig = 0;
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
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
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 {
00335 PC pc; 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
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 ierr = KSPSetConvergenceTest
00359 (solver,&fivestepplan,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00360
00361
00362
00363
00364 {
00365 PetscScalar one=1.; KSPConvergedReason convergereason;
00366 ierr = VecSet(x,one); CHKERRQ(ierr);
00367 ierr = PetscPushErrorHandler(&PetscReturnErrorHandler,NULL); CHKERRQ(ierr);
00368 ierr = KSPSetUp(solver);
00369 if (!ierr) {ierr = KSPSolve(solver,x,b);}
00370 ierr = PetscPopErrorHandler(); CHKERRQ(ierr);
00371 if (ierr) {
00372 *success = PETSC_FALSE; *reason = "KSPSetup/Solve breakdown";
00373 goto exit;}
00374 ierr = KSPGetConvergedReason(solver,&convergereason); CHKERRQ(ierr);
00375 if (convergereason < 0 && convergereason != KSP_DIVERGED_ITS) {
00376 *success = PETSC_FALSE; *reason = "KSPSolve divergence";
00377 goto exit;
00378 }
00379 }
00380
00381
00382
00383
00384 {
00385 PetscReal emax,emin, *real_part,*im_part;
00386 emax = 10; emin = 1;
00387 ierr = KSPComputeExtremeSingularValues(solver,&emax,&emin); CHKERRQ(ierr);
00388 *rx = emax; *rm = emin;
00389
00390 ierr = PetscMalloc((maxit+1)*sizeof(PetscReal),&real_part); CHKERRQ(ierr);
00391 ierr = PetscMalloc((maxit+1)*sizeof(PetscReal),&im_part); CHKERRQ(ierr);
00392 *neig = 0;
00393 ierr = KSPComputeEigenvalues
00394 (solver,maxit,real_part,im_part,neig); CHKERRQ(ierr);
00395
00396 if (trace_spectrum) {
00397 int i;
00398 printf("Spectrum calculated:\n");
00399 for (i=0; i<*neig; i++)
00400 printf("%d: %e+%e i\n",i,real_part[i],im_part[i]);
00401 }
00402 *r = real_part;
00403 *c = im_part;
00404
00405 #ifdef PETSC_INTERNALS
00406 if (trace_hessenberg) {
00407 Mat H;
00408 void *solverdata = solver->data;
00409 KSP_GMRES *gmres = (KSP_GMRES*)solverdata;
00410 PetscInt N = gmres->max_k + 2;
00411
00412 ierr = MatCreateSeqDense
00413 (MPI_COMM_SELF,N,N,gmres->hh_origin,&H); CHKERRQ(ierr);
00414 if (trace_hessenberg) {
00415 ierr = MatView(H,0); CHKERRQ(ierr);
00416 }
00417 ierr = MatDestroy(H); CHKERRQ(ierr);
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428 }
00429 #endif
00430 }
00431
00432 exit:
00433 if (x) {ierr = VecDestroy(x); CHKERRQ(ierr);}
00434 if (b) {ierr = VecDestroy(b); CHKERRQ(ierr);}
00435 ierr = KSPDestroy(solver); CHKERRQ(ierr);
00436
00437 PetscFunctionReturn(0);
00438 }
00439 #endif
00440
00441 #undef __FUNCT__
00442 #define __FUNCT__ "compute_eigenvalues"
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456 static PetscErrorCode compute_eigenvalues
00457 (Mat A,PetscReal **r,PetscReal **c,int *neig,PetscTruth *success,char **reason)
00458 {
00459 #if !defined(HAVE_SLEPC)
00460 PetscReal emax,emin;
00461 #endif
00462 int id; PetscTruth has; PetscErrorCode ierr;
00463 PetscFunctionBegin;
00464 #if defined(HAVE_SLEPC)
00465 ierr = compute_eigenvalues_slepc(A,r,c,neig,success,reason); CHKERRQ(ierr);
00466 #else
00467 ierr = compute_eigenvalues_petsc
00468 (A,r,c,neig,&emax,&emin,success,reason); CHKERRQ(ierr);
00469 #endif
00470
00471 if (*success) {
00472 ierr = GetDataID("spectrum","n-ritz-values",&id,&has); CHKERRQ(ierr);
00473 HASTOEXIST(has);
00474 ierr = PetscObjectComposedDataSetInt
00475 ((PetscObject)A,id,*neig); CHKERRQ(ierr);
00476
00477 ierr = GetDataID("spectrum","ritz-values-r",&id,&has); CHKERRQ(ierr);
00478 HASTOEXIST(has);
00479 ierr = PetscObjectComposedDataSetRealstar
00480 ((PetscObject)A,id,*r); CHKERRQ(ierr);
00481
00482 ierr = GetDataID("spectrum","ritz-values-c",&id,&has); CHKERRQ(ierr);
00483 HASTOEXIST(has);
00484 ierr = PetscObjectComposedDataSetRealstar
00485 ((PetscObject)A,id,*c); CHKERRQ(ierr);
00486
00487 }
00488
00489 PetscFunctionReturn(0);
00490 }
00491
00492 #undef __FUNCT__
00493 #define __FUNCT__ "compute_singularvalues"
00494
00495
00496
00497
00498
00499
00500
00501 static PetscErrorCode compute_singularvalues
00502 (Mat A,PetscTruth *success,char **reason)
00503 {
00504 PetscReal emax,emin; int id; PetscTruth has; PetscErrorCode ierr;
00505 PetscFunctionBegin;
00506 #if defined(HAVE_SLEPC)
00507 ierr = compute_sigmas_slepc(A,&emax,&emin,success,reason); CHKERRQ(ierr);
00508 #else
00509 {
00510 PetscReal *r,*c; int neig;
00511 ierr = compute_eigenvalues_petsc
00512 (A,&r,&c,&neig,&emax,&emin,success,reason); CHKERRQ(ierr);
00513 }
00514 #endif
00515
00516 if (*success) {
00517 ierr = GetDataID("spectrum","sigma-max",&id,&has); CHKERRQ(ierr);
00518 HASTOEXIST(has);
00519 ierr = PetscObjectComposedDataSetReal
00520 ((PetscObject)A,id,emax); CHKERRQ(ierr);
00521
00522 ierr = GetDataID("spectrum","sigma-min",&id,&has); CHKERRQ(ierr);
00523 HASTOEXIST(has);
00524 ierr = PetscObjectComposedDataSetReal
00525 ((PetscObject)A,id,emin); CHKERRQ(ierr);
00526 }
00527
00528 PetscFunctionReturn(0);
00529 }
00530
00531 #undef __FUNCT__
00532 #define __FUNCT__ "compute_ellipse_from_Ritz_values"
00533 static PetscErrorCode compute_ellipse_from_Ritz_values
00534 (Mat A,PetscReal *r,PetscReal *c,int neig)
00535 {
00536 PetscReal np,ax,ay,cx,cy,kappa;
00537 PetscReal mx,Mx,my,My, lambda_min,lambda_max;
00538 int npos,nneg, i,id; PetscTruth has; PetscErrorCode ierr;
00539
00540 PetscFunctionBegin;
00541
00542 nneg = 0;
00543 for (i=0; i<neig; i++) if (r[i]<0) nneg++;
00544 npos = neig-nneg;
00545 if (nneg==0) np = 1; else np = npos/(1.*neig);
00546 ierr = GetDataID("spectrum","positive-fraction",&id,&has); CHKERRQ(ierr);
00547 HASTOEXIST(has);
00548 ierr = PetscObjectComposedDataSetReal
00549 ((PetscObject)A,id,np); CHKERRQ(ierr);
00550
00551
00552 Mx = mx = My = my = 0.; lambda_max = lambda_min = 0;
00553 for (i=0; i<neig; i++) {
00554 double lambda = sqrt(r[i]*r[i]+c[i]*c[i]);
00555 if (lambda_min==0 || lambda<lambda_min) lambda_min = lambda;
00556 if (lambda>lambda_max) lambda_max = lambda;
00557 if (r[i]>Mx) {
00558 if (mx==Mx && mx==0) mx=Mx=r[i]; else Mx=r[i];
00559 }
00560 if (r[i]<mx) {
00561 if (mx==Mx && mx==0) Mx=mx=r[i]; else mx=r[i];
00562 }
00563 if (c[i]>My) {
00564 if (my==My && my==0) my=My=c[i]; else My=c[i];
00565 }
00566 if (c[i]<my) {
00567 if (my==My && my==0) My=my=c[i]; else my=c[i];
00568 }
00569 }
00570
00571 ax = (Mx-mx)/2;
00572 ierr = GetDataID("spectrum","ellipse-ax",&id,&has); CHKERRQ(ierr);
00573 HASTOEXIST(has);
00574 ierr = PetscObjectComposedDataSetReal
00575 ((PetscObject)A,id,ax); CHKERRQ(ierr);
00576
00577 ay = (My-my)/2;
00578 ierr = GetDataID("spectrum","ellipse-ay",&id,&has); CHKERRQ(ierr);
00579 HASTOEXIST(has);
00580 ierr = PetscObjectComposedDataSetReal
00581 ((PetscObject)A,id,ay); CHKERRQ(ierr);
00582
00583 cx = (Mx+mx)/2;
00584 ierr = GetDataID("spectrum","ellipse-cx",&id,&has); CHKERRQ(ierr);
00585 HASTOEXIST(has);
00586
00587 ierr = PetscObjectComposedDataSetReal
00588 ((PetscObject)A,id,cx); CHKERRQ(ierr);
00589
00590 cy = (My+my)/2;
00591 ierr = GetDataID("spectrum","ellipse-cy",&id,&has); CHKERRQ(ierr);
00592 HASTOEXIST(has);
00593 ierr = PetscObjectComposedDataSetReal
00594 ((PetscObject)A,id,cy); CHKERRQ(ierr);
00595
00596 kappa = lambda_max / lambda_min;
00597 ierr = GetDataID("spectrum","kappa",&id,&has); CHKERRQ(ierr);
00598 HASTOEXIST(has);
00599 ierr = PetscObjectComposedDataSetReal
00600 ((PetscObject)A,id,kappa); CHKERRQ(ierr);
00601
00602 PetscFunctionReturn(0);
00603 }
00604
00605 #undef __FUNCT__
00606 #define __FUNCT__ "NRitzValues"
00607 static PetscErrorCode NRitzValues
00608 (AnaModNumericalProblem prob,AnalysisItem *iv,int *lv,PetscTruth *flg)
00609 {
00610 Mat A = (Mat)prob;
00611 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00612 PetscFunctionBegin;
00613 ierr = GetDataID("spectrum","n-ritz-values",&id,&has); CHKERRQ(ierr);
00614 HASTOEXIST(has);
00615 ierr = PetscObjectComposedDataGetInt
00616 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00617 if (*flg) iv->i = v;
00618 else {
00619 PetscReal *r,*c; int neig; PetscTruth success; char *reason;
00620 ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00621 if (!success) {
00622 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00623 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00624 if (neig==0) {
00625 printf("Spectrum computation failed, reason unknown\n");
00626 *flg = PETSC_FALSE;
00627 } else {
00628 ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00629 ierr = PetscObjectComposedDataGetInt
00630 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00631 if (*flg) iv->i = v;
00632 }
00633 }
00634
00635 PetscFunctionReturn(0);
00636 }
00637
00638 #undef __FUNCT__
00639 #define __FUNCT__ "RitzValuesR"
00640 static PetscErrorCode RitzValuesR
00641 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00642 {
00643 Mat A = (Mat)prob;
00644 PetscReal *v = NULL; PetscTruth has; int id,id2; PetscErrorCode ierr;
00645 PetscFunctionBegin;
00646
00647 ierr = GetDataID("spectrum","ritz-values-r",&id,&has); CHKERRQ(ierr);
00648 HASTOEXIST(has);
00649 ierr = PetscObjectComposedDataGetRealstar
00650 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00651
00652 ierr = GetDataID("spectrum","n-ritz-values",&id2,&has); CHKERRQ(ierr);
00653 HASTOEXIST(has);
00654 ierr = PetscObjectComposedDataGetInt
00655 ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00656 if (*flg) rv->rr = v;
00657 else {
00658 PetscReal *r,*c; int neig; PetscTruth success; char *reason;
00659 ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00660 if (!success) {
00661 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00662 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00663 if (neig==0) {
00664 printf("Spectrum computation failed, reason unknown\n");
00665 *flg = PETSC_FALSE;
00666 } else {
00667 ierr = compute_ellipse_from_Ritz_values(A,r+1,c+1,neig); CHKERRQ(ierr);
00668 ierr = PetscObjectComposedDataGetRealstar
00669 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00670 ierr = PetscObjectComposedDataGetInt
00671 ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00672 if (*flg) rv->rr = v;
00673 }
00674 }
00675 PetscFunctionReturn(0);
00676 }
00677
00678 #undef __FUNCT__
00679 #define __FUNCT__ "RitzValuesC"
00680 static PetscErrorCode RitzValuesC
00681 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00682 {
00683 Mat A = (Mat)prob;
00684 PetscReal *v = NULL; PetscTruth has; int id,id2; PetscErrorCode ierr;
00685 PetscFunctionBegin;
00686 ierr = GetDataID("spectrum","ritz-values-c",&id,&has); CHKERRQ(ierr);
00687 HASTOEXIST(has);
00688 ierr = PetscObjectComposedDataGetRealstar
00689 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00690 ierr = GetDataID("spectrum","n-ritz-values",&id2,&has); CHKERRQ(ierr);
00691 HASTOEXIST(has);
00692 ierr = PetscObjectComposedDataGetInt
00693 ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00694 if (*flg) rv->rr = v;
00695 else {
00696 PetscReal *r,*c; int neig; PetscTruth success; char *reason;
00697 ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00698 if (!success) {
00699 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00700 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00701 if (neig==0) {
00702 printf("Spectrum computation failed, reason unknown\n");
00703 *flg = PETSC_FALSE;
00704 } else {
00705 ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00706 ierr = PetscObjectComposedDataGetRealstar
00707 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00708 ierr = PetscObjectComposedDataGetInt
00709 ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00710 if (*flg) rv->rr = v;
00711 }
00712 }
00713 PetscFunctionReturn(0);
00714 }
00715
00716 #if 0
00717 #undef __FUNCT__
00718 #define __FUNCT__ "Hessenberg"
00719 static PetscErrorCode Hessenberg
00720 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00721 {
00722 Mat A = (Mat)prob;
00723 PetscReal *v = NULL; PetscTruth has; int id,id2; PetscErrorCode ierr;
00724 PetscFunctionBegin;
00725 ierr = GetDataID("spectrum","hessenberg",&id,&has); CHKERRQ(ierr);
00726 HASTOEXIST(has);
00727 ierr = PetscObjectComposedDataGetRealstar
00728 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00729 if (*flg) rv->rr = v;
00730 else {
00731 PetscReal *r,*c; int neig; PetscTruth success; char *reason;
00732 ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00733 if (!success) {
00734 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00735 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00736 if (neig==0) {
00737 printf("Spectrum computation failed, reason unknown\n");
00738 *flg = PETSC_FALSE;
00739 } else {
00740 ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00741 ierr = PetscObjectComposedDataGetRealstar
00742 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00743 if (*flg) rv->rr = v;
00744 }
00745 }
00746 PetscFunctionReturn(0);
00747 }
00748 #endif
00749
00750 #undef __FUNCT__
00751 #define __FUNCT__ "SpectrumAX"
00752 static PetscErrorCode SpectrumAX
00753 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00754 {
00755 Mat A = (Mat)prob;
00756 PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00757 PetscFunctionBegin;
00758 ierr = GetDataID("spectrum","ellipse-ax",&id,&has); CHKERRQ(ierr);
00759 HASTOEXIST(has);
00760 ierr = PetscObjectComposedDataGetReal
00761 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00762 if (*flg) rv->r = v;
00763 else {
00764 PetscReal *r,*c; int neig; PetscTruth success; char *reason;
00765 ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00766 if (!success) {
00767 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00768 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00769 if (neig==0) {
00770 printf("Spectrum computation failed, reason unknown\n");
00771 *flg = PETSC_FALSE;
00772 } else {
00773 ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00774 ierr = PetscObjectComposedDataGetReal
00775 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00776 if (*flg) rv->r = v;
00777 }
00778 }
00779 PetscFunctionReturn(0);
00780 }
00781
00782 #undef __FUNCT__
00783 #define __FUNCT__ "SpectrumAY"
00784 static PetscErrorCode SpectrumAY
00785 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00786 {
00787 Mat A = (Mat)prob;
00788 PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00789 PetscFunctionBegin;
00790 ierr = GetDataID("spectrum","ellipse-ay",&id,&has); CHKERRQ(ierr);
00791 HASTOEXIST(has);
00792 ierr = PetscObjectComposedDataGetReal
00793 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00794 if (*flg) rv->r = v;
00795 else {
00796 PetscReal *r,*c; int neig; PetscTruth success; char *reason;
00797 ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00798 if (!success) {
00799 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00800 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00801 if (neig==0) {
00802 printf("Spectrum computation failed, reason unknown\n");
00803 *flg = PETSC_FALSE;
00804 } else {
00805 ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00806 ierr = PetscObjectComposedDataGetReal
00807 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00808 if (*flg) rv->r = v;
00809 }
00810 }
00811 PetscFunctionReturn(0);
00812 }
00813
00814 #undef __FUNCT__
00815 #define __FUNCT__ "SpectrumCX"
00816 static PetscErrorCode SpectrumCX
00817 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00818 {
00819 Mat A = (Mat)prob;
00820 PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00821 PetscFunctionBegin;
00822 ierr = GetDataID("spectrum","ellipse-cx",&id,&has); CHKERRQ(ierr);
00823 HASTOEXIST(has);
00824 ierr = PetscObjectComposedDataGetReal
00825 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00826
00827 if (*flg) rv->r = v;
00828 else {
00829 PetscReal *r,*c; int neig; PetscTruth success; char *reason;
00830 ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00831 if (!success) {
00832 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00833 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00834 if (neig==0) {
00835 printf("Spectrum computation failed, reason unknown\n");
00836 *flg = PETSC_FALSE;
00837 } else {
00838 ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00839 ierr = PetscObjectComposedDataGetReal
00840 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00841 if (*flg) rv->r = v;
00842 }
00843 }
00844 PetscFunctionReturn(0);
00845 }
00846
00847 #undef __FUNCT__
00848 #define __FUNCT__ "SpectrumCY"
00849 static PetscErrorCode SpectrumCY
00850 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00851 {
00852 Mat A = (Mat)prob;
00853 PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00854 PetscFunctionBegin;
00855 ierr = GetDataID("spectrum","ellipse-cy",&id,&has); CHKERRQ(ierr);
00856 HASTOEXIST(has);
00857 ierr = PetscObjectComposedDataGetReal
00858 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00859 if (*flg) rv->r = v;
00860 else {
00861 PetscReal *r,*c; int neig; PetscTruth success; char *reason;
00862 ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00863 if (!success) {
00864 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00865 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00866 if (neig==0) {
00867 printf("Spectrum computation failed, reason unknown\n");
00868 *flg = PETSC_FALSE;
00869 } else {
00870 ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00871 ierr = PetscObjectComposedDataGetReal
00872 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00873 if (*flg) rv->r = v;
00874 }
00875 }
00876 PetscFunctionReturn(0);
00877 }
00878
00879 #undef __FUNCT__
00880 #define __FUNCT__ "Kappa"
00881 static PetscErrorCode Kappa
00882 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00883 {
00884 Mat A = (Mat)prob;
00885 PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00886 PetscFunctionBegin;
00887 ierr = GetDataID("spectrum","kappa",&id,&has); CHKERRQ(ierr);
00888 HASTOEXIST(has);
00889 ierr = PetscObjectComposedDataGetReal
00890 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00891 if (*flg) rv->r = v;
00892 else {
00893 PetscReal *r,*c; int neig; PetscTruth success; char *reason;
00894 ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00895 if (!success) {
00896 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00897 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00898 if (neig==0) {
00899 printf("Spectrum computation failed, reason unknown\n");
00900 *flg = PETSC_FALSE;
00901 } else {
00902 ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00903 ierr = PetscObjectComposedDataGetReal
00904 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00905 if (*flg) rv->r = v;
00906 }
00907 }
00908 PetscFunctionReturn(0);
00909 }
00910
00911 #undef __FUNCT__
00912 #define __FUNCT__ "PosFraction"
00913 static PetscErrorCode PosFraction
00914 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00915 {
00916 Mat A = (Mat)prob;
00917 PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00918 PetscFunctionBegin;
00919 ierr = GetDataID("spectrum","positive-fraction",&id,&has); CHKERRQ(ierr);
00920 HASTOEXIST(has);
00921 ierr = PetscObjectComposedDataGetReal
00922 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00923 if (*flg) rv->r = v;
00924 else {
00925 PetscReal *r,*c; int neig; PetscTruth success; char *reason;
00926 ierr = compute_eigenvalues(A,&r,&c,&neig,&success,&reason); CHKERRQ(ierr);
00927 if (!success) {
00928 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00929 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00930 if (neig==0) {
00931 printf("Spectrum computation failed, reason unknown\n");
00932 *flg = PETSC_FALSE;
00933 } else {
00934 ierr = compute_ellipse_from_Ritz_values(A,r,c,neig); CHKERRQ(ierr);
00935 ierr = PetscObjectComposedDataGetReal
00936 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00937 if (*flg) rv->r = v;
00938 }
00939 }
00940 PetscFunctionReturn(0);
00941 }
00942
00943 #undef __FUNCT__
00944 #define __FUNCT__ "SigmaMax"
00945 static PetscErrorCode SigmaMax
00946 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00947 {
00948 Mat A = (Mat)prob;
00949 PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00950 PetscFunctionBegin;
00951 ierr = GetDataID("spectrum","sigma-max",&id,&has); CHKERRQ(ierr);
00952 HASTOEXIST(has);
00953 ierr = PetscObjectComposedDataGetReal
00954 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00955 if (*flg) rv->r = v;
00956 else {
00957 PetscTruth success; char *reason;
00958 ierr = compute_singularvalues(A,&success,&reason); CHKERRQ(ierr);
00959 if (!success) {
00960 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00961 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00962 ierr = PetscObjectComposedDataGetReal
00963 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00964 if (*flg) rv->r = v;
00965 }
00966 PetscFunctionReturn(0);
00967 }
00968
00969 #undef __FUNCT__
00970 #define __FUNCT__ "SigmaMin"
00971 static PetscErrorCode SigmaMin
00972 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00973 {
00974 Mat A = (Mat)prob;
00975 PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00976 PetscFunctionBegin;
00977 ierr = GetDataID("spectrum","sigma-min",&id,&has); CHKERRQ(ierr);
00978 HASTOEXIST(has);
00979 ierr = PetscObjectComposedDataGetReal
00980 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00981 if (*flg) rv->r = v;
00982 else {
00983 PetscTruth success; char *reason;
00984 ierr = compute_singularvalues(A,&success,&reason); CHKERRQ(ierr);
00985 if (!success) {
00986 printf("Spectrum computation failure in <%s>: <%s>\n",__FUNCT__,reason);
00987 *flg = PETSC_FALSE; PetscFunctionReturn(0);}
00988 ierr = PetscObjectComposedDataGetReal
00989 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00990 if (*flg) rv->r = v;
00991 }
00992 PetscFunctionReturn(0);
00993 }
00994
00995 #undef __FUNCT__
00996 #define __FUNCT__ "MaxEVbyMagRe"
00997 static PetscErrorCode MaxEVbyMagRe
00998 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00999 {
01000 Mat A = (Mat)prob;
01001 PetscTruth has; int id; PetscErrorCode ierr; PetscReal v = 0;
01002 PetscFunctionBegin;
01003 ierr = GetDataID
01004 ("spectrum","lambda-max-by-magnitude-re",&id,&has); CHKERRQ(ierr);
01005 HASTOEXIST(has);
01006 ierr = PetscObjectComposedDataGetReal
01007 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01008 if (*flg) rv->r = v;
01009 else {
01010 AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01011
01012 ierr = ComputeQuantity
01013 (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01014 if (!*flg) PetscFunctionReturn(0);
01015 re = q.rr;
01016 ierr = ComputeQuantity
01017 (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01018 if (!*flg) PetscFunctionReturn(0);
01019 im = q.rr;
01020 *lv = n;
01021
01022 N = 0; mag = 0;
01023 for (i=0; i<n; i++) {
01024 PetscReal m = sqrt(re[i]*re[i]+im[i]*im[i]);
01025 if (m>mag) {mag = m; N = i;}
01026 }
01027 v = re[N];
01028 ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01029 rv->r = v;
01030 }
01031 PetscFunctionReturn(0);
01032 }
01033
01034 #undef __FUNCT__
01035 #define __FUNCT__ "MaxEVbyMagIm"
01036 static PetscErrorCode MaxEVbyMagIm
01037 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
01038 {
01039 Mat A = (Mat)prob;
01040 PetscTruth has; int id; PetscErrorCode ierr; PetscReal v = 0;
01041 PetscFunctionBegin;
01042 ierr = GetDataID
01043 ("spectrum","lambda-max-by-magnitude-im",&id,&has); CHKERRQ(ierr);
01044 HASTOEXIST(has);
01045 ierr = PetscObjectComposedDataGetReal
01046 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01047 if (*flg) rv->r = v;
01048 else {
01049 AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01050
01051 ierr = ComputeQuantity
01052 (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01053 if (!*flg) PetscFunctionReturn(0);
01054 re = q.rr;
01055 ierr = ComputeQuantity
01056 (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01057 if (!*flg) PetscFunctionReturn(0);
01058 im = q.rr;
01059 *lv = n;
01060
01061 N = 0; mag = 0;
01062 for (i=0; i<n; i++) {
01063 PetscReal m = sqrt(re[i]*re[i]+im[i]*im[i]);
01064 if (m>mag) {mag = m; N = i;}
01065 }
01066 v = im[N];
01067 ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01068 rv->r = v;
01069 }
01070 PetscFunctionReturn(0);
01071 }
01072
01073 #undef __FUNCT__
01074 #define __FUNCT__ "MinEVbyMagRe"
01075 static PetscErrorCode MinEVbyMagRe
01076 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
01077 {
01078 Mat A = (Mat)prob;
01079 PetscTruth has; int id; PetscErrorCode ierr; PetscReal v = 0;
01080 PetscFunctionBegin;
01081 ierr = GetDataID
01082 ("spectrum","lambda-min-by-magnitude-re",&id,&has); CHKERRQ(ierr);
01083 HASTOEXIST(has);
01084 ierr = PetscObjectComposedDataGetReal
01085 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01086 if (*flg) rv->r = v;
01087 else {
01088 AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01089
01090 ierr = ComputeQuantity
01091 (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01092 if (!*flg) PetscFunctionReturn(0);
01093 re = q.rr;
01094 ierr = ComputeQuantity
01095 (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01096 if (!*flg) PetscFunctionReturn(0);
01097 im = q.rr;
01098 *lv = n;
01099
01100 N = 0; mag = 0;
01101 for (i=0; i<n; i++) {
01102 PetscReal m = sqrt(re[i]*re[i]+im[i]*im[i]);
01103 if (mag==0 || m<mag) {mag = m; N = i;}
01104 }
01105 v = re[N];
01106 ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01107 rv->r = v;
01108 }
01109 PetscFunctionReturn(0);
01110 }
01111
01112 #undef __FUNCT__
01113 #define __FUNCT__ "MinEVbyMagIm"
01114 static PetscErrorCode MinEVbyMagIm
01115 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
01116 {
01117 Mat A = (Mat)prob;
01118 PetscTruth has; int id; PetscErrorCode ierr; PetscReal v = 0;
01119 PetscFunctionBegin;
01120 ierr = GetDataID
01121 ("spectrum","lambda-min-by-magnitude-im",&id,&has); CHKERRQ(ierr);
01122 HASTOEXIST(has);
01123 ierr = PetscObjectComposedDataGetReal
01124 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01125 if (*flg) rv->r = v;
01126 else {
01127 AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01128
01129 ierr = ComputeQuantity
01130 (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01131 if (!*flg) PetscFunctionReturn(0);
01132 re = q.rr;
01133 ierr = ComputeQuantity
01134 (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01135 if (!*flg) PetscFunctionReturn(0);
01136 im = q.rr;
01137 *lv = n;
01138
01139 N = 0; mag = 0;
01140 for (i=0; i<n; i++) {
01141 PetscReal m = sqrt(re[i]*re[i]+im[i]*im[i]);
01142 if (mag==0 || m<mag) {mag = m; N = i;}
01143 }
01144 v = im[N];
01145 ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01146 rv->r = v;
01147 }
01148 PetscFunctionReturn(0);
01149 }
01150
01151 #undef __FUNCT__
01152 #define __FUNCT__ "MaxEVbyRealRe"
01153 static PetscErrorCode MaxEVbyRealRe
01154 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
01155 {
01156 Mat A = (Mat)prob;
01157 PetscTruth has; int id; PetscErrorCode ierr; PetscReal v = 0;
01158 PetscFunctionBegin;
01159 ierr = GetDataID
01160 ("spectrum","lambda-max-by-real-part-re",&id,&has); CHKERRQ(ierr);
01161 HASTOEXIST(has);
01162 ierr = PetscObjectComposedDataGetReal
01163 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01164 if (*flg) rv->r = v;
01165 else {
01166 AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01167
01168 ierr = ComputeQuantity
01169 (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01170 if (!*flg) PetscFunctionReturn(0);
01171 re = q.rr;
01172 ierr = ComputeQuantity
01173 (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01174 if (!*flg) PetscFunctionReturn(0);
01175 im = q.rr;
01176 *lv = n;
01177
01178 N = 0; mag = 0;
01179 for (i=0; i<n; i++) {
01180 if (PetscAbsScalar(re[i])>mag) {mag = PetscAbsScalar(re[i]); N = i;}
01181 }
01182 v = re[N];
01183 ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01184 rv->r = v;
01185 }
01186 PetscFunctionReturn(0);
01187 }
01188
01189 #undef __FUNCT__
01190 #define __FUNCT__ "MaxEVbyRealIm"
01191 static PetscErrorCode MaxEVbyRealIm
01192 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
01193 {
01194 Mat A = (Mat)prob;
01195 PetscTruth has; int id; PetscErrorCode ierr; PetscReal v = 0;
01196 PetscFunctionBegin;
01197 ierr = GetDataID
01198 ("spectrum","lambda-max-by-real-part-im",&id,&has); CHKERRQ(ierr);
01199 HASTOEXIST(has);
01200 ierr = PetscObjectComposedDataGetReal
01201 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01202 if (*flg) rv->r = v;
01203 else {
01204 AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01205
01206 ierr = ComputeQuantity
01207 (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01208 if (!*flg) PetscFunctionReturn(0);
01209 re = q.rr;
01210 ierr = ComputeQuantity
01211 (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01212 if (!*flg) PetscFunctionReturn(0);
01213 im = q.rr;
01214 *lv = n;
01215
01216 N = 0; mag = 0;
01217 for (i=0; i<n; i++) {
01218 if (PetscAbsScalar(re[i])>mag) {mag = PetscAbsScalar(re[i]); N = i;}
01219 }
01220 v = im[N];
01221 ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01222 rv->r = v;
01223 }
01224 PetscFunctionReturn(0);
01225 }
01226
01227 #undef __FUNCT__
01228 #define __FUNCT__ "MaxEVbyImRe"
01229 static PetscErrorCode MaxEVbyImRe
01230 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
01231 {
01232 Mat A = (Mat)prob;
01233 PetscTruth has; int id; PetscErrorCode ierr; PetscReal v = 0;
01234 PetscFunctionBegin;
01235 ierr = GetDataID
01236 ("spectrum","lambda-max-by-im-part-re",&id,&has); CHKERRQ(ierr);
01237 HASTOEXIST(has);
01238 ierr = PetscObjectComposedDataGetReal
01239 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01240 if (*flg) rv->r = v;
01241 else {
01242 AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01243
01244 ierr = ComputeQuantity
01245 (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01246 if (!*flg) PetscFunctionReturn(0);
01247 re = q.rr;
01248 ierr = ComputeQuantity
01249 (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01250 if (!*flg) PetscFunctionReturn(0);
01251 im = q.rr;
01252 *lv = n;
01253
01254 N = 0; mag = 0;
01255 for (i=0; i<n; i++) {
01256 if (PetscAbsScalar(im[i])>mag) {mag = PetscAbsScalar(im[i]); N = i;}
01257 }
01258 v = re[N];
01259 ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01260 rv->r = v;
01261 }
01262 PetscFunctionReturn(0);
01263 }
01264
01265 #undef __FUNCT__
01266 #define __FUNCT__ "MaxEVbyImIm"
01267 static PetscErrorCode MaxEVbyImIm
01268 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
01269 {
01270 Mat A = (Mat)prob;
01271 PetscTruth has; int id; PetscErrorCode ierr; PetscReal v = 0;
01272 PetscFunctionBegin;
01273 ierr = GetDataID
01274 ("spectrum","lambda-max-by-im-part-im",&id,&has); CHKERRQ(ierr);
01275 HASTOEXIST(has);
01276 ierr = PetscObjectComposedDataGetReal
01277 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
01278 if (*flg) rv->r = v;
01279 else {
01280 AnalysisItem q; int n,i,N; PetscReal *re,*im,mag;
01281
01282 ierr = ComputeQuantity
01283 (prob,"spectrum","ritz-values-r",&q,&n,flg); CHKERRQ(ierr);
01284 if (!*flg) PetscFunctionReturn(0);
01285 re = q.rr;
01286 ierr = ComputeQuantity
01287 (prob,"spectrum","ritz-values-c",&q,&n,flg); CHKERRQ(ierr);
01288 if (!*flg) PetscFunctionReturn(0);
01289 im = q.rr;
01290 *lv = n;
01291
01292 N = 0; mag = 0;
01293 for (i=0; i<n; i++) {
01294 if (PetscAbsScalar(im[i])>mag) {mag = PetscAbsScalar(im[i]); N = i;}
01295 }
01296 v = im[N];
01297 ierr = PetscObjectComposedDataSetReal((PetscObject)A,id,v); CHKERRQ(ierr);
01298 rv->r = v;
01299 }
01300 PetscFunctionReturn(0);
01301 }
01302
01303 #undef __FUNCT__
01304 #define __FUNCT__ "SpectrumOptions"
01305
01306
01307
01308
01309
01310
01311
01312
01313 static PetscErrorCode SpectrumOptions(char *opt)
01314 {
01315 PetscTruth flg; PetscErrorCode ierr;
01316 PetscFunctionBegin;
01317 ierr = PetscStrcmp(opt,"trace_spectrum",&flg); CHKERRQ(ierr);
01318 if (flg)
01319 trace_spectrum = PETSC_TRUE;
01320 ierr = PetscStrcmp(opt,"trace_hessenberg",&flg); CHKERRQ(ierr);
01321 if (flg)
01322 trace_hessenberg = PETSC_TRUE;
01323 PetscFunctionReturn(0);
01324 }
01325
01326
01327 #undef __FUNCT__
01328 #define __FUNCT__ "RegisterSpectrumModules"
01329 PetscErrorCode RegisterSpectrumModules()
01330 {
01331 PetscErrorCode ierr;
01332 PetscFunctionBegin;
01333
01334 #if defined(HAVE_SLEPC)
01335 ierr = SlepcInitialize(0,0,0,0); CHKERRQ(ierr);
01336 #endif
01337
01338 ierr = RegisterModule
01339 ("spectrum","n-ritz-values",ANALYSISINTEGER,&NRitzValues); CHKERRQ(ierr);
01340 ierr = RegisterModule
01341 ("spectrum","ritz-values-r",ANALYSISDBLARRAY,&RitzValuesR); CHKERRQ(ierr);
01342 ierr = RegisterModule
01343 ("spectrum","ritz-values-c",ANALYSISDBLARRAY,&RitzValuesC); CHKERRQ(ierr);
01344
01345
01346
01347
01348
01349 ierr = RegisterModule
01350 ("spectrum","ellipse-ax",ANALYSISDOUBLE,&SpectrumAX); CHKERRQ(ierr);
01351 ierr = RegisterModule
01352 ("spectrum","ellipse-ay",ANALYSISDOUBLE,&SpectrumAY); CHKERRQ(ierr);
01353 ierr = RegisterModule
01354 ("spectrum","ellipse-cx",ANALYSISDOUBLE,&SpectrumCX); CHKERRQ(ierr);
01355 ierr = RegisterModule
01356 ("spectrum","ellipse-cy",ANALYSISDOUBLE,&SpectrumCY); CHKERRQ(ierr);
01357
01358 ierr = RegisterModule
01359 ("spectrum","kappa",ANALYSISDOUBLE,&Kappa); CHKERRQ(ierr);
01360 ierr = RegisterModule
01361 ("spectrum","positive-fraction",ANALYSISDOUBLE,&PosFraction); CHKERRQ(ierr);
01362
01363 ierr = RegisterModule
01364 ("spectrum","sigma-max",ANALYSISDOUBLE,&SigmaMax); CHKERRQ(ierr);
01365 ierr = RegisterModule
01366 ("spectrum","sigma-min",ANALYSISDOUBLE,&SigmaMin); CHKERRQ(ierr);
01367
01368 ierr = RegisterModule
01369 ("spectrum","lambda-max-by-magnitude-re",
01370 ANALYSISDOUBLE,&MaxEVbyMagRe); CHKERRQ(ierr);
01371 ierr = RegisterModule
01372 ("spectrum","lambda-max-by-magnitude-im",
01373 ANALYSISDOUBLE,&MaxEVbyMagIm); CHKERRQ(ierr);
01374 ierr = RegisterModule
01375 ("spectrum","lambda-min-by-magnitude-re",
01376 ANALYSISDOUBLE,&MinEVbyMagRe); CHKERRQ(ierr);
01377 ierr = RegisterModule
01378 ("spectrum","lambda-min-by-magnitude-im",
01379 ANALYSISDOUBLE,&MinEVbyMagIm); CHKERRQ(ierr);
01380 ierr = RegisterModule
01381 ("spectrum","lambda-max-by-real-part-re",
01382 ANALYSISDOUBLE,&MaxEVbyRealRe); CHKERRQ(ierr);
01383 ierr = RegisterModule
01384 ("spectrum","lambda-max-by-real-part-im",
01385 ANALYSISDOUBLE,&MaxEVbyRealIm); CHKERRQ(ierr);
01386 ierr = RegisterModule
01387 ("spectrum","lambda-max-by-im-part-re",
01388 ANALYSISDOUBLE,&MaxEVbyImRe); CHKERRQ(ierr);
01389 ierr = RegisterModule
01390 ("spectrum","lambda-max-by-im-part-im",
01391 ANALYSISDOUBLE,&MaxEVbyImIm); CHKERRQ(ierr);
01392
01393 ierr = DeclareCategoryOptionFunction
01394 ("spectrum",&SpectrumOptions); CHKERRQ(ierr);
01395
01396 PetscFunctionReturn(0);
01397 }
01398
01399 #undef __FUNCT__
01400 #define __FUNCT__ "DeregisterSpectrumModules"
01401 PetscErrorCode DeregisterSpectrumModules()
01402 {
01403 PetscErrorCode ierr;
01404 PetscFunctionBegin;
01405 #if defined(HAVE_SLEPC)
01406 ierr = SlepcFinalize(); CHKERRQ(ierr);
01407 #endif
01408 PetscFunctionReturn(0);
01409 }
01410
01411