150 FactoryMonitor m(*
this,
"Setup blocked SIMPLE Smoother", currentLevel);
153 this->GetOStream(
Warnings0) <<
"MueLu::SimpleSmoother::Setup(): Setup() has already been called";
156 A_ = Factory::Get<RCP<Matrix> > (currentLevel,
"A");
158 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
159 TEUCHOS_TEST_FOR_EXCEPTION(bA == Teuchos::null,
Exceptions::BadCast,
"MueLu::SimpleSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
162 rangeMapExtractor_ = bA->getRangeMapExtractor();
163 domainMapExtractor_ = bA->getDomainMapExtractor();
166 F_ = bA->getMatrix(0, 0);
167 G_ = bA->getMatrix(0, 1);
168 D_ = bA->getMatrix(1, 0);
169 Z_ = bA->getMatrix(1, 1);
172 bool bSIMPLEC = pL.get<
bool>(
"UseSIMPLEC");
176 RCP<Vector> diagFVector = VectorFactory::Build(F_->getRowMap());
178 F_->getLocalDiagCopy(*diagFVector);
199 RCP<BlockedVector> bdiagFinv = Teuchos::rcp_dynamic_cast<BlockedVector>(diagFinv_);
200 if(bdiagFinv.is_null() ==
false && bdiagFinv->getBlockedMap()->getNumMaps() == 1) {
201 RCP<Vector> nestedVec = bdiagFinv->getMultiVector(0,bdiagFinv->getBlockedMap()->getThyraMode())->getVectorNonConst(0);
202 diagFinv_.swap(nestedVec);
208 RCP<const FactoryManagerBase> velpredictFactManager = FactManager_.at(0);
210 velPredictSmoo_ = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother",velpredictFactManager->GetFactory(
"Smoother").get());
213 RCP<const FactoryManagerBase> schurFactManager = FactManager_.at(1);
215 schurCompSmoo_ = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother", schurFactManager->GetFactory(
"Smoother").get());
225 "MueLu::SimpleSmoother::Apply(): Setup() has not been called");
228 RCP<MultiVector> rcpDebugX = Teuchos::rcpFromRef(X);
229 RCP<const MultiVector> rcpDebugB = Teuchos::rcpFromRef(B);
230 RCP<BlockedMultiVector> rcpBDebugX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpDebugX);
231 RCP<const BlockedMultiVector> rcpBDebugB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpDebugB);
232 if(rcpBDebugB.is_null() ==
false) {
237 if(rcpBDebugX.is_null() ==
false) {
244 const SC zero = Teuchos::ScalarTraits<SC>::zero();
245 const SC one = Teuchos::ScalarTraits<SC>::one();
253 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
254 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
259 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
260 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
263 bool bCopyResultX =
false;
264 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
266 "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix");
267 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
268 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
270 if(bX.is_null() ==
true) {
271 RCP<MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
276 if(bB.is_null() ==
true) {
277 RCP<const MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
282 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
283 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
286 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA =
287 Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
288 if(rbA.is_null() ==
false) {
290 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
293 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
295 Teuchos::RCP<MultiVector> test = buildReorderedBlockedMultiVector(brm, bX);
298 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
300 Teuchos::RCP<const MultiVector> test = buildReorderedBlockedMultiVector(brm, bB);
309 RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
310 RCP<BlockedMultiVector> bresidual = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(residual);
311 RCP<MultiVector> r1 = bresidual->getMultiVector(0,bRangeThyraMode);
312 RCP<MultiVector> r2 = bresidual->getMultiVector(1,bRangeThyraMode);
315 RCP<MultiVector> xtilde = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
316 RCP<BlockedMultiVector> bxtilde = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xtilde);
317 RCP<MultiVector> xtilde1 = bxtilde->getMultiVector(0, bDomainThyraMode);
318 RCP<MultiVector> xtilde2 = bxtilde->getMultiVector(1, bDomainThyraMode);
321 RCP<MultiVector> xhat = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
322 RCP<BlockedMultiVector> bxhat = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(xhat);
323 RCP<MultiVector> xhat1 = bxhat->getMultiVector(0, bDomainThyraMode);
324 RCP<MultiVector> xhat2 = bxhat->getMultiVector(1, bDomainThyraMode);
330 residual->update(one, *rcpB, zero);
331 if(InitialGuessIsZero ==
false || run > 0)
332 A_->apply(*rcpX, *residual, Teuchos::NO_TRANS, -one, one);
336 xtilde1->putScalar(zero);
337 xtilde2->putScalar(zero);
338 velPredictSmoo_->Apply(*xtilde1, *r1);
342 RCP<MultiVector> schurCompRHS = rangeMapExtractor_->getVector(1, rcpB->getNumVectors(), bRangeThyraMode);
343 D_->apply(*xtilde1, *schurCompRHS);
345 schurCompRHS->update(one, *r2, -one);
349 schurCompSmoo_->Apply(*xtilde2, *schurCompRHS);
353 xhat2->update(omega, *xtilde2, zero);
356 RCP<MultiVector> xhat1_temp = domainMapExtractor_->getVector(0, rcpX->getNumVectors(), bDomainThyraMode);
357 G_->apply(*xhat2, *xhat1_temp);
359 xhat1->elementWiseMultiply(one, *diagFinv_, *xhat1_temp, zero);
360 xhat1->update(one, *xtilde1, -one);
362 rcpX->update(one, *bxhat, one);
365 if (bCopyResultX ==
true) {
366 RCP<const MultiVector> Xmerged = bX->Merge();
367 X.update(one, *Xmerged, zero);