Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AnasaziMinres.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Anasazi: Block Eigensolvers Package
4 //
5 // Copyright 2004 NTESS and the Anasazi contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
14 #ifndef ANASAZI_MINRES_HPP
15 #define ANASAZI_MINRES_HPP
16 
17 #include "AnasaziConfigDefs.hpp"
18 #include "Teuchos_TimeMonitor.hpp"
19 
20 namespace Anasazi {
21 namespace Experimental {
22 
23 template <class Scalar, class MV, class OP>
24 class PseudoBlockMinres
25 {
26  typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
27  typedef Anasazi::OperatorTraits<Scalar,MV,OP> OPT;
28  const Scalar ONE;
29  const Scalar ZERO;
30 
31 public:
32  // Constructor
33  PseudoBlockMinres(Teuchos::RCP<OP> A, Teuchos::RCP<OP> Prec = Teuchos::null);
34 
35  // Set tolerance and maximum iterations
36  void setTol(const std::vector<Scalar>& tolerances) { tolerances_ = tolerances; };
37  void setMaxIter(const int maxIt) { maxIt_ = maxIt; };
38 
39  // Set solution and RHS
40  void setSol(Teuchos::RCP<MV> X) { X_ = X; };
41  void setRHS(Teuchos::RCP<const MV> B) { B_ = B; };
42 
43  // Solve the linear system
44  void solve();
45 
46 private:
47  Teuchos::RCP<OP> A_;
48  Teuchos::RCP<OP> Prec_;
49  Teuchos::RCP<MV> X_;
50  Teuchos::RCP<const MV> B_;
51  std::vector<Scalar> tolerances_;
52  int maxIt_;
53 
54  void symOrtho(Scalar a, Scalar b, Scalar *c, Scalar *s, Scalar *r);
55 
56 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
57  Teuchos::RCP<Teuchos::Time> AddTime_, ApplyOpTime_, ApplyPrecTime_, AssignTime_, DotTime_, LockTime_, NormTime_, ScaleTime_, TotalTime_;
58 #endif
59 };
60 
61 
62 
63 template <class Scalar, class MV, class OP>
64 PseudoBlockMinres<Scalar,MV,OP>::PseudoBlockMinres(Teuchos::RCP<OP> A, Teuchos::RCP<OP> Prec) :
65  ONE(Teuchos::ScalarTraits<Scalar>::one()),
66  ZERO(Teuchos::ScalarTraits<Scalar>::zero())
67 {
68 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
69  AddTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Add Multivectors");
70  ApplyOpTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Apply Operator");
71  ApplyPrecTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Apply Preconditioner");
72  AssignTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Assignment (no locking)");
73  DotTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Compute Dot Product");
74  LockTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Lock Converged Vectors");
75  NormTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Compute Norm");
76  ScaleTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Scale Multivector");
77  TotalTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Total Time");
78 #endif
79 
80  A_ = A;
81  Prec_ = Prec;
82  maxIt_ = 0;
83 }
84 
85 
86 
87 template <class Scalar, class MV, class OP>
88 void PseudoBlockMinres<Scalar,MV,OP>::solve()
89 {
90  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
91  Teuchos::TimeMonitor outertimer( *TotalTime_ );
92  #endif
93 
94  // Get number of vectors
95  int ncols = MVT::GetNumberVecs(*B_);
96  int newNumConverged;
97 
98  // Declare some variables
99  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);
100  std::vector<int> indicesToRemove, newlyConverged, unconvergedIndices(ncols);
101  Teuchos::RCP<MV> V, Y, R1, R2, W, W1, W2, locX, scaleHelper, swapHelper;
102 
103  // Allocate space for multivectors
104  V = MVT::Clone(*B_, ncols);
105  Y = MVT::Clone(*B_, ncols);
106  R1 = MVT::Clone(*B_, ncols);
107  R2 = MVT::Clone(*B_, ncols);
108  W = MVT::Clone(*B_, ncols);
109  W1 = MVT::Clone(*B_, ncols);
110  W2 = MVT::Clone(*B_, ncols);
111  scaleHelper = MVT::Clone(*B_, ncols);
112 
113  // Reserve space for arrays
114  indicesToRemove.reserve(ncols);
115  newlyConverged.reserve(ncols);
116 
117  // Initialize array of unconverged indices
118  for(int i=0; i<ncols; i++)
119  unconvergedIndices[i] = i;
120 
121  // Get a local copy of X
122  // We want the vectors to remain contiguous even as things converge
123  locX = MVT::CloneCopy(*X_);
124 
125  // Initialize residuals
126  // R1 = B - AX
127  {
128  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
129  Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
130  #endif
131  OPT::Apply(*A_,*locX,*R1);
132  }
133  {
134  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
135  Teuchos::TimeMonitor lcltimer( *AddTime_ );
136  #endif
137  MVT::MvAddMv(ONE, *B_, -ONE, *R1, *R1);
138  }
139 
140  // R2 = R1
141  {
142  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
143  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
144  #endif
145  MVT::Assign(*R1,*R2);
146  }
147 
148  // Initialize the W's to 0.
149  MVT::MvInit (*W);
150  MVT::MvInit (*W2);
151 
152  // Y = M\R1 (preconditioned residual)
153  if(Prec_ != Teuchos::null)
154  {
155  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
156  Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
157  #endif
158 
159  OPT::Apply(*Prec_,*R1,*Y);
160  }
161  else
162  {
163  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
164  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
165  #endif
166  MVT::Assign(*R1,*Y);
167  }
168 
169  // beta1 = sqrt(Y'*R1)
170  {
171  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
172  Teuchos::TimeMonitor lcltimer( *DotTime_ );
173  #endif
174  MVT::MvDot(*Y,*R1, beta1);
175 
176  for(size_t i=0; i<beta1.size(); i++)
177  beta1[i] = sqrt(beta1[i]);
178  }
179 
180  // beta = beta1
181  beta = beta1;
182 
183  // phibar = beta1
184  phibar = beta1;
185 
187  // Begin Lanczos iterations
188  for(int iter=1; iter <= maxIt_; iter++)
189  {
190  // Test convergence
191  indicesToRemove.clear();
192  for(int i=0; i<ncols; i++)
193  {
194  Scalar relres = phibar[i]/beta1[i];
195 // std::cout << "relative residual[" << unconvergedIndices[i] << "] at iteration " << iter << " = " << relres << std::endl;
196 
197  // If the vector converged, mark it for termination
198  // Make sure to add it to
199  if(relres < tolerances_[i])
200  {
201  indicesToRemove.push_back(i);
202  }
203  }
204 
205  // Check whether anything converged
206  newNumConverged = indicesToRemove.size();
207  if(newNumConverged > 0)
208  {
209  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
210  Teuchos::TimeMonitor lcltimer( *LockTime_ );
211  #endif
212 
213  // If something converged, stick the converged vectors in X
214  newlyConverged.resize(newNumConverged);
215  for(int i=0; i<newNumConverged; i++)
216  newlyConverged[i] = unconvergedIndices[indicesToRemove[i]];
217 
218  Teuchos::RCP<MV> helperLocX = MVT::CloneCopy(*locX,indicesToRemove);
219 
220  MVT::SetBlock(*helperLocX,newlyConverged,*X_);
221 
222  // If everything has converged, we are done
223  if(newNumConverged == ncols)
224  return;
225 
226  // Update unconverged indices
227  std::vector<int> helperVec(ncols - newNumConverged);
228 
229  std::set_difference(unconvergedIndices.begin(), unconvergedIndices.end(), newlyConverged.begin(), newlyConverged.end(), helperVec.begin());
230  unconvergedIndices = helperVec;
231 
232  // Get indices of things we want to keep
233  std::vector<int> thingsToKeep(ncols - newNumConverged);
234  helperVec.resize(ncols);
235  for(int i=0; i<ncols; i++)
236  helperVec[i] = i;
237  ncols = ncols - newNumConverged;
238 
239  std::set_difference(helperVec.begin(), helperVec.end(), indicesToRemove.begin(), indicesToRemove.end(), thingsToKeep.begin());
240 
241  // Shrink the multivectors
242  Teuchos::RCP<MV> helperMV;
243  helperMV = MVT::CloneCopy(*V,thingsToKeep);
244  V = helperMV;
245  helperMV = MVT::CloneCopy(*Y,thingsToKeep);
246  Y = helperMV;
247  helperMV = MVT::CloneCopy(*R1,thingsToKeep);
248  R1 = helperMV;
249  helperMV = MVT::CloneCopy(*R2,thingsToKeep);
250  R2 = helperMV;
251  helperMV = MVT::CloneCopy(*W,thingsToKeep);
252  W = helperMV;
253  helperMV = MVT::CloneCopy(*W1,thingsToKeep);
254  W1 = helperMV;
255  helperMV = MVT::CloneCopy(*W2,thingsToKeep);
256  W2 = helperMV;
257  helperMV = MVT::CloneCopy(*locX,thingsToKeep);
258  locX = helperMV;
259  scaleHelper = MVT::Clone(*V,ncols);
260 
261  // Shrink the arrays
262  alpha.resize(ncols);
263  oldeps.resize(ncols);
264  delta.resize(ncols);
265  gbar.resize(ncols);
266  phi.resize(ncols);
267  gamma.resize(ncols);
268  tmpvec.resize(ncols);
269  std::vector<Scalar> helperVecS(ncols);
270  for(int i=0; i<ncols; i++)
271  helperVecS[i] = beta[thingsToKeep[i]];
272  beta = helperVecS;
273  for(int i=0; i<ncols; i++)
274  helperVecS[i] = beta1[thingsToKeep[i]];
275  beta1 = helperVecS;
276  for(int i=0; i<ncols; i++)
277  helperVecS[i] = phibar[thingsToKeep[i]];
278  phibar = helperVecS;
279  for(int i=0; i<ncols; i++)
280  helperVecS[i] = oldBeta[thingsToKeep[i]];
281  oldBeta = helperVecS;
282  for(int i=0; i<ncols; i++)
283  helperVecS[i] = epsln[thingsToKeep[i]];
284  epsln = helperVecS;
285  for(int i=0; i<ncols; i++)
286  helperVecS[i] = cs[thingsToKeep[i]];
287  cs = helperVecS;
288  for(int i=0; i<ncols; i++)
289  helperVecS[i] = sn[thingsToKeep[i]];
290  sn = helperVecS;
291  for(int i=0; i<ncols; i++)
292  helperVecS[i] = dbar[thingsToKeep[i]];
293  dbar = helperVecS;
294 
295  // Tell operator about the removed indices
296  A_->removeIndices(indicesToRemove);
297  }
298 
299  // Normalize previous vector
300  // V = Y / beta
301  {
302  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
303  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
304  #endif
305  MVT::Assign(*Y,*V);
306  }
307  for(int i=0; i<ncols; i++)
308  tmpvec[i] = ONE/beta[i];
309  {
310  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
311  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
312  #endif
313  MVT::MvScale (*V, tmpvec);
314  }
315 
316  // Apply operator
317  // Y = AV
318  {
319  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
320  Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
321  #endif
322  OPT::Apply(*A_, *V, *Y);
323  }
324 
325  if(iter > 1)
326  {
327  // Y = Y - beta/oldBeta R1
328  for(int i=0; i<ncols; i++)
329  tmpvec[i] = beta[i]/oldBeta[i];
330  {
331  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
332  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
333  #endif
334  MVT::Assign(*R1,*scaleHelper);
335  }
336  {
337  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
338  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
339  #endif
340  MVT::MvScale(*scaleHelper,tmpvec);
341  }
342  {
343  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
344  Teuchos::TimeMonitor lcltimer( *AddTime_ );
345  #endif
346  MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
347  }
348  }
349 
350  // alpha = V'*Y
351  {
352  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
353  Teuchos::TimeMonitor lcltimer( *DotTime_ );
354  #endif
355  MVT::MvDot(*V, *Y, alpha);
356  }
357 
358  // Y = Y - alpha/beta R2
359  for(int i=0; i<ncols; i++)
360  tmpvec[i] = alpha[i]/beta[i];
361  {
362  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
363  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
364  #endif
365  MVT::Assign(*R2,*scaleHelper);
366  }
367  {
368  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
369  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
370  #endif
371  MVT::MvScale(*scaleHelper,tmpvec);
372  }
373  {
374  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
375  Teuchos::TimeMonitor lcltimer( *AddTime_ );
376  #endif
377  MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
378  }
379 
380  // R1 = R2
381  // R2 = Y
382  swapHelper = R1;
383  R1 = R2;
384  R2 = Y;
385  Y = swapHelper;
386 
387  // Y = M\R2
388  if(Prec_ != Teuchos::null)
389  {
390  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
391  Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
392  #endif
393 
394  OPT::Apply(*Prec_,*R2,*Y);
395  }
396  else
397  {
398  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
399  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
400  #endif
401  MVT::Assign(*R2,*Y);
402  }
403 
404  // Get new beta
405  // beta = sqrt(R2'*Y)
406  oldBeta = beta;
407  {
408  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
409  Teuchos::TimeMonitor lcltimer( *DotTime_ );
410  #endif
411  MVT::MvDot(*Y,*R2, beta);
412 
413  for(int i=0; i<ncols; i++)
414  beta[i] = sqrt(beta[i]);
415  }
416 
417  // Apply previous rotation
418  oldeps = epsln;
419  for(int i=0; i<ncols; i++)
420  {
421  delta[i] = cs[i]*dbar[i] + sn[i]*alpha[i];
422  gbar[i] = sn[i]*dbar[i] - cs[i]*alpha[i];
423  epsln[i] = sn[i]*beta[i];
424  dbar[i] = - cs[i]*beta[i];
425  }
426 
427  // Compute the next plane rotation
428  for(int i=0; i<ncols; i++)
429  {
430  symOrtho(gbar[i], beta[i], &cs[i], &sn[i], &gamma[i]);
431 
432  phi[i] = cs[i]*phibar[i];
433  phibar[i] = sn[i]*phibar[i];
434  }
435 
436  // w1 = w2
437  // w2 = w
438  swapHelper = W1;
439  W1 = W2;
440  W2 = W;
441  W = swapHelper;
442 
443  // W = (V - oldeps*W1 - delta*W2) / gamma
444  {
445  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
446  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
447  #endif
448  MVT::Assign(*W1,*scaleHelper);
449  }
450  {
451  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
452  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
453  #endif
454  MVT::MvScale(*scaleHelper,oldeps);
455  }
456  {
457  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
458  Teuchos::TimeMonitor lcltimer( *AddTime_ );
459  #endif
460  MVT::MvAddMv(ONE, *V, -ONE, *scaleHelper, *W);
461  }
462  {
463  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
464  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
465  #endif
466  MVT::Assign(*W2,*scaleHelper);
467  }
468  {
469  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
470  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
471  #endif
472  MVT::MvScale(*scaleHelper,delta);
473  }
474  {
475  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
476  Teuchos::TimeMonitor lcltimer( *AddTime_ );
477  #endif
478  MVT::MvAddMv(ONE, *W, -ONE, *scaleHelper, *W);
479  }
480  for(int i=0; i<ncols; i++)
481  tmpvec[i] = ONE/gamma[i];
482  {
483  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
484  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
485  #endif
486  MVT::MvScale(*W,tmpvec);
487  }
488 
489  // Update X
490  // X = X + phi*W
491  {
492  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
493  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
494  #endif
495  MVT::Assign(*W,*scaleHelper);
496  }
497  {
498  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
499  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
500  #endif
501  MVT::MvScale(*scaleHelper,phi);
502  }
503  {
504  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
505  Teuchos::TimeMonitor lcltimer( *AddTime_ );
506  #endif
507  MVT::MvAddMv(ONE, *locX, ONE, *scaleHelper, *locX);
508  }
509  }
510 
511  // Stick unconverged vectors in X
512  {
513  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
514  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
515  #endif
516  MVT::SetBlock(*locX,unconvergedIndices,*X_);
517  }
518 }
519 
520 template <class Scalar, class MV, class OP>
521 void PseudoBlockMinres<Scalar,MV,OP>::symOrtho(Scalar a, Scalar b, Scalar *c, Scalar *s, Scalar *r)
522 {
523  const Scalar absA = std::abs(a);
524  const Scalar absB = std::abs(b);
525  if ( absB == ZERO ) {
526  *s = ZERO;
527  *r = absA;
528  if ( absA == ZERO )
529  *c = ONE;
530  else
531  *c = a / absA;
532  } else if ( absA == ZERO ) {
533  *c = ZERO;
534  *s = b / absB;
535  *r = absB;
536  } else if ( absB >= absA ) { // && a!=0 && b!=0
537  Scalar tau = a / b;
538  if ( b < ZERO )
539  *s = -ONE / sqrt( ONE+tau*tau );
540  else
541  *s = ONE / sqrt( ONE+tau*tau );
542  *c = *s * tau;
543  *r = b / *s;
544  } else { // (absA > absB) && a!=0 && b!=0
545  Scalar tau = b / a;
546  if ( a < ZERO )
547  *c = -ONE / sqrt( ONE+tau*tau );
548  else
549  *c = ONE / sqrt( ONE+tau*tau );
550  *s = *c * tau;
551  *r = a / *c;
552  }
553 }
554 
555 }} // End of namespace
556 
557 #endif
Anasazi::OperatorTraits
Virtual base class which defines basic traits for the operator type.
Definition: AnasaziOperatorTraits.hpp:52
Teuchos::TimeMonitor::getNewTimer
static RCP< Time > getNewTimer(const std::string &name)
Anasazi::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition: AnasaziMultiVecTraits.hpp:95
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Teuchos_TimeMonitor.hpp
Generated on Thu Nov 21 2024 09:18:50 for Anasazi by   doxygen 1.8.5