00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include <stdlib.h>
00030 #include "anamod.h"
00031 #include "anamodsalsamodules.h"
00032 #include "petscerror.h"
00033 #include "petscmat.h"
00034
00035
00036
00037
00038
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 PetscTruth 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
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
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 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
00095
00096
00097 MPI_Allreduce(cmax,cmaxx,N,MPIU_SCALAR,MPI_MAX,comm);
00098 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
00119
00120 static PetscErrorCode RowVariability
00121 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00122 {
00123 Mat A = (Mat)prob;
00124 PetscReal v=0.; PetscTruth 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
00144
00145 static PetscErrorCode ColVariability
00146 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00147 {
00148 Mat A = (Mat)prob;
00149 PetscReal v=0.; PetscTruth 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
00168
00169 #undef __FUNCT__
00170 #define __FUNCT__ "ComputeDiagonal"
00171
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 PetscTruth 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 MPI_Allreduce((void*)&nn,(void*)&gnn,1,MPI_INT,MPI_SUM,comm);
00195 MPI_Allreduce((void*)&nz,(void*)&gnz,1,MPI_INT,MPI_SUM,comm);
00196 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 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 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
00235
00236 static PetscErrorCode DiagonalAverage
00237 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00238 {
00239 Mat A = (Mat)prob;
00240 PetscReal v=0.; PetscTruth 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
00261
00262 static PetscErrorCode DiagonalVariance
00263 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00264 {
00265 Mat A = (Mat)prob;
00266 PetscReal v=0.; PetscTruth 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
00287
00288 static PetscErrorCode DiagonalSign
00289 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00290 {
00291 Mat A = (Mat)prob;
00292 int v = 0; PetscTruth 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 }