SALSA Analysis Modules
simple.c
Go to the documentation of this file.
00001 /*! \file simple.c \ingroup categories 
00002   \brief Norm-like properties
00003 
00004   \section simple Simple (normlike) quantities
00005 
00006   The "simple" module contains simple, normlike quantities,
00007   In general,
00008   these are properties that can be calculated in time proportional to the
00009   number of nonzeros of the matrix. Note that this excludes the 2-norm.
00010 
00011   \section Usage
00012 
00013   Activate this module with
00014 
00015   PetscErrorCode RegisterSimpleModules();
00016 
00017   Compute these elements with
00018 
00019   ComputeQuantity("simple",<element>,A,(AnalysisItem*)&res,&flg);
00020 
00021   Available elements are:
00022   - "trace" : Trace, sum of diagonal elements; see Trace().
00023   - "trace-abs" : Trace Absolute, sum of absolute values of diagonal elements;
00024     see TraceAbs().
00025   - "norm1" : 1-Norm, maximum column sum of absolute element sizes;
00026   see norm1().
00027   - "normInf" : Infinity-Norm, maximum row sum of absolute element sizes;
00028   see normInf().
00029   - "normF" : Frobenius-norm, sqrt of sum of elements squared;
00030   see normF().
00031   - "diagonal-dominance" : least positive or most negative value of
00032     diagonal element minus sum of absolute off-diagonal elements;
00033     see diagonaldominance().
00034   - "symmetry-snorm" : infinity norm of symmetric part; see SymmetrySNorm().
00035   This element and all following are sequential; see \ref options;
00036   - "symmetry-anorm" : infinity norm of anti-symmetric part; see SymmetryANorm().
00037   - "symmetry-fsnorm" : Frobenius norm of symmetric part; see
00038   SymmetryFSNorm().
00039   - "symmetry-fanorm" : Frobenius norm of anti-symmetric part;
00040   see SymmetryFANorm().
00041   - "n-struct-unsymm" : number of structurally unsymmetric elements;
00042   see NUnstruct().
00043     (this is calculated here, but declared as an element of the \ref structure
00044     category.)
00045  */
00046 
00047 #include <stdlib.h>
00048 #include "anamod.h"
00049 #include "anamodsalsamodules.h"
00050 #include "petscerror.h"
00051 #include "petscmat.h"
00052 #include "src/mat/impls/aij/seq/aij.h"
00053 
00054 #undef __FUNCT__
00055 #define __FUNCT__ "computetrace"
00056 /*!
00057   This is the computational routine for trace calculation.
00058   It computes both the trace and the absolute trace, that is, using
00059   absolute values of matrix elements.
00060 
00061   A Petsc vector is created and destroyed.
00062 */
00063 static PetscErrorCode computetrace(Mat A)
00064 {
00065   MPI_Comm comm; Vec D; PetscScalar *d;
00066   PetscBool has; PetscScalar vsum,gvsum; PetscReal asum,gasum;
00067   int local_size,i,id;
00068   PetscErrorCode ierr;
00069 
00070   PetscFunctionBegin;
00071 
00072   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00073   ierr = MatGetLocalSize(A,&local_size,PETSC_NULL); CHKERRQ(ierr);
00074   ierr = VecCreateMPI(comm,local_size,PETSC_DECIDE,&D); CHKERRQ(ierr);
00075   ierr = MatGetDiagonal(A,D); CHKERRQ(ierr);
00076   ierr = VecGetArray(D,&d); CHKERRQ(ierr);
00077   vsum = 0; asum = 0;
00078   for (i=0; i<local_size; i++) {
00079     vsum += d[i]; asum += PetscAbsScalar(d[i]);
00080   }
00081   ierr = VecRestoreArray(D,&d); CHKERRQ(ierr);
00082   ierr = VecDestroy(&D); CHKERRQ(ierr);
00083   ierr = MPI_Allreduce
00084     ((void*)&vsum,(void*)&gvsum,1,MPIU_SCALAR,MPI_SUM,comm); CHKERRQ(ierr);
00085   ierr = MPI_Allreduce
00086     ((void*)&asum,(void*)&gasum,1,MPI_DOUBLE,MPI_SUM,comm); CHKERRQ(ierr);
00087 
00088   ierr = GetDataID("simple","trace",&id,&has); CHKERRQ(ierr);
00089   HASTOEXIST(has);
00090   ierr = PetscObjectComposedDataSetScalar
00091     ((PetscObject)A,id,gvsum); CHKERRQ(ierr);
00092   ierr = GetDataID("simple","trace-abs",&id,&has); CHKERRQ(ierr);
00093   HASTOEXIST(has);
00094   ierr = PetscObjectComposedDataSetReal
00095     ((PetscObject)A,id,gasum); CHKERRQ(ierr);
00096 
00097   PetscFunctionReturn(0);
00098 }
00099 
00100 #undef __FUNCT__
00101 #define __FUNCT__ "Trace"
00102 /*! Compute and store the trace of a matrix; 
00103   the computation is done in computetrace().
00104 */
00105 static PetscErrorCode Trace
00106 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00107 {
00108   Mat A = (Mat)prob;
00109   PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr;
00110   PetscFunctionBegin;
00111   ierr = GetDataID("simple","trace",&id,&has); CHKERRQ(ierr);
00112   HASTOEXIST(has);
00113   ierr = PetscObjectComposedDataGetScalar
00114     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00115   if (*flg) {
00116     rv->r = v;
00117   } else {
00118     ierr = computetrace(A); CHKERRQ(ierr);
00119     ierr = PetscObjectComposedDataGetScalar
00120       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00121     if (*flg) {
00122       rv->r = v;
00123     } else SETERRQ(MPI_COMM_WORLD,1,"Should have computed trace by now");
00124   }
00125   
00126   PetscFunctionReturn(0);
00127 }
00128 
00129 #undef __FUNCT__
00130 #define __FUNCT__ "TraceAbs"
00131 /*! Compute and store the absolute trace of a matrix; 
00132   the computation is done in computetrace().
00133 */
00134 static PetscErrorCode TraceAbs
00135 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00136 {
00137   Mat A = (Mat)prob;
00138   PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr;
00139   PetscFunctionBegin;
00140   ierr = GetDataID("simple","trace-abs",&id,&has); CHKERRQ(ierr);
00141   HASTOEXIST(has);
00142   ierr = PetscObjectComposedDataGetReal
00143     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00144   if (*flg) rv->r = v;
00145   else {
00146     ierr = computetrace(A); CHKERRQ(ierr);
00147     ierr = PetscObjectComposedDataGetReal
00148       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00149     if (*flg) rv->r = v;
00150   }
00151   
00152   PetscFunctionReturn(0);
00153 }
00154 
00155 #undef __FUNCT__
00156 #define __FUNCT__ "norm1"
00157 /*! Compute the 1-norm of a matrix. */
00158 static PetscErrorCode norm1
00159 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00160 {
00161   Mat A = (Mat)prob;
00162   PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr;
00163   PetscFunctionBegin;
00164 
00165   ierr = GetDataID("simple","norm1",&id,&has); CHKERRQ(ierr);
00166   HASTOEXIST(has);
00167   ierr = PetscObjectComposedDataGetReal
00168     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00169   if (*flg) rv->r = v;
00170   else {
00171     ierr = MatNorm(A,NORM_1,&v); CHKERRQ(ierr);
00172     ierr = PetscObjectComposedDataSetReal
00173       ((PetscObject)A,id,v); CHKERRQ(ierr);
00174     *flg = PETSC_TRUE;
00175     rv->r = v;
00176   }
00177   
00178   *flg = PETSC_TRUE;
00179   PetscFunctionReturn(0);
00180 }
00181 
00182 #undef __FUNCT__
00183 #define __FUNCT__ "normF"
00184 /*! Compute the Frobenius norm of a matrix. */
00185 static PetscErrorCode normF
00186 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00187 {
00188   Mat A = (Mat)prob;
00189   PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr;
00190   PetscFunctionBegin;
00191 
00192   ierr = GetDataID("simple","normF",&id,&has); CHKERRQ(ierr);
00193   HASTOEXIST(has);
00194   ierr = PetscObjectComposedDataGetReal
00195     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00196   if (*flg) rv->r = v;
00197   else {
00198     ierr = MatNorm(A,NORM_FROBENIUS,&v); CHKERRQ(ierr);
00199     ierr = PetscObjectComposedDataSetReal
00200       ((PetscObject)A,id,v); CHKERRQ(ierr);
00201     *flg = PETSC_TRUE;
00202     rv->r = v;
00203   }
00204   
00205   *flg = PETSC_TRUE;
00206   PetscFunctionReturn(0);
00207 }
00208 
00209 #undef __FUNCT__
00210 #define __FUNCT__ "normInf"
00211 /*! Compute the infinity-norm of a matrix. */
00212 static PetscErrorCode normInf
00213 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00214 {
00215   Mat A = (Mat)prob;
00216   PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr;
00217   PetscFunctionBegin;
00218 
00219   ierr = GetDataID("simple","normInf",&id,&has); CHKERRQ(ierr);
00220   HASTOEXIST(has);
00221   ierr = PetscObjectComposedDataGetReal
00222     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00223   if (*flg) rv->r = v;
00224   else {
00225     ierr = MatNorm(A,NORM_INFINITY,&v); CHKERRQ(ierr);
00226     ierr = PetscObjectComposedDataSetReal
00227       ((PetscObject)A,id,v); CHKERRQ(ierr);
00228     *flg = PETSC_TRUE;
00229     rv->r = v;
00230   }
00231   
00232   *flg = PETSC_TRUE;
00233   PetscFunctionReturn(0);
00234 }
00235 
00236 #undef __FUNCT__
00237 #define __FUNCT__ "compute_dd"
00238 /*!
00239   Compute the diagonal dominance of a matrix:
00240   minimum value of diagonal element (not absolute!) minus off-diagonal
00241   elements absolute.
00242   \latexonly
00243   \[ \min_i a_{ii}-\sum_{j\not=i}|a_{ij}| \]
00244   \endlatexonly
00245  */
00246 static PetscErrorCode compute_dd(Mat A,PetscBool *flg)
00247 {
00248   PetscReal r,gr;
00249   MPI_Comm comm;
00250   int first,last,m,row; PetscErrorCode ierr;
00251 
00252   PetscFunctionBegin;
00253   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00254   ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00255   m = last-first;
00256   for (row=first; row<last; row++) {
00257     int col,ncols; const int *cols; 
00258     const PetscScalar *vals; PetscReal rr=0.;
00259     ierr = MatGetRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00260     for (col=0; col<ncols; col++) {
00261       PetscReal v = vals[col], vr = PetscRealPart(v), va = PetscAbsScalar(v);
00262       if (cols[col]==row) rr +=vr; else rr -= va;
00263     }
00264     if (row==first) r = rr; else r = PetscMin(r,rr);
00265     ierr = MatRestoreRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00266   }
00267   ierr = MPI_Allreduce((void*)&r,(void*)&gr,1,MPIU_REAL,MPI_MIN,comm);
00268   {
00269     int id; PetscBool has;
00270     ierr = GetDataID("simple","diagonal-dominance",&id,&has); CHKERRQ(ierr);
00271     HASTOEXIST(has);
00272     ierr = PetscObjectComposedDataSetReal
00273       ((PetscObject)A,id,gr); CHKERRQ(ierr);
00274   }
00275   *flg = PETSC_TRUE;
00276   PetscFunctionReturn(0);
00277 }
00278 
00279 #undef __FUNCT__
00280 #define __FUNCT__ "DiagonalDominance"
00281 static PetscErrorCode DiagonalDominance
00282 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00283 {
00284   Mat A = (Mat)prob;
00285   PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr;
00286   PetscFunctionBegin;
00287 
00288   ierr = GetDataID("simple","diagonal-dominance",&id,&has); CHKERRQ(ierr);
00289   HASTOEXIST(has);
00290   ierr = PetscObjectComposedDataGetReal
00291     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00292   if (*flg) rv->r = v;
00293   else {
00294     ierr = compute_dd(A,flg); CHKERRQ(ierr);
00295     if (*flg) {
00296       ierr = PetscObjectComposedDataGetReal
00297         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00298       if (*flg) rv->r = v;
00299     }
00300   }
00301   
00302   PetscFunctionReturn(0);
00303 }
00304 
00305 /*
00306  * Norms of symmetric / unsymmetric part
00307  */
00308 #undef __FUNCT__
00309 #define __FUNCT__ "MatSymmPartNormInf_seq"
00310 /*!
00311   This is where the real work of MatSymmPartNormInf() is done
00312 
00313   A few temporary arrays of the 
00314   size of a vector are allocated and freed. Otherwise this is just a whole
00315   bunch of pointer chasing.
00316 */
00317 static PetscErrorCode MatSymmPartNormInf_seq
00318 (Mat A,PetscReal *sn_r,PetscReal *an_r,PetscReal *srn_r,PetscReal *arn_r,
00319  int *nun_r,PetscBool *done)
00320 {
00321   Mat_SeqAIJ *aij;  
00322   int *adx,*bdx,*aii,*bii,*aptr,*bptr,ma,na,i,nun; PetscScalar *va,*vb;
00323   PetscReal snorm,anorm,*sronorm,*aronorm;
00324   PetscErrorCode ierr;
00325   PetscFunctionBegin;
00326 
00327   aij = (Mat_SeqAIJ *) A->data;
00328   ierr = MatGetLocalSize(A,&ma,&na); CHKERRQ(ierr);
00329 
00330   /* two arrays of pointers into the matrix */
00331   aii = bii = aij->i;
00332   adx = bdx = aij->j;
00333   va = vb = aij->a;
00334   ierr = PetscMalloc(ma*sizeof(int),&aptr); CHKERRQ(ierr);
00335   ierr = PetscMalloc(ma*sizeof(int),&bptr); CHKERRQ(ierr);
00336   for (i=0; i<ma; i++) aptr[i] = bptr[i] = aii[i];
00337 
00338   /* temp storage for the inf norm */
00339   ierr = PetscMalloc(ma*sizeof(PetscReal),&sronorm); CHKERRQ(ierr);
00340   ierr = PetscMalloc(ma*sizeof(PetscReal),&aronorm); CHKERRQ(ierr);
00341   ierr = PetscMemzero(sronorm,ma*sizeof(PetscReal)); CHKERRQ(ierr);
00342   ierr = PetscMemzero(aronorm,ma*sizeof(PetscReal)); CHKERRQ(ierr);
00343 
00344   /* update row & frobenius norms with a(i,j)=v and a(j,i)=w */
00345 #if defined(ANAMODDEBUG)
00346 #define RC_NORM(i,j,v,w) \
00347   printf("s/a norm calculation: (i,j)=(%d,%d), (vij,vji)=(%e,%e)\n",i,j,v,w); \
00348   {                                                                      \
00349     PetscReal sum = PetscAbsScalar(v+w)/2, dif = PetscAbsScalar(v-w)/2;  \
00350     ASSERT3(!(i==j && v!=w),"same elements on diagonal",i,j);            \
00351     sronorm[i] += sum; if (i!=j) sronorm[j] += sum;                      \
00352     aronorm[i] += dif; if (i!=j) aronorm[j] += dif;                      \
00353     if (i==j) snorm += sum*sum;                                          \
00354     else {snorm += 2*sum*sum; anorm += 2*dif*dif;}                       \
00355   }
00356 #else
00357 #define RC_NORM(i,j,v,w) \
00358   {                                                                     \
00359     PetscReal sum = PetscAbsScalar(v+w)/2, dif = PetscAbsScalar(v-w)/2; \
00360     ASSERT3(!(i==j && v!=w),"same elements on diagonal",i,j);            \
00361     sronorm[i] += sum; if (i!=j) sronorm[j] += sum;                     \
00362     aronorm[i] += dif; if (i!=j) aronorm[j] += dif;                     \
00363     if (i==j) snorm += sum*sum;                                         \
00364     else {snorm += 2*sum*sum; anorm += 2*dif*dif;}                      \
00365   }
00366 #endif
00367 #define ASSERT3(v,s,i,j) \
00368   if (!(v)) SETERRQ3(MPI_COMM_WORLD,1,"Violation of <%s> @ (%d,%d): %d",s,i,j);
00369 #define ASSERT(v,s,i,j,t) \
00370   if (!(v)) SETERRQ4(MPI_COMM_WORLD,1,"Violation of <%s> @ (%d,%d): %d",s,i,j,t);
00371 
00372   snorm = anorm = 0; nun = 0;
00373   for (i=0; i<ma; i++) {
00374     int idc,idr; PetscScalar vc,vr;
00375     /* get the a-pointer into the upper triangle */
00376     while (aptr[i]<aii[i+1] && adx[aptr[i]]<i) aptr[i]++;
00377     ASSERT3(aptr[i]==aii[i+1]||adx[aptr[i]]>=i,"A in upper0",i,adx[aptr[i]]);
00378     /* incorporate L elements that are not structurally symmetric */
00379     while (bptr[i]<bii[i+1] && bdx[bptr[i]]<i) {
00380       vr = PetscAbsScalar(vb[bptr[i]]); nun++;
00381       RC_NORM(i,bdx[bptr[i]],vr,0.);
00382       bptr[i]++;
00383     }
00384     ASSERT3(bptr[i]==bii[i+1] || bdx[bptr[i]]>=i,"B part done",i,bdx[bptr[i]]);
00385     while (aptr[i]<aii[i+1]) {
00386       /* in row i we are at element j = idc @ adx( aptr[i] );
00387          vc = A( i, j ) */
00388       idc = adx[aptr[i]]; vc = va[aptr[i]];
00389       ASSERT(idc>=i,"A in upper",i,idc,i-idc);
00390       /* in row j=idc look at elements i' = idr @ bdx(bptr[idc]); vr = B(j,i')
00391          first catch up with structurally unsymmetric elements */
00392       idr = bdx[bptr[idc]]; vr = vb[bptr[idc]];
00393       while (idr<i) {
00394         RC_NORM(idc,idr,vr,0.);
00395         bptr[idc]++; idr = bdx[bptr[idc]]; vr = vb[bptr[idc]]; nun++;
00396       }     
00397       /* now see if (i,j) == (j,i') */
00398       if (idr==i) {
00399         /* structurally symmetric; process these elements of L and U */
00400         RC_NORM(i,idc,vc,vr);
00401         bptr[idc]++;
00402       } else {
00403         /* element (j,i') in L corresponds to a later row in U,
00404            so we have a structurally unsymmetric element in U */
00405         nun++; RC_NORM(i,idc,vc,0.);
00406       }
00407       aptr[i]++;
00408     }
00409     ASSERT(aptr[i]==aii[i+1],"Aptr has traversed row i",i+1,adx[aptr[i]],aii[i+1]);
00410   }
00411 
00412   ierr = PetscFree(aptr); CHKERRQ(ierr);
00413   ierr = PetscFree(bptr); CHKERRQ(ierr);
00414 
00415   snorm = sqrt(snorm); anorm = sqrt(anorm);
00416 
00417   for (i=1; i<ma; i++) {
00418     if (sronorm[i]>sronorm[0]) sronorm[0] = sronorm[i];
00419     if (aronorm[i]>aronorm[0]) aronorm[0] = aronorm[i];
00420   }
00421 
00422   *sn_r = snorm; *an_r = anorm; *srn_r = sronorm[0]; *arn_r = aronorm[0];
00423   *nun_r = nun;
00424   *done = PETSC_TRUE;
00425 
00426   ierr = PetscFree(sronorm); CHKERRQ(ierr);
00427   ierr = PetscFree(aronorm); CHKERRQ(ierr);
00428 
00429   PetscFunctionReturn(0);
00430 }
00431 
00432 #undef __FUNCT__
00433 #define __FUNCT__ "MatSymmPartNormInf"
00434 /*!
00435   This is the computational routine for all quantities relating to
00436   symmetric and anti-symmetric parts of the matrix.
00437 
00438   This computation can only be done in single-processor mode.
00439   If the commandline
00440   option (see \ref options) \c "-anamod-force sequential" is specified,
00441   the matrix is gathered on processor zero to compute this quantity.
00442 
00443   See MatSymmPartNormInf_seq() for the actual computation.
00444 */
00445 static PetscErrorCode MatSymmPartNormInf(Mat A,PetscBool *done)
00446 {
00447   MPI_Comm comm; Mat Awork; int id,nun;
00448   PetscReal snorm,anorm,sronorm,aronorm;
00449   PetscBool has,alloc_work,do_work,any_work; PetscErrorCode ierr;
00450 
00451   PetscFunctionBegin;
00452 
00453   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00454   ierr = AnaModGetSequentialMatrix
00455     (A,&Awork,&alloc_work,&do_work,&any_work); CHKERRQ(ierr);
00456   if (!any_work) {
00457     *done = PETSC_FALSE; PetscFunctionReturn(0);
00458   }
00459   if (do_work) {
00460     ierr = MatSymmPartNormInf_seq
00461       (Awork,&snorm,&anorm,&sronorm,&aronorm,&nun,done); CHKERRQ(ierr);
00462   }
00463   if (alloc_work) {
00464     ierr = MatDestroy(&Awork); CHKERRQ(ierr);
00465   }
00466 
00467   ierr = MPI_Bcast((void*)done,1,MPI_INT,0,comm);
00468   ierr = MPI_Bcast((void*)&snorm,1,MPIU_REAL,0,comm);
00469   ierr = MPI_Bcast((void*)&anorm,1,MPIU_REAL,0,comm);
00470   ierr = MPI_Bcast((void*)&sronorm,1,MPIU_REAL,0,comm);
00471   ierr = MPI_Bcast((void*)&aronorm,1,MPIU_REAL,0,comm);
00472   ierr = MPI_Bcast((void*)&nun,1,MPI_INT,0,comm);
00473 
00474   /* infinity norms of symmetric / anti-symmetric part */
00475   ierr = GetDataID("simple","symmetry-snorm",&id,&has); CHKERRQ(ierr);
00476   HASTOEXIST(has);
00477   ierr = PetscObjectComposedDataSetReal
00478     ((PetscObject)A,id,sronorm); CHKERRQ(ierr);
00479 
00480   ierr = GetDataID("simple","symmetry-anorm",&id,&has); CHKERRQ(ierr);
00481   HASTOEXIST(has);
00482   ierr = PetscObjectComposedDataSetReal
00483     ((PetscObject)A,id,aronorm); CHKERRQ(ierr);
00484 
00485   /* frobenius norms of symmetric / anti-symmetric part */
00486   ierr = GetDataID("simple","symmetry-fsnorm",&id,&has); CHKERRQ(ierr);
00487   HASTOEXIST(has);
00488   ierr = PetscObjectComposedDataSetReal
00489     ((PetscObject)A,id,snorm); CHKERRQ(ierr);
00490 
00491   ierr = GetDataID("simple","symmetry-fanorm",&id,&has); CHKERRQ(ierr);
00492   HASTOEXIST(has);
00493   ierr = PetscObjectComposedDataSetReal
00494     ((PetscObject)A,id,anorm); CHKERRQ(ierr);
00495 
00496   ierr = GetDataID("structure","n-struct-unsymm",&id,&has); CHKERRQ(ierr);
00497   HASTOEXIST(has);
00498   ierr = PetscObjectComposedDataSetInt
00499     ((PetscObject)A,id,nun); CHKERRQ(ierr);
00500 
00501   PetscFunctionReturn(0);
00502 }
00503 
00504 #undef __FUNCT__
00505 #define __FUNCT__ "SymmetrySNorm"
00506 /*! Compute and store the infinity norm of the symmetric part of the matrix.
00507   The actual computation is done in MatSymmPartNormInf();
00508   see that routine for remarks on parallelism.
00509 */
00510 static PetscErrorCode SymmetrySNorm
00511 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00512 {
00513   Mat A = (Mat)prob;
00514   PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr;
00515   PetscFunctionBegin;
00516   ierr = GetDataID("simple","symmetry-snorm",&id,&has); CHKERRQ(ierr);
00517   HASTOEXIST(has);
00518   ierr = PetscObjectComposedDataGetReal
00519     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00520   if (*flg) rv->r = v;
00521   else {
00522     ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr);
00523     if (*flg) {
00524       ierr = PetscObjectComposedDataGetReal
00525         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00526       if (*flg) rv->r = v;
00527     }
00528   }
00529   
00530   PetscFunctionReturn(0);
00531 }
00532 
00533 #undef __FUNCT__
00534 #define __FUNCT__ "SymmetryANorm"
00535 /*! Compute and store the infinity norm of the antisymmetric part of the matrix.
00536   The actual computation is done in MatSymmPartNormInf();
00537   see that routine for remarks on parallelism.
00538 */
00539 static PetscErrorCode SymmetryANorm
00540 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00541 {
00542     Mat A = (Mat)prob;
00543   PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr;
00544   PetscFunctionBegin;
00545   ierr = GetDataID("simple","symmetry-anorm",&id,&has); CHKERRQ(ierr);
00546   HASTOEXIST(has);
00547   ierr = PetscObjectComposedDataGetReal
00548     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00549   if (*flg) rv->r = v;
00550   else {
00551     ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr);
00552     if (*flg) {
00553       ierr = PetscObjectComposedDataGetReal
00554         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00555       if (*flg) rv->r = v;
00556     }
00557   }
00558   
00559   PetscFunctionReturn(0);
00560 }
00561 
00562 #undef __FUNCT__
00563 #define __FUNCT__ "SymmetryFSNorm"
00564 /*! Compute and store the Frobenius norm of the symmetric part of the matrix.
00565   The actual computation is done in MatSymmPartNormInf();
00566   see that routine for remarks on parallelism.
00567 */
00568 static PetscErrorCode SymmetryFSNorm
00569 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00570 {
00571     Mat A = (Mat)prob;
00572   PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr;
00573   PetscFunctionBegin;
00574   ierr = GetDataID("simple","symmetry-fsnorm",&id,&has); CHKERRQ(ierr);
00575   HASTOEXIST(has);
00576   ierr = PetscObjectComposedDataGetReal
00577     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00578   if (*flg) rv->r = v;
00579   else {
00580     ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr);
00581     if (*flg) {
00582       ierr = PetscObjectComposedDataGetReal
00583         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00584       if (*flg) rv->r = v;
00585     }
00586   }
00587   
00588   PetscFunctionReturn(0);
00589 }
00590 
00591 #undef __FUNCT__
00592 #define __FUNCT__ "SymmetryFANorm"
00593 /*! Compute and store the Frobenius norm of the antisymmetric part of the matrix.
00594   The actual computation is done in MatSymmPartNormInf();
00595   see that routine for remarks on parallelism.
00596 */
00597 static PetscErrorCode SymmetryFANorm
00598 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00599 {
00600     Mat A = (Mat)prob;
00601   PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr;
00602   PetscFunctionBegin;
00603   ierr = GetDataID("simple","symmetry-fanorm",&id,&has); CHKERRQ(ierr);
00604   HASTOEXIST(has);
00605   ierr = PetscObjectComposedDataGetReal
00606     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00607   if (*flg) rv->r = v;
00608   else {
00609     ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr);
00610     if (*flg) {
00611       ierr = PetscObjectComposedDataGetReal
00612         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00613       if (*flg) rv->r = v;
00614     }
00615   }
00616   
00617   PetscFunctionReturn(0);
00618 }
00619 
00620 #undef __FUNCT__
00621 #define __FUNCT__ "NUnstruct"
00622 /*! Compute and store the number of structurally unsymmetric elements of a matrix.
00623   The actual computation is done in MatSymmPartNormInf();
00624   see that routine for remarks on parallelism.
00625 
00626   The "n-struct-unsymm" feature is declared to be in category "structure",
00627   rather than simple.
00628 */
00629 static PetscErrorCode NUnstruct
00630 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00631 {
00632   Mat A = (Mat)prob;
00633   PetscReal v=0; PetscBool has; int id; PetscErrorCode ierr;
00634   PetscFunctionBegin;
00635   ierr = GetDataID("structure","n-struct-unsymm",&id,&has); CHKERRQ(ierr);
00636   HASTOEXIST(has);
00637   ierr = PetscObjectComposedDataGetInt
00638     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00639   if (*flg) rv->i = v;
00640   else {
00641     ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr);
00642     if (*flg) {
00643       ierr = PetscObjectComposedDataGetInt
00644         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00645       if (*flg) rv->i = v;
00646     }
00647   }
00648   
00649   PetscFunctionReturn(0);
00650 }
00651 
00652 #undef __FUNCT__
00653 #define __FUNCT__ "RegisterSimpleModules"
00654 /*!
00655   Declare norm-like modules that can be performed with simple calculations.
00656 */
00657 PetscErrorCode RegisterSimpleModules()
00658 {
00659   PetscErrorCode ierr;
00660   PetscFunctionBegin;
00661 
00662   ierr = RegisterModule
00663     ("simple","trace",ANALYSISDOUBLE,&Trace); CHKERRQ(ierr);
00664   ierr = RegisterModule
00665     ("simple","trace-abs",ANALYSISDOUBLE,&TraceAbs); CHKERRQ(ierr);
00666 
00667   ierr = RegisterModule
00668     ("simple","norm1",ANALYSISDOUBLE,&norm1); CHKERRQ(ierr);
00669   ierr = RegisterModule
00670     ("simple","normInf",ANALYSISDOUBLE,&normInf); CHKERRQ(ierr);
00671   ierr = RegisterModule
00672     ("simple","normF",ANALYSISDOUBLE,&normF); CHKERRQ(ierr);
00673 
00674   ierr = RegisterModule
00675     ("simple","diagonal-dominance",ANALYSISDOUBLE,&DiagonalDominance); CHKERRQ(ierr);
00676 
00677   ierr = RegisterModule
00678     ("simple","symmetry-snorm",ANALYSISDOUBLE,&SymmetrySNorm); CHKERRQ(ierr);
00679   ierr = RegisterModule
00680     ("simple","symmetry-anorm",ANALYSISDOUBLE,&SymmetryANorm); CHKERRQ(ierr);
00681   ierr = RegisterModule
00682     ("simple","symmetry-fsnorm",ANALYSISDOUBLE,&SymmetryFSNorm); CHKERRQ(ierr);
00683   ierr = RegisterModule
00684     ("simple","symmetry-fanorm",ANALYSISDOUBLE,&SymmetryFANorm); CHKERRQ(ierr);
00685 
00686   ierr = RegisterModule
00687     ("structure","n-struct-unsymm",ANALYSISINTEGER,&NUnstruct); CHKERRQ(ierr);
00688 
00689   PetscFunctionReturn(0);
00690 }
00691 
00692 #undef __FUNCT__
00693 #define __FUNCT__ "DeRegisterSimpleModules"
00694 PetscErrorCode DeRegisterSimpleModules(void)
00695 {
00696   PetscFunctionBegin;
00697   PetscFunctionReturn(0);
00698 }