SALSA Analysis Modules
|
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 }