normal.c File Reference

Various estimates of the departure from normality. More...

#include <stdlib.h>
#include "anamod.h"
#include "anamodsalsamodules.h"
#include "anamatrix.h"
#include "petscerror.h"
#include "petscksp.h"
#include "petscmat.h"

Include dependency graph for normal.c:

Go to the source code of this file.

Defines

#define min(a, b)   (a<b?a:b)
#define max(a, b)   (a>b?a:b)

Functions

static PetscErrorCode SparseVecProd (int na, const int *ia, const PetscScalar *va, int nb, const int *ib, const PetscScalar *vb, PetscScalar *result)
static PetscErrorCode MatCenter (Mat A)
static PetscErrorCode Lee95bounds (Mat A)
static PetscErrorCode DepartureLee95 (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscTruth *flg)
static PetscErrorCode compute_tracea2_seq (Mat A, PetscScalar *trace)
static PetscErrorCode compute_tracea2 (Mat A, PetscTruth *done)
static PetscErrorCode TraceA2 (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscTruth *flg)
PetscErrorCode CommutatorNormFAllowSqrtTimes (int n)
static PetscErrorCode MatCommutatorNormF_seq (Mat A, PetscReal *result, PetscTruth *done)
static PetscErrorCode MatCommutatorNormF (Mat A, PetscTruth *done)
static PetscErrorCode Commutator (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscTruth *flg)
static PetscErrorCode Lee96bounds (Mat A)
static PetscErrorCode DepartureLee96U (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscTruth *flg)
static PetscErrorCode DepartureLee96L (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscTruth *flg)
static PetscErrorCode DepartureRuhe75 (AnaModNumericalProblem prob, AnalysisItem *rv, int *lv, PetscTruth *flg)
PetscErrorCode RegisterNormalityModules (void)

Variables

static int sqrt_times = 50


Detailed Description

Various estimates of the departure from normality.

Estimates for the departure from normality

The departure from normality is a very hard to calculate quantity. This file provides several estimates.

Usage

Activate this module with:

PetscErrorCode RegisterNormalityModules();

Compute these elements with

ComputeQuantity("normal",<element>,A,&res,&reslen,&flg);

The auxiliary quantity of the Frobenius norm of the commutator can be very expensive to calculate. If the bandwidth of the matrix is more then a specified times (default: 4) the square root of the matrix size, this routine returns with failure. Set this tolerance with

PetscErrorCode CommutatorNormFAllowSqrtTimes(int n);

If this quantity is unavailable, so are the Lee96 bounds.

Available elements are:

References

@article{ElPa:nonnormality,
author = {L. Elsner and M.H.C. Paardekoper},
title = {On Measures of Non-normality for Matrices},
journal = LAA,
volume = {307/308},
year = {1979},
pages = {107--124}
}

@article{Lee:upper-bound,
author = {Steven L. Lee},
title = {A Practical Upper Bound for Departure from Normality},
journal = SIMAT,
volume = {16},
issue = {2},
year = {1995},
pages = {462--468},
abstract = {Upper bound expressed in terms of Hermitian and Anti-Hermitian
    part; hence computable in $O(n^2)$ ops, as opposed to bounds by
    Henrici~\cite{Henrici:nonnormal-bounds}, which are based on $A^HA$ and
    take $O(n^3)$ ops.}
}

@article{Lee:available-bound,
author = {Steven L. Lee},
title = {Best Available Bounds for Departure from Normality},
journal = SIMAT,
volume = {17},
issue = {4},
year = {1996},
pages = {984--991}
}

Definition in file normal.c.


Define Documentation

#define max ( a,
 )     (a>b?a:b)

Referenced by MatCommutatorNormF_seq().

#define min ( a,
 )     (a<b?a:b)

Referenced by MatCommutatorNormF_seq().


Function Documentation

static PetscErrorCode Commutator ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  lv,
PetscTruth *  flg 
) [static]

Definition at line 415 of file normal.c.

References GetDataID(), HASTOEXIST, id, MatCommutatorNormF(), and AnalysisItem::r.

Referenced by Lee96bounds(), and RegisterNormalityModules().

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 }

Here is the call graph for this function:

PetscErrorCode CommutatorNormFAllowSqrtTimes ( int  n  ) 

Definition at line 281 of file normal.c.

Referenced by AnaModRegisterSalsaModules().

00282 {
00283   sqrt_times = n;
00284   return 0;
00285 }

static PetscErrorCode compute_tracea2 ( Mat  A,
PetscTruth *  done 
) [static]

Definition at line 221 of file normal.c.

References AnaModGetSequentialMatrix(), compute_tracea2_seq(), GetDataID(), HASTOEXIST, and id.

Referenced by TraceA2().

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 }

Here is the call graph for this function:

static PetscErrorCode compute_tracea2_seq ( Mat  A,
PetscScalar *  trace 
) [static]

Definition at line 195 of file normal.c.

References SparseVecProd().

Referenced by compute_tracea2().

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 }

Here is the call graph for this function:

static PetscErrorCode DepartureLee95 ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  lv,
PetscTruth *  flg 
) [static]

Definition at line 172 of file normal.c.

References GetDataID(), HASTOEXIST, id, Lee95bounds(), and AnalysisItem::r.

Referenced by RegisterNormalityModules().

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 }

Here is the call graph for this function:

static PetscErrorCode DepartureLee96L ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  lv,
PetscTruth *  flg 
) [static]

Definition at line 506 of file normal.c.

References GetDataID(), HASTOEXIST, id, Lee96bounds(), and AnalysisItem::r.

Referenced by RegisterNormalityModules().

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 }

Here is the call graph for this function:

static PetscErrorCode DepartureLee96U ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  lv,
PetscTruth *  flg 
) [static]

Definition at line 482 of file normal.c.

References GetDataID(), HASTOEXIST, id, Lee96bounds(), and AnalysisItem::r.

Referenced by RegisterNormalityModules().

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 }

Here is the call graph for this function:

static PetscErrorCode DepartureRuhe75 ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  lv,
PetscTruth *  flg 
) [static]

Definition at line 530 of file normal.c.

References ComputeQuantity(), GetDataID(), HASTOEXIST, id, and AnalysisItem::r.

Referenced by RegisterNormalityModules().

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 }

Here is the call graph for this function:

static PetscErrorCode Lee95bounds ( Mat  A  )  [static]

Definition at line 139 of file normal.c.

References GetDataID(), HASTOEXIST, id, MatCenter(), MatrixComputeQuantity(), and AnalysisItem::r.

Referenced by DepartureLee95().

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 }

Here is the call graph for this function:

static PetscErrorCode Lee96bounds ( Mat  A  )  [static]

Definition at line 440 of file normal.c.

References Commutator(), GetDataID(), HASTOEXIST, id, MatCenter(), and TraceA2().

Referenced by DepartureLee96L(), and DepartureLee96U().

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 }

Here is the call graph for this function:

static PetscErrorCode MatCenter ( Mat  A  )  [static]

Definition at line 105 of file normal.c.

References MatrixComputeQuantity(), and AnalysisItem::r.

Referenced by Lee95bounds(), and Lee96bounds().

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 }

Here is the call graph for this function:

static PetscErrorCode MatCommutatorNormF ( Mat  A,
PetscTruth *  done 
) [static]

Compute the Frobenius norm of the commutator

This computation can only be done in single-processor mode. If the commandline option (see Commandline options) "-anamod-force sequential" is specified, the matrix is gathered on processor zero to compute this quantity.

This is an expensive operation which can only be optimized for banded matrices. The maximum allowed bandwidth is set through CommutatorNormFAllowSqrtTimes().

See MatCommutatorNormF_seq() for the actual computation.

Definition at line 379 of file normal.c.

References AnaModGetSequentialMatrix(), GetDataID(), HASTOEXIST, id, and MatCommutatorNormF_seq().

Referenced by Commutator().

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 }

Here is the call graph for this function:

static PetscErrorCode MatCommutatorNormF_seq ( Mat  A,
PetscReal *  result,
PetscTruth *  done 
) [static]

Compute the Frobenius norm of the commutator of a sequential matrix.

This routine is accessed through MatCommutatorNormF().

Definition at line 297 of file normal.c.

References AnaModHasForcedExpensiveComputation(), MatrixComputeQuantity(), max, min, and SparseVecProd().

Referenced by MatCommutatorNormF().

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 }

Here is the call graph for this function:

PetscErrorCode RegisterNormalityModules ( void   ) 

Definition at line 581 of file normal.c.

References ANALYSISDOUBLE, Commutator(), DepartureLee95(), DepartureLee96L(), DepartureLee96U(), DepartureRuhe75(), RegisterModule(), and TraceA2().

Referenced by AnaModRegisterSalsaModules(), and AnaModRegisterStandardModules().

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 }

Here is the call graph for this function:

static PetscErrorCode SparseVecProd ( int  na,
const int *  ia,
const PetscScalar *  va,
int  nb,
const int *  ib,
const PetscScalar *  vb,
PetscScalar *  result 
) [static]

Definition at line 85 of file normal.c.

Referenced by compute_tracea2_seq(), and MatCommutatorNormF_seq().

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 }

static PetscErrorCode TraceA2 ( AnaModNumericalProblem  prob,
AnalysisItem rv,
int *  lv,
PetscTruth *  flg 
) [static]

Definition at line 251 of file normal.c.

References compute_tracea2(), GetDataID(), HASTOEXIST, id, and AnalysisItem::r.

Referenced by Lee96bounds(), and RegisterNormalityModules().

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 }

Here is the call graph for this function:


Variable Documentation

int sqrt_times = 50 [static]

Definition at line 274 of file normal.c.


Generated on Sun Oct 4 04:01:18 2009 for SALSA Analysis Modules by  doxygen 1.5.9