SALSA Analysis Modules
|
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,PetscBool *flg) 00055 { 00056 Mat A = (Mat)prob; 00057 int v = 0; PetscBool 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,PetscBool *flg) 00079 { 00080 Mat A = (Mat)prob; 00081 int v = 0; PetscBool 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 PetscBool 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 PetscBool 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 ierr = 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 ierr = 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 ierr = 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 ierr = 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 ierr = 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,PetscBool *flg) 00181 { 00182 Mat A = (Mat)prob; 00183 int v = 0; PetscBool 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,PetscBool *flg) 00204 { 00205 Mat A = (Mat)prob; 00206 int v = 0; PetscBool 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,PetscBool *flg) 00227 { 00228 Mat A = (Mat)prob; 00229 int v = 0; PetscBool 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,PetscBool *flg) 00250 { 00251 Mat A = (Mat)prob; 00252 int v = 0; PetscBool 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,PetscBool *flg) 00273 { 00274 Mat A = (Mat)prob; 00275 int v = 0; PetscBool 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,PetscBool *flg) 00296 { 00297 Mat A = (Mat)prob; int N; 00298 int *v = NULL; PetscBool 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,PetscBool *flg) 00323 { 00324 Mat A = (Mat)prob; int N; 00325 int *v = NULL; PetscBool 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 PetscBool 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 ierr = MPI_Comm_size(comm,&ntids); 00386 ierr = PetscMalloc(ntids*sizeof(int),&sizes); CHKERRQ(ierr); 00387 ierr = PetscMalloc(ntids*sizeof(int),&offsets); CHKERRQ(ierr); 00388 ierr = 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,PetscBool *flg) 00420 { 00421 Mat A = (Mat)prob; 00422 int v = 0; PetscBool 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,PetscBool *flg) 00443 { 00444 Mat A = (Mat)prob; 00445 int v = 0; PetscBool 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,PetscBool *flg) 00466 { 00467 Mat A = (Mat)prob; 00468 int *v = NULL; PetscBool 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 PetscBool 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 ierr = 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 ierr = 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,PetscBool *flg) 00544 { 00545 Mat A = (Mat)prob; 00546 int v = 0; PetscBool 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,PetscBool *flg) 00567 { 00568 Mat A = (Mat)prob; 00569 int v = 0; PetscBool 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; PetscBool 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 PetscBool 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 ierr = 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 PetscBool 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 ierr = 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 #undef __FUNCT__ 00644 #define __FUNCT__ "BlockSize" 00645 static PetscErrorCode BlockSize 00646 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscBool *flg) 00647 { 00648 Mat A = (Mat)prob; 00649 int v = 0; PetscBool has; int id; PetscErrorCode ierr; 00650 PetscFunctionBegin; 00651 ierr = GetDataID("structure","blocksize",&id,&has); CHKERRQ(ierr); 00652 HASTOEXIST(has); 00653 ierr = PetscObjectComposedDataGetInt 00654 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00655 if (*flg) rv->i = v; 00656 else { 00657 ierr = regularblocks(A); CHKERRQ(ierr); 00658 ierr = PetscObjectComposedDataGetInt 00659 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr); 00660 if (*flg) rv->i = v; 00661 } 00662 00663 PetscFunctionReturn(0); 00664 } 00665 00666 /* 00667 * Registration 00668 */ 00669 #undef __FUNCT__ 00670 #define __FUNCT__ "RegisterStructureModules" 00671 PetscErrorCode RegisterStructureModules() 00672 { 00673 PetscErrorCode ierr; 00674 PetscFunctionBegin; 00675 00676 ierr = RegisterModule 00677 ("structure","nrows",ANALYSISINTEGER,&nRows); CHKERRQ(ierr); 00678 ierr = RegisterModule 00679 ("structure","symmetry",ANALYSISINTEGER,&Symmetry); CHKERRQ(ierr); 00680 00681 ierr = RegisterModule 00682 ("structure","nnzeros",ANALYSISINTEGER,&NNonZeros); CHKERRQ(ierr); 00683 ierr = RegisterModule 00684 ("structure","max-nnzeros-per-row",ANALYSISINTEGER, 00685 &MaxNNonZerosPerRow); CHKERRQ(ierr); 00686 ierr = RegisterModule 00687 ("structure","min-nnzeros-per-row",ANALYSISINTEGER, 00688 &MinNNonZerosPerRow); CHKERRQ(ierr); 00689 00690 ierr = RegisterModule 00691 ("structure","left-bandwidth",ANALYSISINTEGER,&LBandWidth); CHKERRQ(ierr); 00692 ierr = RegisterModule 00693 ("structure","right-bandwidth",ANALYSISINTEGER,&RBandWidth); CHKERRQ(ierr); 00694 00695 ierr = RegisterModule 00696 ("structure","left-skyline",ANALYSISINTARRAY,&LeftSkyline); CHKERRQ(ierr); 00697 ierr = RegisterModule 00698 ("structure","right-skyline",ANALYSISINTARRAY,&RightSkyline); CHKERRQ(ierr); 00699 00700 ierr = RegisterModule 00701 ("structure","n-dummy-rows",ANALYSISINTEGER,&NDummyRows); CHKERRQ(ierr); 00702 ierr = RegisterModule 00703 ("structure","dummy-rows-kind",ANALYSISINTEGER,&DummyRowsKind); CHKERRQ(ierr); 00704 ierr = RegisterModule 00705 ("structure","dummy-rows",ANALYSISINTARRAY,&DummyRows); CHKERRQ(ierr); 00706 00707 ierr = RegisterModule 00708 ("structure","diag-zerostart",ANALYSISINTEGER,&DiagZeroStart); CHKERRQ(ierr); 00709 ierr = RegisterModule 00710 ("structure","diag-definite",ANALYSISINTEGER,&DiagDefinite); CHKERRQ(ierr); 00711 00712 ierr = RegisterModule 00713 ("structure","blocksize",ANALYSISINTEGER,&BlockSize); CHKERRQ(ierr); 00714 00715 PetscFunctionReturn(0); 00716 } 00717 00718 #undef __FUNCT__ 00719 #define __FUNCT__ "DeRegisterStructureModules" 00720 PetscErrorCode DeRegisterStructureModules(void) 00721 { 00722 PetscFunctionBegin; 00723 PetscFunctionReturn(0); 00724 }