SALSA Analysis Modules
structure.c
Go to the documentation of this file.
00001 /*! \file structure.c \ingroup categories 
00002   \brief Structural properties of a matrix
00003 
00004   This file defines the routines for computing structural properties 
00005   of a matrix, as well as the routine RegisterStructureModules() that 
00006   registers them.
00007 
00008   \section structure Structural properties
00009 
00010   This category computes matrix properties that are only a function of
00011   the nonzero structure. Thus, these properties will likely stay invariant
00012   during a nonlinear solve process, or while time-stepping a system
00013   of equations.
00014 
00015   \section Usage
00016 
00017   Activate this module with
00018 
00019   PetscErrorCode RegisterStructureModules();
00020 
00021   Compute these elements with
00022 
00023   ComputeQuantity("structure",<element>,A,(void*)&res,&flg);
00024 
00025   Available elements are:
00026   - "nrows" : number of rows in the matrix
00027   - "nnzeros" : number of nonzero elements in the matrix
00028   - "max-nnzeros-per-row" : maximum number of nonzeros per row
00029   - "min-nnzeros-per-row" : minimum number of nonzeros per row
00030   - "left-bandwidth" : \f$\max_i\{i-j\colon a_{ij}\not=0\}\f$
00031   - "right-bandwidth" : \f$\max_i\{j-i\colon a_{ij}\not=0\}\f$
00032   - "n-dummy-rows" : number of rows with only one element
00033   - "dummy-rows-kind" : 0 if all dummy rows have a one on the diagonal,
00034     1 if the value is not one, 2 if not on the diagonal
00035   - "diag-zerostart" : \f$\min\{i\colon \forall_{j>i}a_{jj}=0\}\f$
00036   - "diag-definite" : 1 if diagonal positive, 0 otherwise
00037   - "block-size" : integer size of blocks that comprise matrix block structure,
00038     1 in the general case
00039   - "symmetry" : 1 for symmetric matrix, 0 otherwise
00040   - "n-struct-unsymm" : number of structurally unsymmetric elements
00041     (this is actually calculated in the simple.c)
00042 
00043 */
00044 #include <stdlib.h>
00045 #include "anamod.h"
00046 #include "anamodutils.h"
00047 #include "anamodsalsamodules.h"
00048 #include "petscerror.h"
00049 #include "petscmat.h"
00050 
00051 #undef __FUNCT__
00052 #define __FUNCT__ "nRows"
00053 static PetscErrorCode nRows
00054 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00055 {
00056   Mat A = (Mat)prob;
00057   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00058   PetscFunctionBegin;
00059   ierr = GetDataID("structure","nrows",&id,&has); CHKERRQ(ierr);
00060   HASTOEXIST(has);
00061   ierr = PetscObjectComposedDataGetInt
00062     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00063   if (*flg) rv->i = v;
00064   else {
00065     int idum;
00066     ierr = MatGetSize(A,&v,&idum); CHKERRQ(ierr);
00067     *flg = PETSC_TRUE;
00068     ierr = PetscObjectComposedDataSetInt
00069       ((PetscObject)A,id,v); CHKERRQ(ierr);
00070       rv->i = v;
00071   }
00072   PetscFunctionReturn(0);
00073 }
00074 
00075 #undef __FUNCT__
00076 #define __FUNCT__ "Symmetry"
00077 static PetscErrorCode Symmetry
00078 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00079 {
00080   Mat A = (Mat)prob;
00081   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00082   PetscFunctionBegin;
00083   ierr = GetDataID("structure","symmetry",&id,&has); CHKERRQ(ierr);
00084   HASTOEXIST(has);
00085   ierr = PetscObjectComposedDataGetInt
00086     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00087   if (*flg) rv->i = v;
00088   else {
00089     PetscTruth f; const MatType type;
00090     ierr = MatGetType(A,&type); CHKERRQ(ierr);
00091     ierr = PetscStrcmp(type,MATSEQAIJ,flg); CHKERRQ(ierr);
00092     if (*flg) {
00093       ierr = MatIsSymmetric(A,0.0,&f); CHKERRQ(ierr);
00094       v = f;
00095       ierr = PetscObjectComposedDataSetInt
00096         ((PetscObject)A,id,v); CHKERRQ(ierr);
00097       rv->i = v;
00098     }
00099   }
00100   
00101   PetscFunctionReturn(0);
00102 }
00103 
00104 /*
00105  * Nonzero structure, bandwidth
00106  */
00107 #undef __FUNCT__
00108 #define __FUNCT__ "compute_nnz_structure"
00109 static PetscErrorCode compute_nnz_structure(AnaModNumericalProblem prob)
00110 {
00111   Mat A = (Mat)prob;
00112   int nnz,p,q,minnz,maxnz,minnnz,maxnnz;
00113   MPI_Comm comm;
00114   int N,first,last,local_size,row,nz,*Bl,*Br,bl,br, id; 
00115   PetscTruth has; PetscErrorCode ierr;
00116 
00117   PetscFunctionBegin;
00118   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00119   ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr);
00120   ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00121   local_size = last-first;
00122   ierr = PetscMalloc(local_size*sizeof(int),&Bl); CHKERRQ(ierr);
00123   ierr = PetscMalloc(local_size*sizeof(int),&Br); CHKERRQ(ierr);
00124   bl = 0; br = 0; nz = 0; minnz = maxnz = 0;
00125   for (row=first; row<last; row++) {
00126     int icol,ncols,localnz=0; const int *cols; const PetscScalar *val;
00127     ierr = MatGetRow(A,row,&ncols,&cols,&val); CHKERRQ(ierr);
00128     for (icol=0; icol<ncols; icol++) 
00129       if (val[icol]!=0.) localnz++;
00130     nz += localnz;
00131     if (row==first) minnz = localnz;
00132     if (localnz<minnz) minnz = localnz;
00133     if (localnz>maxnz) maxnz = localnz;
00134     if (ncols) {
00135       Bl[row-first] = bl = PetscMax(bl,row-cols[0]);
00136       Br[row-first] = br = PetscMax(br,cols[ncols-1]-row);
00137     }
00138     ierr = MatRestoreRow(A,row,&ncols,&cols,&val); CHKERRQ(ierr);
00139   }
00140 
00141   MPI_Allreduce((void*)&nz,(void*)&nnz,1,MPI_INT,MPI_SUM,comm);
00142   ierr = GetDataID("structure","nnzeros",&id,&has); CHKERRQ(ierr);
00143   HASTOEXIST(has);
00144   ierr = PetscObjectComposedDataSetInt
00145     ((PetscObject)A,id,nnz); CHKERRQ(ierr);
00146 
00147   MPI_Allreduce((void*)&maxnz,(void*)&maxnnz,1,MPI_INT,MPI_MAX,comm);
00148   ierr = GetDataID("structure","max-nnzeros-per-row",&id,&has); CHKERRQ(ierr);
00149   HASTOEXIST(has);
00150   ierr = PetscObjectComposedDataSetInt
00151     ((PetscObject)A,id,maxnnz); CHKERRQ(ierr);
00152 
00153   MPI_Allreduce((void*)&minnz,(void*)&minnnz,1,MPI_INT,MPI_MIN,comm);
00154   ierr = GetDataID("structure","min-nnzeros-per-row",&id,&has); CHKERRQ(ierr);
00155   HASTOEXIST(has);
00156   ierr = PetscObjectComposedDataSetInt
00157     ((PetscObject)A,id,minnnz); CHKERRQ(ierr);
00158 
00159   MPI_Allreduce((void*)&bl,(void*)&p,1,MPI_INT,MPI_MAX,comm);
00160   ierr = GetDataID("structure","left-bandwidth",&id,&has); CHKERRQ(ierr);
00161   HASTOEXIST(has);
00162   ierr = PetscObjectComposedDataSetInt
00163     ((PetscObject)A,id,p); CHKERRQ(ierr);
00164 
00165   MPI_Allreduce((void*)&br,(void*)&q,1,MPI_INT,MPI_MAX,comm);
00166   ierr = GetDataID("structure","right-bandwidth",&id,&has); CHKERRQ(ierr);
00167   HASTOEXIST(has);
00168   ierr = PetscObjectComposedDataSetInt
00169     ((PetscObject)A,id,q); CHKERRQ(ierr);
00170 
00171   ierr = PetscFree(Bl); CHKERRQ(ierr);
00172   ierr = PetscFree(Br); CHKERRQ(ierr);
00173 
00174   PetscFunctionReturn(0);
00175 }
00176 
00177 #undef __FUNCT__
00178 #define __FUNCT__ "NNonZeros"
00179 static PetscErrorCode NNonZeros
00180 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00181 {
00182   Mat A = (Mat)prob;
00183   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00184   PetscFunctionBegin;
00185   ierr = GetDataID("structure","nnzeros",&id,&has); CHKERRQ(ierr);
00186   HASTOEXIST(has);
00187   ierr = PetscObjectComposedDataGetInt
00188     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00189   if (*flg) rv->i = v;
00190   else {
00191     ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00192     ierr = PetscObjectComposedDataGetInt
00193       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00194     if (*flg) rv->i = v;
00195   }
00196   
00197   PetscFunctionReturn(0);
00198 }
00199 
00200 #undef __FUNCT__
00201 #define __FUNCT__ "MaxNNonZerosPerRow"
00202 static PetscErrorCode MaxNNonZerosPerRow
00203 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00204 {
00205   Mat A = (Mat)prob;
00206   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00207   PetscFunctionBegin;
00208   ierr = GetDataID("structure","max-nnzeros-per-row",&id,&has); CHKERRQ(ierr);
00209   HASTOEXIST(has);
00210   ierr = PetscObjectComposedDataGetInt
00211     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00212   if (*flg) rv->i = v;
00213   else {
00214     ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00215     ierr = PetscObjectComposedDataGetInt
00216       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00217     if (*flg) rv->i = v;
00218   }
00219   
00220   PetscFunctionReturn(0);
00221 }
00222 
00223 #undef __FUNCT__
00224 #define __FUNCT__ "MinNNonZerosPerRow"
00225 static PetscErrorCode MinNNonZerosPerRow
00226 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00227 {
00228   Mat A = (Mat)prob;
00229   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00230   PetscFunctionBegin;
00231   ierr = GetDataID("structure","min-nnzeros-per-row",&id,&has); CHKERRQ(ierr);
00232   HASTOEXIST(has);
00233   ierr = PetscObjectComposedDataGetInt
00234     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00235   if (*flg) rv->i = v;
00236   else {
00237     ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00238     ierr = PetscObjectComposedDataGetInt
00239       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00240     if (*flg) rv->i = v;
00241   }
00242   
00243   PetscFunctionReturn(0);
00244 }
00245 
00246 #undef __FUNCT__
00247 #define __FUNCT__ "LBandWidth"
00248 static PetscErrorCode LBandWidth
00249 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00250 {
00251   Mat A = (Mat)prob;
00252   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00253   PetscFunctionBegin;
00254   ierr = GetDataID("structure","left-bandwidth",&id,&has); CHKERRQ(ierr);
00255   HASTOEXIST(has);
00256   ierr = PetscObjectComposedDataGetInt
00257     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00258   if (*flg) rv->i = v;
00259   else {
00260     ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00261     ierr = PetscObjectComposedDataGetInt
00262       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00263     if (*flg) rv->i = v;
00264   }
00265   
00266   PetscFunctionReturn(0);
00267 }
00268 
00269 #undef __FUNCT__
00270 #define __FUNCT__ "RBandWidth"
00271 static PetscErrorCode RBandWidth
00272 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00273 {
00274   Mat A = (Mat)prob;
00275   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00276   PetscFunctionBegin;
00277   ierr = GetDataID("structure","right-bandwidth",&id,&has); CHKERRQ(ierr);
00278   HASTOEXIST(has);
00279   ierr = PetscObjectComposedDataGetInt
00280     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00281   if (*flg) rv->i = v;
00282   else {
00283     ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00284     ierr = PetscObjectComposedDataGetInt
00285       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00286     if (*flg) rv->i = v;
00287   }
00288   
00289   PetscFunctionReturn(0);
00290 }
00291 
00292 #undef __FUNCT__
00293 #define __FUNCT__ "LeftSkyline"
00294 static PetscErrorCode LeftSkyline
00295 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00296 {
00297   Mat A = (Mat)prob; int N;
00298   int *v = NULL; PetscTruth has; int id; PetscErrorCode ierr;
00299   PetscFunctionBegin;
00300   ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr);
00301   ierr = GetDataID("structure","left-skyline",&id,&has); CHKERRQ(ierr);
00302   HASTOEXIST(has);
00303   ierr = PetscObjectComposedDataGetIntstar
00304     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00305   if (*flg) {
00306     rv->ii = v; *lv = N;
00307   } else {
00308     ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00309     ierr = PetscObjectComposedDataGetIntstar
00310       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00311     if (*flg) {
00312       rv->ii = v; *lv = N;
00313     }
00314   }
00315   
00316   PetscFunctionReturn(0);
00317 }
00318 
00319 #undef __FUNCT__
00320 #define __FUNCT__ "RightSkyline"
00321 static PetscErrorCode RightSkyline
00322 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00323 {
00324   Mat A = (Mat)prob; int N;
00325   int *v = NULL; PetscTruth has; int id; PetscErrorCode ierr;
00326   PetscFunctionBegin;
00327   ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr);
00328   ierr = GetDataID("structure","right-skyline",&id,&has); CHKERRQ(ierr);
00329   HASTOEXIST(has);
00330   ierr = PetscObjectComposedDataGetIntstar
00331     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00332   if (*flg) {
00333     rv->ii = v; *lv = N;
00334   } else {
00335     ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00336     ierr = PetscObjectComposedDataGetIntstar
00337       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00338     if (*flg) {
00339       rv->ii = v; *lv = N;
00340     }
00341   }
00342   
00343   PetscFunctionReturn(0);
00344 }
00345 
00346 /*
00347  * Dummy rows
00348  */
00349 #undef __FUNCT__
00350 #define __FUNCT__ "compute_dummy_rows"
00351 static PetscErrorCode compute_dummy_rows(Mat A)
00352 {
00353   int gkind;
00354   MPI_Comm comm;
00355   int first,last,row,ncols,ndummies,Ndummies,*dummies,*Dummies,kind,id; 
00356   PetscTruth has; PetscErrorCode ierr;
00357   const int *cols; const PetscScalar *vals;
00358 
00359   PetscFunctionBegin;
00360   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00361   ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00362   ierr = PetscMalloc((last-first+1)*sizeof(int),&dummies); CHKERRQ(ierr);
00363 
00364   ndummies = 0; kind = 0;
00365   for (row=first; row<last; row++) {
00366     ierr = MatGetRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00367     if (ncols==1) {
00368       dummies[ndummies++] = row;
00369       if (vals[0]!=1. && !kind ) kind = 1;
00370       if (cols[0]!=row) kind = 2;
00371     }
00372     ierr = MatRestoreRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00373   }
00374 
00375   ierr = MPI_Allreduce
00376     ((void*)&kind,(void*)&gkind,1,MPI_INT,MPI_MAX,comm); CHKERRQ(ierr);
00377   ierr = GetDataID("structure","dummy-rows-kind",&id,&has); CHKERRQ(ierr);
00378   HASTOEXIST(has);
00379   ierr = PetscObjectComposedDataSetInt
00380     ((PetscObject)A,id,gkind); CHKERRQ(ierr);
00381 
00382   {
00383     int *sizes,*offsets,ntids,i,s;
00384 
00385     MPI_Comm_size(comm,&ntids);
00386     ierr = PetscMalloc(ntids*sizeof(int),&sizes); CHKERRQ(ierr);
00387     ierr = PetscMalloc(ntids*sizeof(int),&offsets); CHKERRQ(ierr);
00388     MPI_Allgather(&ndummies,1,MPI_INT,sizes,1,MPI_INT,comm);
00389     s=0; for (i=0; i<ntids; i++) {offsets[i]=s; s+=sizes[i];} Ndummies = s;
00390     if (Ndummies==0) 
00391       Dummies = NULL;
00392     else {
00393       ierr = PetscMalloc(Ndummies*sizeof(int),&Dummies); CHKERRQ(ierr);
00394       ierr = MPI_Allgatherv
00395         (dummies,ndummies,MPI_INT,Dummies,sizes,offsets,MPI_INT,
00396          comm); CHKERRQ(ierr);
00397     }
00398     ierr = PetscFree(dummies); CHKERRQ(ierr);
00399     ierr = PetscFree(offsets); CHKERRQ(ierr);
00400     ierr = PetscFree(sizes); CHKERRQ(ierr);
00401   }
00402 
00403   ierr = GetDataID("structure","n-dummy-rows",&id,&has); CHKERRQ(ierr);
00404   HASTOEXIST(has);
00405   ierr = PetscObjectComposedDataSetInt
00406     ((PetscObject)A,id,Ndummies); CHKERRQ(ierr);
00407 
00408   ierr = GetDataID("structure","dummy-rows",&id,&has); CHKERRQ(ierr);
00409   HASTOEXIST(has);
00410   ierr = PetscObjectComposedDataSetIntstar
00411     ((PetscObject)A,id,Dummies); CHKERRQ(ierr);
00412   
00413   PetscFunctionReturn(0);
00414 }
00415 
00416 #undef __FUNCT__
00417 #define __FUNCT__ "NDummyRows"
00418 static PetscErrorCode NDummyRows
00419 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00420 {
00421   Mat A = (Mat)prob;
00422   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00423   PetscFunctionBegin;
00424   ierr = GetDataID("structure","n-dummy-rows",&id,&has); CHKERRQ(ierr);
00425   HASTOEXIST(has);
00426   ierr = PetscObjectComposedDataGetInt
00427     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00428   if (*flg) rv->i = v;
00429   else {
00430     ierr = compute_dummy_rows(A); CHKERRQ(ierr);
00431     ierr = PetscObjectComposedDataGetInt
00432       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00433     if (*flg) rv->i = v;
00434   }
00435   
00436   PetscFunctionReturn(0);
00437 }
00438 
00439 #undef __FUNCT__
00440 #define __FUNCT__ "DummyRowsKind"
00441 static PetscErrorCode DummyRowsKind
00442 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00443 {
00444   Mat A = (Mat)prob;
00445   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00446   PetscFunctionBegin;
00447   ierr = GetDataID("structure","dummy-rows-kind",&id,&has); CHKERRQ(ierr);
00448   HASTOEXIST(has);
00449   ierr = PetscObjectComposedDataGetInt
00450     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00451   if (*flg) rv->i = v;
00452   else {
00453     ierr = compute_dummy_rows(A); CHKERRQ(ierr);
00454     ierr = PetscObjectComposedDataGetInt
00455       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00456     if (*flg) rv->i = v;
00457   }
00458   
00459   PetscFunctionReturn(0);
00460 }
00461 
00462 #undef __FUNCT__
00463 #define __FUNCT__ "DummyRows"
00464 static PetscErrorCode DummyRows
00465 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00466 {
00467   Mat A = (Mat)prob;
00468   int *v = NULL; PetscTruth has; int id,id2; PetscErrorCode ierr;
00469   PetscFunctionBegin;
00470   ierr = GetDataID("structure","dummy-rows",&id,&has); CHKERRQ(ierr);
00471   HASTOEXIST(has);
00472   ierr = PetscObjectComposedDataGetIntstar
00473     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00474   ierr = GetDataID("structure","n-dummy-rows",&id2,&has); CHKERRQ(ierr);
00475   HASTOEXIST(has);
00476   ierr = PetscObjectComposedDataGetInt
00477     ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00478   if (*flg) {
00479     rv->ii = v;
00480   } else {
00481     ierr = compute_dummy_rows(A); CHKERRQ(ierr);
00482     ierr = PetscObjectComposedDataGetIntstar
00483       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00484     ierr = PetscObjectComposedDataGetInt
00485       ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00486     if (*flg) rv->ii = v;
00487   }
00488   
00489   PetscFunctionReturn(0);
00490 }
00491 
00492 /*
00493  * Diagonal analysis
00494  */
00495 #undef __FUNCT__
00496 #define __FUNCT__ "compute_posdiag"
00497 static PetscErrorCode compute_posdiag
00498 (AnaModNumericalProblem prob)
00499 {
00500   Mat A = (Mat)prob;
00501   MPI_Comm comm;
00502   Vec D;
00503   PetscScalar *d;
00504   int i,first,last,m,M,di,df,Di,Df,pd,id; 
00505   PetscTruth has; PetscErrorCode ierr;
00506 
00507   PetscFunctionBegin;
00508   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00509   ierr = MatGetSize(A,&M,PETSC_NULL); CHKERRQ(ierr);
00510   ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00511   m = last-first;
00512   ierr = VecCreateMPI(comm,m,PETSC_DECIDE,&D); CHKERRQ(ierr);
00513   ierr = MatGetDiagonal(A,D); CHKERRQ(ierr);
00514   ierr = VecGetArray(D,&d); CHKERRQ(ierr);
00515   di = 1; df = first;
00516   for (i=0; i<m; i++) {
00517     if (d[i]<=0) di = 0;
00518     if (d[i]!=0) df = i+1;
00519   }
00520   MPI_Allreduce((void*)&di,(void*)&Di,1,MPI_INT,MPI_PROD,comm);
00521   ierr = GetDataID("structure","diag-definite",&id,&has); CHKERRQ(ierr);
00522   HASTOEXIST(has);
00523   ierr = PetscObjectComposedDataSetInt
00524     ((PetscObject)A,id,Di); CHKERRQ(ierr);
00525 
00526   MPI_Allreduce((void*)&df,(void*)&Df,1,MPI_INT,MPI_MAX,comm);
00527   ierr = GetDataID("structure","diag-zerostart",&id,&has); CHKERRQ(ierr);
00528   HASTOEXIST(has);
00529   ierr = PetscObjectComposedDataSetInt
00530     ((PetscObject)A,id,Df); CHKERRQ(ierr);
00531 
00532 
00533   if (Df<M) pd = 2+Df; else pd = Di;
00534 
00535   ierr = VecRestoreArray(D,&d); CHKERRQ(ierr);
00536   ierr = VecDestroy(D); CHKERRQ(ierr);
00537   PetscFunctionReturn(0);
00538 }
00539 
00540 #undef __FUNCT__
00541 #define __FUNCT__ "DiagZeroStart"
00542 static PetscErrorCode DiagZeroStart
00543 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00544 {
00545   Mat A = (Mat)prob;
00546   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00547   PetscFunctionBegin;
00548   ierr = GetDataID("structure","diag-zerostart",&id,&has); CHKERRQ(ierr);
00549   HASTOEXIST(has);
00550   ierr = PetscObjectComposedDataGetInt
00551     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00552   if (*flg) rv->i = v;
00553   else {
00554     ierr = compute_posdiag(A); CHKERRQ(ierr);
00555     ierr = PetscObjectComposedDataGetInt
00556       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00557     if (*flg) rv->i = v;
00558   }
00559   
00560   PetscFunctionReturn(0);
00561 }
00562 
00563 #undef __FUNCT__
00564 #define __FUNCT__ "DiagDefinite"
00565 static PetscErrorCode DiagDefinite
00566 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00567 {
00568   Mat A = (Mat)prob;
00569   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00570   PetscFunctionBegin;
00571   ierr = GetDataID("structure","diag-definite",&id,&has); CHKERRQ(ierr);
00572   HASTOEXIST(has);
00573   ierr = PetscObjectComposedDataGetInt
00574     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00575   if (*flg) rv->i = v;
00576   else {
00577     ierr = compute_posdiag(A); CHKERRQ(ierr);
00578     ierr = PetscObjectComposedDataGetInt
00579       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00580     if (*flg) rv->i = v;
00581   }
00582   
00583   PetscFunctionReturn(0);
00584 }
00585 
00586 /*
00587  * Simple block analysis
00588  */
00589 #undef __FUNCT__
00590 #define __FUNCT__ "regularblocks"
00591 static PetscErrorCode regularblocks
00592 (AnaModNumericalProblem prob)
00593 {
00594   Mat A = (Mat)prob;
00595   int first,last,local_size,row, bs,gbs,id,mytid;
00596   MPI_Comm comm; PetscTruth has; PetscErrorCode ierr;
00597   PetscFunctionBegin;
00598   ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00599   ierr = MPI_Comm_rank(comm,&mytid); CHKERRQ(ierr);
00600   ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00601   local_size = last-first;
00602   if (mytid==0) {
00603     int row;
00604     bs = 1;
00605     for (row=first; row<last; row++) {
00606       const int *cols; const PetscScalar *vals; int ncols,icol=0;
00607       PetscTruth found = PETSC_FALSE;
00608       ierr = MatGetRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00609       while (icol<ncols && cols[icol]<row) icol++;
00610       if (cols[icol]>row+1 ||
00611           (cols[icol]==row && icol+1<ncols && 
00612            (cols[icol+1]>row+1 || vals[icol+1]==0.))) {
00613         bs = row+1; found = PETSC_TRUE;
00614       }
00615       ierr = MatRestoreRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00616       if (found) break;
00617     }
00618   }
00619   MPI_Bcast((void*)&bs,1,MPI_INT,0,comm);
00620   for (row=first; row<last; row++) {
00621     const int *cols; const PetscScalar *vals; int ncols,icol=0;
00622       PetscTruth found = PETSC_FALSE;
00623     if ((row+1)%bs!=0) continue;
00624     ierr = MatGetRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00625     while (icol<ncols && cols[icol]<row) icol++;
00626     if ( (cols[icol]==row+1 && vals[icol]!=0.) ||
00627          (cols[icol]==row && 
00628           (icol+1<ncols && cols[icol+1]==row+1 && vals[icol+1]!=0.)) ) {
00629       bs = row+1; found = PETSC_TRUE;
00630     }
00631     ierr = MatRestoreRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00632     if (found) break;
00633   }
00634   MPI_Allreduce((void*)&bs,(void*)&gbs,1,MPI_INT,MPI_MIN,comm);
00635 
00636   ierr = GetDataID("structure","blocksize",&id,&has); CHKERRQ(ierr);
00637   HASTOEXIST(has);
00638   ierr = PetscObjectComposedDataSetInt((PetscObject)A,id,gbs); CHKERRQ(ierr);
00639 
00640   PetscFunctionReturn(0);
00641 }
00642 
00643 static PetscErrorCode BlockSize
00644 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00645 {
00646   Mat A = (Mat)prob;
00647   int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00648   PetscFunctionBegin;
00649   ierr = GetDataID("structure","blocksize",&id,&has); CHKERRQ(ierr);
00650   HASTOEXIST(has);
00651   ierr = PetscObjectComposedDataGetInt
00652     ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00653   if (*flg) rv->i = v;
00654   else {
00655     ierr = regularblocks(A); CHKERRQ(ierr);
00656     ierr = PetscObjectComposedDataGetInt
00657       ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00658     if (*flg) rv->i = v;
00659   }
00660   
00661   PetscFunctionReturn(0);
00662 }
00663 
00664 /*
00665  * Registration
00666  */
00667 #undef __FUNCT__
00668 #define __FUNCT__ "RegisterStructureModules"
00669 PetscErrorCode RegisterStructureModules()
00670 {
00671   PetscErrorCode ierr;
00672   PetscFunctionBegin;
00673 
00674   ierr = RegisterModule
00675     ("structure","nrows",ANALYSISINTEGER,&nRows); CHKERRQ(ierr);
00676   ierr = RegisterModule
00677     ("structure","symmetry",ANALYSISINTEGER,&Symmetry); CHKERRQ(ierr);
00678 
00679   ierr = RegisterModule
00680     ("structure","nnzeros",ANALYSISINTEGER,&NNonZeros); CHKERRQ(ierr);
00681   ierr = RegisterModule
00682     ("structure","max-nnzeros-per-row",ANALYSISINTEGER,
00683      &MaxNNonZerosPerRow); CHKERRQ(ierr);
00684   ierr = RegisterModule
00685     ("structure","min-nnzeros-per-row",ANALYSISINTEGER,
00686      &MinNNonZerosPerRow); CHKERRQ(ierr);
00687 
00688   ierr = RegisterModule
00689     ("structure","left-bandwidth",ANALYSISINTEGER,&LBandWidth); CHKERRQ(ierr);
00690   ierr = RegisterModule
00691     ("structure","right-bandwidth",ANALYSISINTEGER,&RBandWidth); CHKERRQ(ierr);
00692 
00693   ierr = RegisterModule
00694     ("structure","left-skyline",ANALYSISINTARRAY,&LeftSkyline); CHKERRQ(ierr);
00695   ierr = RegisterModule
00696     ("structure","right-skyline",ANALYSISINTARRAY,&RightSkyline); CHKERRQ(ierr);
00697 
00698   ierr = RegisterModule
00699     ("structure","n-dummy-rows",ANALYSISINTEGER,&NDummyRows); CHKERRQ(ierr);
00700   ierr = RegisterModule
00701     ("structure","dummy-rows-kind",ANALYSISINTEGER,&DummyRowsKind); CHKERRQ(ierr);
00702   ierr = RegisterModule
00703     ("structure","dummy-rows",ANALYSISINTARRAY,&DummyRows); CHKERRQ(ierr);
00704 
00705   ierr = RegisterModule
00706     ("structure","diag-zerostart",ANALYSISINTEGER,&DiagZeroStart); CHKERRQ(ierr);
00707   ierr = RegisterModule
00708     ("structure","diag-definite",ANALYSISINTEGER,&DiagDefinite); CHKERRQ(ierr);
00709 
00710   ierr = RegisterModule
00711     ("structure","blocksize",ANALYSISINTEGER,&BlockSize); CHKERRQ(ierr);
00712 
00713   PetscFunctionReturn(0);
00714 }
00715 
00716 #undef __FUNCT__
00717 #define __FUNCT__ "DeRegisterStructureModules"
00718 PetscErrorCode DeRegisterStructureModules(void)
00719 {
00720   PetscFunctionBegin;
00721   PetscFunctionReturn(0);
00722 }