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 //
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:
80  Teuchos::RCP<OP> Prec_;
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.
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...