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
45
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{
58 typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
59 typedef Anasazi::OperatorTraits<Scalar,MV,OP> OPT;
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
AnasaziConfigDefs.hpp
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Anasazi
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
Generated by doxygen 1.13.2