SALSA Analysis Modules
variance.c
Go to the documentation of this file.
00001 /*! \file variance.c \ingroup categories 
00002   \brief Various estimates of how `wild' a matrix is
00003 
00004   \section variance Measurements of element variance
00005 
00006   This module contains various estimates of how different elements in a
00007   matrix are. These measures are completely heuristic, and do not relate
00008   to any known mathematical theory.
00009 
00010   \section Usage
00011 
00012   Activate this module with
00013 
00014   PetscErrorCode RegisterVarianceModules();
00015 
00016   Compute these elements with
00017 
00018   ComputeQuantity("variance",<element>,A,(AnalysisItem*)&res,&flg);
00019 
00020   Available elements are:
00021   - "row-variability" : \f$\max_i \log_{10} {\max_j|a_{ij}|\over\min_j|a_{ij}|} \f$, computed by RowVariability()
00022   - "col-variability" : \f$\max_j \log_{10} {\max_i|a_{ij}|\over\min_i|a_{ij}|} \f$, computed by ColVariability()
00023   - "diagonal-average" : average value of absolute diagonal elements, by DiagonalAverage()
00024   - "diagonal-variance" : standard deviation of diagonal average, by DiagonalVariance()
00025   - "diagonal-sign" : indicator of diagonal sign pattern, by DiagonalSign():
00026     -10 all zero, -2 all negative, -1 nonpositive, 0 indefinite, 1 nonnegative,
00027     2 all positive
00028 */
00029 #include <stdlib.h>
00030 #include "anamod.h"
00031 #include "anamodsalsamodules.h"
00032 #include "petscerror.h"
00033 #include "petscmat.h"
00034 
00035 /*!
00036   Variability in rows and columns.
00037 
00038   This is a computational routine.
00039  */
00040 #undef __FUNCT__
00041 #define __FUNCT__ "ComputeVariability"
00042 static PetscErrorCode ComputeVariability(Mat A)
00043 {
00044   MPI_Comm comm;
00045   PetscReal rr,cr;
00046   PetscReal *rmax,*rmin,*cmax,*cmin,*cmaxx,*cminn,rr_local;
00047   int m,n,M,N,first,last,i,id,row; 
00048   PetscBool has; PetscErrorCode ierr;
00049 
00050   PetscFunctionBegin;
00051   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00052   ierr = MatGetSize(A,&M,&N); CHKERRQ(ierr);
00053   ierr = MatGetLocalSize(A,&m,&n); CHKERRQ(ierr);
00054   ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00055   ierr = PetscMalloc(m*sizeof(PetscReal),&rmax); CHKERRQ(ierr);
00056   ierr = PetscMalloc(m*sizeof(PetscReal),&rmin); CHKERRQ(ierr);
00057   ierr = PetscMalloc(N*sizeof(PetscReal),&cmax); CHKERRQ(ierr);
00058   ierr = PetscMalloc(N*sizeof(PetscReal),&cmin); CHKERRQ(ierr);
00059   ierr = PetscMalloc(N*sizeof(PetscReal),&cmaxx); CHKERRQ(ierr);
00060   ierr = PetscMalloc(N*sizeof(PetscReal),&cminn); CHKERRQ(ierr);
00061 
00062   /* init; we take logs, so 1.e+5 is beyond infinity */
00063   for (i=0; i<m; i++) {rmax[i] = 0; rmin[i] = 1.e+5;}
00064   for (i=0; i<N; i++) {cmax[i] = 0; cmin[i] = 1.e+5;}
00065   for (row=first; row<last; row++) {
00066     int irow=row-first,ncols,icol; const int *cols; const PetscScalar *vals;
00067     ierr = MatGetRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00068     for (icol=0; icol<ncols; icol++) {
00069       int col=cols[icol];
00070       if (vals[icol]) {
00071         PetscReal logval = (PetscReal) log10(PetscAbsScalar(vals[icol]));
00072         if (logval>rmax[irow]) rmax[irow] = logval;
00073         if (logval<rmin[irow]) rmin[irow] = logval;
00074         if (logval>cmax[col]) cmax[col] = logval;
00075         if (logval<cmin[col]) cmin[col] = logval;
00076       }
00077     }
00078     ierr = MatRestoreRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00079   }
00080 
00081   /*
00082    * now get global row spread
00083    */
00084   for (i=0; i<m; i++) rmax[i] -= rmin[i];
00085   for (i=1; i<m; i++) if (rmax[i]>rmax[0]) rmax[0] = rmax[i];
00086   rr_local = rmax[0];
00087   ierr = MPI_Allreduce(&rr_local,&rr,1,MPIU_SCALAR,MPI_MAX,comm);
00088   ierr = GetDataID("variance","row-variability",&id,&has); CHKERRQ(ierr);
00089   HASTOEXIST(has);
00090   ierr = PetscObjectComposedDataSetReal
00091     ((PetscObject)A,id,rr); CHKERRQ(ierr);
00092 
00093   /*
00094    * global column spread: reduce the full length row of values,
00095    * then local max is global
00096    */
00097   ierr = MPI_Allreduce(cmax,cmaxx,N,MPIU_SCALAR,MPI_MAX,comm);
00098   ierr = MPI_Allreduce(cmin,cminn,N,MPIU_SCALAR,MPI_MIN,comm);
00099   for (i=0; i<N; i++) cmaxx[i] -= cminn[i];
00100   for (i=1; i<N; i++) if (cmaxx[i]>cmaxx[0]) cmaxx[0] = cmaxx[i];
00101   cr = cmaxx[0];
00102   ierr = GetDataID("variance","col-variability",&id,&has); CHKERRQ(ierr);
00103   HASTOEXIST(has);
00104   ierr = PetscObjectComposedDataSetReal
00105     ((PetscObject)A,id,cr); CHKERRQ(ierr);
00106 
00107   ierr = PetscFree(rmax); CHKERRQ(ierr);
00108   ierr = PetscFree(cmax); CHKERRQ(ierr);
00109   ierr = PetscFree(cmaxx); CHKERRQ(ierr);
00110   ierr = PetscFree(rmin); CHKERRQ(ierr);
00111   ierr = PetscFree(cmin); CHKERRQ(ierr);
00112   ierr = PetscFree(cminn); CHKERRQ(ierr);
00113   PetscFunctionReturn(0);
00114 }
00115 
00116 #undef __FUNCT__
00117 #define __FUNCT__ "RowVariability"
00118 /*! This routine relies on a call to ComputeVariability()
00119  */
00120 static PetscErrorCode RowVariability
00121 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00122 {
00123   Mat A = (Mat)prob;
00124   PetscReal v=0.; PetscBool has; int id; PetscErrorCode ierr;
00125   PetscFunctionBegin;
00126   ierr = GetDataID("variance","row-variability",&id,&has); CHKERRQ(ierr);
00127   HASTOEXIST(has);
00128   ierr = PetscObjectComposedDataGetReal
00129     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00130   if (*flg) rv->r = v;
00131   else {
00132     ierr = ComputeVariability(A); CHKERRQ(ierr);
00133     ierr = PetscObjectComposedDataGetReal
00134       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00135     if (*flg) rv->r = v;
00136   }
00137   
00138   PetscFunctionReturn(0);
00139 }
00140 
00141 #undef __FUNCT__
00142 #define __FUNCT__ "ColVariability"
00143 /*! This routine relies on a call to ComputeVariability()
00144  */
00145 static PetscErrorCode ColVariability
00146 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00147 {
00148   Mat A = (Mat)prob;
00149   PetscReal v=0.; PetscBool has; int id; PetscErrorCode ierr;
00150   PetscFunctionBegin;
00151   ierr = GetDataID("variance","col-variability",&id,&has); CHKERRQ(ierr);
00152   HASTOEXIST(has);
00153   ierr = PetscObjectComposedDataGetReal
00154     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00155   if (*flg) rv->r = v;
00156   else {
00157     ierr = ComputeVariability(A); CHKERRQ(ierr);
00158     ierr = PetscObjectComposedDataGetReal
00159       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00160     if (*flg) rv->r = v;
00161   }
00162   
00163   PetscFunctionReturn(0);
00164 }
00165 
00166 /*
00167  * Analysis of the matrix diagonal
00168  */
00169 #undef __FUNCT__
00170 #define __FUNCT__ "ComputeDiagonal"
00171 /*! This is a computational routine; it computes the value of several modules
00172  */
00173 static PetscErrorCode ComputeDiagonal(Mat A)
00174 {
00175   int sign;
00176   Vec V; PetscScalar *v; MPI_Comm comm;
00177   PetscReal average,g_average,variance,g_variance;
00178   int local,nn,nz,np,gnn,gnz,gnp,row,id, N;
00179   PetscBool has; PetscErrorCode ierr;
00180 
00181   PetscFunctionBegin;
00182   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00183   ierr = MatGetLocalSize(A,&local,PETSC_NULL); CHKERRQ(ierr);
00184   ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr);
00185   ierr = VecCreateMPI(comm,local,PETSC_DECIDE,&V); CHKERRQ(ierr);
00186   ierr = MatGetDiagonal(A,V); CHKERRQ(ierr);
00187   ierr = VecGetArray(V,&v); CHKERRQ(ierr);
00188 
00189   average = 0; nn = np = nz = 0;
00190   for (row=0; row<local; row++) {
00191     if (v[row]>0) np++; else if (v[row]<0) nn++; else nz++;
00192     average += PetscAbsScalar(v[row]);
00193   }
00194   ierr = MPI_Allreduce((void*)&nn,(void*)&gnn,1,MPI_INT,MPI_SUM,comm);
00195   ierr = MPI_Allreduce((void*)&nz,(void*)&gnz,1,MPI_INT,MPI_SUM,comm);
00196   ierr = MPI_Allreduce((void*)&np,(void*)&gnp,1,MPI_INT,MPI_SUM,comm);
00197   if (np==N) sign = DIAGONAL_POSITIVE;
00198   else if (nn==0) {
00199     if (np==0) sign = DIAGONAL_ZERO; else sign = DIAGONAL_NONNEGATIVE;}
00200   else if (nn==N) sign = DIAGONAL_NEGATIVE;
00201   else if (np==0) sign = DIAGONAL_NONPOSITIVE;
00202   else sign = DIAGONAL_INDEFINITE;
00203   ierr = MPI_Allreduce((void*)&average,(void*)&g_average,1,MPIU_REAL,MPI_SUM,comm);
00204   average = g_average / N; variance = 0;
00205   for (row=0; row<local; row++) {
00206     PetscReal d = PetscAbsScalar(v[row])-average;
00207     variance += d*d;
00208   }
00209   ierr = MPI_Allreduce((void*)&variance,(void*)&g_variance,1,MPIU_REAL,MPI_SUM,comm);
00210   variance = sqrt(g_variance / N);
00211   ierr = GetDataID("variance","diagonal-average",&id,&has); CHKERRQ(ierr);
00212   HASTOEXIST(has);
00213   ierr = PetscObjectComposedDataSetReal
00214     ((PetscObject)A,id,average); CHKERRQ(ierr);
00215 
00216   ierr = GetDataID("variance","diagonal-variance",&id,&has); CHKERRQ(ierr);
00217   HASTOEXIST(has);
00218   ierr = PetscObjectComposedDataSetReal
00219     ((PetscObject)A,id,variance); CHKERRQ(ierr);
00220 
00221   ierr = GetDataID("variance","diagonal-sign",&id,&has); CHKERRQ(ierr);
00222   HASTOEXIST(has);
00223   ierr = PetscObjectComposedDataSetInt
00224     ((PetscObject)A,id,sign); CHKERRQ(ierr);
00225 
00226   ierr = VecRestoreArray(V,&v); CHKERRQ(ierr);
00227   ierr = VecDestroy(&V); CHKERRQ(ierr);
00228   PetscFunctionReturn(0);
00229 }
00230 
00231 #undef __FUNCT__
00232 #define __FUNCT__ "DiagonalAverage"
00233 /*!
00234   This routine relies on a call to ComputeDiagonal()
00235 */
00236 static PetscErrorCode DiagonalAverage
00237 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00238 {
00239   Mat A = (Mat)prob;
00240   PetscReal v=0.; PetscBool has; int id; PetscErrorCode ierr;
00241   PetscFunctionBegin;
00242   ierr = GetDataID("variance","diagonal-average",&id,&has); CHKERRQ(ierr);
00243   HASTOEXIST(has);
00244   ierr = PetscObjectComposedDataGetReal
00245     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00246   if (*flg) rv->r = v;
00247   else {
00248     ierr = ComputeDiagonal(A); CHKERRQ(ierr);
00249     ierr = PetscObjectComposedDataGetReal
00250       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00251     if (*flg) rv->r = v;
00252   }
00253   
00254   PetscFunctionReturn(0);
00255 }
00256 
00257 #undef __FUNCT__
00258 #define __FUNCT__ "DiagonalVariance"
00259 /*!
00260   This routine relies on a call to ComputeDiagonal()
00261 */
00262 static PetscErrorCode DiagonalVariance
00263 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00264 {
00265   Mat A = (Mat)prob;
00266   PetscReal v=0.; PetscBool has; int id; PetscErrorCode ierr;
00267   PetscFunctionBegin;
00268   ierr = GetDataID("variance","diagonal-variance",&id,&has); CHKERRQ(ierr);
00269   HASTOEXIST(has);
00270   ierr = PetscObjectComposedDataGetReal
00271     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00272   if (*flg) rv->r = v;
00273   else {
00274     ierr = ComputeDiagonal(A); CHKERRQ(ierr);
00275     ierr = PetscObjectComposedDataGetReal
00276       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00277     if (*flg) rv->r = v;
00278   }
00279   
00280   PetscFunctionReturn(0);
00281 }
00282 
00283 #undef __FUNCT__
00284 #define __FUNCT__ "DiagonalSign"
00285 /*!
00286   This routine relies on a call to ComputeDiagonal()
00287 */
00288 static PetscErrorCode DiagonalSign
00289 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg)
00290 {
00291   Mat A = (Mat)prob;
00292   int v = 0; PetscBool has; int id; PetscErrorCode ierr;
00293   PetscFunctionBegin;
00294   ierr = GetDataID("variance","diagonal-sign",&id,&has); CHKERRQ(ierr);
00295   HASTOEXIST(has);
00296   ierr = PetscObjectComposedDataGetInt
00297     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00298   if (*flg) rv->i = v;
00299   else {
00300     ierr = ComputeDiagonal(A); CHKERRQ(ierr);
00301     ierr = PetscObjectComposedDataGetInt
00302       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00303     if (*flg) rv->i = v;
00304   }
00305   
00306   PetscFunctionReturn(0);
00307 }
00308 
00309 #undef __FUNCT__
00310 #define __FUNCT__ "RegisterVarianceModules"
00311 PetscErrorCode RegisterVarianceModules()
00312 {
00313   PetscErrorCode ierr;
00314   PetscFunctionBegin;
00315 
00316   ierr = RegisterModule
00317     ("variance","row-variability",ANALYSISDOUBLE,&RowVariability); CHKERRQ(ierr);
00318   ierr = RegisterModule
00319     ("variance","col-variability",ANALYSISDOUBLE,&ColVariability); CHKERRQ(ierr);
00320 
00321   ierr = RegisterModule
00322     ("variance","diagonal-average",ANALYSISDOUBLE,&DiagonalAverage); CHKERRQ(ierr);
00323   ierr = RegisterModule
00324     ("variance","diagonal-variance",ANALYSISDOUBLE,&DiagonalVariance); CHKERRQ(ierr);
00325   ierr = RegisterModule
00326     ("variance","diagonal-sign",ANALYSISINTEGER,&DiagonalSign); CHKERRQ(ierr);
00327 
00328   PetscFunctionReturn(0);
00329 }
00330 
00331 #undef __FUNCT__
00332 #define __FUNCT__ "DeRegisterVarianceModules"
00333 PetscErrorCode DeRegisterVarianceModules(void)
00334 {
00335   PetscFunctionBegin;
00336   PetscFunctionReturn(0);
00337 }