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 );