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   PetscTruth 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,PetscTruth *flg)
00107 {
00108   Mat A = (Mat)prob;
00109   PetscReal v=0; PetscTruth 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(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,PetscTruth *flg)
00136 {
00137   Mat A = (Mat)prob;
00138   PetscReal v=0; PetscTruth 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,PetscTruth *flg)
00160 {
00161   Mat A = (Mat)prob;
00162   PetscReal v=0; PetscTruth 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,PetscTruth *flg)
00187 {
00188   Mat A = (Mat)prob;
00189   PetscReal v=0; PetscTruth 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,PetscTruth *flg)
00214 {
00215   Mat A = (Mat)prob;
00216   PetscReal v=0; PetscTruth 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,PetscTruth *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   MPI_Allreduce((void*)&r,(void*)&gr,1,MPIU_REAL,MPI_MIN,comm);
00268   {
00269     int id; PetscTruth 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 static PetscErrorCode DiagonalDominance
00280 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00281 {
00282   Mat A = (Mat)prob;
00283   PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr;
00284   PetscFunctionBegin;
00285 
00286   ierr = GetDataID("simple","diagonal-dominance",&id,&has); CHKERRQ(ierr);
00287   HASTOEXIST(has);
00288   ierr = PetscObjectComposedDataGetReal
00289     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00290   if (*flg) rv->r = v;
00291   else {
00292     ierr = compute_dd(A,flg); CHKERRQ(ierr);
00293     if (*flg) {
00294       ierr = PetscObjectComposedDataGetReal
00295         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00296       if (*flg) rv->r = v;
00297     }
00298   }
00299   
00300   PetscFunctionReturn(0);
00301 }
00302 
00303 /*
00304  * Norms of symmetric / unsymmetric part
00305  */
00306 #undef __FUNCT__
00307 #define __FUNCT__ "MatSymmPartNormInf_seq"
00308 /*!
00309   This is where the real work of MatSymmPartNormInf() is done
00310 
00311   A few temporary arrays of the 
00312   size of a vector are allocated and freed. Otherwise this is just a whole
00313   bunch of pointer chasing.
00314 */
00315 static PetscErrorCode MatSymmPartNormInf_seq
00316 (Mat A,PetscReal *sn_r,PetscReal *an_r,PetscReal *srn_r,PetscReal *arn_r,
00317  int *nun_r,PetscTruth *done)
00318 {
00319   Mat_SeqAIJ *aij;  
00320   int *adx,*bdx,*aii,*bii,*aptr,*bptr,ma,na,i,nun; PetscScalar *va,*vb;
00321   PetscReal snorm,anorm,*sronorm,*aronorm;
00322   PetscErrorCode ierr;
00323   PetscFunctionBegin;
00324 
00325   aij = (Mat_SeqAIJ *) A->data;
00326   ierr = MatGetLocalSize(A,&ma,&na); CHKERRQ(ierr);
00327 
00328   /* two arrays of pointers into the matrix */
00329   aii = bii = aij->i;
00330   adx = bdx = aij->j;
00331   va = vb = aij->a;
00332   ierr = PetscMalloc(ma*sizeof(int),&aptr); CHKERRQ(ierr);
00333   ierr = PetscMalloc(ma*sizeof(int),&bptr); CHKERRQ(ierr);
00334   for (i=0; i<ma; i++) aptr[i] = bptr[i] = aii[i];
00335 
00336   /* temp storage for the inf norm */
00337   ierr = PetscMalloc(ma*sizeof(PetscReal),&sronorm); CHKERRQ(ierr);
00338   ierr = PetscMalloc(ma*sizeof(PetscReal),&aronorm); CHKERRQ(ierr);
00339   ierr = PetscMemzero(sronorm,ma*sizeof(PetscReal)); CHKERRQ(ierr);
00340   ierr = PetscMemzero(aronorm,ma*sizeof(PetscReal)); CHKERRQ(ierr);
00341 
00342   /* update row & frobenius norms with a(i,j)=v and a(j,i)=w */
00343 #if defined(ANAMODDEBUG)
00344 #define RC_NORM(i,j,v,w) \
00345   printf("s/a norm calculation: (i,j)=(%d,%d), (vij,vji)=(%e,%e)\n",i,j,v,w); \
00346   {                                                                      \
00347     PetscReal sum = PetscAbsScalar(v+w)/2, dif = PetscAbsScalar(v-w)/2;  \
00348     ASSERT3(!(i==j && v!=w),"same elements on diagonal",i,j);            \
00349     sronorm[i] += sum; if (i!=j) sronorm[j] += sum;                      \
00350     aronorm[i] += dif; if (i!=j) aronorm[j] += dif;                      \
00351     if (i==j) snorm += sum*sum;                                          \
00352     else {snorm += 2*sum*sum; anorm += 2*dif*dif;}                       \
00353   }
00354 #else
00355 #define RC_NORM(i,j,v,w) \
00356   {                                                                     \
00357     PetscReal sum = PetscAbsScalar(v+w)/2, dif = PetscAbsScalar(v-w)/2; \
00358     ASSERT3(!(i==j && v!=w),"same elements on diagonal",i,j);            \
00359     sronorm[i] += sum; if (i!=j) sronorm[j] += sum;                     \
00360     aronorm[i] += dif; if (i!=j) aronorm[j] += dif;                     \
00361     if (i==j) snorm += sum*sum;                                         \
00362     else {snorm += 2*sum*sum; anorm += 2*dif*dif;}                      \
00363   }
00364 #endif
00365 #define ASSERT3(v,s,i,j) \
00366     if (!(v)) SETERRQ3(1,"Violation of <%s> @ (%d,%d): %d",s,i,j);
00367 #define ASSERT(v,s,i,j,t) \
00368     if (!(v)) SETERRQ4(1,"Violation of <%s> @ (%d,%d): %d",s,i,j,t);
00369 
00370   snorm = anorm = 0; nun = 0;
00371   for (i=0; i<ma; i++) {
00372     int idc,idr; PetscScalar vc,vr;
00373     /* get the a-pointer into the upper triangle */
00374     while (aptr[i]<aii[i+1] && adx[aptr[i]]<i) aptr[i]++;
00375     ASSERT3(aptr[i]==aii[i+1]||adx[aptr[i]]>=i,"A in upper0",i,adx[aptr[i]]);
00376     /* incorporate L elements that are not structurally symmetric */
00377     while (bptr[i]<bii[i+1] && bdx[bptr[i]]<i) {
00378       vr = PetscAbsScalar(vb[bptr[i]]); nun++;
00379       RC_NORM(i,bdx[bptr[i]],vr,0.);
00380       bptr[i]++;
00381     }
00382     ASSERT3(bptr[i]==bii[i+1] || bdx[bptr[i]]>=i,"B part done",i,bdx[bptr[i]]);
00383     while (aptr[i]<aii[i+1]) {
00384       /* in row i we are at element j = idc @ adx( aptr[i] );
00385          vc = A( i, j ) */
00386       idc = adx[aptr[i]]; vc = va[aptr[i]];
00387       ASSERT(idc>=i,"A in upper",i,idc,i-idc);
00388       /* in row j=idc look at elements i' = idr @ bdx(bptr[idc]); vr = B(j,i')
00389          first catch up with structurally unsymmetric elements */
00390       idr = bdx[bptr[idc]]; vr = vb[bptr[idc]];
00391       while (idr<i) {
00392         RC_NORM(idc,idr,vr,0.);
00393         bptr[idc]++; idr = bdx[bptr[idc]]; vr = vb[bptr[idc]]; nun++;
00394       }     
00395       /* now see if (i,j) == (j,i') */
00396       if (idr==i) {
00397         /* structurally symmetric; process these elements of L and U */
00398         RC_NORM(i,idc,vc,vr);
00399         bptr[idc]++;
00400       } else {
00401         /* element (j,i') in L corresponds to a later row in U,
00402            so we have a structurally unsymmetric element in U */
00403         nun++; RC_NORM(i,idc,vc,0.);
00404       }
00405       aptr[i]++;
00406     }
00407     ASSERT(aptr[i]==aii[i+1],"Aptr has traversed row i",i+1,adx[aptr[i]],aii[i+1]);
00408   }
00409 
00410   ierr = PetscFree(aptr); CHKERRQ(ierr);
00411   ierr = PetscFree(bptr); CHKERRQ(ierr);
00412 
00413   snorm = sqrt(snorm); anorm = sqrt(anorm);
00414 
00415   for (i=1; i<ma; i++) {
00416     if (sronorm[i]>sronorm[0]) sronorm[0] = sronorm[i];
00417     if (aronorm[i]>aronorm[0]) aronorm[0] = aronorm[i];
00418   }
00419 
00420   *sn_r = snorm; *an_r = anorm; *srn_r = sronorm[0]; *arn_r = aronorm[0];
00421   *nun_r = nun;
00422   *done = PETSC_TRUE;
00423 
00424   ierr = PetscFree(sronorm); CHKERRQ(ierr);
00425   ierr = PetscFree(aronorm); CHKERRQ(ierr);
00426 
00427   PetscFunctionReturn(0);
00428 }
00429 
00430 #undef __FUNCT__
00431 #define __FUNCT__ "MatSymmPartNormInf"
00432 /*!
00433   This is the computational routine for all quantities relating to
00434   symmetric and anti-symmetric parts of the matrix.
00435 
00436   This computation can only be done in single-processor mode.
00437   If the commandline
00438   option (see \ref options) \c "-anamod-force sequential" is specified,
00439   the matrix is gathered on processor zero to compute this quantity.
00440 
00441   See MatSymmPartNormInf_seq() for the actual computation.
00442 */
00443 static PetscErrorCode MatSymmPartNormInf(Mat A,PetscTruth *done)
00444 {
00445   MPI_Comm comm; Mat Awork; int id,nun;
00446   PetscReal snorm,anorm,sronorm,aronorm;
00447   PetscTruth has,alloc_work,do_work,any_work; PetscErrorCode ierr;
00448 
00449   PetscFunctionBegin;
00450 
00451   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00452   ierr = AnaModGetSequentialMatrix
00453     (A,&Awork,&alloc_work,&do_work,&any_work); CHKERRQ(ierr);
00454   if (!any_work) {
00455     *done = PETSC_FALSE; PetscFunctionReturn(0);
00456   }
00457   if (do_work) {
00458     ierr = MatSymmPartNormInf_seq
00459       (Awork,&snorm,&anorm,&sronorm,&aronorm,&nun,done); CHKERRQ(ierr);
00460   }
00461   if (alloc_work) {
00462     ierr = MatDestroy(Awork); CHKERRQ(ierr);
00463   }
00464 
00465   MPI_Bcast((void*)done,1,MPI_INT,0,comm);
00466   MPI_Bcast((void*)&snorm,1,MPIU_REAL,0,comm);
00467   MPI_Bcast((void*)&anorm,1,MPIU_REAL,0,comm);
00468   MPI_Bcast((void*)&sronorm,1,MPIU_REAL,0,comm);
00469   MPI_Bcast((void*)&aronorm,1,MPIU_REAL,0,comm);
00470   MPI_Bcast((void*)&nun,1,MPI_INT,0,comm);
00471 
00472   /* infinity norms of symmetric / anti-symmetric part */
00473   ierr = GetDataID("simple","symmetry-snorm",&id,&has); CHKERRQ(ierr);
00474   HASTOEXIST(has);
00475   ierr = PetscObjectComposedDataSetReal
00476     ((PetscObject)A,id,sronorm); CHKERRQ(ierr);
00477 
00478   ierr = GetDataID("simple","symmetry-anorm",&id,&has); CHKERRQ(ierr);
00479   HASTOEXIST(has);
00480   ierr = PetscObjectComposedDataSetReal
00481     ((PetscObject)A,id,aronorm); CHKERRQ(ierr);
00482 
00483   /* frobenius norms of symmetric / anti-symmetric part */
00484   ierr = GetDataID("simple","symmetry-fsnorm",&id,&has); CHKERRQ(ierr);
00485   HASTOEXIST(has);
00486   ierr = PetscObjectComposedDataSetReal
00487     ((PetscObject)A,id,snorm); CHKERRQ(ierr);
00488 
00489   ierr = GetDataID("simple","symmetry-fanorm",&id,&has); CHKERRQ(ierr);
00490   HASTOEXIST(has);
00491   ierr = PetscObjectComposedDataSetReal
00492     ((PetscObject)A,id,anorm); CHKERRQ(ierr);
00493 
00494   ierr = GetDataID("structure","n-struct-unsymm",&id,&has); CHKERRQ(ierr);
00495   HASTOEXIST(has);
00496   ierr = PetscObjectComposedDataSetInt
00497     ((PetscObject)A,id,nun); CHKERRQ(ierr);
00498 
00499   PetscFunctionReturn(0);
00500 }
00501 
00502 #undef __FUNCT__
00503 #define __FUNCT__ "SymmetrySNorm"
00504 /*! Compute and store the infinity norm of the symmetric part of the matrix.
00505   The actual computation is done in MatSymmPartNormInf();
00506   see that routine for remarks on parallelism.
00507 */
00508 static PetscErrorCode SymmetrySNorm
00509 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00510 {
00511   Mat A = (Mat)prob;
00512   PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr;
00513   PetscFunctionBegin;
00514   ierr = GetDataID("simple","symmetry-snorm",&id,&has); CHKERRQ(ierr);
00515   HASTOEXIST(has);
00516   ierr = PetscObjectComposedDataGetReal
00517     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00518   if (*flg) rv->r = v;
00519   else {
00520     ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr);
00521     if (*flg) {
00522       ierr = PetscObjectComposedDataGetReal
00523         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00524       if (*flg) rv->r = v;
00525     }
00526   }
00527   
00528   PetscFunctionReturn(0);
00529 }
00530 
00531 #undef __FUNCT__
00532 #define __FUNCT__ "SymmetryANorm"
00533 /*! Compute and store the infinity norm of the antisymmetric part of the matrix.
00534   The actual computation is done in MatSymmPartNormInf();
00535   see that routine for remarks on parallelism.
00536 */
00537 static PetscErrorCode SymmetryANorm
00538 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00539 {
00540     Mat A = (Mat)prob;
00541   PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr;
00542   PetscFunctionBegin;
00543   ierr = GetDataID("simple","symmetry-anorm",&id,&has); CHKERRQ(ierr);
00544   HASTOEXIST(has);
00545   ierr = PetscObjectComposedDataGetReal
00546     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00547   if (*flg) rv->r = v;
00548   else {
00549     ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr);
00550     if (*flg) {
00551       ierr = PetscObjectComposedDataGetReal
00552         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00553       if (*flg) rv->r = v;
00554     }
00555   }
00556   
00557   PetscFunctionReturn(0);
00558 }
00559 
00560 #undef __FUNCT__
00561 #define __FUNCT__ "SymmetryFSNorm"
00562 /*! Compute and store the Frobenius norm of the symmetric part of the matrix.
00563   The actual computation is done in MatSymmPartNormInf();
00564   see that routine for remarks on parallelism.
00565 */
00566 static PetscErrorCode SymmetryFSNorm
00567 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00568 {
00569     Mat A = (Mat)prob;
00570   PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr;
00571   PetscFunctionBegin;
00572   ierr = GetDataID("simple","symmetry-fsnorm",&id,&has); CHKERRQ(ierr);
00573   HASTOEXIST(has);
00574   ierr = PetscObjectComposedDataGetReal
00575     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00576   if (*flg) rv->r = v;
00577   else {
00578     ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr);
00579     if (*flg) {
00580       ierr = PetscObjectComposedDataGetReal
00581         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00582       if (*flg) rv->r = v;
00583     }
00584   }
00585   
00586   PetscFunctionReturn(0);
00587 }
00588 
00589 #undef __FUNCT__
00590 #define __FUNCT__ "SymmetryFANorm"
00591 /*! Compute and store the Frobenius norm of the antisymmetric part of the matrix.
00592   The actual computation is done in MatSymmPartNormInf();
00593   see that routine for remarks on parallelism.
00594 */
00595 static PetscErrorCode SymmetryFANorm
00596 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00597 {
00598     Mat A = (Mat)prob;
00599   PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr;
00600   PetscFunctionBegin;
00601   ierr = GetDataID("simple","symmetry-fanorm",&id,&has); CHKERRQ(ierr);
00602   HASTOEXIST(has);
00603   ierr = PetscObjectComposedDataGetReal
00604     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00605   if (*flg) rv->r = v;
00606   else {
00607     ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr);
00608     if (*flg) {
00609       ierr = PetscObjectComposedDataGetReal
00610         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00611       if (*flg) rv->r = v;
00612     }
00613   }
00614   
00615   PetscFunctionReturn(0);
00616 }
00617 
00618 #undef __FUNCT__
00619 #define __FUNCT__ "NUnstruct"
00620 /*! Compute and store the number of structurally unsymmetric elements of a matrix.
00621   The actual computation is done in MatSymmPartNormInf();
00622   see that routine for remarks on parallelism.
00623 
00624   The "n-struct-unsymm" feature is declared to be in category "structure",
00625   rather than simple.
00626 */
00627 static PetscErrorCode NUnstruct
00628 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00629 {
00630   Mat A = (Mat)prob;
00631   PetscReal v=0; PetscTruth has; int id; PetscErrorCode ierr;
00632   PetscFunctionBegin;
00633   ierr = GetDataID("structure","n-struct-unsymm",&id,&has); CHKERRQ(ierr);
00634   HASTOEXIST(has);
00635   ierr = PetscObjectComposedDataGetInt
00636     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00637   if (*flg) rv->i = v;
00638   else {
00639     ierr = MatSymmPartNormInf(A,flg); CHKERRQ(ierr);
00640     if (*flg) {
00641       ierr = PetscObjectComposedDataGetInt
00642         ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00643       if (*flg) rv->i = v;
00644     }
00645   }
00646   
00647   PetscFunctionReturn(0);
00648 }
00649 
00650 #undef __FUNCT__
00651 #define __FUNCT__ "RegisterSimpleModules"
00652 /*!
00653   Declare norm-like modules that can be performed with simple calculations.
00654 */
00655 PetscErrorCode RegisterSimpleModules()
00656 {
00657   PetscErrorCode ierr;
00658   PetscFunctionBegin;
00659 
00660   ierr = RegisterModule
00661     ("simple","trace",ANALYSISDOUBLE,&Trace); CHKERRQ(ierr);
00662   ierr = RegisterModule
00663     ("simple","trace-abs",ANALYSISDOUBLE,&TraceAbs); CHKERRQ(ierr);
00664 
00665   ierr = RegisterModule
00666     ("simple","norm1",ANALYSISDOUBLE,&norm1); CHKERRQ(ierr);
00667   ierr = RegisterModule
00668     ("simple","normInf",ANALYSISDOUBLE,&normInf); CHKERRQ(ierr);
00669   ierr = RegisterModule
00670     ("simple","normF",ANALYSISDOUBLE,&normF); CHKERRQ(ierr);
00671 
00672   ierr = RegisterModule
00673     ("simple","diagonal-dominance",ANALYSISDOUBLE,&DiagonalDominance); CHKERRQ(ierr);
00674 
00675   ierr = RegisterModule
00676     ("simple","symmetry-snorm",ANALYSISDOUBLE,&SymmetrySNorm); CHKERRQ(ierr);
00677   ierr = RegisterModule
00678     ("simple","symmetry-anorm",ANALYSISDOUBLE,&SymmetryANorm); CHKERRQ(ierr);
00679   ierr = RegisterModule
00680     ("simple","symmetry-fsnorm",ANALYSISDOUBLE,&SymmetryFSNorm); CHKERRQ(ierr);
00681   ierr = RegisterModule
00682     ("simple","symmetry-fanorm",ANALYSISDOUBLE,&SymmetryFANorm); CHKERRQ(ierr);
00683 
00684   ierr = RegisterModule
00685     ("structure","n-struct-unsymm",ANALYSISINTEGER,&NUnstruct); CHKERRQ(ierr);
00686 
00687   PetscFunctionReturn(0);
00688 }
00689 
00690 #undef __FUNCT__
00691 #define __FUNCT__ "DeRegisterSimpleModules"
00692 PetscErrorCode DeRegisterSimpleModules(void)
00693 {
00694   PetscFunctionBegin;
00695   PetscFunctionReturn(0);
00696 }