94 TEUCHOS_TEST_FOR_EXCEPTION(pos < 0,
Exceptions::RuntimeError,
"MueLu::BlockedGaussSeidelSmoother::AddFactoryManager: parameter \'pos\' must not be negative! error.");
96 size_t myPos = Teuchos::as<size_t>(pos);
98 if (myPos < FactManager_.size()) {
100 FactManager_.at(myPos) = FactManager;
101 }
else if(myPos == FactManager_.size()) {
103 FactManager_.push_back(FactManager);
105 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
106 *out <<
"Warning: cannot add new FactoryManager at proper position " << pos <<
". The FactoryManager is just appended to the end. Check this!" << std::endl;
109 FactManager_.push_back(FactManager);
135 RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
137 FactoryMonitor m(*
this,
"Setup blocked Gauss-Seidel Smoother", currentLevel);
141 A_ = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
142 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
143 TEUCHOS_TEST_FOR_EXCEPTION(bA==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: input matrix A is not of type BlockedCrsMatrix! error.");
146 TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != FactManager_.size(),
Exceptions::RuntimeError,
"MueLu::BlockedGaussSeidelSmoother::Setup: number of block rows of A is " << bA->Rows() <<
" and does not match number of SubFactoryManagers " << FactManager_.size() <<
". error.");
147 TEUCHOS_TEST_FOR_EXCEPTION(bA->Cols() != FactManager_.size(),
Exceptions::RuntimeError,
"MueLu::BlockedGaussSeidelSmoother::Setup: number of block cols of A is " << bA->Cols() <<
" and does not match number of SubFactoryManagers " << FactManager_.size() <<
". error.");
150 rangeMapExtractor_ = bA->getRangeMapExtractor();
151 domainMapExtractor_ = bA->getDomainMapExtractor();
154 std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
155 for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
159 RCP<const SmootherBase> Smoo = currentLevel.
Get< RCP<SmootherBase> >(
"PreSmoother",(*it)->GetFactory(
"Smoother").get());
160 Inverse_.push_back(Smoo);
163 RCP<Matrix> Aii = currentLevel.
Get< RCP<Matrix> >(
"A",(*it)->GetFactory(
"A").get());
164 bIsBlockedOperator_.push_back(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Aii)!=Teuchos::null);
177 RCP<MultiVector> rcpDebugX = Teuchos::rcpFromRef(X);
178 RCP<const MultiVector> rcpDebugB = Teuchos::rcpFromRef(B);
179 RCP<BlockedMultiVector> rcpBDebugX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpDebugX);
180 RCP<const BlockedMultiVector> rcpBDebugB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpDebugB);
182 if(rcpBDebugB.is_null() ==
false) {
189 if(rcpBDebugX.is_null() ==
false) {
197 SC zero = Teuchos::ScalarTraits<SC>::zero(), one = Teuchos::ScalarTraits<SC>::one();
200 RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
201 RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
204 bool bCopyResultX =
false;
205 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
207 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
208 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
210 if(bX.is_null() ==
true) {
211 RCP<MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
216 if(bB.is_null() ==
true) {
217 RCP<const MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
222 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
223 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
226 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
227 if(rbA.is_null() ==
false) {
229 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
232 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
234 Teuchos::RCP<MultiVector> test =
235 buildReorderedBlockedMultiVector(brm, bX);
238 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
240 Teuchos::RCP<const MultiVector> test =
241 buildReorderedBlockedMultiVector(brm, bB);
286 RCP<MultiVector> residual = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
287 RCP<MultiVector> tempres = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
296 if( InitialGuessIsZero==
true )
297 rcpX->putScalar(Teuchos::ScalarTraits<Scalar>::zero());
304 for(
size_t i = 0; i<Inverse_.size(); i++) {
308 residual->update(1.0,*rcpB,0.0);
309 if(InitialGuessIsZero ==
false || i > 0 || run > 0)
310 bA->bgs_apply(*rcpX, *residual, i, Teuchos::NO_TRANS, -1.0, 1.0);
313 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
314 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
315 Teuchos::RCP<MultiVector> Xi = domainMapExtractor_->ExtractVector(rcpX, i, bDomainThyraMode);
316 Teuchos::RCP<MultiVector> ri = rangeMapExtractor_->ExtractVector(residual, i, bRangeThyraMode);
317 Teuchos::RCP<MultiVector> tXi = domainMapExtractor_->getVector(i, X.getNumVectors(), bDomainThyraMode);
320 Inverse_.at(i)->Apply(*tXi, *ri,
false);
323 Xi->update(omega,*tXi,1.0);
328 if (bCopyResultX ==
true) {
329 RCP<MultiVector> Xmerged = bX->Merge();
330 X.update(one, *Xmerged, zero);