#include <stdlib.h>
#include "anamod.h"
#include "anamodsalsamodules.h"
#include "anamatrix.h"
#include "petscerror.h"
#include "petscksp.h"
#include "petscmat.h"
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 |
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:
@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 max | ( | a, | |||
b | ) | (a>b?a:b) |
Referenced by MatCommutatorNormF_seq().
#define min | ( | a, | |||
b | ) | (a<b?a:b) |
Referenced by MatCommutatorNormF_seq().
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 }
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 }
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 }
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 }
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 }
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 }
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 }
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 }
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 }
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 }
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 }
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 }
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 }
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 }
int sqrt_times = 50 [static] |