120 this->GetOStream(
Warnings0) <<
"MueLu::BreaessSarazinSmoother::Setup(): Setup() has already been called";
123 A_ = Factory::Get<RCP<Matrix> > (currentLevel,
"A");
124 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
126 "MueLu::BraessSarazinSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
129 rangeMapExtractor_ = bA->getRangeMapExtractor();
130 domainMapExtractor_ = bA->getDomainMapExtractor();
133 A00_ = bA->getMatrix(0,0);
134 A01_ = bA->getMatrix(0,1);
135 A10_ = bA->getMatrix(1,0);
136 A11_ = bA->getMatrix(1,1);
139 SC omega = pL.get<SC>(
"Damping factor");
143 RCP<Vector> diagFVector = VectorFactory::Build(A00_->getRowMap());
144 if (pL.get<
bool>(
"lumping") ==
false) {
145 A00_->getLocalDiagCopy(*diagFVector);
149 diagFVector->scale(omega);
153 RCP<BlockedVector> bD = Teuchos::rcp_dynamic_cast<BlockedVector>(D_);
154 if(bD.is_null() ==
false && bD->getBlockedMap()->getNumMaps() == 1) {
155 RCP<Vector> nestedVec = bD->getMultiVector(0,bD->getBlockedMap()->getThyraMode())->getVectorNonConst(0);
163 smoo_ = currentLevel.
Get<RCP<SmootherBase> >(
"PreSmoother", FactManager_->GetFactory(
"Smoother").get());
164 S_ = currentLevel.
Get<RCP<Matrix> > (
"A", FactManager_->GetFactory(
"A").get());
173 "MueLu::BraessSarazinSmoother::Apply(): Setup() has not been called");
175 RCP<MultiVector> rcpX = rcpFromRef(X);
176 RCP<const MultiVector> rcpB = rcpFromRef(B);
179 bool bCopyResultX =
false;
180 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
182 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
183 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
185 if(bX.is_null() ==
true) {
186 RCP<MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
191 if(bB.is_null() ==
true) {
192 RCP<const MultiVector> test = Teuchos::rcp(
new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
197 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
198 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
201 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
202 if(rbA.is_null() ==
false) {
204 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
207 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
209 Teuchos::RCP<MultiVector> test =
210 buildReorderedBlockedMultiVector(brm, bX);
213 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
215 Teuchos::RCP<const MultiVector> test =
216 buildReorderedBlockedMultiVector(brm, bB);
224 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
225 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
227 RCP<MultiVector> deltaX = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
228 RCP<BlockedMultiVector> bdeltaX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(deltaX);
229 RCP<MultiVector> deltaX0 = bdeltaX->getMultiVector(0,bDomainThyraMode);
230 RCP<MultiVector> deltaX1 = bdeltaX->getMultiVector(1,bDomainThyraMode);
232 RCP<MultiVector> Rtmp = rangeMapExtractor_->getVector(1, rcpB->getNumVectors(), bRangeThyraMode);
235 typedef Teuchos::ScalarTraits<SC> STS;
236 SC one = STS::one(), zero = STS::zero();
240 LO nSweeps = pL.get<LO>(
"Sweeps");
243 if (InitialGuessIsZero) {
244 R = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
245 R->update(one, *rcpB, zero);
251 RCP<Vector> diagSVector = VectorFactory::Build(S_->getRowMap());
252 S_->getLocalDiagCopy(*diagSVector);
254 for (LO run = 0; run < nSweeps; ++run) {
258 RCP<MultiVector> R0 = rangeMapExtractor_ ->ExtractVector(R, 0, bRangeThyraMode);
259 RCP<MultiVector> R1 = rangeMapExtractor_ ->ExtractVector(R, 1, bRangeThyraMode);
262 deltaX0->putScalar(zero);
263 deltaX0->elementWiseMultiply(one, *D_, *R0, zero);
264 A10_->apply(*deltaX0, *Rtmp);
265 Rtmp->update(one, *R1, -one);
267 if (!pL.get<
bool>(
"q2q1 mode")) {
268 deltaX1->putScalar(zero);
271 if(Teuchos::rcp_dynamic_cast<BlockedVector>(diagSVector) == Teuchos::null) {
272 ArrayRCP<SC> Sdiag = diagSVector->getDataNonConst(0);
273 ArrayRCP<SC> deltaX1data = deltaX1->getDataNonConst(0);
274 ArrayRCP<SC> Rtmpdata = Rtmp->getDataNonConst(0);
275 for (GO row = 0; row < deltaX1data.size(); row++)
276 deltaX1data[row] = Teuchos::as<SC>(1.1)*Rtmpdata[row] / Sdiag[row];
282 smoo_->Apply(*deltaX1,*Rtmp);
285 deltaX0->putScalar(zero);
286 A01_->apply(*deltaX1, *deltaX0);
287 deltaX0->update(one, *R0, -one);
289 deltaX0->elementWiseMultiply(one, *D_, *R0, zero);
291 RCP<MultiVector> X0 = domainMapExtractor_->ExtractVector(rcpX, 0, bDomainThyraMode);
292 RCP<MultiVector> X1 = domainMapExtractor_->ExtractVector(rcpX, 1, bDomainThyraMode);
295 X0->update(one, *deltaX0, one);
296 X1->update(one, *deltaX1, one);
298 domainMapExtractor_->InsertVector(X0, 0, rcpX, bDomainThyraMode);
299 domainMapExtractor_->InsertVector(X1, 1, rcpX, bDomainThyraMode);
301 if (run < nSweeps-1) {
307 if (bCopyResultX ==
true) {
308 RCP<MultiVector> Xmerged = bX->Merge();
309 X.update(one, *Xmerged, zero);