46#ifndef ANASAZI_MINRES_HPP
47#define ANASAZI_MINRES_HPP
50#include "Teuchos_TimeMonitor.hpp"
55template <
class Scalar,
class MV,
class OP>
56class 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_;
95template <
class Scalar,
class MV,
class OP>
96PseudoBlockMinres<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");
119template <
class Scalar,
class MV,
class OP>
120void 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_);
552template <
class Scalar,
class MV,
class OP>
553void 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 );
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Traits class which defines basic operations on multivectors.
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.