33 #ifndef ANASAZI_MINRES_HPP 34 #define ANASAZI_MINRES_HPP 37 #include "Teuchos_TimeMonitor.hpp" 42 template <
class Scalar,
class MV,
class OP>
43 class PseudoBlockMinres
52 PseudoBlockMinres(Teuchos::RCP<OP> A, Teuchos::RCP<OP> Prec = Teuchos::null);
55 void setTol(
const std::vector<Scalar>& tolerances) { tolerances_ = tolerances; };
56 void setMaxIter(
const int maxIt) { maxIt_ = maxIt; };
59 void setSol(Teuchos::RCP<MV> X) { X_ = X; };
60 void setRHS(Teuchos::RCP<const MV> B) { B_ = B; };
67 Teuchos::RCP<OP> Prec_;
69 Teuchos::RCP<const MV> B_;
70 std::vector<Scalar> tolerances_;
73 void symOrtho(Scalar a, Scalar b, Scalar *c, Scalar *s, Scalar *r);
75 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 76 Teuchos::RCP<Teuchos::Time> AddTime_, ApplyOpTime_, ApplyPrecTime_, AssignTime_, DotTime_, LockTime_, NormTime_, ScaleTime_, TotalTime_;
82 template <
class Scalar,
class MV,
class OP>
83 PseudoBlockMinres<Scalar,MV,OP>::PseudoBlockMinres(Teuchos::RCP<OP> A, Teuchos::RCP<OP> Prec) :
84 ONE(Teuchos::ScalarTraits<Scalar>::one()),
85 ZERO(Teuchos::ScalarTraits<Scalar>::zero())
87 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 88 AddTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Add Multivectors");
89 ApplyOpTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Apply Operator");
90 ApplyPrecTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Apply Preconditioner");
91 AssignTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Assignment (no locking)");
92 DotTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Compute Dot Product");
93 LockTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Lock Converged Vectors");
94 NormTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Compute Norm");
95 ScaleTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Scale Multivector");
96 TotalTime_ = Teuchos::TimeMonitor::getNewTimer(
"Anasazi: PseudoBlockMinres::Total Time");
106 template <
class Scalar,
class MV,
class OP>
107 void PseudoBlockMinres<Scalar,MV,OP>::solve()
109 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 110 Teuchos::TimeMonitor outertimer( *TotalTime_ );
114 int ncols = MVT::GetNumberVecs(*B_);
118 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);
119 std::vector<int> indicesToRemove, newlyConverged, unconvergedIndices(ncols);
120 Teuchos::RCP<MV> V, Y, R1, R2, W, W1, W2, locX, scaleHelper, swapHelper;
123 V = MVT::Clone(*B_, ncols);
124 Y = MVT::Clone(*B_, ncols);
125 R1 = MVT::Clone(*B_, ncols);
126 R2 = MVT::Clone(*B_, ncols);
127 W = MVT::Clone(*B_, ncols);
128 W1 = MVT::Clone(*B_, ncols);
129 W2 = MVT::Clone(*B_, ncols);
130 scaleHelper = MVT::Clone(*B_, ncols);
133 indicesToRemove.reserve(ncols);
134 newlyConverged.reserve(ncols);
137 for(
int i=0; i<ncols; i++)
138 unconvergedIndices[i] = i;
142 locX = MVT::CloneCopy(*X_);
147 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 148 Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
150 OPT::Apply(*A_,*locX,*R1);
153 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 154 Teuchos::TimeMonitor lcltimer( *AddTime_ );
156 MVT::MvAddMv(ONE, *B_, -ONE, *R1, *R1);
161 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 162 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
164 MVT::Assign(*R1,*R2);
172 if(Prec_ != Teuchos::null)
174 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 175 Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
178 OPT::Apply(*Prec_,*R1,*Y);
182 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 183 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
190 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 191 Teuchos::TimeMonitor lcltimer( *DotTime_ );
193 MVT::MvDot(*Y,*R1, beta1);
195 for(
size_t i=0; i<beta1.size(); i++)
196 beta1[i] = sqrt(beta1[i]);
207 for(
int iter=1; iter <= maxIt_; iter++)
210 indicesToRemove.clear();
211 for(
int i=0; i<ncols; i++)
213 Scalar relres = phibar[i]/beta1[i];
218 if(relres < tolerances_[i])
220 indicesToRemove.push_back(i);
225 newNumConverged = indicesToRemove.size();
226 if(newNumConverged > 0)
228 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 229 Teuchos::TimeMonitor lcltimer( *LockTime_ );
233 newlyConverged.resize(newNumConverged);
234 for(
int i=0; i<newNumConverged; i++)
235 newlyConverged[i] = unconvergedIndices[indicesToRemove[i]];
237 Teuchos::RCP<MV> helperLocX = MVT::CloneCopy(*locX,indicesToRemove);
239 MVT::SetBlock(*helperLocX,newlyConverged,*X_);
242 if(newNumConverged == ncols)
246 std::vector<int> helperVec(ncols - newNumConverged);
248 std::set_difference(unconvergedIndices.begin(), unconvergedIndices.end(), newlyConverged.begin(), newlyConverged.end(), helperVec.begin());
249 unconvergedIndices = helperVec;
252 std::vector<int> thingsToKeep(ncols - newNumConverged);
253 helperVec.resize(ncols);
254 for(
int i=0; i<ncols; i++)
256 ncols = ncols - newNumConverged;
258 std::set_difference(helperVec.begin(), helperVec.end(), indicesToRemove.begin(), indicesToRemove.end(), thingsToKeep.begin());
261 Teuchos::RCP<MV> helperMV;
262 helperMV = MVT::CloneCopy(*V,thingsToKeep);
264 helperMV = MVT::CloneCopy(*Y,thingsToKeep);
266 helperMV = MVT::CloneCopy(*R1,thingsToKeep);
268 helperMV = MVT::CloneCopy(*R2,thingsToKeep);
270 helperMV = MVT::CloneCopy(*W,thingsToKeep);
272 helperMV = MVT::CloneCopy(*W1,thingsToKeep);
274 helperMV = MVT::CloneCopy(*W2,thingsToKeep);
276 helperMV = MVT::CloneCopy(*locX,thingsToKeep);
278 scaleHelper = MVT::Clone(*V,ncols);
282 oldeps.resize(ncols);
287 tmpvec.resize(ncols);
288 std::vector<Scalar> helperVecS(ncols);
289 for(
int i=0; i<ncols; i++)
290 helperVecS[i] = beta[thingsToKeep[i]];
292 for(
int i=0; i<ncols; i++)
293 helperVecS[i] = beta1[thingsToKeep[i]];
295 for(
int i=0; i<ncols; i++)
296 helperVecS[i] = phibar[thingsToKeep[i]];
298 for(
int i=0; i<ncols; i++)
299 helperVecS[i] = oldBeta[thingsToKeep[i]];
300 oldBeta = helperVecS;
301 for(
int i=0; i<ncols; i++)
302 helperVecS[i] = epsln[thingsToKeep[i]];
304 for(
int i=0; i<ncols; i++)
305 helperVecS[i] = cs[thingsToKeep[i]];
307 for(
int i=0; i<ncols; i++)
308 helperVecS[i] = sn[thingsToKeep[i]];
310 for(
int i=0; i<ncols; i++)
311 helperVecS[i] = dbar[thingsToKeep[i]];
315 A_->removeIndices(indicesToRemove);
321 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 322 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
326 for(
int i=0; i<ncols; i++)
327 tmpvec[i] = ONE/beta[i];
329 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 330 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
332 MVT::MvScale (*V, tmpvec);
338 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 339 Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
341 OPT::Apply(*A_, *V, *Y);
347 for(
int i=0; i<ncols; i++)
348 tmpvec[i] = beta[i]/oldBeta[i];
350 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 351 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
353 MVT::Assign(*R1,*scaleHelper);
356 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 357 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
359 MVT::MvScale(*scaleHelper,tmpvec);
362 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 363 Teuchos::TimeMonitor lcltimer( *AddTime_ );
365 MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
371 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 372 Teuchos::TimeMonitor lcltimer( *DotTime_ );
374 MVT::MvDot(*V, *Y, alpha);
378 for(
int i=0; i<ncols; i++)
379 tmpvec[i] = alpha[i]/beta[i];
381 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 382 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
384 MVT::Assign(*R2,*scaleHelper);
387 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 388 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
390 MVT::MvScale(*scaleHelper,tmpvec);
393 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 394 Teuchos::TimeMonitor lcltimer( *AddTime_ );
396 MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
407 if(Prec_ != Teuchos::null)
409 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 410 Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
413 OPT::Apply(*Prec_,*R2,*Y);
417 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 418 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
427 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 428 Teuchos::TimeMonitor lcltimer( *DotTime_ );
430 MVT::MvDot(*Y,*R2, beta);
432 for(
int i=0; i<ncols; i++)
433 beta[i] = sqrt(beta[i]);
438 for(
int i=0; i<ncols; i++)
440 delta[i] = cs[i]*dbar[i] + sn[i]*alpha[i];
441 gbar[i] = sn[i]*dbar[i] - cs[i]*alpha[i];
442 epsln[i] = sn[i]*beta[i];
443 dbar[i] = - cs[i]*beta[i];
447 for(
int i=0; i<ncols; i++)
449 symOrtho(gbar[i], beta[i], &cs[i], &sn[i], &gamma[i]);
451 phi[i] = cs[i]*phibar[i];
452 phibar[i] = sn[i]*phibar[i];
464 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 465 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
467 MVT::Assign(*W1,*scaleHelper);
470 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 471 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
473 MVT::MvScale(*scaleHelper,oldeps);
476 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 477 Teuchos::TimeMonitor lcltimer( *AddTime_ );
479 MVT::MvAddMv(ONE, *V, -ONE, *scaleHelper, *W);
482 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 483 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
485 MVT::Assign(*W2,*scaleHelper);
488 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 489 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
491 MVT::MvScale(*scaleHelper,delta);
494 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 495 Teuchos::TimeMonitor lcltimer( *AddTime_ );
497 MVT::MvAddMv(ONE, *W, -ONE, *scaleHelper, *W);
499 for(
int i=0; i<ncols; i++)
500 tmpvec[i] = ONE/gamma[i];
502 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 503 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
505 MVT::MvScale(*W,tmpvec);
511 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 512 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
514 MVT::Assign(*W,*scaleHelper);
517 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 518 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
520 MVT::MvScale(*scaleHelper,phi);
523 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 524 Teuchos::TimeMonitor lcltimer( *AddTime_ );
526 MVT::MvAddMv(ONE, *locX, ONE, *scaleHelper, *locX);
532 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 533 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
535 MVT::SetBlock(*locX,unconvergedIndices,*X_);
539 template <
class Scalar,
class MV,
class OP>
540 void PseudoBlockMinres<Scalar,MV,OP>::symOrtho(Scalar a, Scalar b, Scalar *c, Scalar *s, Scalar *r)
542 const Scalar absA = std::abs(a);
543 const Scalar absB = std::abs(b);
544 if ( absB == ZERO ) {
551 }
else if ( absA == ZERO ) {
555 }
else if ( absB >= absA ) {
558 *s = -ONE / sqrt( ONE+tau*tau );
560 *s = ONE / sqrt( ONE+tau*tau );
566 *c = -ONE / sqrt( ONE+tau*tau );
568 *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...