33 #ifndef ANASAZI_MINRES_HPP 34 #define ANASAZI_MINRES_HPP 42 template <
class Scalar,
class MV,
class OP>
43 class PseudoBlockMinres
55 void setTol(
const std::vector<Scalar>& tolerances) { tolerances_ = tolerances; };
56 void setMaxIter(
const int maxIt) { maxIt_ = maxIt; };
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>
84 ONE(
Teuchos::ScalarTraits<Scalar>::one()),
85 ZERO(
Teuchos::ScalarTraits<Scalar>::zero())
87 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 106 template <
class Scalar,
class MV,
class OP>
107 void PseudoBlockMinres<Scalar,MV,OP>::solve()
109 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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);
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 150 OPT::Apply(*A_,*locX,*R1);
153 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 156 MVT::MvAddMv(ONE, *B_, -ONE, *R1, *R1);
161 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 164 MVT::Assign(*R1,*R2);
172 if(Prec_ != Teuchos::null)
174 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 178 OPT::Apply(*Prec_,*R1,*Y);
182 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 190 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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 233 newlyConverged.resize(newNumConverged);
234 for(
int i=0; i<newNumConverged; i++)
235 newlyConverged[i] = unconvergedIndices[indicesToRemove[i]];
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());
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 326 for(
int i=0; i<ncols; i++)
327 tmpvec[i] = ONE/beta[i];
329 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 332 MVT::MvScale (*V, tmpvec);
338 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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 353 MVT::Assign(*R1,*scaleHelper);
356 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 359 MVT::MvScale(*scaleHelper,tmpvec);
362 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 365 MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
371 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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 384 MVT::Assign(*R2,*scaleHelper);
387 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 390 MVT::MvScale(*scaleHelper,tmpvec);
393 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 396 MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
407 if(Prec_ != Teuchos::null)
409 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 413 OPT::Apply(*Prec_,*R2,*Y);
417 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 427 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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 467 MVT::Assign(*W1,*scaleHelper);
470 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 473 MVT::MvScale(*scaleHelper,oldeps);
476 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 479 MVT::MvAddMv(ONE, *V, -ONE, *scaleHelper, *W);
482 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 485 MVT::Assign(*W2,*scaleHelper);
488 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 491 MVT::MvScale(*scaleHelper,delta);
494 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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 505 MVT::MvScale(*W,tmpvec);
511 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 514 MVT::Assign(*W,*scaleHelper);
517 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 520 MVT::MvScale(*scaleHelper,phi);
523 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 526 MVT::MvAddMv(ONE, *locX, ONE, *scaleHelper, *locX);
532 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 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...
static RCP< Time > getNewTimer(const std::string &name)
Traits class which defines basic operations on multivectors.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...