Anasazi  Version of the Day
AnasaziMinres.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
33 #ifndef ANASAZI_MINRES_HPP
34 #define ANASAZI_MINRES_HPP
35 
36 #include "AnasaziConfigDefs.hpp"
37 #include "Teuchos_TimeMonitor.hpp"
38 
39 namespace Anasazi {
40 namespace Experimental {
41 
42 template <class Scalar, class MV, class OP>
43 class PseudoBlockMinres
44 {
47  const Scalar ONE;
48  const Scalar ZERO;
49 
50 public:
51  // Constructor
52  PseudoBlockMinres(Teuchos::RCP<OP> A, Teuchos::RCP<OP> Prec = Teuchos::null);
53 
54  // Set tolerance and maximum iterations
55  void setTol(const std::vector<Scalar>& tolerances) { tolerances_ = tolerances; };
56  void setMaxIter(const int maxIt) { maxIt_ = maxIt; };
57 
58  // Set solution and RHS
59  void setSol(Teuchos::RCP<MV> X) { X_ = X; };
60  void setRHS(Teuchos::RCP<const MV> B) { B_ = B; };
61 
62  // Solve the linear system
63  void solve();
64 
65 private:
67  Teuchos::RCP<OP> Prec_;
70  std::vector<Scalar> tolerances_;
71  int maxIt_;
72 
73  void symOrtho(Scalar a, Scalar b, Scalar *c, Scalar *s, Scalar *r);
74 
75 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
76  Teuchos::RCP<Teuchos::Time> AddTime_, ApplyOpTime_, ApplyPrecTime_, AssignTime_, DotTime_, LockTime_, NormTime_, ScaleTime_, TotalTime_;
77 #endif
78 };
79 
80 
81 
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())
86 {
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");
97 #endif
98 
99  A_ = A;
100  Prec_ = Prec;
101  maxIt_ = 0;
102 }
103 
104 
105 
106 template <class Scalar, class MV, class OP>
107 void PseudoBlockMinres<Scalar,MV,OP>::solve()
108 {
109  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
110  Teuchos::TimeMonitor outertimer( *TotalTime_ );
111  #endif
112 
113  // Get number of vectors
114  int ncols = MVT::GetNumberVecs(*B_);
115  int newNumConverged;
116 
117  // Declare some variables
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;
121 
122  // Allocate space for multivectors
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);
131 
132  // Reserve space for arrays
133  indicesToRemove.reserve(ncols);
134  newlyConverged.reserve(ncols);
135 
136  // Initialize array of unconverged indices
137  for(int i=0; i<ncols; i++)
138  unconvergedIndices[i] = i;
139 
140  // Get a local copy of X
141  // We want the vectors to remain contiguous even as things converge
142  locX = MVT::CloneCopy(*X_);
143 
144  // Initialize residuals
145  // R1 = B - AX
146  {
147  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
148  Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
149  #endif
150  OPT::Apply(*A_,*locX,*R1);
151  }
152  {
153  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
154  Teuchos::TimeMonitor lcltimer( *AddTime_ );
155  #endif
156  MVT::MvAddMv(ONE, *B_, -ONE, *R1, *R1);
157  }
158 
159  // R2 = R1
160  {
161  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
162  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
163  #endif
164  MVT::Assign(*R1,*R2);
165  }
166 
167  // Initialize the W's to 0.
168  MVT::MvInit (*W);
169  MVT::MvInit (*W2);
170 
171  // Y = M\R1 (preconditioned residual)
172  if(Prec_ != Teuchos::null)
173  {
174  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
175  Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
176  #endif
177 
178  OPT::Apply(*Prec_,*R1,*Y);
179  }
180  else
181  {
182  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
183  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
184  #endif
185  MVT::Assign(*R1,*Y);
186  }
187 
188  // beta1 = sqrt(Y'*R1)
189  {
190  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
191  Teuchos::TimeMonitor lcltimer( *DotTime_ );
192  #endif
193  MVT::MvDot(*Y,*R1, beta1);
194 
195  for(size_t i=0; i<beta1.size(); i++)
196  beta1[i] = sqrt(beta1[i]);
197  }
198 
199  // beta = beta1
200  beta = beta1;
201 
202  // phibar = beta1
203  phibar = beta1;
204 
206  // Begin Lanczos iterations
207  for(int iter=1; iter <= maxIt_; iter++)
208  {
209  // Test convergence
210  indicesToRemove.clear();
211  for(int i=0; i<ncols; i++)
212  {
213  Scalar relres = phibar[i]/beta1[i];
214 // std::cout << "relative residual[" << unconvergedIndices[i] << "] at iteration " << iter << " = " << relres << std::endl;
215 
216  // If the vector converged, mark it for termination
217  // Make sure to add it to
218  if(relres < tolerances_[i])
219  {
220  indicesToRemove.push_back(i);
221  }
222  }
223 
224  // Check whether anything converged
225  newNumConverged = indicesToRemove.size();
226  if(newNumConverged > 0)
227  {
228  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
229  Teuchos::TimeMonitor lcltimer( *LockTime_ );
230  #endif
231 
232  // If something converged, stick the converged vectors in X
233  newlyConverged.resize(newNumConverged);
234  for(int i=0; i<newNumConverged; i++)
235  newlyConverged[i] = unconvergedIndices[indicesToRemove[i]];
236 
237  Teuchos::RCP<MV> helperLocX = MVT::CloneCopy(*locX,indicesToRemove);
238 
239  MVT::SetBlock(*helperLocX,newlyConverged,*X_);
240 
241  // If everything has converged, we are done
242  if(newNumConverged == ncols)
243  return;
244 
245  // Update unconverged indices
246  std::vector<int> helperVec(ncols - newNumConverged);
247 
248  std::set_difference(unconvergedIndices.begin(), unconvergedIndices.end(), newlyConverged.begin(), newlyConverged.end(), helperVec.begin());
249  unconvergedIndices = helperVec;
250 
251  // Get indices of things we want to keep
252  std::vector<int> thingsToKeep(ncols - newNumConverged);
253  helperVec.resize(ncols);
254  for(int i=0; i<ncols; i++)
255  helperVec[i] = i;
256  ncols = ncols - newNumConverged;
257 
258  std::set_difference(helperVec.begin(), helperVec.end(), indicesToRemove.begin(), indicesToRemove.end(), thingsToKeep.begin());
259 
260  // Shrink the multivectors
261  Teuchos::RCP<MV> helperMV;
262  helperMV = MVT::CloneCopy(*V,thingsToKeep);
263  V = helperMV;
264  helperMV = MVT::CloneCopy(*Y,thingsToKeep);
265  Y = helperMV;
266  helperMV = MVT::CloneCopy(*R1,thingsToKeep);
267  R1 = helperMV;
268  helperMV = MVT::CloneCopy(*R2,thingsToKeep);
269  R2 = helperMV;
270  helperMV = MVT::CloneCopy(*W,thingsToKeep);
271  W = helperMV;
272  helperMV = MVT::CloneCopy(*W1,thingsToKeep);
273  W1 = helperMV;
274  helperMV = MVT::CloneCopy(*W2,thingsToKeep);
275  W2 = helperMV;
276  helperMV = MVT::CloneCopy(*locX,thingsToKeep);
277  locX = helperMV;
278  scaleHelper = MVT::Clone(*V,ncols);
279 
280  // Shrink the arrays
281  alpha.resize(ncols);
282  oldeps.resize(ncols);
283  delta.resize(ncols);
284  gbar.resize(ncols);
285  phi.resize(ncols);
286  gamma.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]];
291  beta = helperVecS;
292  for(int i=0; i<ncols; i++)
293  helperVecS[i] = beta1[thingsToKeep[i]];
294  beta1 = helperVecS;
295  for(int i=0; i<ncols; i++)
296  helperVecS[i] = phibar[thingsToKeep[i]];
297  phibar = helperVecS;
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]];
303  epsln = helperVecS;
304  for(int i=0; i<ncols; i++)
305  helperVecS[i] = cs[thingsToKeep[i]];
306  cs = helperVecS;
307  for(int i=0; i<ncols; i++)
308  helperVecS[i] = sn[thingsToKeep[i]];
309  sn = helperVecS;
310  for(int i=0; i<ncols; i++)
311  helperVecS[i] = dbar[thingsToKeep[i]];
312  dbar = helperVecS;
313 
314  // Tell operator about the removed indices
315  A_->removeIndices(indicesToRemove);
316  }
317 
318  // Normalize previous vector
319  // V = Y / beta
320  {
321  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
322  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
323  #endif
324  MVT::Assign(*Y,*V);
325  }
326  for(int i=0; i<ncols; i++)
327  tmpvec[i] = ONE/beta[i];
328  {
329  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
330  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
331  #endif
332  MVT::MvScale (*V, tmpvec);
333  }
334 
335  // Apply operator
336  // Y = AV
337  {
338  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
339  Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
340  #endif
341  OPT::Apply(*A_, *V, *Y);
342  }
343 
344  if(iter > 1)
345  {
346  // Y = Y - beta/oldBeta R1
347  for(int i=0; i<ncols; i++)
348  tmpvec[i] = beta[i]/oldBeta[i];
349  {
350  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
351  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
352  #endif
353  MVT::Assign(*R1,*scaleHelper);
354  }
355  {
356  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
357  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
358  #endif
359  MVT::MvScale(*scaleHelper,tmpvec);
360  }
361  {
362  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
363  Teuchos::TimeMonitor lcltimer( *AddTime_ );
364  #endif
365  MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
366  }
367  }
368 
369  // alpha = V'*Y
370  {
371  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
372  Teuchos::TimeMonitor lcltimer( *DotTime_ );
373  #endif
374  MVT::MvDot(*V, *Y, alpha);
375  }
376 
377  // Y = Y - alpha/beta R2
378  for(int i=0; i<ncols; i++)
379  tmpvec[i] = alpha[i]/beta[i];
380  {
381  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
382  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
383  #endif
384  MVT::Assign(*R2,*scaleHelper);
385  }
386  {
387  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
388  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
389  #endif
390  MVT::MvScale(*scaleHelper,tmpvec);
391  }
392  {
393  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
394  Teuchos::TimeMonitor lcltimer( *AddTime_ );
395  #endif
396  MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
397  }
398 
399  // R1 = R2
400  // R2 = Y
401  swapHelper = R1;
402  R1 = R2;
403  R2 = Y;
404  Y = swapHelper;
405 
406  // Y = M\R2
407  if(Prec_ != Teuchos::null)
408  {
409  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
410  Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
411  #endif
412 
413  OPT::Apply(*Prec_,*R2,*Y);
414  }
415  else
416  {
417  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
418  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
419  #endif
420  MVT::Assign(*R2,*Y);
421  }
422 
423  // Get new beta
424  // beta = sqrt(R2'*Y)
425  oldBeta = beta;
426  {
427  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
428  Teuchos::TimeMonitor lcltimer( *DotTime_ );
429  #endif
430  MVT::MvDot(*Y,*R2, beta);
431 
432  for(int i=0; i<ncols; i++)
433  beta[i] = sqrt(beta[i]);
434  }
435 
436  // Apply previous rotation
437  oldeps = epsln;
438  for(int i=0; i<ncols; i++)
439  {
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];
444  }
445 
446  // Compute the next plane rotation
447  for(int i=0; i<ncols; i++)
448  {
449  symOrtho(gbar[i], beta[i], &cs[i], &sn[i], &gamma[i]);
450 
451  phi[i] = cs[i]*phibar[i];
452  phibar[i] = sn[i]*phibar[i];
453  }
454 
455  // w1 = w2
456  // w2 = w
457  swapHelper = W1;
458  W1 = W2;
459  W2 = W;
460  W = swapHelper;
461 
462  // W = (V - oldeps*W1 - delta*W2) / gamma
463  {
464  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
465  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
466  #endif
467  MVT::Assign(*W1,*scaleHelper);
468  }
469  {
470  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
471  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
472  #endif
473  MVT::MvScale(*scaleHelper,oldeps);
474  }
475  {
476  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
477  Teuchos::TimeMonitor lcltimer( *AddTime_ );
478  #endif
479  MVT::MvAddMv(ONE, *V, -ONE, *scaleHelper, *W);
480  }
481  {
482  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
483  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
484  #endif
485  MVT::Assign(*W2,*scaleHelper);
486  }
487  {
488  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
489  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
490  #endif
491  MVT::MvScale(*scaleHelper,delta);
492  }
493  {
494  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
495  Teuchos::TimeMonitor lcltimer( *AddTime_ );
496  #endif
497  MVT::MvAddMv(ONE, *W, -ONE, *scaleHelper, *W);
498  }
499  for(int i=0; i<ncols; i++)
500  tmpvec[i] = ONE/gamma[i];
501  {
502  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
503  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
504  #endif
505  MVT::MvScale(*W,tmpvec);
506  }
507 
508  // Update X
509  // X = X + phi*W
510  {
511  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
512  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
513  #endif
514  MVT::Assign(*W,*scaleHelper);
515  }
516  {
517  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
518  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
519  #endif
520  MVT::MvScale(*scaleHelper,phi);
521  }
522  {
523  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
524  Teuchos::TimeMonitor lcltimer( *AddTime_ );
525  #endif
526  MVT::MvAddMv(ONE, *locX, ONE, *scaleHelper, *locX);
527  }
528  }
529 
530  // Stick unconverged vectors in X
531  {
532  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
533  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
534  #endif
535  MVT::SetBlock(*locX,unconvergedIndices,*X_);
536  }
537 }
538 
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)
541 {
542  const Scalar absA = std::abs(a);
543  const Scalar absB = std::abs(b);
544  if ( absB == ZERO ) {
545  *s = ZERO;
546  *r = absA;
547  if ( absA == ZERO )
548  *c = ONE;
549  else
550  *c = a / absA;
551  } else if ( absA == ZERO ) {
552  *c = ZERO;
553  *s = b / absB;
554  *r = absB;
555  } else if ( absB >= absA ) { // && a!=0 && b!=0
556  Scalar tau = a / b;
557  if ( b < ZERO )
558  *s = -ONE / sqrt( ONE+tau*tau );
559  else
560  *s = ONE / sqrt( ONE+tau*tau );
561  *c = *s * tau;
562  *r = b / *s;
563  } else { // (absA > absB) && a!=0 && b!=0
564  Scalar tau = b / a;
565  if ( a < ZERO )
566  *c = -ONE / sqrt( ONE+tau*tau );
567  else
568  *c = ONE / sqrt( ONE+tau*tau );
569  *s = *c * tau;
570  *r = a / *c;
571  }
572 }
573 
574 }} // End of namespace
575 
576 #endif
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...