SALSA Analysis Modules
Defines | Functions | Variables
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().

{
  Mat A = (Mat)prob;
  PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = GetDataID("normal","commutator-normF",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetReal
    ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  if (lv) *lv = 1;
  if (*flg) rv->r = v;
  else {
    ierr = MatCommutatorNormF(A,flg); CHKERRQ(ierr);
    if (*flg) {
      ierr = PetscObjectComposedDataGetReal
        ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
      if (*flg) rv->r = v;
    }
  }
  
  PetscFunctionReturn(0);
}

Here is the call graph for this function:

PetscErrorCode CommutatorNormFAllowSqrtTimes ( int  n)

Definition at line 281 of file normal.c.

Referenced by AnaModRegisterSalsaModules().

{
  sqrt_times = n;
  return 0;
}
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().

{
  MPI_Comm comm; Mat Awork; PetscScalar trace; int id;
  PetscTruth has,alloc_work,do_work,any_work; PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
  ierr = AnaModGetSequentialMatrix
    (A,&Awork,&alloc_work,&do_work,&any_work); CHKERRQ(ierr);
  if (!any_work) {
    *done = PETSC_FALSE; PetscFunctionReturn(0);
  }
  if (do_work) {
    ierr = compute_tracea2_seq(Awork,&trace); CHKERRQ(ierr);
  }
  if (alloc_work) {
    ierr = MatDestroy(Awork); CHKERRQ(ierr);
  }
  MPI_Bcast((void*)&trace,1,MPIU_REAL,0,comm);
  ierr = GetDataID("normal","trace-asquared",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataSetScalar
    ((PetscObject)A,id,trace); CHKERRQ(ierr);
  *done = PETSC_TRUE;
  PetscFunctionReturn(0);
}

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().

{
  MPI_Comm comm;
  Mat B; int first,last,row,nacols,nbcols; const int *acols,*bcols;
  const PetscScalar *avals,*bvals; PetscScalar vsum,v;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
  ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
  ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&B); CHKERRQ(ierr);
  vsum = 0;
  for (row=first; row<last; row++) {
    ierr = MatGetRow(A,row,&nacols,&acols,&avals); CHKERRQ(ierr);
    ierr = MatGetRow(B,row,&nbcols,&bcols,&bvals); CHKERRQ(ierr);
    ierr = SparseVecProd
      (nacols,acols,avals,nbcols,bcols,bvals,&v); CHKERRQ(ierr);
    vsum += v;
    ierr = MatRestoreRow(A,row,&nacols,&acols,&avals); CHKERRQ(ierr);
    ierr = MatRestoreRow(B,row,&nbcols,&bcols,&bvals); CHKERRQ(ierr);
  }
  ierr = MatDestroy(B); CHKERRQ(ierr);
  *trace = vsum;
  PetscFunctionReturn(0);
}

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().

{
  Mat A = (Mat)prob;
  PetscReal v = 0; PetscTruth has; int id; PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = GetDataID("normal","lee95-bound",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetReal
    ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  if (lv) *lv = 1;
  if (*flg) rv->r = v;
  else {
    ierr = Lee95bounds(A); CHKERRQ(ierr);
    ierr = PetscObjectComposedDataGetReal
      ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
    if (*flg) rv->r = v;
  }
  
  PetscFunctionReturn(0);
}

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().

{
  Mat A = (Mat)prob;
  PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = GetDataID("normal","lee96-lbound",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetReal
    ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  if (lv) *lv = 1;
  if (*flg) rv->r = v;
  else {
    ierr = Lee96bounds(A); CHKERRQ(ierr);
    ierr = PetscObjectComposedDataGetReal
      ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
    if (*flg) rv->r = v;
  }
  
  PetscFunctionReturn(0);
}

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().

{
  Mat A = (Mat)prob;
  PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = GetDataID("normal","lee96-ubound",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetReal
    ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  if (lv) *lv = 1;
  if (*flg) rv->r = v;
  else {
    ierr = Lee96bounds(A); CHKERRQ(ierr);
    ierr = PetscObjectComposedDataGetReal
      ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
    if (*flg) rv->r = v;
  }
  
  PetscFunctionReturn(0);
}

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().

{
  Mat A = (Mat)prob;
  PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = GetDataID("normal","ruhe75-bound",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetReal
    ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  if (lv) *lv = 1;
  if (*flg) rv->r = v;
  else {
    AnalysisItem q; PetscReal s,er,ec,dep=0; int ql;

    ierr = ComputeQuantity
      (prob,"spectrum","sigma-min",(AnalysisItem*)&s,&ql,flg); CHKERRQ(ierr);
    if (!*flg) PetscFunctionReturn(0);
    ierr = ComputeQuantity
      (prob,"spectrum","lambda-min-by-magnitude-re",&q,&ql,flg); CHKERRQ(ierr);
    if (!*flg) PetscFunctionReturn(0);
    er = q.r;
    ierr = ComputeQuantity
      (prob,"spectrum","lambda-min-by-magnitude-im",&q,&ql,flg); CHKERRQ(ierr);
    if (!*flg) PetscFunctionReturn(0);
    ec = q.r;
    dep = PetscAbsScalar(s-sqrt(er*er+ec*ec));

    ierr = ComputeQuantity
      (prob,"spectrum","sigma-max",&q,&ql,flg); CHKERRQ(ierr);
    if (!*flg) PetscFunctionReturn(0);
    s = q.r;
    ierr = ComputeQuantity
      (prob,"spectrum","lambda-max-by-magnitude-re",&q,&ql,flg); CHKERRQ(ierr);
    if (!*flg) PetscFunctionReturn(0);
    er = q.r;
    ierr = ComputeQuantity
      (prob,"spectrum","lambda-max-by-magnitude-im",&q,&ql,flg); CHKERRQ(ierr);
    if (!*flg) PetscFunctionReturn(0);
    ec = q.r;
    dep = PetscMax(dep,PetscAbsScalar(s-sqrt(er*er+ec*ec)));

    ierr = PetscObjectComposedDataSetReal
      ((PetscObject)A,id,dep); CHKERRQ(ierr);
    rv->r = dep;
  }
  
  PetscFunctionReturn(0);
}

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().

{
  Mat AA; AnalysisItem q;
  PetscReal mnorm,nnorm,bound; int id; 
  PetscTruth f,has; PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = MatDuplicate(A,MAT_COPY_VALUES,&AA); CHKERRQ(ierr);
  ierr = MatCenter(AA); CHKERRQ(ierr);

  ierr = MatrixComputeQuantity
    (AA,"simple","symmetry-fsnorm",&q,PETSC_NULL,&f); CHKERRQ(ierr);
  if (!f) goto exit;
  mnorm = q.r;
  ierr = MatrixComputeQuantity
    (A,"simple","symmetry-fanorm",&q,PETSC_NULL,&f); CHKERRQ(ierr);
  if (!f) goto exit;
  nnorm = q.r;
  bound = sqrt(2) * PetscMin(mnorm,nnorm);

  ierr = GetDataID("normal","lee95-bound",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataSetReal
    ((PetscObject)A,id,bound); CHKERRQ(ierr);

 exit:
  ierr = MatDestroy(AA); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

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().

{
  Mat AA; PetscReal upper,lower;
  PetscScalar tracea2; PetscReal cnorm,fnorm; int id,ql; 
  PetscTruth f,has; PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = MatDuplicate(A,MAT_COPY_VALUES,&AA); CHKERRQ(ierr);
  ierr = MatCenter(AA); CHKERRQ(ierr);
  ierr = TraceA2(AA,(AnalysisItem*)&tracea2,&ql,&f); CHKERRQ(ierr);
  if (!f) {
    /*printf(">>warning: did not compute trace of A squared\n");*/
    goto exit;
  }
  ierr = MatNorm(AA,NORM_FROBENIUS,&fnorm); CHKERRQ(ierr);

  upper = fnorm*fnorm-PetscAbsScalar(tracea2);
  ierr = Commutator(A,(AnalysisItem*)&cnorm,&ql,&f); CHKERRQ(ierr);
  if (!f) {
    /*printf(">>warning: did not compute commutator Frobenius norm\n");*/
    goto exit;
  }
  fnorm *= fnorm;
  lower = sqrt(fnorm - sqrt(fnorm*fnorm-cnorm*cnorm/2));

  ierr = GetDataID("normal","lee96-ubound",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataSetReal
    ((PetscObject)A,id,upper); CHKERRQ(ierr);
  ierr = GetDataID("normal","lee96-lbound",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataSetReal
    ((PetscObject)A,id,lower); CHKERRQ(ierr);

 exit:
  ierr = MatDestroy(AA); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

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().

{
  MPI_Comm comm; AnalysisItem q; PetscTruth f;
  Vec D; PetscScalar trace,*d; int N,n,i; PetscErrorCode ierr;

  PetscFunctionBegin;
  ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
  ierr = MatGetLocalSize(A,&n,PETSC_NULL); CHKERRQ(ierr);
  ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr);

  ierr = MatrixComputeQuantity
    (A,"simple","trace",&q,PETSC_NULL,&f); CHKERRQ(ierr);
  if (!f) 
    SETERRQ(1,"Failed to compute trace");
  trace = q.r;

  ierr = VecCreateMPI(comm,n,PETSC_DECIDE,&D); CHKERRQ(ierr);
  ierr = MatGetDiagonal(A,D); CHKERRQ(ierr);
  ierr = VecGetArray(D,&d); CHKERRQ(ierr);
  for (i=0; i<n; i++) d[i] -= trace/N;
  ierr = VecRestoreArray(D,&d); CHKERRQ(ierr);
  ierr = MatDiagonalSet(A,D,INSERT_VALUES); CHKERRQ(ierr);
  ierr = VecDestroy(D); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

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().

{
  MPI_Comm comm; Mat Awork; PetscReal result;
  PetscTruth alloc_work,do_work,any_work,has;
  int id; PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
  ierr = AnaModGetSequentialMatrix
    (A,&Awork,&alloc_work,&do_work,&any_work); CHKERRQ(ierr);
  if (!any_work) {
    *done = PETSC_FALSE; PetscFunctionReturn(0);
  }
  if (do_work) {
    ierr = MatCommutatorNormF_seq(Awork,&result,done); CHKERRQ(ierr);
    /*printf("result %d:%e\n",*done,result);*/
  }
  if (alloc_work) {
    ierr = MatDestroy(Awork); CHKERRQ(ierr);
  }

  MPI_Bcast((void*)done,1,MPI_INT,0,comm);
  if (*done) {
    MPI_Bcast((void*)&result,1,MPIU_REAL,0,comm);
  } else PetscFunctionReturn(0);

  ierr = GetDataID("normal","commutator-normF",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataSetReal
    ((PetscObject)A,id,result); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

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().

{
  Mat B,E,F;
  PetscReal vsum; PetscTruth force,flg;
  int afirst,alast,bfirst,blast,p,q,N;
  PetscErrorCode ierr;

  PetscFunctionBegin;

  ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr);

  ierr = AnaModHasForcedExpensiveComputation(&force); CHKERRQ(ierr);
  ierr = MatrixComputeQuantity
    (A,"structure","left-bandwidth",(AnalysisItem*)&p,PETSC_NULL,&flg); CHKERRQ(ierr);
  if (!flg) p = N;
  ierr = MatrixComputeQuantity
    (A,"structure","right-bandwidth",(AnalysisItem*)&q,PETSC_NULL,&flg); CHKERRQ(ierr);
  if (!flg) q = N;
  if (p+q>sqrt_times*sqrt(N) && !force) {
    /*printf("too expensive\n");*/
    *done = PETSC_FALSE;
    PetscFunctionReturn(0);
  }

  /* set up matrix and transpose */
  ierr = MatDuplicate(A,MAT_COPY_VALUES,&B); CHKERRQ(ierr);
  afirst = 0; alast = N; bfirst = 0; blast = N;
  ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&E); CHKERRQ(ierr);
  ierr = MatDuplicate(E,MAT_COPY_VALUES,&F); CHKERRQ(ierr);
#define min(a,b) (a<b?a:b)
#define max(a,b) (a>b?a:b)
  {
    const PetscReal *avals,*bvals,*evals,*fvals; PetscScalar v1,v2;
    int a,b, nacols,nbcols, necols,nfcols;
    const int *acols,*bcols,*ecols,*fcols;
    vsum = 0;
    for (a=afirst; a<alast; a++) {
      int bmin=max(bfirst,a-(p+q)),bmax=min(blast,a+(p+q));
      ierr = MatGetRow(A,a,&nacols,&acols,&avals); CHKERRQ(ierr);
      ierr = MatGetRow(E,a,&necols,&ecols,&evals); CHKERRQ(ierr);
      for (b=bmin; b<bmax; b++) {
        ierr = MatGetRow(B,b,&nbcols,&bcols,&bvals); CHKERRQ(ierr);
        ierr = MatGetRow(F,b,&nfcols,&fcols,&fvals); CHKERRQ(ierr);
        ierr = SparseVecProd
          (nacols,acols,avals,nbcols,bcols,bvals,&v1); CHKERRQ(ierr);
        ierr = SparseVecProd
          (necols,ecols,evals,nfcols,fcols,fvals,&v2); CHKERRQ(ierr);
        v1 -= v2;
        v1 = PetscAbsScalar(v1);
        vsum += v1*v1;
        ierr = MatRestoreRow(B,b,&nbcols,&bcols,&bvals); CHKERRQ(ierr);
        ierr = MatRestoreRow(F,b,&nfcols,&fcols,&fvals); CHKERRQ(ierr);
      }
      ierr = MatRestoreRow(A,a,&nacols,&acols,&avals); CHKERRQ(ierr);
      ierr = MatRestoreRow(E,a,&necols,&ecols,&evals); CHKERRQ(ierr);
    }
  }
  ierr = MatDestroy(B); CHKERRQ(ierr);
  ierr = MatDestroy(E); CHKERRQ(ierr);
  ierr = MatDestroy(F); CHKERRQ(ierr);
  *result = sqrt(vsum); *done = PETSC_TRUE;

  PetscFunctionReturn(0);
}

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().

{
  PetscErrorCode ierr;
  PetscFunctionBegin;

  ierr = RegisterModule
    ("normal","trace-asquared",ANALYSISDOUBLE,&TraceA2); CHKERRQ(ierr);
  ierr = RegisterModule
    ("normal","commutator-normF",ANALYSISDOUBLE,&Commutator); CHKERRQ(ierr);

  ierr = RegisterModule
    ("normal","ruhe75-bound",ANALYSISDOUBLE,&DepartureRuhe75); CHKERRQ(ierr);

  ierr = RegisterModule
    ("normal","lee95-bound",ANALYSISDOUBLE,&DepartureLee95); CHKERRQ(ierr);

  ierr = RegisterModule
    ("normal","lee96-ubound",ANALYSISDOUBLE,&DepartureLee96U); CHKERRQ(ierr);
  ierr = RegisterModule
    ("normal","lee96-lbound",ANALYSISDOUBLE,&DepartureLee96L); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

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().

{
  int pa = 0, pb = 0; PetscScalar v = 0;
  PetscFunctionBegin;

  while (pa<na && pb<nb) {
    if (ia[pa]==ib[pb]) {
      v += va[pa]*vb[pb]; pa++; pb++;
    } else if (ia[pa]<ib[pb]) pa++;
    else pb++;
  }
  *result = v;

  PetscFunctionReturn(0);
}
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().

{
  Mat A = (Mat)prob;
  PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr;
  PetscFunctionBegin;
  ierr = GetDataID("normal","trace-asquared",&id,&has); CHKERRQ(ierr);
  HASTOEXIST(has);
  ierr = PetscObjectComposedDataGetScalar
    ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
  if (lv) *lv = 1;
  if (*flg) rv->r = v;
  else {
    ierr = compute_tracea2(A,flg); CHKERRQ(ierr);
    if (flg) {
      ierr = PetscObjectComposedDataGetScalar
        ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
      if (*flg) rv->r = v;
    }
  }
  
  PetscFunctionReturn(0);
}

Here is the call graph for this function:


Variable Documentation

int sqrt_times = 50 [static]

Definition at line 274 of file normal.c.