SALSA Analysis Modules
|
00001 /*! \file normal.c \ingroup categories 00002 \brief Various estimates of the departure from normality 00003 00004 \section normal Estimates for the departure from normality 00005 00006 The departure from normality is a very hard to calculate quantity. 00007 This file provides several estimates. 00008 00009 \section Usage 00010 00011 Activate this module with: 00012 00013 PetscErrorCode RegisterNormalityModules(); 00014 00015 Compute these elements with 00016 00017 ComputeQuantity("normal",<element>,A,&res,&reslen,&flg); 00018 00019 The auxiliary quantity of the Frobenius norm of the commutator 00020 can be very expensive to calculate. If the bandwidth of the matrix is more 00021 then a specified times (default: 4) the square root of the matrix size, 00022 this routine returns with failure. Set this tolerance with 00023 00024 PetscErrorCode CommutatorNormFAllowSqrtTimes(int n); 00025 00026 If this quantity is unavailable, so are the Lee96 bounds. 00027 00028 Available elements are: 00029 - "trace-asquared" : an auxiliary quantity 00030 - "commutator-normF" : the Frobenius norm of the commutator \f$AA^t-A^tA\f$ 00031 - "ruhe75-bound" : the bound from Ruhe[1975] 00032 - "lee95-bound" : the bound from Lee[1995] 00033 - "lee96-ubound" : the upper bound from Lee[1996] 00034 - "lee96-lbound" : the lower bound from Lee[1996] 00035 00036 \section References 00037 00038 \verbatim 00039 @article{ElPa:nonnormality, 00040 author = {L. Elsner and M.H.C. Paardekoper}, 00041 title = {On Measures of Non-normality for Matrices}, 00042 journal = LAA, 00043 volume = {307/308}, 00044 year = {1979}, 00045 pages = {107--124} 00046 } 00047 00048 @article{Lee:upper-bound, 00049 author = {Steven L. Lee}, 00050 title = {A Practical Upper Bound for Departure from Normality}, 00051 journal = SIMAT, 00052 volume = {16}, 00053 issue = {2}, 00054 year = {1995}, 00055 pages = {462--468}, 00056 abstract = {Upper bound expressed in terms of Hermitian and Anti-Hermitian 00057 part; hence computable in $O(n^2)$ ops, as opposed to bounds by 00058 Henrici~\cite{Henrici:nonnormal-bounds}, which are based on $A^HA$ and 00059 take $O(n^3)$ ops.} 00060 } 00061 00062 @article{Lee:available-bound, 00063 author = {Steven L. Lee}, 00064 title = {Best Available Bounds for Departure from Normality}, 00065 journal = SIMAT, 00066 volume = {17}, 00067 issue = {4}, 00068 year = {1996}, 00069 pages = {984--991} 00070 \endverbatim 00071 } 00072 00073 */ 00074 #include <stdlib.h> 00075 #include "anamod.h" 00076 #include "anamodsalsamodules.h" 00077 #include "anamatrix.h" 00078 #include "petscerror.h" 00079 #include "petscksp.h" 00080 #include "petscmat.h" 00081 00082 #undef __FUNCT__ 00083 #define __FUNCT__ "SparseVecProd" 00084 static PetscErrorCode SparseVecProd 00085 (int na,const int *ia,const PetscScalar *va, 00086 int nb,const int *ib,const PetscScalar *vb, 00087 PetscScalar *result) 00088 { 00089 int pa = 0, pb = 0; PetscScalar v = 0; 00090 PetscFunctionBegin; 00091 00092 while (pa<na && pb<nb) { 00093 if (ia[pa]==ib[pb]) { 00094 v += va[pa]*vb[pb]; pa++; pb++; 00095 } else if (ia[pa]<ib[pb]) pa++; 00096 else pb++; 00097 } 00098 *result = v; 00099 00100 PetscFunctionReturn(0); 00101 } 00102 00103 #undef __FUNCT__ 00104 #define __FUNCT__ "MatCenter" 00105 static PetscErrorCode MatCenter(Mat A) 00106 { 00107 MPI_Comm comm; AnalysisItem q; PetscTruth f; 00108 Vec D; PetscScalar trace,*d; int N,n,i; PetscErrorCode ierr; 00109 00110 PetscFunctionBegin; 00111 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00112 ierr = MatGetLocalSize(A,&n,PETSC_NULL); CHKERRQ(ierr); 00113 ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr); 00114 00115 ierr = MatrixComputeQuantity 00116 (A,"simple","trace",&q,PETSC_NULL,&f); CHKERRQ(ierr); 00117 if (!f) 00118 SETERRQ(1,"Failed to compute trace"); 00119 trace = q.r; 00120 00121 ierr = VecCreateMPI(comm,n,PETSC_DECIDE,&D); CHKERRQ(ierr); 00122 ierr = MatGetDiagonal(A,D); CHKERRQ(ierr); 00123 ierr = VecGetArray(D,&d); CHKERRQ(ierr); 00124 for (i=0; i<n; i++) d[i] -= trace/N; 00125 ierr = VecRestoreArray(D,&d); CHKERRQ(ierr); 00126 ierr = MatDiagonalSet(A,D,INSERT_VALUES); CHKERRQ(ierr); 00127 ierr = VecDestroy(D); CHKERRQ(ierr); 00128 00129 PetscFunctionReturn(0); 00130 } 00131 00132 /* 00133 * this routine needs to be rewritten for the complex case: 00134 * AA needs to be centred with only the Re of the trace, 00135 * A needs to be centred by the Im of the trace 00136 */ 00137 #undef __FUNCT__ 00138 #define __FUNCT__ "Lee95bounds" 00139 static PetscErrorCode Lee95bounds(Mat A) 00140 { 00141 Mat AA; AnalysisItem q; 00142 PetscReal mnorm,nnorm,bound; int id; 00143 PetscTruth f,has; PetscErrorCode ierr; 00144 PetscFunctionBegin; 00145 00146 ierr = MatDuplicate(A,MAT_COPY_VALUES,&AA); CHKERRQ(ierr); 00147 ierr = MatCenter(AA); CHKERRQ(ierr); 00148 00149 ierr = MatrixComputeQuantity 00150 (AA,"simple","symmetry-fsnorm",&q,PETSC_NULL,&f); CHKERRQ(ierr); 00151 if (!f) goto exit; 00152 mnorm = q.r; 00153 ierr = MatrixComputeQuantity 00154 (A,"simple","symmetry-fanorm",&q,PETSC_NULL,&f); CHKERRQ(ierr); 00155 if (!f) goto exit; 00156 nnorm = q.r; 00157 bound = sqrt(2) * PetscMin(mnorm,nnorm); 00158 00159 ierr = GetDataID("normal","lee95-bound",&id,&has); CHKERRQ(ierr); 00160 HASTOEXIST(has); 00161 ierr = PetscObjectComposedDataSetReal 00162 ((PetscObject)A,id,bound); CHKERRQ(ierr); 00163 00164 exit: 00165 ierr = MatDestroy(AA); CHKERRQ(ierr); 00166 PetscFunctionReturn(0); 00167 } 00168 00169 #undef __FUNCT__ 00170 #define __FUNCT__ "DepartureLee95" 00171 static PetscErrorCode DepartureLee95 00172 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00173 { 00174 Mat A = (Mat)prob; 00175 PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr; 00176 PetscFunctionBegin; 00177 ierr = GetDataID("normal","lee95-bound",&id,&has); CHKERRQ(ierr); 00178 HASTOEXIST(has); 00179 ierr = PetscObjectComposedDataGetReal 00180 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00181 if (lv) *lv = 1; 00182 if (*flg) rv->r = v; 00183 else { 00184 ierr = Lee95bounds(A); CHKERRQ(ierr); 00185 ierr = PetscObjectComposedDataGetReal 00186 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00187 if (*flg) rv->r = v; 00188 } 00189 00190 PetscFunctionReturn(0); 00191 } 00192 00193 #undef __FUNCT__ 00194 #define __FUNCT__ "compute_tracea2_seq" 00195 static PetscErrorCode compute_tracea2_seq(Mat A,PetscScalar *trace) 00196 { 00197 MPI_Comm comm; 00198 Mat B; int first,last,row,nacols,nbcols; const int *acols,*bcols; 00199 const PetscScalar *avals,*bvals; PetscScalar vsum,v; 00200 PetscErrorCode ierr; 00201 00202 PetscFunctionBegin; 00203 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00204 ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr); 00205 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&B); CHKERRQ(ierr); 00206 vsum = 0; 00207 for (row=first; row<last; row++) { 00208 ierr = MatGetRow(A,row,&nacols,&acols,&avals); CHKERRQ(ierr); 00209 ierr = MatGetRow(B,row,&nbcols,&bcols,&bvals); CHKERRQ(ierr); 00210 ierr = SparseVecProd 00211 (nacols,acols,avals,nbcols,bcols,bvals,&v); CHKERRQ(ierr); 00212 vsum += v; 00213 ierr = MatRestoreRow(A,row,&nacols,&acols,&avals); CHKERRQ(ierr); 00214 ierr = MatRestoreRow(B,row,&nbcols,&bcols,&bvals); CHKERRQ(ierr); 00215 } 00216 ierr = MatDestroy(B); CHKERRQ(ierr); 00217 *trace = vsum; 00218 PetscFunctionReturn(0); 00219 } 00220 00221 static PetscErrorCode compute_tracea2(Mat A,PetscTruth *done) 00222 { 00223 MPI_Comm comm; Mat Awork; PetscScalar trace; int id; 00224 PetscTruth has,alloc_work,do_work,any_work; PetscErrorCode ierr; 00225 00226 PetscFunctionBegin; 00227 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00228 ierr = AnaModGetSequentialMatrix 00229 (A,&Awork,&alloc_work,&do_work,&any_work); CHKERRQ(ierr); 00230 if (!any_work) { 00231 *done = PETSC_FALSE; PetscFunctionReturn(0); 00232 } 00233 if (do_work) { 00234 ierr = compute_tracea2_seq(Awork,&trace); CHKERRQ(ierr); 00235 } 00236 if (alloc_work) { 00237 ierr = MatDestroy(Awork); CHKERRQ(ierr); 00238 } 00239 MPI_Bcast((void*)&trace,1,MPIU_REAL,0,comm); 00240 ierr = GetDataID("normal","trace-asquared",&id,&has); CHKERRQ(ierr); 00241 HASTOEXIST(has); 00242 ierr = PetscObjectComposedDataSetScalar 00243 ((PetscObject)A,id,trace); CHKERRQ(ierr); 00244 *done = PETSC_TRUE; 00245 PetscFunctionReturn(0); 00246 } 00247 00248 #undef __FUNCT__ 00249 #define __FUNCT__ "TraceA2" 00250 static PetscErrorCode TraceA2 00251 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00252 { 00253 Mat A = (Mat)prob; 00254 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00255 PetscFunctionBegin; 00256 ierr = GetDataID("normal","trace-asquared",&id,&has); CHKERRQ(ierr); 00257 HASTOEXIST(has); 00258 ierr = PetscObjectComposedDataGetScalar 00259 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00260 if (lv) *lv = 1; 00261 if (*flg) rv->r = v; 00262 else { 00263 ierr = compute_tracea2(A,flg); CHKERRQ(ierr); 00264 if (flg) { 00265 ierr = PetscObjectComposedDataGetScalar 00266 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00267 if (*flg) rv->r = v; 00268 } 00269 } 00270 00271 PetscFunctionReturn(0); 00272 } 00273 00274 static int sqrt_times = 50; 00275 /* This sets a bound on the bandwidth of a matrix; matrices wider than 00276 this bound will not have their Commutator norm computed. 00277 00278 The bandwidth is specified as the number of multiples of the square root 00279 of the matrix size. 00280 */ 00281 PetscErrorCode CommutatorNormFAllowSqrtTimes(int n) 00282 { 00283 sqrt_times = n; 00284 return 0; 00285 } 00286 00287 #undef __FUNCT__ 00288 #define __FUNCT__ "MatCommutatorNormF_seq" 00289 /*! 00290 Compute the Frobenius norm of the commutator 00291 \latexonly $AA^t-A^t$ \endlatexonly 00292 of a sequential matrix. 00293 00294 This routine is accessed through MatCommutatorNormF(). 00295 */ 00296 static PetscErrorCode MatCommutatorNormF_seq 00297 (Mat A,PetscReal *result,PetscTruth *done) 00298 { 00299 Mat B,E,F; 00300 PetscReal vsum; PetscTruth force,flg; 00301 int afirst,alast,bfirst,blast,p,q,N; 00302 PetscErrorCode ierr; 00303 00304 PetscFunctionBegin; 00305 00306 ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr); 00307 00308 ierr = AnaModHasForcedExpensiveComputation(&force); CHKERRQ(ierr); 00309 ierr = MatrixComputeQuantity 00310 (A,"structure","left-bandwidth",(AnalysisItem*)&p,PETSC_NULL,&flg); CHKERRQ(ierr); 00311 if (!flg) p = N; 00312 ierr = MatrixComputeQuantity 00313 (A,"structure","right-bandwidth",(AnalysisItem*)&q,PETSC_NULL,&flg); CHKERRQ(ierr); 00314 if (!flg) q = N; 00315 if (p+q>sqrt_times*sqrt(N) && !force) { 00316 /*printf("too expensive\n");*/ 00317 *done = PETSC_FALSE; 00318 PetscFunctionReturn(0); 00319 } 00320 00321 /* set up matrix and transpose */ 00322 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B); CHKERRQ(ierr); 00323 afirst = 0; alast = N; bfirst = 0; blast = N; 00324 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&E); CHKERRQ(ierr); 00325 ierr = MatDuplicate(E,MAT_COPY_VALUES,&F); CHKERRQ(ierr); 00326 #define min(a,b) (a<b?a:b) 00327 #define max(a,b) (a>b?a:b) 00328 { 00329 const PetscReal *avals,*bvals,*evals,*fvals; PetscScalar v1,v2; 00330 int a,b, nacols,nbcols, necols,nfcols; 00331 const int *acols,*bcols,*ecols,*fcols; 00332 vsum = 0; 00333 for (a=afirst; a<alast; a++) { 00334 int bmin=max(bfirst,a-(p+q)),bmax=min(blast,a+(p+q)); 00335 ierr = MatGetRow(A,a,&nacols,&acols,&avals); CHKERRQ(ierr); 00336 ierr = MatGetRow(E,a,&necols,&ecols,&evals); CHKERRQ(ierr); 00337 for (b=bmin; b<bmax; b++) { 00338 ierr = MatGetRow(B,b,&nbcols,&bcols,&bvals); CHKERRQ(ierr); 00339 ierr = MatGetRow(F,b,&nfcols,&fcols,&fvals); CHKERRQ(ierr); 00340 ierr = SparseVecProd 00341 (nacols,acols,avals,nbcols,bcols,bvals,&v1); CHKERRQ(ierr); 00342 ierr = SparseVecProd 00343 (necols,ecols,evals,nfcols,fcols,fvals,&v2); CHKERRQ(ierr); 00344 v1 -= v2; 00345 v1 = PetscAbsScalar(v1); 00346 vsum += v1*v1; 00347 ierr = MatRestoreRow(B,b,&nbcols,&bcols,&bvals); CHKERRQ(ierr); 00348 ierr = MatRestoreRow(F,b,&nfcols,&fcols,&fvals); CHKERRQ(ierr); 00349 } 00350 ierr = MatRestoreRow(A,a,&nacols,&acols,&avals); CHKERRQ(ierr); 00351 ierr = MatRestoreRow(E,a,&necols,&ecols,&evals); CHKERRQ(ierr); 00352 } 00353 } 00354 ierr = MatDestroy(B); CHKERRQ(ierr); 00355 ierr = MatDestroy(E); CHKERRQ(ierr); 00356 ierr = MatDestroy(F); CHKERRQ(ierr); 00357 *result = sqrt(vsum); *done = PETSC_TRUE; 00358 00359 PetscFunctionReturn(0); 00360 } 00361 00362 #undef __FUNCT__ 00363 #define __FUNCT__ "MatCommutatorNormF" 00364 /*! 00365 Compute the Frobenius norm of the commutator 00366 \latexonly $AA^t-A^t$ \endlatexonly 00367 00368 This computation can only be done in single-processor mode. 00369 If the commandline 00370 option (see \ref options) \c "-anamod-force sequential" is specified, 00371 the matrix is gathered on processor zero to compute this quantity. 00372 00373 This is an expensive operation \latexonly ($O(N^3)$) \endlatexonly 00374 which can only be optimized for banded matrices. The maximum allowed 00375 bandwidth is set through CommutatorNormFAllowSqrtTimes(). 00376 00377 See MatCommutatorNormF_seq() for the actual computation. 00378 */ 00379 static PetscErrorCode MatCommutatorNormF(Mat A,PetscTruth *done) 00380 { 00381 MPI_Comm comm; Mat Awork; PetscReal result; 00382 PetscTruth alloc_work,do_work,any_work,has; 00383 int id; PetscErrorCode ierr; 00384 PetscFunctionBegin; 00385 00386 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00387 ierr = AnaModGetSequentialMatrix 00388 (A,&Awork,&alloc_work,&do_work,&any_work); CHKERRQ(ierr); 00389 if (!any_work) { 00390 *done = PETSC_FALSE; PetscFunctionReturn(0); 00391 } 00392 if (do_work) { 00393 ierr = MatCommutatorNormF_seq(Awork,&result,done); CHKERRQ(ierr); 00394 /*printf("result %d:%e\n",*done,result);*/ 00395 } 00396 if (alloc_work) { 00397 ierr = MatDestroy(Awork); CHKERRQ(ierr); 00398 } 00399 00400 MPI_Bcast((void*)done,1,MPI_INT,0,comm); 00401 if (*done) { 00402 MPI_Bcast((void*)&result,1,MPIU_REAL,0,comm); 00403 } else PetscFunctionReturn(0); 00404 00405 ierr = GetDataID("normal","commutator-normF",&id,&has); CHKERRQ(ierr); 00406 HASTOEXIST(has); 00407 ierr = PetscObjectComposedDataSetReal 00408 ((PetscObject)A,id,result); CHKERRQ(ierr); 00409 PetscFunctionReturn(0); 00410 } 00411 00412 #undef __FUNCT__ 00413 #define __FUNCT__ "Commutator" 00414 static PetscErrorCode Commutator 00415 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00416 { 00417 Mat A = (Mat)prob; 00418 PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr; 00419 PetscFunctionBegin; 00420 ierr = GetDataID("normal","commutator-normF",&id,&has); CHKERRQ(ierr); 00421 HASTOEXIST(has); 00422 ierr = PetscObjectComposedDataGetReal 00423 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00424 if (lv) *lv = 1; 00425 if (*flg) rv->r = v; 00426 else { 00427 ierr = MatCommutatorNormF(A,flg); CHKERRQ(ierr); 00428 if (*flg) { 00429 ierr = PetscObjectComposedDataGetReal 00430 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00431 if (*flg) rv->r = v; 00432 } 00433 } 00434 00435 PetscFunctionReturn(0); 00436 } 00437 00438 #undef __FUNCT__ 00439 #define __FUNCT__ "Lee96bounds" 00440 static PetscErrorCode Lee96bounds(Mat A) 00441 { 00442 Mat AA; PetscReal upper,lower; 00443 PetscScalar tracea2; PetscReal cnorm,fnorm; int id,ql; 00444 PetscTruth f,has; PetscErrorCode ierr; 00445 PetscFunctionBegin; 00446 00447 ierr = MatDuplicate(A,MAT_COPY_VALUES,&AA); CHKERRQ(ierr); 00448 ierr = MatCenter(AA); CHKERRQ(ierr); 00449 ierr = TraceA2(AA,(AnalysisItem*)&tracea2,&ql,&f); CHKERRQ(ierr); 00450 if (!f) { 00451 /*printf(">>warning: did not compute trace of A squared\n");*/ 00452 goto exit; 00453 } 00454 ierr = MatNorm(AA,NORM_FROBENIUS,&fnorm); CHKERRQ(ierr); 00455 00456 upper = fnorm*fnorm-PetscAbsScalar(tracea2); 00457 ierr = Commutator(A,(AnalysisItem*)&cnorm,&ql,&f); CHKERRQ(ierr); 00458 if (!f) { 00459 /*printf(">>warning: did not compute commutator Frobenius norm\n");*/ 00460 goto exit; 00461 } 00462 fnorm *= fnorm; 00463 lower = sqrt(fnorm - sqrt(fnorm*fnorm-cnorm*cnorm/2)); 00464 00465 ierr = GetDataID("normal","lee96-ubound",&id,&has); CHKERRQ(ierr); 00466 HASTOEXIST(has); 00467 ierr = PetscObjectComposedDataSetReal 00468 ((PetscObject)A,id,upper); CHKERRQ(ierr); 00469 ierr = GetDataID("normal","lee96-lbound",&id,&has); CHKERRQ(ierr); 00470 HASTOEXIST(has); 00471 ierr = PetscObjectComposedDataSetReal 00472 ((PetscObject)A,id,lower); CHKERRQ(ierr); 00473 00474 exit: 00475 ierr = MatDestroy(AA); CHKERRQ(ierr); 00476 PetscFunctionReturn(0); 00477 } 00478 00479 #undef __FUNCT__ 00480 #define __FUNCT__ "DepartureLee96U" 00481 static PetscErrorCode DepartureLee96U 00482 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00483 { 00484 Mat A = (Mat)prob; 00485 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00486 PetscFunctionBegin; 00487 ierr = GetDataID("normal","lee96-ubound",&id,&has); CHKERRQ(ierr); 00488 HASTOEXIST(has); 00489 ierr = PetscObjectComposedDataGetReal 00490 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00491 if (lv) *lv = 1; 00492 if (*flg) rv->r = v; 00493 else { 00494 ierr = Lee96bounds(A); CHKERRQ(ierr); 00495 ierr = PetscObjectComposedDataGetReal 00496 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00497 if (*flg) rv->r = v; 00498 } 00499 00500 PetscFunctionReturn(0); 00501 } 00502 00503 #undef __FUNCT__ 00504 #define __FUNCT__ "DepartureLee96L" 00505 static PetscErrorCode DepartureLee96L 00506 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00507 { 00508 Mat A = (Mat)prob; 00509 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00510 PetscFunctionBegin; 00511 ierr = GetDataID("normal","lee96-lbound",&id,&has); CHKERRQ(ierr); 00512 HASTOEXIST(has); 00513 ierr = PetscObjectComposedDataGetReal 00514 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00515 if (lv) *lv = 1; 00516 if (*flg) rv->r = v; 00517 else { 00518 ierr = Lee96bounds(A); CHKERRQ(ierr); 00519 ierr = PetscObjectComposedDataGetReal 00520 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00521 if (*flg) rv->r = v; 00522 } 00523 00524 PetscFunctionReturn(0); 00525 } 00526 00527 #undef __FUNCT__ 00528 #define __FUNCT__ "DepartureRuhe75" 00529 static PetscErrorCode DepartureRuhe75 00530 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg) 00531 { 00532 Mat A = (Mat)prob; 00533 PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr; 00534 PetscFunctionBegin; 00535 ierr = GetDataID("normal","ruhe75-bound",&id,&has); CHKERRQ(ierr); 00536 HASTOEXIST(has); 00537 ierr = PetscObjectComposedDataGetReal 00538 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00539 if (lv) *lv = 1; 00540 if (*flg) rv->r = v; 00541 else { 00542 AnalysisItem q; PetscReal s,er,ec,dep=0; int ql; 00543 00544 ierr = ComputeQuantity 00545 (prob,"spectrum","sigma-min",(AnalysisItem*)&s,&ql,flg); CHKERRQ(ierr); 00546 if (!*flg) PetscFunctionReturn(0); 00547 ierr = ComputeQuantity 00548 (prob,"spectrum","lambda-min-by-magnitude-re",&q,&ql,flg); CHKERRQ(ierr); 00549 if (!*flg) PetscFunctionReturn(0); 00550 er = q.r; 00551 ierr = ComputeQuantity 00552 (prob,"spectrum","lambda-min-by-magnitude-im",&q,&ql,flg); CHKERRQ(ierr); 00553 if (!*flg) PetscFunctionReturn(0); 00554 ec = q.r; 00555 dep = PetscAbsScalar(s-sqrt(er*er+ec*ec)); 00556 00557 ierr = ComputeQuantity 00558 (prob,"spectrum","sigma-max",&q,&ql,flg); CHKERRQ(ierr); 00559 if (!*flg) PetscFunctionReturn(0); 00560 s = q.r; 00561 ierr = ComputeQuantity 00562 (prob,"spectrum","lambda-max-by-magnitude-re",&q,&ql,flg); CHKERRQ(ierr); 00563 if (!*flg) PetscFunctionReturn(0); 00564 er = q.r; 00565 ierr = ComputeQuantity 00566 (prob,"spectrum","lambda-max-by-magnitude-im",&q,&ql,flg); CHKERRQ(ierr); 00567 if (!*flg) PetscFunctionReturn(0); 00568 ec = q.r; 00569 dep = PetscMax(dep,PetscAbsScalar(s-sqrt(er*er+ec*ec))); 00570 00571 ierr = PetscObjectComposedDataSetReal 00572 ((PetscObject)A,id,dep); CHKERRQ(ierr); 00573 rv->r = dep; 00574 } 00575 00576 PetscFunctionReturn(0); 00577 } 00578 00579 #undef __FUNCT__ 00580 #define __FUNCT__ "RegisterNormalityModules" 00581 PetscErrorCode RegisterNormalityModules(void) 00582 { 00583 PetscErrorCode ierr; 00584 PetscFunctionBegin; 00585 00586 ierr = RegisterModule 00587 ("normal","trace-asquared",ANALYSISDOUBLE,&TraceA2); CHKERRQ(ierr); 00588 ierr = RegisterModule 00589 ("normal","commutator-normF",ANALYSISDOUBLE,&Commutator); CHKERRQ(ierr); 00590 00591 ierr = RegisterModule 00592 ("normal","ruhe75-bound",ANALYSISDOUBLE,&DepartureRuhe75); CHKERRQ(ierr); 00593 00594 ierr = RegisterModule 00595 ("normal","lee95-bound",ANALYSISDOUBLE,&DepartureLee95); CHKERRQ(ierr); 00596 00597 ierr = RegisterModule 00598 ("normal","lee96-ubound",ANALYSISDOUBLE,&DepartureLee96U); CHKERRQ(ierr); 00599 ierr = RegisterModule 00600 ("normal","lee96-lbound",ANALYSISDOUBLE,&DepartureLee96L); CHKERRQ(ierr); 00601 00602 PetscFunctionReturn(0); 00603 }