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 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
46 #ifndef ANASAZI_MINRES_HPP
47 #define ANASAZI_MINRES_HPP
48 
49 #include "AnasaziConfigDefs.hpp"
50 #include "Teuchos_TimeMonitor.hpp"
51 
52 namespace Anasazi {
53 namespace Experimental {
54 
55 template <class Scalar, class MV, class OP>
56 class PseudoBlockMinres
57 {
60  const Scalar ONE;
61  const Scalar ZERO;
62 
63 public:
64  // Constructor
65  PseudoBlockMinres(Teuchos::RCP<OP> A, Teuchos::RCP<OP> Prec = Teuchos::null);
66 
67  // Set tolerance and maximum iterations
68  void setTol(const std::vector<Scalar>& tolerances) { tolerances_ = tolerances; };
69  void setMaxIter(const int maxIt) { maxIt_ = maxIt; };
70 
71  // Set solution and RHS
72  void setSol(Teuchos::RCP<MV> X) { X_ = X; };
73  void setRHS(Teuchos::RCP<const MV> B) { B_ = B; };
74 
75  // Solve the linear system
76  void solve();
77 
78 private:
79  Teuchos::RCP<OP> A_;
80  Teuchos::RCP<OP> Prec_;
81  Teuchos::RCP<MV> X_;
82  Teuchos::RCP<const MV> B_;
83  std::vector<Scalar> tolerances_;
84  int maxIt_;
85 
86  void symOrtho(Scalar a, Scalar b, Scalar *c, Scalar *s, Scalar *r);
87 
88 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
89  Teuchos::RCP<Teuchos::Time> AddTime_, ApplyOpTime_, ApplyPrecTime_, AssignTime_, DotTime_, LockTime_, NormTime_, ScaleTime_, TotalTime_;
90 #endif
91 };
92 
93 
94 
95 template <class Scalar, class MV, class OP>
96 PseudoBlockMinres<Scalar,MV,OP>::PseudoBlockMinres(Teuchos::RCP<OP> A, Teuchos::RCP<OP> Prec) :
97  ONE(Teuchos::ScalarTraits<Scalar>::one()),
98  ZERO(Teuchos::ScalarTraits<Scalar>::zero())
99 {
100 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
101  AddTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Add Multivectors");
102  ApplyOpTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Apply Operator");
103  ApplyPrecTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Apply Preconditioner");
104  AssignTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Assignment (no locking)");
105  DotTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Compute Dot Product");
106  LockTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Lock Converged Vectors");
107  NormTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Compute Norm");
108  ScaleTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Scale Multivector");
109  TotalTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Total Time");
110 #endif
111 
112  A_ = A;
113  Prec_ = Prec;
114  maxIt_ = 0;
115 }
116 
117 
118 
119 template <class Scalar, class MV, class OP>
120 void PseudoBlockMinres<Scalar,MV,OP>::solve()
121 {
122  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
123  Teuchos::TimeMonitor outertimer( *TotalTime_ );
124  #endif
125 
126  // Get number of vectors
127  int ncols = MVT::GetNumberVecs(*B_);
128  int newNumConverged;
129 
130  // Declare some variables
131  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);
132  std::vector<int> indicesToRemove, newlyConverged, unconvergedIndices(ncols);
133  Teuchos::RCP<MV> V, Y, R1, R2, W, W1, W2, locX, scaleHelper, swapHelper;
134 
135  // Allocate space for multivectors
136  V = MVT::Clone(*B_, ncols);
137  Y = MVT::Clone(*B_, ncols);
138  R1 = MVT::Clone(*B_, ncols);
139  R2 = MVT::Clone(*B_, ncols);
140  W = MVT::Clone(*B_, ncols);
141  W1 = MVT::Clone(*B_, ncols);
142  W2 = MVT::Clone(*B_, ncols);
143  scaleHelper = MVT::Clone(*B_, ncols);
144 
145  // Reserve space for arrays
146  indicesToRemove.reserve(ncols);
147  newlyConverged.reserve(ncols);
148 
149  // Initialize array of unconverged indices
150  for(int i=0; i<ncols; i++)
151  unconvergedIndices[i] = i;
152 
153  // Get a local copy of X
154  // We want the vectors to remain contiguous even as things converge
155  locX = MVT::CloneCopy(*X_);
156 
157  // Initialize residuals
158  // R1 = B - AX
159  {
160  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
161  Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
162  #endif
163  OPT::Apply(*A_,*locX,*R1);
164  }
165  {
166  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
167  Teuchos::TimeMonitor lcltimer( *AddTime_ );
168  #endif
169  MVT::MvAddMv(ONE, *B_, -ONE, *R1, *R1);
170  }
171 
172  // R2 = R1
173  {
174  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
175  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
176  #endif
177  MVT::Assign(*R1,*R2);
178  }
179 
180  // Initialize the W's to 0.
181  MVT::MvInit (*W);
182  MVT::MvInit (*W2);
183 
184  // Y = M\R1 (preconditioned residual)
185  if(Prec_ != Teuchos::null)
186  {
187  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
188  Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
189  #endif
190 
191  OPT::Apply(*Prec_,*R1,*Y);
192  }
193  else
194  {
195  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
196  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
197  #endif
198  MVT::Assign(*R1,*Y);
199  }
200 
201  // beta1 = sqrt(Y'*R1)
202  {
203  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
204  Teuchos::TimeMonitor lcltimer( *DotTime_ );
205  #endif
206  MVT::MvDot(*Y,*R1, beta1);
207 
208  for(size_t i=0; i<beta1.size(); i++)
209  beta1[i] = sqrt(beta1[i]);
210  }
211 
212  // beta = beta1
213  beta = beta1;
214 
215  // phibar = beta1
216  phibar = beta1;
217 
219  // Begin Lanczos iterations
220  for(int iter=1; iter <= maxIt_; iter++)
221  {
222  // Test convergence
223  indicesToRemove.clear();
224  for(int i=0; i<ncols; i++)
225  {
226  Scalar relres = phibar[i]/beta1[i];
227 // std::cout << "relative residual[" << unconvergedIndices[i] << "] at iteration " << iter << " = " << relres << std::endl;
228 
229  // If the vector converged, mark it for termination
230  // Make sure to add it to
231  if(relres < tolerances_[i])
232  {
233  indicesToRemove.push_back(i);
234  }
235  }
236 
237  // Check whether anything converged
238  newNumConverged = indicesToRemove.size();
239  if(newNumConverged > 0)
240  {
241  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
242  Teuchos::TimeMonitor lcltimer( *LockTime_ );
243  #endif
244 
245  // If something converged, stick the converged vectors in X
246  newlyConverged.resize(newNumConverged);
247  for(int i=0; i<newNumConverged; i++)
248  newlyConverged[i] = unconvergedIndices[indicesToRemove[i]];
249 
250  Teuchos::RCP<MV> helperLocX = MVT::CloneCopy(*locX,indicesToRemove);
251 
252  MVT::SetBlock(*helperLocX,newlyConverged,*X_);
253 
254  // If everything has converged, we are done
255  if(newNumConverged == ncols)
256  return;
257 
258  // Update unconverged indices
259  std::vector<int> helperVec(ncols - newNumConverged);
260 
261  std::set_difference(unconvergedIndices.begin(), unconvergedIndices.end(), newlyConverged.begin(), newlyConverged.end(), helperVec.begin());
262  unconvergedIndices = helperVec;
263 
264  // Get indices of things we want to keep
265  std::vector<int> thingsToKeep(ncols - newNumConverged);
266  helperVec.resize(ncols);
267  for(int i=0; i<ncols; i++)
268  helperVec[i] = i;
269  ncols = ncols - newNumConverged;
270 
271  std::set_difference(helperVec.begin(), helperVec.end(), indicesToRemove.begin(), indicesToRemove.end(), thingsToKeep.begin());
272 
273  // Shrink the multivectors
274  Teuchos::RCP<MV> helperMV;
275  helperMV = MVT::CloneCopy(*V,thingsToKeep);
276  V = helperMV;
277  helperMV = MVT::CloneCopy(*Y,thingsToKeep);
278  Y = helperMV;
279  helperMV = MVT::CloneCopy(*R1,thingsToKeep);
280  R1 = helperMV;
281  helperMV = MVT::CloneCopy(*R2,thingsToKeep);
282  R2 = helperMV;
283  helperMV = MVT::CloneCopy(*W,thingsToKeep);
284  W = helperMV;
285  helperMV = MVT::CloneCopy(*W1,thingsToKeep);
286  W1 = helperMV;
287  helperMV = MVT::CloneCopy(*W2,thingsToKeep);
288  W2 = helperMV;
289  helperMV = MVT::CloneCopy(*locX,thingsToKeep);
290  locX = helperMV;
291  scaleHelper = MVT::Clone(*V,ncols);
292 
293  // Shrink the arrays
294  alpha.resize(ncols);
295  oldeps.resize(ncols);
296  delta.resize(ncols);
297  gbar.resize(ncols);
298  phi.resize(ncols);
299  gamma.resize(ncols);
300  tmpvec.resize(ncols);
301  std::vector<Scalar> helperVecS(ncols);
302  for(int i=0; i<ncols; i++)
303  helperVecS[i] = beta[thingsToKeep[i]];
304  beta = helperVecS;
305  for(int i=0; i<ncols; i++)
306  helperVecS[i] = beta1[thingsToKeep[i]];
307  beta1 = helperVecS;
308  for(int i=0; i<ncols; i++)
309  helperVecS[i] = phibar[thingsToKeep[i]];
310  phibar = helperVecS;
311  for(int i=0; i<ncols; i++)
312  helperVecS[i] = oldBeta[thingsToKeep[i]];
313  oldBeta = helperVecS;
314  for(int i=0; i<ncols; i++)
315  helperVecS[i] = epsln[thingsToKeep[i]];
316  epsln = helperVecS;
317  for(int i=0; i<ncols; i++)
318  helperVecS[i] = cs[thingsToKeep[i]];
319  cs = helperVecS;
320  for(int i=0; i<ncols; i++)
321  helperVecS[i] = sn[thingsToKeep[i]];
322  sn = helperVecS;
323  for(int i=0; i<ncols; i++)
324  helperVecS[i] = dbar[thingsToKeep[i]];
325  dbar = helperVecS;
326 
327  // Tell operator about the removed indices
328  A_->removeIndices(indicesToRemove);
329  }
330 
331  // Normalize previous vector
332  // V = Y / beta
333  {
334  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
335  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
336  #endif
337  MVT::Assign(*Y,*V);
338  }
339  for(int i=0; i<ncols; i++)
340  tmpvec[i] = ONE/beta[i];
341  {
342  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
343  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
344  #endif
345  MVT::MvScale (*V, tmpvec);
346  }
347 
348  // Apply operator
349  // Y = AV
350  {
351  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
352  Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
353  #endif
354  OPT::Apply(*A_, *V, *Y);
355  }
356 
357  if(iter > 1)
358  {
359  // Y = Y - beta/oldBeta R1
360  for(int i=0; i<ncols; i++)
361  tmpvec[i] = beta[i]/oldBeta[i];
362  {
363  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
364  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
365  #endif
366  MVT::Assign(*R1,*scaleHelper);
367  }
368  {
369  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
370  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
371  #endif
372  MVT::MvScale(*scaleHelper,tmpvec);
373  }
374  {
375  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
376  Teuchos::TimeMonitor lcltimer( *AddTime_ );
377  #endif
378  MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
379  }
380  }
381 
382  // alpha = V'*Y
383  {
384  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
385  Teuchos::TimeMonitor lcltimer( *DotTime_ );
386  #endif
387  MVT::MvDot(*V, *Y, alpha);
388  }
389 
390  // Y = Y - alpha/beta R2
391  for(int i=0; i<ncols; i++)
392  tmpvec[i] = alpha[i]/beta[i];
393  {
394  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
395  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
396  #endif
397  MVT::Assign(*R2,*scaleHelper);
398  }
399  {
400  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
401  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
402  #endif
403  MVT::MvScale(*scaleHelper,tmpvec);
404  }
405  {
406  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
407  Teuchos::TimeMonitor lcltimer( *AddTime_ );
408  #endif
409  MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
410  }
411 
412  // R1 = R2
413  // R2 = Y
414  swapHelper = R1;
415  R1 = R2;
416  R2 = Y;
417  Y = swapHelper;
418 
419  // Y = M\R2
420  if(Prec_ != Teuchos::null)
421  {
422  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
423  Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
424  #endif
425 
426  OPT::Apply(*Prec_,*R2,*Y);
427  }
428  else
429  {
430  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
431  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
432  #endif
433  MVT::Assign(*R2,*Y);
434  }
435 
436  // Get new beta
437  // beta = sqrt(R2'*Y)
438  oldBeta = beta;
439  {
440  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
441  Teuchos::TimeMonitor lcltimer( *DotTime_ );
442  #endif
443  MVT::MvDot(*Y,*R2, beta);
444 
445  for(int i=0; i<ncols; i++)
446  beta[i] = sqrt(beta[i]);
447  }
448 
449  // Apply previous rotation
450  oldeps = epsln;
451  for(int i=0; i<ncols; i++)
452  {
453  delta[i] = cs[i]*dbar[i] + sn[i]*alpha[i];
454  gbar[i] = sn[i]*dbar[i] - cs[i]*alpha[i];
455  epsln[i] = sn[i]*beta[i];
456  dbar[i] = - cs[i]*beta[i];
457  }
458 
459  // Compute the next plane rotation
460  for(int i=0; i<ncols; i++)
461  {
462  symOrtho(gbar[i], beta[i], &cs[i], &sn[i], &gamma[i]);
463 
464  phi[i] = cs[i]*phibar[i];
465  phibar[i] = sn[i]*phibar[i];
466  }
467 
468  // w1 = w2
469  // w2 = w
470  swapHelper = W1;
471  W1 = W2;
472  W2 = W;
473  W = swapHelper;
474 
475  // W = (V - oldeps*W1 - delta*W2) / gamma
476  {
477  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
478  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
479  #endif
480  MVT::Assign(*W1,*scaleHelper);
481  }
482  {
483  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
484  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
485  #endif
486  MVT::MvScale(*scaleHelper,oldeps);
487  }
488  {
489  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
490  Teuchos::TimeMonitor lcltimer( *AddTime_ );
491  #endif
492  MVT::MvAddMv(ONE, *V, -ONE, *scaleHelper, *W);
493  }
494  {
495  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
496  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
497  #endif
498  MVT::Assign(*W2,*scaleHelper);
499  }
500  {
501  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
502  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
503  #endif
504  MVT::MvScale(*scaleHelper,delta);
505  }
506  {
507  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
508  Teuchos::TimeMonitor lcltimer( *AddTime_ );
509  #endif
510  MVT::MvAddMv(ONE, *W, -ONE, *scaleHelper, *W);
511  }
512  for(int i=0; i<ncols; i++)
513  tmpvec[i] = ONE/gamma[i];
514  {
515  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
516  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
517  #endif
518  MVT::MvScale(*W,tmpvec);
519  }
520 
521  // Update X
522  // X = X + phi*W
523  {
524  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
525  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
526  #endif
527  MVT::Assign(*W,*scaleHelper);
528  }
529  {
530  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
531  Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
532  #endif
533  MVT::MvScale(*scaleHelper,phi);
534  }
535  {
536  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
537  Teuchos::TimeMonitor lcltimer( *AddTime_ );
538  #endif
539  MVT::MvAddMv(ONE, *locX, ONE, *scaleHelper, *locX);
540  }
541  }
542 
543  // Stick unconverged vectors in X
544  {
545  #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
546  Teuchos::TimeMonitor lcltimer( *AssignTime_ );
547  #endif
548  MVT::SetBlock(*locX,unconvergedIndices,*X_);
549  }
550 }
551 
552 template <class Scalar, class MV, class OP>
553 void PseudoBlockMinres<Scalar,MV,OP>::symOrtho(Scalar a, Scalar b, Scalar *c, Scalar *s, Scalar *r)
554 {
555  const Scalar absA = std::abs(a);
556  const Scalar absB = std::abs(b);
557  if ( absB == ZERO ) {
558  *s = ZERO;
559  *r = absA;
560  if ( absA == ZERO )
561  *c = ONE;
562  else
563  *c = a / absA;
564  } else if ( absA == ZERO ) {
565  *c = ZERO;
566  *s = b / absB;
567  *r = absB;
568  } else if ( absB >= absA ) { // && a!=0 && b!=0
569  Scalar tau = a / b;
570  if ( b < ZERO )
571  *s = -ONE / sqrt( ONE+tau*tau );
572  else
573  *s = ONE / sqrt( ONE+tau*tau );
574  *c = *s * tau;
575  *r = b / *s;
576  } else { // (absA > absB) && a!=0 && b!=0
577  Scalar tau = b / a;
578  if ( a < ZERO )
579  *c = -ONE / sqrt( ONE+tau*tau );
580  else
581  *c = ONE / sqrt( ONE+tau*tau );
582  *s = *c * tau;
583  *r = a / *c;
584  }
585 }
586 
587 }} // End of namespace
588 
589 #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...