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
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
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
00106
00107 #undef __FUNCT__
00108 #define __FUNCT__ "compute_nnz_structure"
00109 static PetscErrorCode compute_nnz_structure
00110 (AnaModNumericalProblem prob)
00111 {
00112 Mat A = (Mat)prob;
00113 int nnz,p,q,minnz,maxnz,minnnz,maxnnz;
00114 MPI_Comm comm;
00115 int N,first,last,local_size,row,nz,*Bl,*Br,bl,br, id;
00116 PetscTruth has; PetscErrorCode ierr;
00117
00118 PetscFunctionBegin;
00119 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00120 ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr);
00121 ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00122 local_size = last-first;
00123 ierr = PetscMalloc(local_size*sizeof(int),&Bl); CHKERRQ(ierr);
00124 ierr = PetscMalloc(local_size*sizeof(int),&Br); CHKERRQ(ierr);
00125 bl = 0; br = 0; nz = 0; minnz = maxnz = 0;
00126 for (row=first; row<last; row++) {
00127 int ncols; const int *cols;
00128 ierr = MatGetRow(A,row,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr);
00129 nz += ncols;
00130 if (row==first) minnz = ncols;
00131 if (ncols<minnz) minnz = ncols;
00132 if (ncols>maxnz) maxnz = ncols;
00133 if (ncols) {
00134 Bl[row-first] = bl = PetscMax(bl,row-cols[0]);
00135 Br[row-first] = br = PetscMax(br,cols[ncols-1]-row);
00136 }
00137 ierr = MatRestoreRow(A,row,&ncols,&cols,PETSC_NULL); CHKERRQ(ierr);
00138 }
00139
00140 MPI_Allreduce((void*)&nz,(void*)&nnz,1,MPI_INT,MPI_SUM,comm);
00141 ierr = GetDataID("structure","nnzeros",&id,&has); CHKERRQ(ierr);
00142 HASTOEXIST(has);
00143 ierr = PetscObjectComposedDataSetInt
00144 ((PetscObject)A,id,nnz); CHKERRQ(ierr);
00145
00146 MPI_Allreduce((void*)&maxnz,(void*)&maxnnz,1,MPI_INT,MPI_MAX,comm);
00147 ierr = GetDataID("structure","max-nnzeros-per-row",&id,&has); CHKERRQ(ierr);
00148 HASTOEXIST(has);
00149 ierr = PetscObjectComposedDataSetInt
00150 ((PetscObject)A,id,maxnnz); CHKERRQ(ierr);
00151
00152 MPI_Allreduce((void*)&minnz,(void*)&minnnz,1,MPI_INT,MPI_MIN,comm);
00153 ierr = GetDataID("structure","min-nnzeros-per-row",&id,&has); CHKERRQ(ierr);
00154 HASTOEXIST(has);
00155 ierr = PetscObjectComposedDataSetInt
00156 ((PetscObject)A,id,minnnz); CHKERRQ(ierr);
00157
00158 MPI_Allreduce((void*)&bl,(void*)&p,1,MPI_INT,MPI_MAX,comm);
00159 ierr = GetDataID("structure","left-bandwidth",&id,&has); CHKERRQ(ierr);
00160 HASTOEXIST(has);
00161 ierr = PetscObjectComposedDataSetInt
00162 ((PetscObject)A,id,p); CHKERRQ(ierr);
00163
00164 MPI_Allreduce((void*)&br,(void*)&q,1,MPI_INT,MPI_MAX,comm);
00165 ierr = GetDataID("structure","right-bandwidth",&id,&has); CHKERRQ(ierr);
00166 HASTOEXIST(has);
00167 ierr = PetscObjectComposedDataSetInt
00168 ((PetscObject)A,id,q); CHKERRQ(ierr);
00169
00170 {
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181 ierr = PetscFree(Bl); CHKERRQ(ierr);
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 ierr = PetscFree(Br); CHKERRQ(ierr);
00192
00193 }
00194 PetscFunctionReturn(0);
00195 }
00196
00197 #undef __FUNCT__
00198 #define __FUNCT__ "NNonZeros"
00199 static PetscErrorCode NNonZeros
00200 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00201 {
00202 Mat A = (Mat)prob;
00203 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00204 PetscFunctionBegin;
00205 ierr = GetDataID("structure","nnzeros",&id,&has); CHKERRQ(ierr);
00206 HASTOEXIST(has);
00207 ierr = PetscObjectComposedDataGetInt
00208 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00209 if (*flg) rv->i = v;
00210 else {
00211 ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00212 ierr = PetscObjectComposedDataGetInt
00213 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00214 if (*flg) rv->i = v;
00215 }
00216
00217 PetscFunctionReturn(0);
00218 }
00219
00220 #undef __FUNCT__
00221 #define __FUNCT__ "MaxNNonZerosPerRow"
00222 static PetscErrorCode MaxNNonZerosPerRow
00223 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00224 {
00225 Mat A = (Mat)prob;
00226 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00227 PetscFunctionBegin;
00228 ierr = GetDataID("structure","max-nnzeros-per-row",&id,&has); CHKERRQ(ierr);
00229 HASTOEXIST(has);
00230 ierr = PetscObjectComposedDataGetInt
00231 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00232 if (*flg) rv->i = v;
00233 else {
00234 ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00235 ierr = PetscObjectComposedDataGetInt
00236 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00237 if (*flg) rv->i = v;
00238 }
00239
00240 PetscFunctionReturn(0);
00241 }
00242
00243 #undef __FUNCT__
00244 #define __FUNCT__ "MinNNonZerosPerRow"
00245 static PetscErrorCode MinNNonZerosPerRow
00246 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00247 {
00248 Mat A = (Mat)prob;
00249 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00250 PetscFunctionBegin;
00251 ierr = GetDataID("structure","min-nnzeros-per-row",&id,&has); CHKERRQ(ierr);
00252 HASTOEXIST(has);
00253 ierr = PetscObjectComposedDataGetInt
00254 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00255 if (*flg) rv->i = v;
00256 else {
00257 ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00258 ierr = PetscObjectComposedDataGetInt
00259 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00260 if (*flg) rv->i = v;
00261 }
00262
00263 PetscFunctionReturn(0);
00264 }
00265
00266 #undef __FUNCT__
00267 #define __FUNCT__ "LBandWidth"
00268 static PetscErrorCode LBandWidth
00269 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00270 {
00271 Mat A = (Mat)prob;
00272 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00273 PetscFunctionBegin;
00274 ierr = GetDataID("structure","left-bandwidth",&id,&has); CHKERRQ(ierr);
00275 HASTOEXIST(has);
00276 ierr = PetscObjectComposedDataGetInt
00277 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00278 if (*flg) rv->i = v;
00279 else {
00280 ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00281 ierr = PetscObjectComposedDataGetInt
00282 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00283 if (*flg) rv->i = v;
00284 }
00285
00286 PetscFunctionReturn(0);
00287 }
00288
00289 #undef __FUNCT__
00290 #define __FUNCT__ "RBandWidth"
00291 static PetscErrorCode RBandWidth
00292 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00293 {
00294 Mat A = (Mat)prob;
00295 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00296 PetscFunctionBegin;
00297 ierr = GetDataID("structure","right-bandwidth",&id,&has); CHKERRQ(ierr);
00298 HASTOEXIST(has);
00299 ierr = PetscObjectComposedDataGetInt
00300 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00301 if (*flg) rv->i = v;
00302 else {
00303 ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00304 ierr = PetscObjectComposedDataGetInt
00305 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00306 if (*flg) rv->i = v;
00307 }
00308
00309 PetscFunctionReturn(0);
00310 }
00311
00312 #undef __FUNCT__
00313 #define __FUNCT__ "LeftSkyline"
00314 static PetscErrorCode LeftSkyline
00315 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00316 {
00317 Mat A = (Mat)prob; int N;
00318 int *v = NULL; PetscTruth has; int id; PetscErrorCode ierr;
00319 PetscFunctionBegin;
00320 ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr);
00321 ierr = GetDataID("structure","left-skyline",&id,&has); CHKERRQ(ierr);
00322 HASTOEXIST(has);
00323 ierr = PetscObjectComposedDataGetIntstar
00324 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00325 if (*flg) {
00326 rv->ii = v; *lv = N;
00327 } else {
00328 ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00329 ierr = PetscObjectComposedDataGetIntstar
00330 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00331 if (*flg) {
00332 rv->ii = v; *lv = N;
00333 }
00334 }
00335
00336 PetscFunctionReturn(0);
00337 }
00338
00339 #undef __FUNCT__
00340 #define __FUNCT__ "RightSkyline"
00341 static PetscErrorCode RightSkyline
00342 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00343 {
00344 Mat A = (Mat)prob; int N;
00345 int *v = NULL; PetscTruth has; int id; PetscErrorCode ierr;
00346 PetscFunctionBegin;
00347 ierr = MatGetSize(A,&N,PETSC_NULL); CHKERRQ(ierr);
00348 ierr = GetDataID("structure","right-skyline",&id,&has); CHKERRQ(ierr);
00349 HASTOEXIST(has);
00350 ierr = PetscObjectComposedDataGetIntstar
00351 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00352 if (*flg) {
00353 rv->ii = v; *lv = N;
00354 } else {
00355 ierr = compute_nnz_structure(A); CHKERRQ(ierr);
00356 ierr = PetscObjectComposedDataGetIntstar
00357 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00358 if (*flg) {
00359 rv->ii = v; *lv = N;
00360 }
00361 }
00362
00363 PetscFunctionReturn(0);
00364 }
00365
00366
00367
00368
00369 #undef __FUNCT__
00370 #define __FUNCT__ "compute_dummy_rows"
00371 static PetscErrorCode compute_dummy_rows(Mat A)
00372 {
00373 int gkind;
00374 MPI_Comm comm;
00375 int first,last,row,ncols,ndummies,Ndummies,*dummies,*Dummies,kind,id;
00376 PetscTruth has; PetscErrorCode ierr;
00377 const int *cols; const PetscScalar *vals;
00378
00379 PetscFunctionBegin;
00380 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00381 ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00382 ierr = PetscMalloc((last-first+1)*sizeof(int),&dummies); CHKERRQ(ierr);
00383
00384 ndummies = 0; kind = 0;
00385 for (row=first; row<last; row++) {
00386 ierr = MatGetRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00387 if (ncols==1) {
00388 dummies[ndummies++] = row;
00389 if (vals[0]!=1. && !kind ) kind = 1;
00390 if (cols[0]!=row) kind = 2;
00391 }
00392 ierr = MatRestoreRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00393 }
00394
00395 ierr = MPI_Allreduce
00396 ((void*)&kind,(void*)&gkind,1,MPI_INT,MPI_MAX,comm); CHKERRQ(ierr);
00397 ierr = GetDataID("structure","dummy-rows-kind",&id,&has); CHKERRQ(ierr);
00398 HASTOEXIST(has);
00399 ierr = PetscObjectComposedDataSetInt
00400 ((PetscObject)A,id,gkind); CHKERRQ(ierr);
00401
00402 {
00403 int *sizes,*offsets,ntids,i,s;
00404
00405 MPI_Comm_size(comm,&ntids);
00406 ierr = PetscMalloc(ntids*sizeof(int),&sizes); CHKERRQ(ierr);
00407 ierr = PetscMalloc(ntids*sizeof(int),&offsets); CHKERRQ(ierr);
00408 MPI_Allgather(&ndummies,1,MPI_INT,sizes,1,MPI_INT,comm);
00409 s=0; for (i=0; i<ntids; i++) {offsets[i]=s; s+=sizes[i];} Ndummies = s;
00410 if (Ndummies==0)
00411 Dummies = NULL;
00412 else {
00413 ierr = PetscMalloc(Ndummies*sizeof(int),&Dummies); CHKERRQ(ierr);
00414 ierr = MPI_Allgatherv
00415 (dummies,ndummies,MPI_INT,Dummies,sizes,offsets,MPI_INT,
00416 comm); CHKERRQ(ierr);
00417 }
00418 ierr = PetscFree(dummies); CHKERRQ(ierr);
00419 ierr = PetscFree(offsets); CHKERRQ(ierr);
00420 ierr = PetscFree(sizes); CHKERRQ(ierr);
00421 }
00422
00423 ierr = GetDataID("structure","n-dummy-rows",&id,&has); CHKERRQ(ierr);
00424 HASTOEXIST(has);
00425 ierr = PetscObjectComposedDataSetInt
00426 ((PetscObject)A,id,Ndummies); CHKERRQ(ierr);
00427
00428 ierr = GetDataID("structure","dummy-rows",&id,&has); CHKERRQ(ierr);
00429 HASTOEXIST(has);
00430 ierr = PetscObjectComposedDataSetIntstar
00431 ((PetscObject)A,id,Dummies); CHKERRQ(ierr);
00432
00433 PetscFunctionReturn(0);
00434 }
00435
00436 #undef __FUNCT__
00437 #define __FUNCT__ "NDummyRows"
00438 static PetscErrorCode NDummyRows
00439 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00440 {
00441 Mat A = (Mat)prob;
00442 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00443 PetscFunctionBegin;
00444 ierr = GetDataID("structure","n-dummy-rows",&id,&has); CHKERRQ(ierr);
00445 HASTOEXIST(has);
00446 ierr = PetscObjectComposedDataGetInt
00447 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00448 if (*flg) rv->i = v;
00449 else {
00450 ierr = compute_dummy_rows(A); CHKERRQ(ierr);
00451 ierr = PetscObjectComposedDataGetInt
00452 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00453 if (*flg) rv->i = v;
00454 }
00455
00456 PetscFunctionReturn(0);
00457 }
00458
00459 #undef __FUNCT__
00460 #define __FUNCT__ "DummyRowsKind"
00461 static PetscErrorCode DummyRowsKind
00462 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00463 {
00464 Mat A = (Mat)prob;
00465 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00466 PetscFunctionBegin;
00467 ierr = GetDataID("structure","dummy-rows-kind",&id,&has); CHKERRQ(ierr);
00468 HASTOEXIST(has);
00469 ierr = PetscObjectComposedDataGetInt
00470 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00471 if (*flg) rv->i = v;
00472 else {
00473 ierr = compute_dummy_rows(A); CHKERRQ(ierr);
00474 ierr = PetscObjectComposedDataGetInt
00475 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00476 if (*flg) rv->i = v;
00477 }
00478
00479 PetscFunctionReturn(0);
00480 }
00481
00482 #undef __FUNCT__
00483 #define __FUNCT__ "DummyRows"
00484 static PetscErrorCode DummyRows
00485 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00486 {
00487 Mat A = (Mat)prob;
00488 int *v = NULL; PetscTruth has; int id,id2; PetscErrorCode ierr;
00489 PetscFunctionBegin;
00490 ierr = GetDataID("structure","dummy-rows",&id,&has); CHKERRQ(ierr);
00491 HASTOEXIST(has);
00492 ierr = PetscObjectComposedDataGetIntstar
00493 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00494 ierr = GetDataID("structure","n-dummy-rows",&id2,&has); CHKERRQ(ierr);
00495 HASTOEXIST(has);
00496 ierr = PetscObjectComposedDataGetInt
00497 ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00498 if (*flg) {
00499 rv->ii = v;
00500 } else {
00501 ierr = compute_dummy_rows(A); CHKERRQ(ierr);
00502 ierr = PetscObjectComposedDataGetIntstar
00503 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00504 ierr = PetscObjectComposedDataGetInt
00505 ((PetscObject)A,id2,*lv,*flg); CHKERRQ(ierr);
00506 if (*flg) rv->ii = v;
00507 }
00508
00509 PetscFunctionReturn(0);
00510 }
00511
00512
00513
00514
00515 #undef __FUNCT__
00516 #define __FUNCT__ "compute_posdiag"
00517 static PetscErrorCode compute_posdiag
00518 (AnaModNumericalProblem prob)
00519 {
00520 Mat A = (Mat)prob;
00521 MPI_Comm comm;
00522 Vec D;
00523 PetscScalar *d;
00524 int i,first,last,m,M,di,df,Di,Df,pd,id;
00525 PetscTruth has; PetscErrorCode ierr;
00526
00527 PetscFunctionBegin;
00528 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00529 ierr = MatGetSize(A,&M,PETSC_NULL); CHKERRQ(ierr);
00530 ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00531 m = last-first;
00532 ierr = VecCreateMPI(comm,m,PETSC_DECIDE,&D); CHKERRQ(ierr);
00533 ierr = MatGetDiagonal(A,D); CHKERRQ(ierr);
00534 ierr = VecGetArray(D,&d); CHKERRQ(ierr);
00535 di = 1; df = first;
00536 for (i=0; i<m; i++) {
00537 if (d[i]<=0) di = 0;
00538 if (d[i]!=0) df = i+1;
00539 }
00540 MPI_Allreduce((void*)&di,(void*)&Di,1,MPI_INT,MPI_PROD,comm);
00541 ierr = GetDataID("structure","diag-definite",&id,&has); CHKERRQ(ierr);
00542 HASTOEXIST(has);
00543 ierr = PetscObjectComposedDataSetInt
00544 ((PetscObject)A,id,Di); CHKERRQ(ierr);
00545
00546 MPI_Allreduce((void*)&df,(void*)&Df,1,MPI_INT,MPI_MAX,comm);
00547 ierr = GetDataID("structure","diag-zerostart",&id,&has); CHKERRQ(ierr);
00548 HASTOEXIST(has);
00549 ierr = PetscObjectComposedDataSetInt
00550 ((PetscObject)A,id,Df); CHKERRQ(ierr);
00551
00552
00553 if (Df<M) pd = 2+Df; else pd = Di;
00554
00555 ierr = VecRestoreArray(D,&d); CHKERRQ(ierr);
00556 ierr = VecDestroy(D); CHKERRQ(ierr);
00557 PetscFunctionReturn(0);
00558 }
00559
00560 #undef __FUNCT__
00561 #define __FUNCT__ "DiagZeroStart"
00562 static PetscErrorCode DiagZeroStart
00563 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00564 {
00565 Mat A = (Mat)prob;
00566 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00567 PetscFunctionBegin;
00568 ierr = GetDataID("structure","diag-zerostart",&id,&has); CHKERRQ(ierr);
00569 HASTOEXIST(has);
00570 ierr = PetscObjectComposedDataGetInt
00571 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00572 if (*flg) rv->i = v;
00573 else {
00574 ierr = compute_posdiag(A); CHKERRQ(ierr);
00575 ierr = PetscObjectComposedDataGetInt
00576 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00577 if (*flg) rv->i = v;
00578 }
00579
00580 PetscFunctionReturn(0);
00581 }
00582
00583 #undef __FUNCT__
00584 #define __FUNCT__ "DiagDefinite"
00585 static PetscErrorCode DiagDefinite
00586 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00587 {
00588 Mat A = (Mat)prob;
00589 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00590 PetscFunctionBegin;
00591 ierr = GetDataID("structure","diag-definite",&id,&has); CHKERRQ(ierr);
00592 HASTOEXIST(has);
00593 ierr = PetscObjectComposedDataGetInt
00594 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00595 if (*flg) rv->i = v;
00596 else {
00597 ierr = compute_posdiag(A); CHKERRQ(ierr);
00598 ierr = PetscObjectComposedDataGetInt
00599 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00600 if (*flg) rv->i = v;
00601 }
00602
00603 PetscFunctionReturn(0);
00604 }
00605
00606
00607
00608
00609 #undef __FUNCT__
00610 #define __FUNCT__ "regularblocks"
00611 static PetscErrorCode regularblocks
00612 (AnaModNumericalProblem prob)
00613 {
00614 Mat A = (Mat)prob;
00615 int first,last,local_size,row, bs,gbs,id,mytid;
00616 MPI_Comm comm; PetscTruth has; PetscErrorCode ierr;
00617 PetscFunctionBegin;
00618 ierr = PetscObjectGetComm((PetscObject)A,&comm); CHKERRQ(ierr);
00619 ierr = MPI_Comm_rank(comm,&mytid); CHKERRQ(ierr);
00620 ierr = MatGetOwnershipRange(A,&first,&last); CHKERRQ(ierr);
00621 local_size = last-first;
00622 if (mytid==0) {
00623 int row;
00624 bs = 1;
00625 for (row=first; row<last; row++) {
00626 const int *cols; const PetscScalar *vals; int ncols,icol=0;
00627 PetscTruth found = PETSC_FALSE;
00628 ierr = MatGetRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00629 while (icol<ncols && cols[icol]<row) icol++;
00630 if (cols[icol]>row+1 ||
00631 (cols[icol]==row && icol+1<ncols &&
00632 (cols[icol+1]>row+1 || vals[icol+1]==0.))) {
00633 bs = row+1; found = PETSC_TRUE;
00634 }
00635 ierr = MatRestoreRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00636 if (found) break;
00637 }
00638 }
00639 MPI_Bcast((void*)&bs,1,MPI_INT,0,comm);
00640 for (row=first; row<last; row++) {
00641 const int *cols; const PetscScalar *vals; int ncols,icol=0;
00642 PetscTruth found = PETSC_FALSE;
00643 if ((row+1)%bs!=0) continue;
00644 ierr = MatGetRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00645 while (icol<ncols && cols[icol]<row) icol++;
00646 if ( (cols[icol]==row+1 && vals[icol]!=0.) ||
00647 (cols[icol]==row &&
00648 (icol+1<ncols && cols[icol+1]==row+1 && vals[icol+1]!=0.)) ) {
00649 bs = row+1; found = PETSC_TRUE;
00650 }
00651 ierr = MatRestoreRow(A,row,&ncols,&cols,&vals); CHKERRQ(ierr);
00652 if (found) break;
00653 }
00654 MPI_Allreduce((void*)&bs,(void*)&gbs,1,MPI_INT,MPI_MIN,comm);
00655
00656 ierr = GetDataID("structure","blocksize",&id,&has); CHKERRQ(ierr);
00657 HASTOEXIST(has);
00658 ierr = PetscObjectComposedDataSetInt((PetscObject)A,id,gbs); CHKERRQ(ierr);
00659
00660 PetscFunctionReturn(0);
00661 }
00662
00663 static PetscErrorCode BlockSize
00664 (AnaModNumericalProblem prob,AnalysisItem *rv,int *lv,PetscTruth *flg)
00665 {
00666 Mat A = (Mat)prob;
00667 int v = 0; PetscTruth has; int id; PetscErrorCode ierr;
00668 PetscFunctionBegin;
00669 ierr = GetDataID("structure","blocksize",&id,&has); CHKERRQ(ierr);
00670 HASTOEXIST(has);
00671 ierr = PetscObjectComposedDataGetInt
00672 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00673 if (*flg) rv->i = v;
00674 else {
00675 ierr = regularblocks(A); CHKERRQ(ierr);
00676 ierr = PetscObjectComposedDataGetInt
00677 ((PetscObject)A,id,v,*flg); CHKERRQ(ierr);
00678 if (*flg) rv->i = v;
00679 }
00680
00681 PetscFunctionReturn(0);
00682 }
00683
00684
00685
00686
00687 #undef __FUNCT__
00688 #define __FUNCT__ "RegisterStructureModules"
00689 PetscErrorCode RegisterStructureModules()
00690 {
00691 PetscErrorCode ierr;
00692 PetscFunctionBegin;
00693
00694 ierr = RegisterModule
00695 ("structure","nrows",ANALYSISINTEGER,&nRows); CHKERRQ(ierr);
00696 ierr = RegisterModule
00697 ("structure","symmetry",ANALYSISINTEGER,&Symmetry); CHKERRQ(ierr);
00698
00699 ierr = RegisterModule
00700 ("structure","nnzeros",ANALYSISINTEGER,&NNonZeros); CHKERRQ(ierr);
00701 ierr = RegisterModule
00702 ("structure","max-nnzeros-per-row",ANALYSISINTEGER,
00703 &MaxNNonZerosPerRow); CHKERRQ(ierr);
00704 ierr = RegisterModule
00705 ("structure","min-nnzeros-per-row",ANALYSISINTEGER,
00706 &MinNNonZerosPerRow); CHKERRQ(ierr);
00707
00708 ierr = RegisterModule
00709 ("structure","left-bandwidth",ANALYSISINTEGER,&LBandWidth); CHKERRQ(ierr);
00710 ierr = RegisterModule
00711 ("structure","right-bandwidth",ANALYSISINTEGER,&RBandWidth); CHKERRQ(ierr);
00712
00713 ierr = RegisterModule
00714 ("structure","left-skyline",ANALYSISINTARRAY,&LeftSkyline); CHKERRQ(ierr);
00715 ierr = RegisterModule
00716 ("structure","right-skyline",ANALYSISINTARRAY,&RightSkyline); CHKERRQ(ierr);
00717
00718 ierr = RegisterModule
00719 ("structure","n-dummy-rows",ANALYSISINTEGER,&NDummyRows); CHKERRQ(ierr);
00720 ierr = RegisterModule
00721 ("structure","dummy-rows-kind",ANALYSISINTEGER,&DummyRowsKind); CHKERRQ(ierr);
00722 ierr = RegisterModule
00723 ("structure","dummy-rows",ANALYSISINTARRAY,&DummyRows); CHKERRQ(ierr);
00724
00725 ierr = RegisterModule
00726 ("structure","diag-zerostart",ANALYSISINTEGER,&DiagZeroStart); CHKERRQ(ierr);
00727 ierr = RegisterModule
00728 ("structure","diag-definite",ANALYSISINTEGER,&DiagDefinite); CHKERRQ(ierr);
00729
00730 ierr = RegisterModule
00731 ("structure","blocksize",ANALYSISINTEGER,&BlockSize); CHKERRQ(ierr);
00732
00733 PetscFunctionReturn(0);
00734 }
00735
00736 #undef __FUNCT__
00737 #define __FUNCT__ "DeRegisterStructureModules"
00738 PetscErrorCode DeRegisterStructureModules(void)
00739 {
00740 PetscFunctionBegin;
00741 PetscFunctionReturn(0);
00742 }