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; 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 }