46 #ifndef ANASAZI_MINRES_HPP 47 #define ANASAZI_MINRES_HPP 50 #include "Teuchos_TimeMonitor.hpp" 55 template <
class Scalar,
class MV,
class OP>
56 class PseudoBlockMinres
65 PseudoBlockMinres(Teuchos::RCP<OP> A, Teuchos::RCP<OP> Prec = Teuchos::null);
68 void setTol(
const std::vector<Scalar>& tolerances) { tolerances_ = tolerances; };
69 void setMaxIter(
const int maxIt) { maxIt_ = maxIt; };
72 void setSol(Teuchos::RCP<MV> X) { X_ = X; };
73 void setRHS(Teuchos::RCP<const MV> B) { B_ = B; };
80 Teuchos::RCP<OP> Prec_;
82 Teuchos::RCP<const MV> B_;
83 std::vector<Scalar> tolerances_;
86 void symOrtho(Scalar a, Scalar b, Scalar *c, Scalar *s, Scalar *r);
88 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 89 Teuchos::RCP<Teuchos::Time> AddTime_, ApplyOpTime_, ApplyPrecTime_, AssignTime_, DotTime_, LockTime_, NormTime_, ScaleTime_, TotalTime_;
95 template <
class Scalar,
class MV,
class OP>
96 PseudoBlockMinres<Scalar,MV,OP>::PseudoBlockMinres(Teuchos::RCP<OP> A, Teuchos::RCP<OP> Prec) :
97 ONE(Teuchos::ScalarTraits<Scalar>::one()),
98 ZERO(Teuchos::ScalarTraits<Scalar>::zero())
100 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 101 AddTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Add Multivectors");
102 ApplyOpTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Apply Operator");
103 ApplyPrecTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Apply Preconditioner");
104 AssignTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Assignment (no locking)");
105 DotTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Compute Dot Product");
106 LockTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Lock Converged Vectors");
107 NormTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Compute Norm");
108 ScaleTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Scale Multivector");
109 TotalTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Total Time");
119 template <
class Scalar,
class MV,
class OP>
120 void PseudoBlockMinres<Scalar,MV,OP>::solve()
122 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 123 Teuchos::TimeMonitor outertimer( *TotalTime_ );
127 int ncols = MVT::GetNumberVecs(*B_);
131 std::vector<Scalar> alpha(ncols), beta, beta1(ncols), phibar, oldBeta(ncols,ZERO), epsln(ncols,ZERO), cs(ncols,-ONE), sn(ncols,ZERO), dbar(ncols,ZERO), oldeps(ncols), delta(ncols), gbar(ncols), phi(ncols), gamma(ncols), tmpvec(ncols);
132 std::vector<int> indicesToRemove, newlyConverged, unconvergedIndices(ncols);
133 Teuchos::RCP<MV> V, Y, R1, R2, W, W1, W2, locX, scaleHelper, swapHelper;
136 V = MVT::Clone(*B_, ncols);
137 Y = MVT::Clone(*B_, ncols);
138 R1 = MVT::Clone(*B_, ncols);
139 R2 = MVT::Clone(*B_, ncols);
140 W = MVT::Clone(*B_, ncols);
141 W1 = MVT::Clone(*B_, ncols);
142 W2 = MVT::Clone(*B_, ncols);
143 scaleHelper = MVT::Clone(*B_, ncols);
146 indicesToRemove.reserve(ncols);
147 newlyConverged.reserve(ncols);
150 for(
int i=0; i<ncols; i++)
151 unconvergedIndices[i] = i;
155 locX = MVT::CloneCopy(*X_);
160 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 161 Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
163 OPT::Apply(*A_,*locX,*R1);
166 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 167 Teuchos::TimeMonitor lcltimer( *AddTime_ );
169 MVT::MvAddMv(ONE, *B_, -ONE, *R1, *R1);
174 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 175 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
177 MVT::Assign(*R1,*R2);
185 if(Prec_ != Teuchos::null)
187 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 188 Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
191 OPT::Apply(*Prec_,*R1,*Y);
195 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 196 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
203 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 204 Teuchos::TimeMonitor lcltimer( *DotTime_ );
206 MVT::MvDot(*Y,*R1, beta1);
208 for(
size_t i=0; i<beta1.size(); i++)
209 beta1[i] = sqrt(beta1[i]);
220 for(
int iter=1; iter <= maxIt_; iter++)
223 indicesToRemove.clear();
224 for(
int i=0; i<ncols; i++)
226 Scalar relres = phibar[i]/beta1[i];
231 if(relres < tolerances_[i])
233 indicesToRemove.push_back(i);
238 newNumConverged = indicesToRemove.size();
239 if(newNumConverged > 0)
241 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 242 Teuchos::TimeMonitor lcltimer( *LockTime_ );
246 newlyConverged.resize(newNumConverged);
247 for(
int i=0; i<newNumConverged; i++)
248 newlyConverged[i] = unconvergedIndices[indicesToRemove[i]];
250 Teuchos::RCP<MV> helperLocX = MVT::CloneCopy(*locX,indicesToRemove);
252 MVT::SetBlock(*helperLocX,newlyConverged,*X_);
255 if(newNumConverged == ncols)
259 std::vector<int> helperVec(ncols - newNumConverged);
261 std::set_difference(unconvergedIndices.begin(), unconvergedIndices.end(), newlyConverged.begin(), newlyConverged.end(), helperVec.begin());
262 unconvergedIndices = helperVec;
265 std::vector<int> thingsToKeep(ncols - newNumConverged);
266 helperVec.resize(ncols);
267 for(
int i=0; i<ncols; i++)
269 ncols = ncols - newNumConverged;
271 std::set_difference(helperVec.begin(), helperVec.end(), indicesToRemove.begin(), indicesToRemove.end(), thingsToKeep.begin());
274 Teuchos::RCP<MV> helperMV;
275 helperMV = MVT::CloneCopy(*V,thingsToKeep);
277 helperMV = MVT::CloneCopy(*Y,thingsToKeep);
279 helperMV = MVT::CloneCopy(*R1,thingsToKeep);
281 helperMV = MVT::CloneCopy(*R2,thingsToKeep);
283 helperMV = MVT::CloneCopy(*W,thingsToKeep);
285 helperMV = MVT::CloneCopy(*W1,thingsToKeep);
287 helperMV = MVT::CloneCopy(*W2,thingsToKeep);
289 helperMV = MVT::CloneCopy(*locX,thingsToKeep);
291 scaleHelper = MVT::Clone(*V,ncols);
295 oldeps.resize(ncols);
300 tmpvec.resize(ncols);
301 std::vector<Scalar> helperVecS(ncols);
302 for(
int i=0; i<ncols; i++)
303 helperVecS[i] = beta[thingsToKeep[i]];
305 for(
int i=0; i<ncols; i++)
306 helperVecS[i] = beta1[thingsToKeep[i]];
308 for(
int i=0; i<ncols; i++)
309 helperVecS[i] = phibar[thingsToKeep[i]];
311 for(
int i=0; i<ncols; i++)
312 helperVecS[i] = oldBeta[thingsToKeep[i]];
313 oldBeta = helperVecS;
314 for(
int i=0; i<ncols; i++)
315 helperVecS[i] = epsln[thingsToKeep[i]];
317 for(
int i=0; i<ncols; i++)
318 helperVecS[i] = cs[thingsToKeep[i]];
320 for(
int i=0; i<ncols; i++)
321 helperVecS[i] = sn[thingsToKeep[i]];
323 for(
int i=0; i<ncols; i++)
324 helperVecS[i] = dbar[thingsToKeep[i]];
328 A_->removeIndices(indicesToRemove);
334 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 335 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
339 for(
int i=0; i<ncols; i++)
340 tmpvec[i] = ONE/beta[i];
342 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 343 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
345 MVT::MvScale (*V, tmpvec);
351 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 352 Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
354 OPT::Apply(*A_, *V, *Y);
360 for(
int i=0; i<ncols; i++)
361 tmpvec[i] = beta[i]/oldBeta[i];
363 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 364 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
366 MVT::Assign(*R1,*scaleHelper);
369 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 370 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
372 MVT::MvScale(*scaleHelper,tmpvec);
375 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 376 Teuchos::TimeMonitor lcltimer( *AddTime_ );
378 MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
384 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 385 Teuchos::TimeMonitor lcltimer( *DotTime_ );
387 MVT::MvDot(*V, *Y, alpha);
391 for(
int i=0; i<ncols; i++)
392 tmpvec[i] = alpha[i]/beta[i];
394 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 395 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
397 MVT::Assign(*R2,*scaleHelper);
400 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 401 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
403 MVT::MvScale(*scaleHelper,tmpvec);
406 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 407 Teuchos::TimeMonitor lcltimer( *AddTime_ );
409 MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
420 if(Prec_ != Teuchos::null)
422 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 423 Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
426 OPT::Apply(*Prec_,*R2,*Y);
430 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 431 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
440 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 441 Teuchos::TimeMonitor lcltimer( *DotTime_ );
443 MVT::MvDot(*Y,*R2, beta);
445 for(
int i=0; i<ncols; i++)
446 beta[i] = sqrt(beta[i]);
451 for(
int i=0; i<ncols; i++)
453 delta[i] = cs[i]*dbar[i] + sn[i]*alpha[i];
454 gbar[i] = sn[i]*dbar[i] - cs[i]*alpha[i];
455 epsln[i] = sn[i]*beta[i];
456 dbar[i] = - cs[i]*beta[i];
460 for(
int i=0; i<ncols; i++)
462 symOrtho(gbar[i], beta[i], &cs[i], &sn[i], &gamma[i]);
464 phi[i] = cs[i]*phibar[i];
465 phibar[i] = sn[i]*phibar[i];
477 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 478 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
480 MVT::Assign(*W1,*scaleHelper);
483 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 484 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
486 MVT::MvScale(*scaleHelper,oldeps);
489 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 490 Teuchos::TimeMonitor lcltimer( *AddTime_ );
492 MVT::MvAddMv(ONE, *V, -ONE, *scaleHelper, *W);
495 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 496 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
498 MVT::Assign(*W2,*scaleHelper);
501 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 502 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
504 MVT::MvScale(*scaleHelper,delta);
507 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 508 Teuchos::TimeMonitor lcltimer( *AddTime_ );
510 MVT::MvAddMv(ONE, *W, -ONE, *scaleHelper, *W);
512 for(
int i=0; i<ncols; i++)
513 tmpvec[i] = ONE/gamma[i];
515 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 516 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
518 MVT::MvScale(*W,tmpvec);
524 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 525 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
527 MVT::Assign(*W,*scaleHelper);
530 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 531 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
533 MVT::MvScale(*scaleHelper,phi);
536 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 537 Teuchos::TimeMonitor lcltimer( *AddTime_ );
539 MVT::MvAddMv(ONE, *locX, ONE, *scaleHelper, *locX);
545 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 546 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
548 MVT::SetBlock(*locX,unconvergedIndices,*X_);
552 template <
class Scalar,
class MV,
class OP>
553 void PseudoBlockMinres<Scalar,MV,OP>::symOrtho(Scalar a, Scalar b, Scalar *c, Scalar *s, Scalar *r)
555 const Scalar absA = std::abs(a);
556 const Scalar absB = std::abs(b);
557 if ( absB == ZERO ) {
564 }
else if ( absA == ZERO ) {
568 }
else if ( absB >= absA ) {
571 *s = -ONE / sqrt( ONE+tau*tau );
573 *s = ONE / sqrt( ONE+tau*tau );
579 *c = -ONE / sqrt( ONE+tau*tau );
581 *c = ONE / sqrt( ONE+tau*tau );
Virtual base class which defines basic traits for the operator type.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Traits class which defines basic operations on multivectors.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...