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:
66  Teuchos::RCP<OP> A_;
67  Teuchos::RCP<OP> Prec_;
68  Teuchos::RCP<MV> X_;
69  Teuchos::RCP<const MV> B_;
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...
Traits class which defines basic operations on multivectors.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...