prometheus.hh Source File

00001 // $Id: prometheus.hh,v 1.34 2004/05/07 20:47:41 adams Exp $ 00002 // Author: Mark F. Adams 00003 // Copyright (c) 2000 by Mark F. Adams 00004 // Filename: prometheus.hh 00005 // 00006 #ifndef __PROMETHEUS_H__ 00007 #define __PROMETHEUS_H__ 00008 00009 #include "prom_grid.hh" 00010 00011 //#define PROM_MEASURE_MPATH 00012 00013 /***************************************************************************** 00014 * Prometheus 00015 ****************************************************************************/ 00016 class PromSolver; 00017 class PromMG; 00019 00022 class Prometheus : public PromContext 00023 { 00024 public: 00025 Prometheus( const PromComm &comm, PromOptions &opts, 00026 PromPerfMonitor &perf ); 00027 ~Prometheus(); 00028 int Init( bool init1 ); 00029 virtual int Archive( FILE *file = stderr, 00030 const PromArchiveType type = PROM_PRINT); 00032 int NewLevel( int ); 00033 private: 00034 bool allDone( PromGrid *lastGrid, const int topMaxEq ) const ; 00035 int MakeNewGrid(bool islast, const int topMaxEq, 00036 PromGrid *const lastGrid, PromGrid *& newGrid ) const; 00037 public: 00038 int FinalizeSymbolic(); 00039 int UpdateGeometry( double * const crds, double * const deltas ); 00040 int SetMatrixWithNewData( PromMatrix *AA ); 00041 int CreateFineGrid(const int npLoc, const int lid_ndof[], const int ndf, 00042 const PromTable<int> &gid_ghostLid, 00043 const PromList<int*> *elemList); 00044 int getNumLevels() const; 00045 int nBlocksWithGlobalDOF( const int gl ) const; 00046 int WritePartitioned(); 00047 int MakeSplitCoarseGrid(PromGrid *lastGrid, const int nloc, 00048 const int nnodes, PromGrid **newGrid, 00049 const int factor ) const; 00050 int ConvertToSymm(PromMatrix *const A, const PromOptions &opts)const; 00051 PromMGType getMGType() const { return mgType_; } 00052 00053 // Data 00054 bool finalized_; 00055 PromMGType mgType_; 00056 PromGrid *top_grid_; 00057 PromGrid *grid0_; 00058 PromVector *D1_2inv_; 00059 PromSolver *solver_; 00060 //PromVector *xVec_; // used in non-fei only! 00061 //PromVector *rVec_; // used in non-fei only! 00062 //PromVector *bVec_; // used in non-fei only! 00063 // 00064 const PromComm &Comm_; 00065 PromOptions &options_; 00066 PromPerfMonitor &perf_mon_; 00067 PromParentBank parBank_; 00068 }; // Prometheus 00069 00070 /***************************************************************************** 00071 * PromProjKKTSolver 00072 ****************************************************************************/ 00074 00078 class PromProjKKTSolver 00079 { 00080 public: 00081 PromProjKKTSolver() : projC_(NULL), CCt_inv_(NULL), workCR_(NULL), 00082 ipiv_(NULL), ksp_mult_work_(NULL), projLEq_geq_(NULL){} 00083 ~PromProjKKTSolver(){ 00084 if( projC_ != NULL ) ClearProjC(); 00085 else assert(workCR_ == NULL); 00086 } 00087 int ClearProjC(){ 00088 assert(projC_ != NULL); 00089 delete projC_; delete workCR_; delete ksp_mult_work_; 00090 PetscFree( ipiv_ ); PetscFree( CCt_inv_ ); 00091 if( projLEq_geq_ != NULL ) PetscFree( projLEq_geq_ ); 00092 // clear other proj data 00093 projC_ = NULL; 00094 projLEq_geq_ = ipiv_ = NULL; 00095 CCt_inv_ = NULL; 00096 ksp_mult_work_ = workCR_ = NULL; 00097 return 0; 00098 } 00099 int SetProjC( const PromMatrix * const C, 00100 PromTable<int> *projCREqs ); 00102 int PreSolve( const PromMatrix * const A, const PromVector* const, 00103 PromVector * const bb, PromVector * const x_0 ) const; 00104 int PostSolve( const PromVector * const x_0, PromVector* const xx ) const; 00105 int ComputeMultiplier( const PromMatrix * const A, 00106 const PromMatrix * const otherC, 00107 const PromVector* const xx, PromVector* const pp, 00108 const PromVector * const bb ) const; 00109 int ApplyProj( const PromVector * const XX, PromVector * const YY ) const; 00110 int Apply_CCInv_C(const PromVector * const XX, PromVector * const YY) const; 00111 int Apply_Ct_CCInv(const PromVector * const XX, PromVector*const YY) const; 00112 int Apply_Ct_add( const PromVector * const GG, PromVector * const XX, 00113 PromVector * const YY ) const; 00114 int ApplyCCInv( PromVector * const XX ) const; 00115 private: 00116 int MapProjCRsToAllCRs(const PromVector*const XX,PromVector*const YY )const; 00117 int MapAllCRsToProjCRs(const PromVector*const GG,PromVector * const YY)const; 00118 public: 00119 bool isConforming( const PromVector * const XX, double &beta ) const; 00120 static int SparseDot(const int nz1, const int ip1[], const double vv1[], 00121 const int nz2, const int ip2[], const double vv2[], 00122 double &dot ); 00123 // data 00124 const PromMatrix *projC_; // C 00125 mutable PromVector *ksp_mult_work_; 00126 private: 00127 double *CCt_inv_; // ( C C' ) ^-1proj 00128 int *ipiv_; 00129 int *projLEq_geq_; 00130 mutable PromVector *workCR_; 00131 }; 00132 00133 /***************************************************************************** 00134 * PromSolver 00135 ****************************************************************************/ 00137 00141 class PromSolver 00142 { 00143 friend class PromPCMG; 00144 friend class PromUzawaSolver; 00145 friend class PromMatrix_base; 00146 friend class Prometheus_LinSysCore; 00147 public: 00149 PromSolver( const PromOptions &opts, PromPerfMonitor &perf, 00150 const PromMap &map, const PromMap * const pmap = NULL ) : 00151 options_(opts), perf_mon_(perf), alias_(NULL), pc_(NULL), 00152 rtol_(1.e-5), dtol_(1.e5), atol_(1.e-29), A_(NULL), 00153 nwork_(0),work_(NULL),projector_(NULL), 00154 verbose_(opts.verbose_), cg_restart_(-300), map_(map), p_map_(pmap), alloced_(TRUE), 00155 ksp_type_(PROMKSPGMRES), maxits_(1),nEig_(0){} 00157 PromSolver( PromSolver * const solver ) : 00158 options_(solver->options_), perf_mon_(solver->perf_mon_), alias_(solver), 00159 pc_(NULL), rtol_(solver->rtol_), dtol_(solver->dtol_), 00160 atol_(solver->atol_), projector_(NULL), 00161 A_(NULL),nwork_(0),work_(NULL),verbose_(solver->options_.verbose_), 00162 alloced_(FALSE), cg_restart_(-300), 00163 map_(solver->map_),p_map_(solver->p_map_), ksp_type_(PROMNUMKSP), 00164 maxits_(solver->maxits_),nEig_(0) 00165 { 00166 assert(solver != NULL && solver->alias_ == NULL); solver->alias_ = this; 00167 } 00169 ~PromSolver(); 00170 // methods 00171 int SetPC( PromPC_base *const pc, const PromMatrix_base *AA ){ 00172 A_ = AA; 00173 assert(pc_==NULL); pc_ = pc; // hack for rho solver 00174 return 0; 00175 } 00176 int PrintInfo( const int nactivep, const int lev = -1 ) const; 00177 int SetUp(); 00178 // 00180 int KSPMult( const PromMatrix_base * const A, 00181 const PromVector_base * const XX, 00182 PromVector_base * const YY ) const; 00183 int SetTolerances(const int maxits, double rtol = -1., 00184 double atol = -1., double dtol = -1. ); 00185 int SetGMRESRestart( int kmax ){ maxits_ = kmax; return 0; } 00186 // KSP 00187 int SetNewPVec( const PromMap &map ); 00188 int SetKSPType( PromKSPType type ); 00189 PromKSPType getKSPType() const { 00190 assert( this != NULL ); 00191 if( !alloced_ ) return alias_->getKSPType(); 00192 return ksp_type_; 00193 } 00194 int SetKSP( const int maxit, const char* const prestring = NULL, 00195 const char* const ksp_string = NULL ); 00196 int SetRTol( double tol ) { 00197 assert(this!=NULL); // !alloc 00198 rtol_ = tol; 00199 return 0; 00200 } 00201 double getRTol() const { 00202 assert(this!=NULL); // !alloc 00203 return rtol_; 00204 } 00205 double getATol() const { 00206 assert(this!=NULL); // !alloc 00207 return atol_; 00208 } 00209 int nEig()const{ return nEig_;} 00210 #undef __FUNCT__ 00211 #define __FUNCT__ "PromSolver::SetNEig" 00212 int SetNEig( int nn ) { 00213 if(this==NULL || nn<=0 )SETERRQ(-1,"this==NULL || nn<=0"); 00214 if( origin_.e != NULL ) { 00215 #if defined(PROM_HAVE_COMPLEX) 00216 PromComplex *pt = dynamic_cast<PromComplex*>(origin_.c); 00217 if( pt != NULL ) delete [] pt; 00218 else { 00219 int ierr = PetscFree( origin_.e ); CHKERRQ(ierr); 00220 } 00221 #else 00222 assert(dd_.e[nEig_+1] == -13.0); 00223 int ierr = PetscFree( origin_.e ); CHKERRQ(ierr); 00224 #endif 00225 origin_.e = d_.e = ee_.e = dd_.e = NULL; 00226 } 00227 nEig_ = nn; 00228 return 0; 00229 } 00230 // 00231 bool computeSVD() const { 00232 return ( (nEig() > 0 && ksp_type_ == PROMKSPCG) || 00233 ksp_type_ == PROMKSPGMRES || ksp_type_ == PROMKSPFGMRES ); 00234 } 00235 int SetComputeSVD(); 00236 int ComputeExtremeSVD( int nits, double &emax, double &emin ); 00237 private: 00238 int ComputeExtremeSVD_cg( int nn, double &emax, double &emin ); 00239 int ComputeExtremeSVD_gmres( int n, double &emax, double &emin ); 00240 public: 00241 // KSP methods 00242 int SolveGMRES(const PromVector_base * const, 00243 PromVector_base* const,int *its, 00244 const bool zeroguess = TRUE, 00245 const bool avoidnorms = FALSE, 00246 const PromPCSide = PROM_RIGHT ); 00247 private: 00248 int ResidualGMRES( const PromPCSide pc_side, const PromVector_base *const xx, 00249 PromVector_base* const vv_0, 00250 const PromVector_base* const binvf, 00251 PromVector_base* const temp ); 00252 int CycleGMRES( const PromPCSide pc_side, const PromVector_base * const bb, 00253 PromVector_base * const xx, int *reason, 00254 int *itcount, double *rnorm0, 00255 const bool avoidnorms ); 00256 int UpdateGMRESHessenberg( const int it, bool hapend, double *res ); 00257 int BuildGMRESSoln( const PromPCSide pc_side, double* nrs, 00258 PromVector_base *vdest, 00259 const int it, const PromVector_base *const* vv_0, 00260 PromVector_base * const temp); 00261 int ModifiedGramSchmidtOrthogonalization(const int it); 00262 int UnmodifiedGramSchmidtOrthogonalization(const int it); 00263 int TwoUnmodifiedGramSchmidtOrthogonalization(const int it); 00264 #if defined(PROM_HAVE_COMPLEX) 00265 int UpdateGMRESHessenberg_complex( const int it, bool hapend, double *res ); 00266 int BuildGMRESSoln_complex( const PromPCSide pc_side, PromComplex* nrs, 00267 PromVector_base *vdest, 00268 const int it, const PromVector_base *const* vv_0, 00269 PromVector_base * const temp); 00270 int ModifiedGramSchmidtOrthogonalization_complex(const int it); 00271 int UnmodifiedGramSchmidtOrthogonalization_complex(const int it); 00272 int TwoUnmodifiedGramSchmidtOrthogonalization_complex(const int it); 00273 #endif 00274 int PromResidual( const PromPCSide pc_side, const bool zerox, 00275 PromVector_base *vsoln, PromVector_base *vt1, 00276 PromVector_base *vt2, PromVector_base *vres, 00277 PromVector_base *vbinvf, 00278 const PromVector_base * const vb ); 00279 int SandwPre( const double omV[], const int deg, PromPC_base * const pc, 00280 const PromMatrix_base * const Amat, PromVector_base *const xx, 00281 PromVector_base *const yy, PromVector_base *const work) const; 00282 int SandwPost( const double omV[], const int deg, PromPC_base * const pc, 00283 const PromMatrix_base * const Amat, PromVector_base *const xx, 00284 PromVector_base *const yy, PromVector_base *const work)const; 00285 int BrezParameters( const PromMatrix_base * const AA, 00286 double &mlsOver, double &mlsCf, 00287 double &mlsOm, double &mlsOm2 ) ; 00288 public: 00289 int SolveRichardson(const PromVector_base* const, 00290 PromVector_base* const, int *its, 00291 const bool zeroguess = TRUE, 00292 const bool avoidnorms = FALSE, 00293 const PromPCSide = PROM_RIGHT ); 00294 int SolveBCGS(const PromVector_base* const, 00295 PromVector_base* const,int *its, 00296 const bool zeroguess = TRUE, 00297 const bool avoidnorms = FALSE, 00298 const PromPCSide = PROM_RIGHT ); 00299 int SolveCR(const PromVector_base*const, 00300 PromVector_base* const,int *its, 00301 const bool zeroguess = TRUE, 00302 const bool avoidnorms = FALSE, 00303 const PromPCSide = PROM_RIGHT ); 00304 int SolveBrez(const PromVector_base* const, 00305 PromVector_base* const, int *its, 00306 const bool zeroguess = TRUE, 00307 const bool avoidnorms = FALSE, 00308 const PromPCSide = PROM_RIGHT ); 00309 int SolveCheb(const PromVector_base*const, 00310 PromVector_base* const, int *its, 00311 const bool zeroguess = TRUE, 00312 const bool avoidnorms = FALSE, 00313 const PromPCSide = PROM_RIGHT ); 00314 int SolveCG(const PromVector_base* const, 00315 PromVector_base* const, int *its, 00316 const bool zeroguess = TRUE, 00317 const bool avoidnorms = FALSE, 00318 const PromPCSide = PROM_RIGHT ); 00319 int SolveTFQMR(const PromVector_base* const, 00320 PromVector_base* const, int *its, 00321 const bool zeroguess = TRUE, 00322 const bool avoidnorms = FALSE, 00323 const PromPCSide = PROM_LEFT ); 00324 int SolveSYMMLQ(const PromVector_base* const, 00325 PromVector_base* const, int *its, 00326 const bool zeroguess = TRUE, 00327 const bool avoidnorms = FALSE, 00328 const PromPCSide = PROM_LEFT ); 00329 int SolveMINRES(const PromVector_base* const, 00330 PromVector_base* const, int *its, 00331 const bool zeroguess = TRUE, 00332 const bool avoidnorms = FALSE, 00333 const PromPCSide = PROM_LEFT ); 00334 // PC 00335 int CreatePC( PromPCType type, const PromVector_base * const v ); 00336 int MakePC( PromGrid * const grid, 00337 int &nblocks, PromIS *&blockISs, bool &kspset, 00338 const PromVector *const xx, 00339 const char* const prestring = NULL, 00340 const char* const pcname = NULL ); 00341 int SetOperator( const PromMatrix_base * const A ); 00342 int SetNewOperator( const PromMatrix_base * const A ); 00343 int SetOperator_private( const PromMatrix_base * const A ); 00345 int PromPCApplyBAorAB( const PromMatrix_base * const Amat, 00346 PromPC_base * const pc, 00347 const PromPCSide pc_side, 00348 const PromVector_base * const xx, 00349 PromVector_base * const dest, 00350 PromVector_base * const zz )const; 00352 int Solve(const PromVector_base* const bb, 00353 PromVector_base* const xx, 00354 int *its, const bool zeroguess = TRUE, 00355 const bool avoidnorms =FALSE); 00356 protected: 00357 int TestSymmMA( const PromMatrix_base * const Amat, 00358 PromPC_base * const pc, PromVector_base *const Z, 00359 PromVector_base *const P, PromVector_base *const R )const; 00360 int Monitor( const int it, const double dp ){ 00361 if(verbose_>1 || (verbose_==1 && it%5==0)) { 00362 return PetscPrintf(A_->MPIComm(),"%d KSP Residual norm %12.7e\n",it,dp); 00363 } 00364 return 0; 00365 } 00366 int converged( double r0, double r ){ 00367 assert(this!=NULL); // !alloc 00368 if( r<r0*rtol_ || r<atol_ ) return 1; 00369 else if ( r > r0*dtol_) return -1; 00370 return 0; 00371 } 00372 public: 00373 // data 00374 static const char * const PromKSPStrings[PROMNUMKSP]; 00375 00376 PromPC_base *pc_; 00377 const PromMatrix_base *A_; 00378 PromVector_base **work_; 00379 const PromMap &map_; 00380 const PromMap *const p_map_; 00381 PromProjKKTSolver *projector_; 00382 protected: 00383 PromSolver *alias_; // !NULL if a smoother 00384 const PromOptions &options_; 00385 PromPerfMonitor &perf_mon_; 00386 float rtol_; // used by aliased 00387 float dtol_; // used by aliased 00388 float atol_; // used by aliased 00389 const bool alloced_; // is 'alias_' a dependent or am I? 00390 public: 00391 mutable short int verbose_; 00392 protected: 00393 short int nwork_; 00394 short int cg_restart_; 00395 struct prom_double_complex{ 00396 prom_double_complex(): e(NULL) 00397 #if defined(PROM_HAVE_COMPLEX) 00398 , c(NULL) 00399 #endif 00400 {}; 00401 double *e; 00402 #if defined(PROM_HAVE_COMPLEX) 00403 PromComplex *c; 00404 #endif 00405 }origin_, d_, hes_, dd_, ee_; 00406 00407 PromKSPType ksp_type_; 00408 short int nEig_; // CG, == nEig_ for GMRES 00409 short int maxits_; // used by aliased 00410 public: 00411 struct prom_history{ 00412 prom_history():helmholtz(0),iters(0),change(-511),level(0),istop(0),isbase(0){} 00413 unsigned int helmholtz : 1; 00414 unsigned int iters : 10; 00415 int change : 10; 00416 unsigned int level : 4; 00417 unsigned int istop : 1; 00418 unsigned int isbase : 1; 00419 unsigned int aux : 5; 00420 }history_; 00421 }; // PromSolver 00422 00423 /***************************************************************************** 00424 * G-S PromGSCont 00425 ****************************************************************************/ 00426 typedef enum BlktypeTAG{PROM_INF=0,PROM_INTER,PROM_INF_SUP,PROM_SUP}PromBlockType; 00427 typedef enum BlkListtypeTAG{PROM_INFL=0,PROM_INTER1,PROM_INTER1B,PROM_INF_SUP_WO_HIGHERIS,PROM_INF_SUPL,PROM_INF_SUP_WO_LOWERIS,PROM_INTER2,PROM_INTER2B,PROM_SUPL,PROM_NBLKLISTTYPES}PromBlockListType; 00428 00429 struct PromGSCont{ 00430 PromGSCont():arr_size_(0),list_(NULL){} 00431 int arr_size_; 00432 union{ 00433 PromList<int> *list_; 00434 int *array_; 00435 }; 00436 }; 00437 00438 /***************************************************************************** 00439 * PromProcXComm 00440 ****************************************************************************/ 00441 class PromProcXComm 00442 { 00443 public: 00444 PromProcXComm( int nn, int neq, int ppart, MPI_Comm cm, const int tg ) : 00445 nNds_(nn), 00446 nEqs_(neq), 00447 part_(ppart), 00448 req_(MPI_REQUEST_NULL), 00449 comm_(cm), 00450 tag_(tg), 00451 lids_(NULL), 00452 xarr_(NULL) 00453 { 00454 assert(nNds_>0); 00455 int isz=sizeof(int), dsz=sizeof(double), xx=nEqs_*(dsz/isz)+nNds_; 00456 PetscMalloc((xx+1)*isz,&lids_); assert(lids_ != NULL); 00457 lids_[xx] = lids_[xx-1] = -7; 00458 xarr_ = (double*)(lids_); 00459 xx = nEqs_; while(xx--) xarr_[xx] = 1.e30; 00460 lids_ += nEqs_ * (dsz/isz); 00461 InfBlocks_.arr_size_ = 0; 00462 InfBlocks_.list_ = new PromList<int>(10); 00463 } 00464 ~PromProcXComm() { 00465 assert(lids_[nNds_-1] != -7); assert(lids_[nNds_] == -7); 00466 PetscFree( xarr_ ); 00467 if(InfBlocks_.arr_size_ != 0) PetscFree(InfBlocks_.array_); 00468 } 00469 // data 00470 const int part_; 00471 const int nNds_; 00472 const int nEqs_; 00473 const int tag_; 00474 PromGSCont InfBlocks_; 00475 int *lids_; 00476 double *xarr_; 00477 MPI_Request req_; 00478 MPI_Comm comm_; 00479 }; 00480 00481 /***************************************************************************** 00482 * PromBlockXComm 00483 *****************************************************************************/ 00484 class PromBlockXComm 00485 { 00486 public: 00487 PromBlockXComm( int nn, int neq, int ppart, const int ggbid ) : 00488 nNds_(nn), 00489 nEqs_(neq), 00490 part_(ppart), 00491 #if defined(PROM_MEASURE_MPATH) 00492 path_(-9999), // debug 00493 #endif 00494 gbid_(ggbid) 00495 { 00496 assert(nNds_>0); int isz = sizeof(int); 00497 PetscMalloc((nNds_+1)*isz,&locBoundIDs_); assert(locBoundIDs_!=NULL); 00498 locBoundIDs_[nNds_] = locBoundIDs_[nNds_-1] = -7; 00499 } 00500 ~PromBlockXComm() { 00501 assert(locBoundIDs_[nNds_-1] != -7 && locBoundIDs_[nNds_] == -7); 00502 PetscFree(locBoundIDs_); 00503 } 00504 // data 00505 int *locBoundIDs_; 00506 const int nEqs_; 00507 const int nNds_; // keep all nodes in agg the same! 00508 const int part_; 00509 const int gbid_; 00510 #if defined(PROM_MEASURE_MPATH) 00511 int path_; // debug 00512 #endif 00513 }; 00514 00515 /************************************************************* 00516 * PromGSBlock 00517 **************************************************************/ 00518 class PromGSBlock 00519 { 00520 public: 00521 PromGSBlock( int iidd, int pproc ) : 00522 #if defined(PROM_MEASURE_MPATH) 00523 path_(-8888), 00524 #endif 00525 proc(pproc) 00526 { 00527 bits.done = FALSE; bits.gid = iidd; 00528 } 00529 virtual ~PromGSBlock(){} 00530 int getProc() const { assert(this!=NULL); return proc; } 00531 int isMarked() const { assert(this!=NULL); return bits.done; } 00532 void mark( int b = 1 ){ assert(this!=NULL); bits.done = (b==0) ? 0 : 1; } 00533 int getGBID() const { assert(this!=NULL); return bits.gid; } 00534 //void setGBID( int b ){ assert(this!=NULL && b>=0); bits.gid = b; } 00535 #if defined(PROM_MEASURE_MPATH) 00536 int path_; // debug 00537 #endif 00538 private: 00539 // data 00540 struct { 00541 unsigned int done : 1; 00542 unsigned int gid : 31; 00543 }bits; 00544 const int proc; 00545 }; 00546 00547 /************************************************************* 00548 * PromInfSupBlock 00549 **************************************************************/ 00550 class PromInfSupBlock : public PromGSBlock 00551 { 00552 public: 00553 PromInfSupBlock( int iidd, int pproc ) : 00554 PromGSBlock(iidd,pproc), 00555 high_ISBlocks_(10), 00556 low_ISBlocks_(10), 00557 ISBlockDeps_(NULL), 00558 lowProc_sendXComm_(10), 00559 highProc_sendXComm_(10), 00560 sendXComm_(NULL) 00561 { 00562 } 00563 ~PromInfSupBlock(){ 00564 PromBlockXComm *commx; 00565 while ( !highProc_sendXComm_.isEmpty() ) { 00566 commx = highProc_sendXComm_.removeHead(); 00567 delete commx; 00568 } 00569 while ( !lowProc_sendXComm_.isEmpty() ) { 00570 commx = lowProc_sendXComm_.removeHead(); 00571 delete commx; 00572 } 00573 } 00574 int collectSends( PromTable<PromList<PromBlockXComm*>*> &proc_list ); 00575 PromList<PromGSBlock*> *ISBlockDeps_; // pointer to low(high)Proc_ISBlocks 00576 PromList<PromBlockXComm*> *sendXComm_; // pointer to low(high)Proc_sendXComm 00577 // lists of dependents, etc. 00578 PromList<PromGSBlock*> high_ISBlocks_; // list of PromGSBlock pointers 00579 PromList<PromGSBlock*> low_ISBlocks_; // list of PromGSBlock pointers 00580 // send communicator if proc == myproc, else recv communicator 00581 PromList<PromBlockXComm*> lowProc_sendXComm_; // list of PromProcXComm; 00582 PromList<PromBlockXComm*> highProc_sendXComm_; // list of PromProcXComm; 00583 }; 00584 00585 template <class T> class PromAList : public PromList<T> 00586 { 00587 public: 00588 PromAList(int sz = 5) : PromList<T>(sz), count_(0) {} 00589 int count_; 00590 }; 00591 00592 /***************************************************************************** 00593 * PromGSNode 00594 ****************************************************************************/ 00595 class PromGSNode 00596 { 00597 public: 00598 PromGSNode(int lid, int gid, unsigned int n=0) : b(lid,n), gid_(gid) {} 00599 int global() const { return gid_; } 00600 int nb()const{return b.multi_;} 00601 // data 00602 struct T{ 00603 T(int l, unsigned int n):lid_(l),multi_(n){} 00604 const unsigned int multi_ : 4; 00605 const int lid_ : 28; 00606 }b; 00607 private: 00608 const int gid_; 00609 }; 00610 00611 /***************************************************************************** 00612 * PromPCGaussSeidel 00613 ****************************************************************************/ 00615 00618 class PromPCGaussSeidel : public PromPC_base 00619 { 00620 friend class PromSolver; 00621 public: 00622 PromPCGaussSeidel( const PromOptions &opt, 00623 PromPerfMonitor &perf ); 00624 int Create( const PromMatrix * const, const PromIS * const domainISs, 00625 const int nblocks ); 00626 ~PromPCGaussSeidel(); 00627 // methods 00628 private: 00629 int Create_private( const PromIS * const domainISs ); 00630 public: 00631 virtual int Apply_private(const PromVector_base* const bb, 00632 PromVector_base* const xx, const bool zerox=TRUE); 00633 virtual int SetUp_private(); 00634 virtual const char * getPCString() const { 00635 return PromPCStrings[PROMPCGAUSSSEIDEL + (symmetric_ ? 1 : 0) ]; 00636 } 00637 virtual const char * getSubPCString() const { 00638 if( kkt_ != NULL ) return kkt_->getPCString(); 00639 else if( pc_ != NULL ) return pc_->getPCString(); 00640 else return PromPCStrings[PROMPCNONE]; 00641 } 00642 virtual int getNumBlocks() const { 00643 if( pc_ != NULL ) return pc_->getNumBlocks(); 00644 else if( kkt_ != NULL ) return kkt_->getNumBlocks(); 00645 else return 0; 00646 } 00647 virtual PromPCType getType() const { 00648 assert( this != NULL ); 00649 return PROMPCGAUSSSEIDEL; 00650 } 00651 virtual int takeNonZeroGuess()const{ return pre_smooth_ ? 1 : -1; } 00652 virtual int SetNewOperator_private( const PromMatrix_base * const Ab ) { 00653 assert(this!=NULL); int ierr; assert(Ab!=NULL); 00654 const PromMatrix *AA = dynamic_cast<const PromMatrix*>(Ab); 00655 assert(AA!=NULL); 00656 if( pc_ != NULL ) { 00657 ierr = pc_->SetNewOperator( AA ); CHKERRQ(ierr); 00658 } 00659 if( A_ != NULL && AA->map_ != A_->map_ )SETERRQ(1,"A->map_ != map_"); 00660 A_ = AA; 00661 return 0; 00662 } 00663 virtual int SetOperator_private( const PromMatrix_base * const Ab ) { 00664 assert(this!=NULL); int ierr; assert(Ab!=NULL); 00665 const PromMatrix *AA = dynamic_cast<const PromMatrix*>(Ab); 00666 assert(AA!=NULL); 00667 if( pc_ != NULL ) { 00668 ierr = pc_->SetOperator( AA ); CHKERRQ(ierr); 00669 } 00670 if( A_ != NULL && AA->map_ != A_->map_ )SETERRQ(1,"A->map_ != map_"); 00671 A_ = AA; 00672 return 0; 00673 } 00674 private: 00675 int Apply_one( const PromVector* const bb, PromVector* const xx, 00676 bool zerog = FALSE, bool send_post = TRUE ); 00677 // setup methods 00678 virtual int Setup_private( PromList<int> ** const block_procs, 00679 PromTable<int> **const block_gid_ghostLid, 00680 PromList<int> **block_nodeLids, 00681 PromList<int> **block_boundLids, 00682 PromList<int> **block_adjacBLids, 00683 PromTable<PromTable<int>*>*const proc_boundNodes, 00684 const PromIS * const domainISs, 00685 PromList<int> *highlist, 00686 PromList<int> *lowlist ); 00687 int ColorProcs(PromList<int> *const highlist, PromList<int> *const lowlist); 00688 int SetBlockTypes( const PromList<int> *const * const block_procs, 00689 const PromTable<PromTable<int>*>*const proc_bNodes, 00690 PromBlockType block_type[], int lid_blockGID[], 00691 const PromIS * const domainISs ); 00692 int MakeLists( PromList<int> * const * const block_procs, 00693 const PromTable<int> * const * const block_ghosts, 00694 const PromBlockType block_type[], const int lid_blockGID[], 00695 const PromList<int> * const * const block_nodeLids, 00696 const PromList<int> * const * const block_boundNodes, 00697 const PromList<int> * const * const block_adajacBlids, 00698 const PromIS*const domainISs); 00699 int SetupNodalGS(); 00700 // solver 00701 int DoBlocks( const int nb, const int blocks[], 00702 const PromVector* const bb, 00703 PromVector *const xx ) const; 00704 #if defined(PROM_USE_PETSC) 00705 int DoBlocks_nodal_1(const int nb, const int block[], 00706 const PromVector* const b, 00707 PromVector* const x) const; 00708 int DoBlocks_nodal_3(const int nb, const int block[], 00709 const PromVector* const bb, 00710 PromVector*const xx) const; 00711 int DoBlocks_nodal_6(const int nb, const int block[], 00712 const PromVector* const bb, 00713 PromVector* const xx)const; 00714 int DoBlocks_nodal(const int nb, const int block[], 00715 const PromVector* const bb,PromVector*const xx, 00716 const int ndf ) const; 00717 int DoBlocks_schwarz_1(const int nb, const int block[], 00718 const PromVector* const bb, 00719 PromVector* const xx, 00720 const PromPC * const pc )const; 00721 int DoBlocks_schwarz_3(const int nb, const int block[], 00722 const PromVector* const bb, 00723 PromVector* const xx, 00724 const PromPC * const pc )const; 00725 int DoBlocks_schwarz_6(const int nb, const int block[], 00726 const PromVector* const bb, 00727 PromVector* const xx, 00728 const PromPC * const pc )const; 00729 int DoBlocks_schwarz(const int nb, const int block[], 00730 const PromVector*const bb,PromVector*const xx, 00731 const PromPC * const pc, 00732 const int ndf ) const; 00733 #else 00734 // petra specific G-S kernals 00735 #endif 00736 00737 // utils 00738 int GetProc( const int gid, int &proc ) const; 00739 int PostRecvs( const PromAList<PromProcXComm*> *commxList ) const; 00740 int SendLocalx(const PromList<PromProcXComm*> *commxList, 00741 const PromVector * const xx)const; 00742 int RecvDoGhostx( PromAList<PromProcXComm*> *const commxList, 00743 const PromVector* const bb = NULL, 00744 PromVector *const xx = NULL ) const; 00745 int WaitForList( const PromList<PromProcXComm*> *commxList ) const; 00746 int SendISblocks( PromTable<PromList<PromBlockXComm*>*> *proc_list, 00747 PromList<MPI_Request*> *smessages, 00748 const PromVector *const xx, const int tag ) const; 00749 int RecvBlockGhostx( double ghostx[], const int nghostEqs, 00750 const int tag, const bool iprobe, 00751 int &ndone ) const; 00752 public: 00753 int ReorderNodes( double wieghts[], int len ); 00754 private: 00755 static int BfsReorder( PromList<int> *blidList, 00756 const PromList<int> *const * const block_adjacBLids, 00757 PromList<int> *list, PromList<int> *nextRoots, 00758 bool buff1[], bool buff2[], const bool print = FALSE); 00759 // data 00760 int maxnGSBlocks_; // used for "globalizing" block IDs 00761 // 'DoBlocks(...)' data 00762 public: 00763 int nGSBlocks_; 00764 int bs_; 00765 bool pre_smooth_; // hack for alternating up/down smoothin 00766 bool symmetric_; 00767 MPI_Comm comm_; 00768 //int *ghostLid_locEQ_; 00769 //int *lid_geq_; 00770 double *NodalFactors_; 00771 double *ghostx_; 00772 int *NodalFactorPointers_; 00773 PromPC *pc_; 00774 const PromMatrix *A_; 00775 private: 00776 PromTable<PromGSNode*> gbid_ghost_; // gid --> local ghost ID (ie, 0:nghost-1) 00777 PromGSCont blocksLists_[PROM_NBLKLISTTYPES];// local block IDs each type 00778 PromList<PromProcXComm*> lowProc_send_; // list of X_COMM 00779 PromAList<PromProcXComm*> lowProc_recv_; // list of X_COMM 00780 PromAList<PromProcXComm*> highProc_recv_; // list of X_COMM 00781 PromList<PromProcXComm*> highProc_send_; // list of X_COMM 00782 PromTable<PromGSBlock*> gbid_ISBlock_; // list of IS_BLOCKS 00783 PromTable<int> proc_color_; // proc --> color 00784 }; // PromPCGaussSeidel 00785 00786 /***************************************************************************** 00787 * PromMG 00788 ****************************************************************************/ 00790 00793 class PromMG 00794 { 00795 friend class PromPCMG; 00796 friend class Prometheus_LinSysCore; 00797 public: 00798 PromMG( PromPerfMonitor &perf ); 00799 ~PromMG(); 00800 // methods 00801 int SetCycles(int nn){ assert(this!=NULL); cycles_ = nn; return 0; } 00802 int getCycles()const{ assert(this!=NULL); return cycles_; } 00803 int Create( PromGrid * const grid0, int npre, int npost, 00804 const PromMap *const p_map = NULL, 00805 const bool use_iter_on_top = FALSE ); 00806 int SetOperator( const PromMatrix_base * const K ) { 00807 A_ = K; int ierr; 00808 if( presmooth_ != NULL ) ierr = presmooth_->SetOperator( K ); 00809 if( postsmooth_ != NULL ) ierr = postsmooth_->SetOperator( K ); 00810 CHKERRQ(ierr); 00811 return 0; 00812 } 00813 int SetNewOperator( const PromMatrix_base * const K ) { 00814 A_ = K; int ierr; 00815 if( presmooth_ != NULL ) ierr = presmooth_->SetNewOperator(K); 00816 if( postsmooth_ != NULL ) ierr = postsmooth_->SetNewOperator(K); 00817 CHKERRQ(ierr); 00818 return 0; 00819 } 00820 int SetNewPVec( const PromMap &map ); 00821 // data 00822 PromSolver *presmooth_; 00823 PromSolver *postsmooth_; 00824 PromVector_base *b_; 00825 PromVector_base *x_; 00826 PromVector_base *r_; 00827 const PromMatrix_base *interpolate_; // with next (up) grid 00828 const PromMatrix_base *restrct_; 00829 const PromMatrix_base *A_; 00830 PromMG *prev_; 00831 PromMG *next_; 00832 int factor_; 00833 protected: 00834 PromPerfMonitor &perf_mon_; 00835 int cycles_; 00836 }; // PromMG 00837 00838 /***************************************************************************** 00839 * PromPCMG 00840 ****************************************************************************/ 00842 00845 class PromPCMG : public PromPC_base 00846 { 00847 public: 00848 PromPCMG( const PromOptions &opts, PromPerfMonitor &perf ) : 00849 PromPC_base( opts, perf ), mgstack_(NULL) {} 00850 virtual ~PromPCMG(); 00851 // methods 00852 PromPCType getType() const { 00853 assert( this != NULL ); 00854 return PROMPCMG; 00855 } 00856 // applies PC part in shell PC 00857 int Apply_private( const PromVector_base* const bb, 00858 PromVector_base* const xx ,const bool zerox = TRUE ); 00859 virtual int SetUp_private(); 00860 const char * getPCString() const { 00861 return PromPCStrings[PROMPCMG]; 00862 } 00863 const char * getSubPCString() const { 00864 if( kkt_ != NULL ) return kkt_->getPCString(); 00865 return PromPCStrings[PROMPCNONE]; 00866 } 00867 virtual int SetOperator_private( const PromMatrix_base * const K ) { 00868 assert(this!=NULL); int ierr; 00869 if( mgstack_ != NULL ) { 00870 ierr = mgstack_->SetOperator( K ); CHKERRQ(ierr); 00871 } 00872 return 0; 00873 } 00874 virtual int SetNewOperator_private( const PromMatrix_base * const K ) { 00875 assert(this!=NULL); int ierr; 00876 if( mgstack_ != NULL ) { 00877 ierr = mgstack_->SetNewOperator( K ); CHKERRQ(ierr); 00878 } 00879 return 0; 00880 } 00881 int PrintLevelInfo() const; 00882 int Create( PromGrid * const grid0, const PromMap *p_map = NULL ); 00883 private: 00884 int ACycle( PromMG * const mg, const bool zerox = TRUE ); 00885 int MCycle( PromMG * const mg, int level, const bool zerox = TRUE ); 00886 int FCycle( PromMG * const mg, const bool zerox = TRUE ); 00887 public: 00888 // data 00889 PromMG* mgstack_; 00890 PromMGAlgo algo_; 00891 }; // PromPCMG 00892 00893 /***************************************************************************** 00894 * PromPCNodalASM 00895 ****************************************************************************/ 00897 00900 class PromPCNodalASM : public PromPC_base 00901 { 00902 public: 00904 PromPCNodalASM( const PromOptions &opts, PromPerfMonitor &perf ) 00905 : PromPC_base(opts,perf),invBlocks_(NULL),A_(NULL),r_(NULL),work_(NULL){} 00907 ~PromPCNodalASM(){ 00908 if( invBlocks_ ) PetscFree( invBlocks_ ); 00909 if( r_ ) delete r_; 00910 if( work_ ) delete work_; //PromVector::Delete( work_ ); 00911 } 00913 virtual PromPCType getType() const { 00914 assert( this != NULL ); 00915 return PROMPCNODALASM; 00916 } 00917 virtual const char * getPCString() const { 00918 return PromPCStrings[PROMPCNODALASM]; 00919 } 00920 virtual const char * getSubPCString() const { 00921 if( kkt_ != NULL ) return kkt_->getPCString(); 00922 return PromPCStrings[PROMPCNONE]; 00923 } 00925 virtual int SetOperator_private( const PromMatrix_base * const AA_b ) { 00926 assert(this!=NULL); 00927 if( AA_b == NULL ) SETERRQ(1,"AA_b == NULL"); 00928 const PromMatrix * const AA = dynamic_cast<const PromMatrix*>(AA_b); 00929 if(AA==NULL)SETERRQ(1,"A not a PromMatrix"); 00930 A_ = AA; 00931 return 0; 00932 } 00934 virtual int SetNewOperator_private( const PromMatrix_base * const AA_b ) { 00935 assert(this!=NULL); 00936 if( AA_b == NULL ) SETERRQ(1,"AA_b == NULL"); 00937 const PromMatrix * const AA = dynamic_cast<const PromMatrix*>(AA_b); 00938 if(AA==NULL)SETERRQ(1,"A not a PromMatrix"); 00939 A_ = AA; 00940 00941 return 0; 00942 } 00944 int Create(); 00945 int SetUp_private(); 00946 virtual int GetOperator( const PromMatrix **K_out, 00947 PromVector **r_out, const PromVector_base*const XX) 00948 { 00949 if( work_ == NULL ) { 00950 work_ = new PromVector( XX ); assert(r_==NULL); 00951 r_ = new PromVector( XX ); 00952 } 00953 *K_out = A_; *r_out = work_; 00954 return 0; 00955 } 00957 virtual int Apply_private(const PromVector_base* const bb, 00958 PromVector_base* const xx,const bool zerox = TRUE); 00959 private: 00960 int Apply_one( const PromVector* const bb, PromVector* const xx ); 00961 // data 00962 double *invBlocks_; 00963 PromVector *r_; 00964 PromVector *work_; 00965 const PromMatrix *A_; 00966 }; // PromPCNodalASM 00967 00968 #endif //__PROMETHEUS_H__