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; PetscBool 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(MPI_COMM_WORLD,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 PetscBool 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,PetscBool *flg) 00173 { 00174 Mat A = (Mat)prob; 00175 PetscReal v = 0; PetscBool 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 #undef __FUNCT__ 00222 #define __FUNCT__ "compute_tracea2" 00223 static PetscErrorCode compute_tracea2(Mat A,PetscBool *done) 00224 { 00225 MPI_Comm comm; Mat Awork; PetscScalar trace; int id; 00226 PetscBool has,alloc_work,do_work,any_work; PetscErrorCode ierr; 00227 00228 PetscFunctionBegin; 00229 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00230 ierr = AnaModGetSequentialMatrix 00231 (A,&Awork,&alloc_work,&do_work,&any_work); CHKERRQ(ierr); 00232 if (!any_work) { 00233 *done = PETSC_FALSE; PetscFunctionReturn(0); 00234 } 00235 if (do_work) { 00236 ierr = compute_tracea2_seq(Awork,&trace); CHKERRQ(ierr); 00237 } 00238 if (alloc_work) { 00239 ierr = MatDestroy(&Awork); CHKERRQ(ierr); 00240 } 00241 ierr = MPI_Bcast((void*)&trace,1,MPIU_REAL,0,comm); 00242 ierr = GetDataID("normal","trace-asquared",&id,&has); CHKERRQ(ierr); 00243 HASTOEXIST(has); 00244 ierr = PetscObjectComposedDataSetScalar 00245 ((PetscObject)A,id,trace); CHKERRQ(ierr); 00246 *done = PETSC_TRUE; 00247 PetscFunctionReturn(0); 00248 } 00249 00250 #undef __FUNCT__ 00251 #define __FUNCT__ "TraceA2" 00252 static PetscErrorCode TraceA2 00253 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00254 { 00255 Mat A = (Mat)prob; 00256 PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; 00257 PetscFunctionBegin; 00258 ierr = GetDataID("normal","trace-asquared",&id,&has); CHKERRQ(ierr); 00259 HASTOEXIST(has); 00260 ierr = PetscObjectComposedDataGetScalar 00261 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00262 if (lv) *lv = 1; 00263 if (*flg) rv->r = v; 00264 else { 00265 ierr = compute_tracea2(A,flg); CHKERRQ(ierr); 00266 if (flg) { 00267 ierr = PetscObjectComposedDataGetScalar 00268 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00269 if (*flg) rv->r = v; 00270 } 00271 } 00272 00273 PetscFunctionReturn(0); 00274 } 00275 00276 static int sqrt_times = 50; 00277 /* This sets a bound on the bandwidth of a matrix; matrices wider than 00278 this bound will not have their Commutator norm computed. 00279 00280 The bandwidth is specified as the number of multiples of the square root 00281 of the matrix size. 00282 */ 00283 PetscErrorCode CommutatorNormFAllowSqrtTimes(int n) 00284 { 00285 sqrt_times = n; 00286 return 0; 00287 } 00288 00289 #undef __FUNCT__ 00290 #define __FUNCT__ "MatCommutatorNormF_seq" 00291 /*! 00292 Compute the Frobenius norm of the commutator 00293 \latexonly $AA^t-A^t$ \endlatexonly 00294 of a sequential matrix. 00295 00296 This routine is accessed through MatCommutatorNormF(). 00297 */ 00298 static PetscErrorCode MatCommutatorNormF_seq 00299 (Mat A,PetscReal *result,PetscBool *done) 00300 { 00301 Mat B,E,F; 00302 PetscReal vsum; PetscBool force,flg; 00303 int afirst,alast,bfirst,blast,p,q,N; 00304 PetscErrorCode ierr; 00305 00306 PetscFunctionBegin; 00307 00308 ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr); 00309 00310 ierr = AnaModHasForcedExpensiveComputation(&force); CHKERRQ(ierr); 00311 ierr = MatrixComputeQuantity 00312 (A,"structure","left-bandwidth",(AnalysisItem*)&p,PETSC_NULL,&flg); CHKERRQ(ierr); 00313 if (!flg) p = N; 00314 ierr = MatrixComputeQuantity 00315 (A,"structure","right-bandwidth",(AnalysisItem*)&q,PETSC_NULL,&flg); CHKERRQ(ierr); 00316 if (!flg) q = N; 00317 if (p+q>sqrt_times*sqrt(N) && !force) { 00318 /*printf("too expensive\n");*/ 00319 *done = PETSC_FALSE; 00320 PetscFunctionReturn(0); 00321 } 00322 00323 /* set up matrix and transpose */ 00324 ierr = MatDuplicate(A,MAT_COPY_VALUES,&B); CHKERRQ(ierr); 00325 afirst = 0; alast = N; bfirst = 0; blast = N; 00326 ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&E); CHKERRQ(ierr); 00327 ierr = MatDuplicate(E,MAT_COPY_VALUES,&F); CHKERRQ(ierr); 00328 #define min(a,b) (a<b?a:b) 00329 #define max(a,b) (a>b?a:b) 00330 { 00331 const PetscReal *avals,*bvals,*evals,*fvals; PetscScalar v1,v2; 00332 int a,b, nacols,nbcols, necols,nfcols; 00333 const int *acols,*bcols,*ecols,*fcols; 00334 vsum = 0; 00335 for (a=afirst; a<alast; a++) { 00336 int bmin=max(bfirst,a-(p+q)),bmax=min(blast,a+(p+q)); 00337 ierr = MatGetRow(A,a,&nacols,&acols,&avals); CHKERRQ(ierr); 00338 ierr = MatGetRow(E,a,&necols,&ecols,&evals); CHKERRQ(ierr); 00339 for (b=bmin; b<bmax; b++) { 00340 ierr = MatGetRow(B,b,&nbcols,&bcols,&bvals); CHKERRQ(ierr); 00341 ierr = MatGetRow(F,b,&nfcols,&fcols,&fvals); CHKERRQ(ierr); 00342 ierr = SparseVecProd 00343 (nacols,acols,avals,nbcols,bcols,bvals,&v1); CHKERRQ(ierr); 00344 ierr = SparseVecProd 00345 (necols,ecols,evals,nfcols,fcols,fvals,&v2); CHKERRQ(ierr); 00346 v1 -= v2; 00347 v1 = PetscAbsScalar(v1); 00348 vsum += v1*v1; 00349 ierr = MatRestoreRow(B,b,&nbcols,&bcols,&bvals); CHKERRQ(ierr); 00350 ierr = MatRestoreRow(F,b,&nfcols,&fcols,&fvals); CHKERRQ(ierr); 00351 } 00352 ierr = MatRestoreRow(A,a,&nacols,&acols,&avals); CHKERRQ(ierr); 00353 ierr = MatRestoreRow(E,a,&necols,&ecols,&evals); CHKERRQ(ierr); 00354 } 00355 } 00356 ierr = MatDestroy(&B); CHKERRQ(ierr); 00357 ierr = MatDestroy(&E); CHKERRQ(ierr); 00358 ierr = MatDestroy(&F); CHKERRQ(ierr); 00359 *result = sqrt(vsum); *done = PETSC_TRUE; 00360 00361 PetscFunctionReturn(0); 00362 } 00363 00364 #undef __FUNCT__ 00365 #define __FUNCT__ "MatCommutatorNormF" 00366 /*! 00367 Compute the Frobenius norm of the commutator 00368 \latexonly $AA^t-A^t$ \endlatexonly 00369 00370 This computation can only be done in single-processor mode. 00371 If the commandline 00372 option (see \ref options) \c "-anamod-force sequential" is specified, 00373 the matrix is gathered on processor zero to compute this quantity. 00374 00375 This is an expensive operation \latexonly ($O(N^3)$) \endlatexonly 00376 which can only be optimized for banded matrices. The maximum allowed 00377 bandwidth is set through CommutatorNormFAllowSqrtTimes(). 00378 00379 See MatCommutatorNormF_seq() for the actual computation. 00380 */ 00381 static PetscErrorCode MatCommutatorNormF(Mat A,PetscBool *done) 00382 { 00383 MPI_Comm comm; Mat Awork; PetscReal result; 00384 PetscBool alloc_work,do_work,any_work,has; 00385 int id; PetscErrorCode ierr; 00386 PetscFunctionBegin; 00387 00388 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr); 00389 ierr = AnaModGetSequentialMatrix 00390 (A,&Awork,&alloc_work,&do_work,&any_work); CHKERRQ(ierr); 00391 if (!any_work) { 00392 *done = PETSC_FALSE; PetscFunctionReturn(0); 00393 } 00394 if (do_work) { 00395 ierr = MatCommutatorNormF_seq(Awork,&result,done); CHKERRQ(ierr); 00396 /*printf("result %d:%e\n",*done,result);*/ 00397 } 00398 if (alloc_work) { 00399 ierr = MatDestroy(&Awork); CHKERRQ(ierr); 00400 } 00401 00402 ierr = MPI_Bcast((void*)done,1,MPI_INT,0,comm); 00403 if (*done) { 00404 ierr = MPI_Bcast((void*)&result,1,MPIU_REAL,0,comm); 00405 } else PetscFunctionReturn(0); 00406 00407 ierr = GetDataID("normal","commutator-normF",&id,&has); CHKERRQ(ierr); 00408 HASTOEXIST(has); 00409 ierr = PetscObjectComposedDataSetReal 00410 ((PetscObject)A,id,result); CHKERRQ(ierr); 00411 PetscFunctionReturn(0); 00412 } 00413 00414 #undef __FUNCT__ 00415 #define __FUNCT__ "Commutator" 00416 static PetscErrorCode Commutator 00417 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00418 { 00419 Mat A = (Mat)prob; 00420 PetscReal v = 0; PetscBool has; int id; PetscErrorCode ierr; 00421 PetscFunctionBegin; 00422 ierr = GetDataID("normal","commutator-normF",&id,&has); CHKERRQ(ierr); 00423 HASTOEXIST(has); 00424 ierr = PetscObjectComposedDataGetReal 00425 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00426 if (lv) *lv = 1; 00427 if (*flg) rv->r = v; 00428 else { 00429 ierr = MatCommutatorNormF(A,flg); CHKERRQ(ierr); 00430 if (*flg) { 00431 ierr = PetscObjectComposedDataGetReal 00432 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00433 if (*flg) rv->r = v; 00434 } 00435 } 00436 00437 PetscFunctionReturn(0); 00438 } 00439 00440 #undef __FUNCT__ 00441 #define __FUNCT__ "Lee96bounds" 00442 static PetscErrorCode Lee96bounds(Mat A) 00443 { 00444 Mat AA; PetscReal upper,lower; 00445 PetscScalar tracea2; PetscReal cnorm,fnorm; int id,ql; 00446 PetscBool f,has; PetscErrorCode ierr; 00447 PetscFunctionBegin; 00448 00449 ierr = MatDuplicate(A,MAT_COPY_VALUES,&AA); CHKERRQ(ierr); 00450 ierr = MatCenter(AA); CHKERRQ(ierr); 00451 ierr = TraceA2(AA,(AnalysisItem*)&tracea2,&ql,&f); CHKERRQ(ierr); 00452 if (!f) { 00453 /*printf(">>warning: did not compute trace of A squared\n");*/ 00454 goto exit; 00455 } 00456 ierr = MatNorm(AA,NORM_FROBENIUS,&fnorm); CHKERRQ(ierr); 00457 00458 upper = fnorm*fnorm-PetscAbsScalar(tracea2); 00459 ierr = Commutator(A,(AnalysisItem*)&cnorm,&ql,&f); CHKERRQ(ierr); 00460 if (!f) { 00461 /*printf(">>warning: did not compute commutator Frobenius norm\n");*/ 00462 goto exit; 00463 } 00464 fnorm *= fnorm; 00465 lower = sqrt(fnorm - sqrt(fnorm*fnorm-cnorm*cnorm/2)); 00466 00467 ierr = GetDataID("normal","lee96-ubound",&id,&has); CHKERRQ(ierr); 00468 HASTOEXIST(has); 00469 ierr = PetscObjectComposedDataSetReal 00470 ((PetscObject)A,id,upper); CHKERRQ(ierr); 00471 ierr = GetDataID("normal","lee96-lbound",&id,&has); CHKERRQ(ierr); 00472 HASTOEXIST(has); 00473 ierr = PetscObjectComposedDataSetReal 00474 ((PetscObject)A,id,lower); CHKERRQ(ierr); 00475 00476 exit: 00477 ierr = MatDestroy(&AA); CHKERRQ(ierr); 00478 PetscFunctionReturn(0); 00479 } 00480 00481 #undef __FUNCT__ 00482 #define __FUNCT__ "DepartureLee96U" 00483 static PetscErrorCode DepartureLee96U 00484 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00485 { 00486 Mat A = (Mat)prob; 00487 PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; 00488 PetscFunctionBegin; 00489 ierr = GetDataID("normal","lee96-ubound",&id,&has); CHKERRQ(ierr); 00490 HASTOEXIST(has); 00491 ierr = PetscObjectComposedDataGetReal 00492 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00493 if (lv) *lv = 1; 00494 if (*flg) rv->r = v; 00495 else { 00496 ierr = Lee96bounds(A); CHKERRQ(ierr); 00497 ierr = PetscObjectComposedDataGetReal 00498 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00499 if (*flg) rv->r = v; 00500 } 00501 00502 PetscFunctionReturn(0); 00503 } 00504 00505 #undef __FUNCT__ 00506 #define __FUNCT__ "DepartureLee96L" 00507 static PetscErrorCode DepartureLee96L 00508 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00509 { 00510 Mat A = (Mat)prob; 00511 PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; 00512 PetscFunctionBegin; 00513 ierr = GetDataID("normal","lee96-lbound",&id,&has); CHKERRQ(ierr); 00514 HASTOEXIST(has); 00515 ierr = PetscObjectComposedDataGetReal 00516 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00517 if (lv) *lv = 1; 00518 if (*flg) rv->r = v; 00519 else { 00520 ierr = Lee96bounds(A); CHKERRQ(ierr); 00521 ierr = PetscObjectComposedDataGetReal 00522 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00523 if (*flg) rv->r = v; 00524 } 00525 00526 PetscFunctionReturn(0); 00527 } 00528 00529 #undef __FUNCT__ 00530 #define __FUNCT__ "DepartureRuhe75" 00531 static PetscErrorCode DepartureRuhe75 00532 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00533 { 00534 Mat A = (Mat)prob; 00535 PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr; 00536 PetscFunctionBegin; 00537 ierr = GetDataID("normal","ruhe75-bound",&id,&has); CHKERRQ(ierr); 00538 HASTOEXIST(has); 00539 ierr = PetscObjectComposedDataGetReal 00540 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00541 if (lv) *lv = 1; 00542 if (*flg) rv->r = v; 00543 else { 00544 AnalysisItem q; PetscReal s,er,ec,dep=0; int ql; 00545 00546 ierr = ComputeQuantity 00547 (prob,"spectrum","sigma-min",(AnalysisItem*)&s,&ql,flg); CHKERRQ(ierr); 00548 if (!*flg) PetscFunctionReturn(0); 00549 ierr = ComputeQuantity 00550 (prob,"spectrum","lambda-min-by-magnitude-re",&q,&ql,flg); CHKERRQ(ierr); 00551 if (!*flg) PetscFunctionReturn(0); 00552 er = q.r; 00553 ierr = ComputeQuantity 00554 (prob,"spectrum","lambda-min-by-magnitude-im",&q,&ql,flg); CHKERRQ(ierr); 00555 if (!*flg) PetscFunctionReturn(0); 00556 ec = q.r; 00557 dep = PetscAbsScalar(s-sqrt(er*er+ec*ec)); 00558 00559 ierr = ComputeQuantity 00560 (prob,"spectrum","sigma-max",&q,&ql,flg); CHKERRQ(ierr); 00561 if (!*flg) PetscFunctionReturn(0); 00562 s = q.r; 00563 ierr = ComputeQuantity 00564 (prob,"spectrum","lambda-max-by-magnitude-re",&q,&ql,flg); CHKERRQ(ierr); 00565 if (!*flg) PetscFunctionReturn(0); 00566 er = q.r; 00567 ierr = ComputeQuantity 00568 (prob,"spectrum","lambda-max-by-magnitude-im",&q,&ql,flg); CHKERRQ(ierr); 00569 if (!*flg) PetscFunctionReturn(0); 00570 ec = q.r; 00571 dep = PetscMax(dep,PetscAbsScalar(s-sqrt(er*er+ec*ec))); 00572 00573 ierr = PetscObjectComposedDataSetReal 00574 ((PetscObject)A,id,dep); CHKERRQ(ierr); 00575 rv->r = dep; 00576 } 00577 00578 PetscFunctionReturn(0); 00579 } 00580 00581 #undef __FUNCT__ 00582 #define __FUNCT__ "RegisterNormalityModules" 00583 PetscErrorCode RegisterNormalityModules(void) 00584 { 00585 PetscErrorCode ierr; 00586 PetscFunctionBegin; 00587 00588 ierr = RegisterModule 00589 ("normal","trace-asquared",ANALYSISDOUBLE,&TraceA2); CHKERRQ(ierr); 00590 ierr = RegisterModule 00591 ("normal","commutator-normF",ANALYSISDOUBLE,&Commutator); CHKERRQ(ierr); 00592 00593 ierr = RegisterModule 00594 ("normal","ruhe75-bound",ANALYSISDOUBLE,&DepartureRuhe75); CHKERRQ(ierr); 00595 00596 ierr = RegisterModule 00597 ("normal","lee95-bound",ANALYSISDOUBLE,&DepartureLee95); CHKERRQ(ierr); 00598 00599 ierr = RegisterModule 00600 ("normal","lee96-ubound",ANALYSISDOUBLE,&DepartureLee96U); CHKERRQ(ierr); 00601 ierr = RegisterModule 00602 ("normal","lee96-lbound",ANALYSISDOUBLE,&DepartureLee96L); CHKERRQ(ierr); 00603 00604 PetscFunctionReturn(0); 00605 }