/usr/include/trilinos/Belos_StatusTest_GenResNorm_MP_Vector.hpp is in libtrilinos-stokhos-dev 12.10.1-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 | //@HEADER
// ************************************************************************
//
// Belos: Block Linear Solvers Package
// Copyright 2004 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER
#ifndef BELOS_STATUS_TEST_GEN_RESNORM_MP_VECTOR_H
#define BELOS_STATUS_TEST_GEN_RESNORM_MP_VECTOR_H
/*!
\file Belos_StatusTest_GenResNorm_MP_Vector.hpp
\brief Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Specialized for Sacado::MP::Vector scalar type to track how many iterations
each component of ensemble takes.
*/
#include "Stokhos_Sacado_Kokkos_MP_Vector.hpp"
#include "BelosStatusTestGenResNorm.hpp"
/*!
\class Belos::StatusTestGenResNorm
\brief An implementation of StatusTestResNorm using a family of residual norms.
StatusTestGenResNorm is an implementation of StatusTestResNorm that allows a user to construct
one of a family of residual tests for use as a status/convergence test for Belos.
The form of the test is
\f[
\frac{\|r_i\|}{\sigma_i} \le \tau
\f]
where
<ul>
<li> \f$r_i\f$ is the i-th residual std::vector, implicitly or explicitly computed (determined by enum ResType),
<li> \f$\|r_i\|\f$ is the i-th residual norm determined by the enum NormType (1-norm, 2-norm or inf-norm),
<li> \f$\sigma_i\f$ is the i-th scale factor that can be passed in as a precomputed number of the templated type,
or can be selected from by the enum ScaleType (norm of RHS, norm of initial residual).
<li> \f$\tau\f$ is the tolerance that is passed in as a number of the templated type to the constructor.
The value of \f$\tau\f$ can be reset using the setTolerance() method.
</ul>
This implementation is the same as the default, but is specialized for
the Sacado::MP::Vector scalar type to keep track of when each component
of the ensemble converges.
*/
namespace Belos {
template <class Storage, class MV, class OP>
class StatusTestGenResNorm<Sacado::MP::Vector<Storage>, MV, OP> :
public StatusTestResNorm<Sacado::MP::Vector<Storage>,MV,OP> {
public:
// Convenience typedefs
typedef Sacado::MP::Vector<Storage> ScalarType;
typedef Teuchos::ScalarTraits<ScalarType> SCT;
typedef typename SCT::magnitudeType MagnitudeType;
typedef MultiVecTraits<ScalarType,MV> MVT;
//! @name Enums.
//@{
/*!
\brief Select how the residual std::vector is produced.
*/
enum ResType {Implicit, /*!< Use the residual std::vector produced by the iterative solver. */
Explicit /*!< Explicitly compute the residual std::vector r = b - A*x using the
linear problem. */
};
//@}
//! @name Constructors/destructors.
//@{
//! Constructor
/*! The constructor takes a single argument specifying the tolerance (\f$\tau\f$).
If none of the form definition methods are called, we use \f$\|r\|_2/\|r^{(0)}\|_2 \le \tau\f$
as the stopping criterion, where \f$\|r\|_2\f$ uses the least costly form of the 2-norm of
residual available from the iterative method and \f$\|r^{(0)}\|_2\f$ is the corresponding norm
of the initial residual. The least costly form of the 2-norm depends on the chosen iterative
method. Most Krylov methods produce the preconditioned residual std::vector in a form that would be
exact in infinite precision arithmetic. This std::vector may be different from the true residual
either because left scaling or preconditioning was used, or because round-off error has
introduced significant error, or both.
You can also state the number of vectors that must pass the convergence criteria before the
status test passes by using the \c quorum argument.
*/
StatusTestGenResNorm( MagnitudeType Tolerance, int quorum = -1, bool showMaxResNormOnly = false );
//! Destructor
virtual ~StatusTestGenResNorm();
//@}
//! @name Form and parameter definition methods.
//@{
//! Define form of the residual, its norm and optional weighting std::vector.
/*! This method defines the form of \f$\|r\|\f$. We specify:
<ul>
<li> Whether the residual std::vector should be explicitly computed, or taken from the iterative method.
<li> The norm to be used on the residual (this may be different than the norm used in
DefineScaleForm()).
</ul>
*/
int defineResForm( ResType TypeOfResidual, NormType TypeOfNorm);
//! Define form of the scaling, its norm, its optional weighting std::vector, or, alternatively, define an explicit value.
/*! This method defines the form of how the residual is scaled (if at all). It operates in two modes:
<ol>
<li> User-provided scaling value:
<ul>
<li> Set argument TypeOfScaling to UserProvided.
<li> Set ScaleValue to a non-zero value that the residual norm will be divided by.
<li> TypeOfNorm argument will be ignored.
<li> Sample use: Define ScaleValue = \f$\|A\|_{\infty}\f$ where \f$ A \f$ is the matrix
of the linear problem.
</ul>
<li> Use a supported Scaling Form:
<ul>
<li> Define TypeOfScaling to be the norm of the right hand side, the initial residual std::vector,
or to none.
<li> Define norm to be used on the scaling std::vector (this may be different than the norm used
in DefineResForm()).
</ul>
</ol>
*/
int defineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
//! Set the value of the tolerance
/*! We allow the tolerance to be reset for cases where, in the process of testing the residual,
we find that the initial tolerance was too tight or too lax.
*/
int setTolerance(MagnitudeType tolerance) {tolerance_ = tolerance; return(0);}
//! Sets the number of residuals that must pass the convergence test before Passed is returned.
//! \note If \c quorum=-1 then all residuals must pass the convergence test before Passed is returned.
int setQuorum(int quorum) {quorum_ = quorum; return(0);}
//! Set whether the only maximum residual norm is displayed when the print() method is called
int setShowMaxResNormOnly(bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly; return(0);}
//@}
//! @name Status methods
//@{
//! Check convergence status: Passed, Failed, or Undefined.
/*! This method checks to see if the convergence criteria are met.
Depending on how the residual test is constructed this method will return
the appropriate status type.
\return StatusType: Passed, Failed, or Undefined.
*/
StatusType checkStatus(Iteration<ScalarType,MV,OP>* iSolver);
//! Return the result of the most recent CheckStatus call.
StatusType getStatus() const {return(status_);};
//@}
//! @name Reset methods
//@{
//! Resets the internal configuration to the initial state.
void reset();
//@}
//! @name Print methods
//@{
//! Output formatted description of stopping test to output stream.
void print(std::ostream& os, int indent = 0) const;
//! Print message for each status specific to this stopping test.
void printStatus(std::ostream& os, StatusType type) const;
//@}
//! @name Methods to access data members.
//@{
//! Returns the current solution estimate that was computed for the most recent residual test.
//! \note This is useful for explicit residual tests, if this test is an implicit residual test
//! a null pointer will be returned.
Teuchos::RCP<MV> getSolution() { if (restype_==Implicit) { return Teuchos::null; } else { return curSoln_; } }
//! Returns the number of residuals that must pass the convergence test before Passed is returned.
//! \note If \c quorum=-1 then all residuals must pass the convergence test before Passed is returned.
int getQuorum() const { return quorum_; }
//! Returns whether the only maximum residual norm is displayed when the print() method is called
bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
//! Returns the std::vector containing the indices of the residuals that passed the test.
std::vector<int> convIndices() { return ind_; }
//! Returns the value of the tolerance, \f$ \tau \f$, set in the constructor.
MagnitudeType getTolerance() const {return(tolerance_);};
//! Returns the test value, \f$ \frac{\|r\|}{\sigma} \f$, computed in most recent call to CheckStatus.
const std::vector<MagnitudeType>* getTestValue() const {return(&testvector_);};
//! Returns the residual norm value, \f$ \|r\| \f$, computed in most recent call to CheckStatus.
const std::vector<MagnitudeType>* getResNormValue() const {return(&resvector_);};
//! Returns the scaled norm value, \f$ \sigma \f$.
const std::vector<MagnitudeType>* getScaledNormValue() const {return(&scalevector_);};
//! Returns a boolean indicating a loss of accuracy has been detected in computing the residual.
//! \note This status test does not check for loss of accuracy, so this method will always return false.
bool getLOADetected() const { return false; }
//! Returns number of ensemble iterations
const std::vector<int> getEnsembleIterations() const { return ensemble_iterations; }
//@}
/** @name Misc. */
//@{
/** \brief Call to setup initial scaling std::vector.
*
* After this function is called <tt>getScaledNormValue()</tt> can be called
* to get the scaling std::vector.
*/
StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);
//@}
/** \name Overridden from Teuchos::Describable */
//@{
/** \brief Method to return description of the maximum iteration status test */
std::string description() const
{
std::ostringstream oss;
oss << "Belos::StatusTestGenResNorm<>: " << resFormStr();
oss << ", tol = " << tolerance_;
return oss.str();
}
//@}
protected:
private:
//! @name Private methods.
//@{
/** \brief Description of current residual form */
std::string resFormStr() const
{
std::ostringstream oss;
oss << "(";
oss << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
oss << ((restype_==Explicit) ? " Exp" : " Imp");
oss << " Res Vec) ";
// If there is no residual scaling, return current string.
if (scaletype_!=None)
{
// Insert division sign.
oss << "/ ";
// Determine output string for scaling, if there is any.
if (scaletype_==UserProvided)
oss << " (User Scale)";
else {
oss << "(";
oss << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
if (scaletype_==NormOfInitRes)
oss << " Res0";
else if (scaletype_==NormOfPrecInitRes)
oss << " Prec Res0";
else
oss << " RHS ";
oss << ")";
}
}
return oss.str();
}
//@}
//! @name Private data members.
//@{
//! Tolerance used to determine convergence
MagnitudeType tolerance_;
//! Number of residuals that must pass the convergence test before Passed is returned.
int quorum_;
//! Determines if the entries for all of the residuals are shown or just the max.
bool showMaxResNormOnly_;
//! Type of residual to use (explicit or implicit)
ResType restype_;
//! Type of norm to use on residual (OneNorm, TwoNorm, or InfNorm).
NormType resnormtype_;
//! Type of scaling to use (Norm of RHS, Norm of Initial Residual, None or User provided)
ScaleType scaletype_;
//! Type of norm to use on the scaling (OneNorm, TwoNorm, or InfNorm)
NormType scalenormtype_;
//! Scaling value.
MagnitudeType scalevalue_;
//! Scaling std::vector.
std::vector<MagnitudeType> scalevector_;
//! Residual norm std::vector.
std::vector<MagnitudeType> resvector_;
//! Test std::vector = resvector_ / scalevector_
std::vector<MagnitudeType> testvector_;
//! Vector containing the indices for the vectors that passed the test.
std::vector<int> ind_;
//! Most recent solution vector used by this status test.
Teuchos::RCP<MV> curSoln_;
//! Status
StatusType status_;
//! The current blocksize of the linear system being solved.
int curBlksz_;
//! The current number of right-hand sides being solved for.
int curNumRHS_;
//! The indices of the current number of right-hand sides being solved for.
std::vector<int> curLSIdx_;
//! The current number of linear systems that have been loaded into the linear problem.
int curLSNum_;
//! The total number of right-hand sides being solved for.
int numrhs_;
//! Is this the first time CheckStatus is called?
bool firstcallCheckStatus_;
//! Is this the first time DefineResForm is called?
bool firstcallDefineResForm_;
//! Is this the first time DefineScaleForm is called?
bool firstcallDefineScaleForm_;
//! Which ensemble components have converged
std::vector<int> ensemble_converged;
//! The number of iterations at which point each ensemble component converges
std::vector<int> ensemble_iterations;
//@}
};
template <class StorageType, class MV, class OP>
StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::
StatusTestGenResNorm (MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly)
: tolerance_(Tolerance),
quorum_(quorum),
showMaxResNormOnly_(showMaxResNormOnly),
restype_(Implicit),
resnormtype_(TwoNorm),
scaletype_(NormOfInitRes),
scalenormtype_(TwoNorm),
scalevalue_(Teuchos::ScalarTraits<MagnitudeType>::one ()),
status_(Undefined),
curBlksz_(0),
curNumRHS_(0),
curLSNum_(0),
numrhs_(0),
firstcallCheckStatus_(true),
firstcallDefineResForm_(true),
firstcallDefineScaleForm_(true),
ensemble_converged(StorageType::static_size, 0),
ensemble_iterations(StorageType::static_size, 0)
{
// This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
// the implicit residual std::vector.
}
template <class StorageType, class MV, class OP>
StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::~StatusTestGenResNorm()
{}
template <class StorageType, class MV, class OP>
void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::reset()
{
status_ = Undefined;
curBlksz_ = 0;
curLSNum_ = 0;
curLSIdx_.resize(0);
numrhs_ = 0;
ind_.resize(0);
firstcallCheckStatus_ = true;
curSoln_ = Teuchos::null;
ensemble_converged = std::vector<int>(StorageType::static_size, 0);
ensemble_iterations = std::vector<int>(StorageType::static_size, 0);
}
template <class StorageType, class MV, class OP>
int StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::defineResForm( ResType TypeOfResidual, NormType TypeOfNorm )
{
TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
"StatusTestGenResNorm::defineResForm(): The residual form has already been defined.");
firstcallDefineResForm_ = false;
restype_ = TypeOfResidual;
resnormtype_ = TypeOfNorm;
return(0);
}
template <class StorageType, class MV, class OP>
int StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
MagnitudeType ScaleValue )
{
TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
"StatusTestGenResNorm::defineScaleForm(): The scaling type has already been defined.");
firstcallDefineScaleForm_ = false;
scaletype_ = TypeOfScaling;
scalenormtype_ = TypeOfNorm;
scalevalue_ = ScaleValue;
return(0);
}
template <class StorageType, class MV, class OP>
StatusType StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::checkStatus( Iteration<ScalarType,MV,OP>* iSolver )
{
MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
// Compute scaling term (done once for each block that's being solved)
if (firstcallCheckStatus_) {
StatusType status = firstCallCheckStatusSetup(iSolver);
if(status==Failed) {
status_ = Failed;
return(status_);
}
}
//
// This section computes the norm of the residual std::vector
//
if ( curLSNum_ != lp.getLSNumber() ) {
//
// We have moved on to the next rhs block
//
curLSNum_ = lp.getLSNumber();
curLSIdx_ = lp.getLSIndex();
curBlksz_ = (int)curLSIdx_.size();
int validLS = 0;
for (int i=0; i<curBlksz_; ++i) {
if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
validLS++;
}
curNumRHS_ = validLS;
curSoln_ = Teuchos::null;
//
} else {
//
// We are in the same rhs block, return if we are converged
//
if (status_==Passed) { return status_; }
}
if (restype_==Implicit) {
//
// get the native residual norms from the solver for this block of right-hand sides.
// If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
// Otherwise the native residual is assumed to be stored in the resvector_.
//
std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
Teuchos::RCP<const MV> residMV = iSolver->getNativeResiduals( &tmp_resvector );
if ( residMV != Teuchos::null ) {
tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
MVT::MvNorm( *residMV, tmp_resvector, resnormtype_ );
typename std::vector<int>::iterator p = curLSIdx_.begin();
for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
// Check if this index is valid
if (*p != -1)
resvector_[*p] = tmp_resvector[i];
}
} else {
typename std::vector<int>::iterator p = curLSIdx_.begin();
for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
// Check if this index is valid
if (*p != -1)
resvector_[*p] = tmp_resvector[i];
}
}
}
else if (restype_==Explicit) {
//
// Request the true residual for this block of right-hand sides.
//
Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
curSoln_ = lp.updateSolution( cur_update );
Teuchos::RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_ ) );
lp.computeCurrResVec( &*cur_res, &*curSoln_ );
std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( *cur_res ) );
MVT::MvNorm( *cur_res, tmp_resvector, resnormtype_ );
typename std::vector<int>::iterator p = curLSIdx_.begin();
for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
// Check if this index is valid
if (*p != -1)
resvector_[*p] = tmp_resvector[i];
}
}
//
// Compute the new linear system residuals for testing.
// (if any of them don't meet the tolerance or are NaN, then we exit with that status)
//
if ( scalevector_.size() > 0 ) {
typename std::vector<int>::iterator p = curLSIdx_.begin();
for (; p<curLSIdx_.end(); ++p) {
// Check if this index is valid
if (*p != -1) {
// Scale the std::vector accordingly
if ( scalevector_[ *p ] != zero ) {
// Don't intentionally divide by zero.
testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
} else {
testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
}
}
}
}
else {
typename std::vector<int>::iterator p = curLSIdx_.begin();
for (; p<curLSIdx_.end(); ++p) {
// Check if this index is valid
if (*p != -1)
testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
}
}
// Check status of new linear system residuals and see if we have the quorum.
int have = 0;
ind_.resize( curLSIdx_.size() );
typename std::vector<int>::iterator p = curLSIdx_.begin();
for (; p<curLSIdx_.end(); ++p) {
// Check if this index is valid
if (*p != -1) {
// Check if any of the residuals are larger than the tolerance.
const int ensemble_size = StorageType::static_size;
bool all_converged = true;
for (int i=0; i<ensemble_size; ++i) {
if (!ensemble_converged[i]) {
if (testvector_[ *p ].coeff(i) > tolerance_) {
++ensemble_iterations[i];
all_converged = false;
}
else if (testvector_[ *p ].coeff(i) <= tolerance_) {
ensemble_converged[i] = 1;
}
else {
// Throw an std::exception if a NaN is found.
status_ = Failed;
TEUCHOS_TEST_FOR_EXCEPTION(true,StatusTestError,"StatusTestGenResNorm::checkStatus(): NaN has been detected.");
}
}
}
if (all_converged) {
ind_[have] = *p;
have++;
}
}
}
ind_.resize(have);
int need = (quorum_ == -1) ? curNumRHS_: quorum_;
status_ = (have >= need) ? Passed : Failed;
// Return the current status
return status_;
}
template <class StorageType, class MV, class OP>
void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::print(std::ostream& os, int indent) const
{
for (int j = 0; j < indent; j ++)
os << ' ';
printStatus(os, status_);
os << resFormStr();
if (status_==Undefined)
os << ", tol = " << tolerance_ << std::endl;
else {
os << std::endl;
if(showMaxResNormOnly_ && curBlksz_ > 1) {
const MagnitudeType maxRelRes = *std::max_element(
testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
);
for (int j = 0; j < indent + 13; j ++)
os << ' ';
os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
<< ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::endl;
}
else {
for ( int i=0; i<numrhs_; i++ ) {
for (int j = 0; j < indent + 13; j ++)
os << ' ';
os << "residual [ " << i << " ] = " << testvector_[ i ];
os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " " ) << tolerance_ << std::endl;
}
}
}
os << std::endl;
}
template <class StorageType, class MV, class OP>
void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::printStatus(std::ostream& os, StatusType type) const
{
os << std::left << std::setw(13) << std::setfill('.');
switch (type) {
case Passed:
os << "Converged";
break;
case Failed:
os << "Unconverged";
break;
case Undefined:
default:
os << "**";
break;
}
os << std::left << std::setfill(' ');
return;
}
template <class StorageType, class MV, class OP>
StatusType StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::firstCallCheckStatusSetup( Iteration<ScalarType,MV,OP>* iSolver )
{
int i;
MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
// Compute scaling term (done once for each block that's being solved)
if (firstcallCheckStatus_) {
//
// Get some current solver information.
//
firstcallCheckStatus_ = false;
if (scaletype_== NormOfRHS) {
Teuchos::RCP<const MV> rhs = lp.getRHS();
numrhs_ = MVT::GetNumberVecs( *rhs );
scalevector_.resize( numrhs_ );
MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
}
else if (scaletype_==NormOfInitRes) {
Teuchos::RCP<const MV> init_res = lp.getInitResVec();
numrhs_ = MVT::GetNumberVecs( *init_res );
scalevector_.resize( numrhs_ );
MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
}
else if (scaletype_==NormOfPrecInitRes) {
Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
numrhs_ = MVT::GetNumberVecs( *init_res );
scalevector_.resize( numrhs_ );
MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
}
else {
numrhs_ = MVT::GetNumberVecs( *(lp.getRHS()) );
}
resvector_.resize( numrhs_ );
testvector_.resize( numrhs_ );
curLSNum_ = lp.getLSNumber();
curLSIdx_ = lp.getLSIndex();
curBlksz_ = (int)curLSIdx_.size();
int validLS = 0;
for (i=0; i<curBlksz_; ++i) {
if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
validLS++;
}
curNumRHS_ = validLS;
//
// Initialize the testvector.
for (i=0; i<numrhs_; i++) { testvector_[i] = one; }
// Return an error if the scaling is zero.
if (scalevalue_ == zero) {
return Failed;
}
}
return Undefined;
}
} // end namespace Belos
#endif /* BELOS_STATUS_TEST_RESNORM_H */
|