SALSA Analysis Modules
normal.c
Go to the documentation of this file.
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 }