Anasazi Version of the Day
Loading...
Searching...
No Matches
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
52namespace Anasazi {
53namespace Experimental {
54
55template <class Scalar, class MV, class OP>
56class PseudoBlockMinres
57{
60 const Scalar ONE;
61 const Scalar ZERO;
62
63public:
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
78private:
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
95template <class Scalar, class MV, class OP>
96PseudoBlockMinres<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
119template <class Scalar, class MV, class OP>
120void 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
552template <class Scalar, class MV, class OP>
553void 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
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Traits class which defines basic operations on multivectors.
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.