Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziMinres.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Anasazi: Block Eigensolvers Package
4//
5// Copyright 2004 NTESS and the Anasazi contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
14#ifndef ANASAZI_MINRES_HPP
15#define ANASAZI_MINRES_HPP
16
17#include "AnasaziConfigDefs.hpp"
18#include "Teuchos_TimeMonitor.hpp"
19
20namespace Anasazi {
21namespace Experimental {
22
23template <class Scalar, class MV, class OP>
24class PseudoBlockMinres
25{
26 typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
27 typedef Anasazi::OperatorTraits<Scalar,MV,OP> OPT;
28 const Scalar ONE;
29 const Scalar ZERO;
30
31public:
32 // Constructor
33 PseudoBlockMinres(Teuchos::RCP<OP> A, Teuchos::RCP<OP> Prec = Teuchos::null);
34
35 // Set tolerance and maximum iterations
36 void setTol(const std::vector<Scalar>& tolerances) { tolerances_ = tolerances; };
37 void setMaxIter(const int maxIt) { maxIt_ = maxIt; };
38
39 // Set solution and RHS
40 void setSol(Teuchos::RCP<MV> X) { X_ = X; };
41 void setRHS(Teuchos::RCP<const MV> B) { B_ = B; };
42
43 // Solve the linear system
44 void solve();
45
46private:
47 Teuchos::RCP<OP> A_;
48 Teuchos::RCP<OP> Prec_;
49 Teuchos::RCP<MV> X_;
50 Teuchos::RCP<const MV> B_;
51 std::vector<Scalar> tolerances_;
52 int maxIt_;
53
54 void symOrtho(Scalar a, Scalar b, Scalar *c, Scalar *s, Scalar *r);
55
56#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
57 Teuchos::RCP<Teuchos::Time> AddTime_, ApplyOpTime_, ApplyPrecTime_, AssignTime_, DotTime_, LockTime_, NormTime_, ScaleTime_, TotalTime_;
58#endif
59};
60
61
62
63template <class Scalar, class MV, class OP>
64PseudoBlockMinres<Scalar,MV,OP>::PseudoBlockMinres(Teuchos::RCP<OP> A, Teuchos::RCP<OP> Prec) :
65 ONE(Teuchos::ScalarTraits<Scalar>::one()),
66 ZERO(Teuchos::ScalarTraits<Scalar>::zero())
67{
68#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
69 AddTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Add Multivectors");
70 ApplyOpTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Apply Operator");
71 ApplyPrecTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Apply Preconditioner");
72 AssignTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Assignment (no locking)");
73 DotTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Compute Dot Product");
74 LockTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Lock Converged Vectors");
75 NormTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Compute Norm");
76 ScaleTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Scale Multivector");
77 TotalTime_ = Teuchos::TimeMonitor::getNewTimer("Anasazi: PseudoBlockMinres::Total Time");
78#endif
79
80 A_ = A;
81 Prec_ = Prec;
82 maxIt_ = 0;
83}
84
85
86
87template <class Scalar, class MV, class OP>
88void PseudoBlockMinres<Scalar,MV,OP>::solve()
89{
90 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
91 Teuchos::TimeMonitor outertimer( *TotalTime_ );
92 #endif
93
94 // Get number of vectors
95 int ncols = MVT::GetNumberVecs(*B_);
96 int newNumConverged;
97
98 // Declare some variables
99 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);
100 std::vector<int> indicesToRemove, newlyConverged, unconvergedIndices(ncols);
101 Teuchos::RCP<MV> V, Y, R1, R2, W, W1, W2, locX, scaleHelper, swapHelper;
102
103 // Allocate space for multivectors
104 V = MVT::Clone(*B_, ncols);
105 Y = MVT::Clone(*B_, ncols);
106 R1 = MVT::Clone(*B_, ncols);
107 R2 = MVT::Clone(*B_, ncols);
108 W = MVT::Clone(*B_, ncols);
109 W1 = MVT::Clone(*B_, ncols);
110 W2 = MVT::Clone(*B_, ncols);
111 scaleHelper = MVT::Clone(*B_, ncols);
112
113 // Reserve space for arrays
114 indicesToRemove.reserve(ncols);
115 newlyConverged.reserve(ncols);
116
117 // Initialize array of unconverged indices
118 for(int i=0; i<ncols; i++)
119 unconvergedIndices[i] = i;
120
121 // Get a local copy of X
122 // We want the vectors to remain contiguous even as things converge
123 locX = MVT::CloneCopy(*X_);
124
125 // Initialize residuals
126 // R1 = B - AX
127 {
128 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
129 Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
130 #endif
131 OPT::Apply(*A_,*locX,*R1);
132 }
133 {
134 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
135 Teuchos::TimeMonitor lcltimer( *AddTime_ );
136 #endif
137 MVT::MvAddMv(ONE, *B_, -ONE, *R1, *R1);
138 }
139
140 // R2 = R1
141 {
142 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
143 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
144 #endif
145 MVT::Assign(*R1,*R2);
146 }
147
148 // Initialize the W's to 0.
149 MVT::MvInit (*W);
150 MVT::MvInit (*W2);
151
152 // Y = M\R1 (preconditioned residual)
153 if(Prec_ != Teuchos::null)
154 {
155 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
156 Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
157 #endif
158
159 OPT::Apply(*Prec_,*R1,*Y);
160 }
161 else
162 {
163 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
164 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
165 #endif
166 MVT::Assign(*R1,*Y);
167 }
168
169 // beta1 = sqrt(Y'*R1)
170 {
171 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
172 Teuchos::TimeMonitor lcltimer( *DotTime_ );
173 #endif
174 MVT::MvDot(*Y,*R1, beta1);
175
176 for(size_t i=0; i<beta1.size(); i++)
177 beta1[i] = sqrt(beta1[i]);
178 }
179
180 // beta = beta1
181 beta = beta1;
182
183 // phibar = beta1
184 phibar = beta1;
185
187 // Begin Lanczos iterations
188 for(int iter=1; iter <= maxIt_; iter++)
189 {
190 // Test convergence
191 indicesToRemove.clear();
192 for(int i=0; i<ncols; i++)
193 {
194 Scalar relres = phibar[i]/beta1[i];
195// std::cout << "relative residual[" << unconvergedIndices[i] << "] at iteration " << iter << " = " << relres << std::endl;
196
197 // If the vector converged, mark it for termination
198 // Make sure to add it to
199 if(relres < tolerances_[i])
200 {
201 indicesToRemove.push_back(i);
202 }
203 }
204
205 // Check whether anything converged
206 newNumConverged = indicesToRemove.size();
207 if(newNumConverged > 0)
208 {
209 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
210 Teuchos::TimeMonitor lcltimer( *LockTime_ );
211 #endif
212
213 // If something converged, stick the converged vectors in X
214 newlyConverged.resize(newNumConverged);
215 for(int i=0; i<newNumConverged; i++)
216 newlyConverged[i] = unconvergedIndices[indicesToRemove[i]];
217
218 Teuchos::RCP<MV> helperLocX = MVT::CloneCopy(*locX,indicesToRemove);
219
220 MVT::SetBlock(*helperLocX,newlyConverged,*X_);
221
222 // If everything has converged, we are done
223 if(newNumConverged == ncols)
224 return;
225
226 // Update unconverged indices
227 std::vector<int> helperVec(ncols - newNumConverged);
228
229 std::set_difference(unconvergedIndices.begin(), unconvergedIndices.end(), newlyConverged.begin(), newlyConverged.end(), helperVec.begin());
230 unconvergedIndices = helperVec;
231
232 // Get indices of things we want to keep
233 std::vector<int> thingsToKeep(ncols - newNumConverged);
234 helperVec.resize(ncols);
235 for(int i=0; i<ncols; i++)
236 helperVec[i] = i;
237 ncols = ncols - newNumConverged;
238
239 std::set_difference(helperVec.begin(), helperVec.end(), indicesToRemove.begin(), indicesToRemove.end(), thingsToKeep.begin());
240
241 // Shrink the multivectors
242 Teuchos::RCP<MV> helperMV;
243 helperMV = MVT::CloneCopy(*V,thingsToKeep);
244 V = helperMV;
245 helperMV = MVT::CloneCopy(*Y,thingsToKeep);
246 Y = helperMV;
247 helperMV = MVT::CloneCopy(*R1,thingsToKeep);
248 R1 = helperMV;
249 helperMV = MVT::CloneCopy(*R2,thingsToKeep);
250 R2 = helperMV;
251 helperMV = MVT::CloneCopy(*W,thingsToKeep);
252 W = helperMV;
253 helperMV = MVT::CloneCopy(*W1,thingsToKeep);
254 W1 = helperMV;
255 helperMV = MVT::CloneCopy(*W2,thingsToKeep);
256 W2 = helperMV;
257 helperMV = MVT::CloneCopy(*locX,thingsToKeep);
258 locX = helperMV;
259 scaleHelper = MVT::Clone(*V,ncols);
260
261 // Shrink the arrays
262 alpha.resize(ncols);
263 oldeps.resize(ncols);
264 delta.resize(ncols);
265 gbar.resize(ncols);
266 phi.resize(ncols);
267 gamma.resize(ncols);
268 tmpvec.resize(ncols);
269 std::vector<Scalar> helperVecS(ncols);
270 for(int i=0; i<ncols; i++)
271 helperVecS[i] = beta[thingsToKeep[i]];
272 beta = helperVecS;
273 for(int i=0; i<ncols; i++)
274 helperVecS[i] = beta1[thingsToKeep[i]];
275 beta1 = helperVecS;
276 for(int i=0; i<ncols; i++)
277 helperVecS[i] = phibar[thingsToKeep[i]];
278 phibar = helperVecS;
279 for(int i=0; i<ncols; i++)
280 helperVecS[i] = oldBeta[thingsToKeep[i]];
281 oldBeta = helperVecS;
282 for(int i=0; i<ncols; i++)
283 helperVecS[i] = epsln[thingsToKeep[i]];
284 epsln = helperVecS;
285 for(int i=0; i<ncols; i++)
286 helperVecS[i] = cs[thingsToKeep[i]];
287 cs = helperVecS;
288 for(int i=0; i<ncols; i++)
289 helperVecS[i] = sn[thingsToKeep[i]];
290 sn = helperVecS;
291 for(int i=0; i<ncols; i++)
292 helperVecS[i] = dbar[thingsToKeep[i]];
293 dbar = helperVecS;
294
295 // Tell operator about the removed indices
296 A_->removeIndices(indicesToRemove);
297 }
298
299 // Normalize previous vector
300 // V = Y / beta
301 {
302 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
303 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
304 #endif
305 MVT::Assign(*Y,*V);
306 }
307 for(int i=0; i<ncols; i++)
308 tmpvec[i] = ONE/beta[i];
309 {
310 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
311 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
312 #endif
313 MVT::MvScale (*V, tmpvec);
314 }
315
316 // Apply operator
317 // Y = AV
318 {
319 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
320 Teuchos::TimeMonitor lcltimer( *ApplyOpTime_ );
321 #endif
322 OPT::Apply(*A_, *V, *Y);
323 }
324
325 if(iter > 1)
326 {
327 // Y = Y - beta/oldBeta R1
328 for(int i=0; i<ncols; i++)
329 tmpvec[i] = beta[i]/oldBeta[i];
330 {
331 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
332 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
333 #endif
334 MVT::Assign(*R1,*scaleHelper);
335 }
336 {
337 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
338 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
339 #endif
340 MVT::MvScale(*scaleHelper,tmpvec);
341 }
342 {
343 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
344 Teuchos::TimeMonitor lcltimer( *AddTime_ );
345 #endif
346 MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
347 }
348 }
349
350 // alpha = V'*Y
351 {
352 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
353 Teuchos::TimeMonitor lcltimer( *DotTime_ );
354 #endif
355 MVT::MvDot(*V, *Y, alpha);
356 }
357
358 // Y = Y - alpha/beta R2
359 for(int i=0; i<ncols; i++)
360 tmpvec[i] = alpha[i]/beta[i];
361 {
362 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
363 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
364 #endif
365 MVT::Assign(*R2,*scaleHelper);
366 }
367 {
368 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
369 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
370 #endif
371 MVT::MvScale(*scaleHelper,tmpvec);
372 }
373 {
374 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
375 Teuchos::TimeMonitor lcltimer( *AddTime_ );
376 #endif
377 MVT::MvAddMv(ONE, *Y, -ONE, *scaleHelper, *Y);
378 }
379
380 // R1 = R2
381 // R2 = Y
382 swapHelper = R1;
383 R1 = R2;
384 R2 = Y;
385 Y = swapHelper;
386
387 // Y = M\R2
388 if(Prec_ != Teuchos::null)
389 {
390 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
391 Teuchos::TimeMonitor lcltimer( *ApplyPrecTime_ );
392 #endif
393
394 OPT::Apply(*Prec_,*R2,*Y);
395 }
396 else
397 {
398 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
399 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
400 #endif
401 MVT::Assign(*R2,*Y);
402 }
403
404 // Get new beta
405 // beta = sqrt(R2'*Y)
406 oldBeta = beta;
407 {
408 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
409 Teuchos::TimeMonitor lcltimer( *DotTime_ );
410 #endif
411 MVT::MvDot(*Y,*R2, beta);
412
413 for(int i=0; i<ncols; i++)
414 beta[i] = sqrt(beta[i]);
415 }
416
417 // Apply previous rotation
418 oldeps = epsln;
419 for(int i=0; i<ncols; i++)
420 {
421 delta[i] = cs[i]*dbar[i] + sn[i]*alpha[i];
422 gbar[i] = sn[i]*dbar[i] - cs[i]*alpha[i];
423 epsln[i] = sn[i]*beta[i];
424 dbar[i] = - cs[i]*beta[i];
425 }
426
427 // Compute the next plane rotation
428 for(int i=0; i<ncols; i++)
429 {
430 symOrtho(gbar[i], beta[i], &cs[i], &sn[i], &gamma[i]);
431
432 phi[i] = cs[i]*phibar[i];
433 phibar[i] = sn[i]*phibar[i];
434 }
435
436 // w1 = w2
437 // w2 = w
438 swapHelper = W1;
439 W1 = W2;
440 W2 = W;
441 W = swapHelper;
442
443 // W = (V - oldeps*W1 - delta*W2) / gamma
444 {
445 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
446 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
447 #endif
448 MVT::Assign(*W1,*scaleHelper);
449 }
450 {
451 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
452 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
453 #endif
454 MVT::MvScale(*scaleHelper,oldeps);
455 }
456 {
457 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
458 Teuchos::TimeMonitor lcltimer( *AddTime_ );
459 #endif
460 MVT::MvAddMv(ONE, *V, -ONE, *scaleHelper, *W);
461 }
462 {
463 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
464 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
465 #endif
466 MVT::Assign(*W2,*scaleHelper);
467 }
468 {
469 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
470 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
471 #endif
472 MVT::MvScale(*scaleHelper,delta);
473 }
474 {
475 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
476 Teuchos::TimeMonitor lcltimer( *AddTime_ );
477 #endif
478 MVT::MvAddMv(ONE, *W, -ONE, *scaleHelper, *W);
479 }
480 for(int i=0; i<ncols; i++)
481 tmpvec[i] = ONE/gamma[i];
482 {
483 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
484 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
485 #endif
486 MVT::MvScale(*W,tmpvec);
487 }
488
489 // Update X
490 // X = X + phi*W
491 {
492 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
493 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
494 #endif
495 MVT::Assign(*W,*scaleHelper);
496 }
497 {
498 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
499 Teuchos::TimeMonitor lcltimer( *ScaleTime_ );
500 #endif
501 MVT::MvScale(*scaleHelper,phi);
502 }
503 {
504 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
505 Teuchos::TimeMonitor lcltimer( *AddTime_ );
506 #endif
507 MVT::MvAddMv(ONE, *locX, ONE, *scaleHelper, *locX);
508 }
509 }
510
511 // Stick unconverged vectors in X
512 {
513 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
514 Teuchos::TimeMonitor lcltimer( *AssignTime_ );
515 #endif
516 MVT::SetBlock(*locX,unconvergedIndices,*X_);
517 }
518}
519
520template <class Scalar, class MV, class OP>
521void PseudoBlockMinres<Scalar,MV,OP>::symOrtho(Scalar a, Scalar b, Scalar *c, Scalar *s, Scalar *r)
522{
523 const Scalar absA = std::abs(a);
524 const Scalar absB = std::abs(b);
525 if ( absB == ZERO ) {
526 *s = ZERO;
527 *r = absA;
528 if ( absA == ZERO )
529 *c = ONE;
530 else
531 *c = a / absA;
532 } else if ( absA == ZERO ) {
533 *c = ZERO;
534 *s = b / absB;
535 *r = absB;
536 } else if ( absB >= absA ) { // && a!=0 && b!=0
537 Scalar tau = a / b;
538 if ( b < ZERO )
539 *s = -ONE / sqrt( ONE+tau*tau );
540 else
541 *s = ONE / sqrt( ONE+tau*tau );
542 *c = *s * tau;
543 *r = b / *s;
544 } else { // (absA > absB) && a!=0 && b!=0
545 Scalar tau = b / a;
546 if ( a < ZERO )
547 *c = -ONE / sqrt( ONE+tau*tau );
548 else
549 *c = ONE / sqrt( ONE+tau*tau );
550 *s = *c * tau;
551 *r = a / *c;
552 }
553}
554
555}} // End of namespace
556
557#endif
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Anasazi::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition AnasaziMultiVecTraits.hpp:95
Anasazi::OperatorTraits
Virtual base class which defines basic traits for the operator type.
Definition AnasaziOperatorTraits.hpp:53
Anasazi
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Generated by doxygen 1.9.8