|
Belos
Version of the Day
|
|
Belos
Version of the Day
|
00001 //@HEADER 00002 // ************************************************************************ 00003 // 00004 // Belos: Block Linear Solvers Package 00005 // Copyright 2004 Sandia Corporation 00006 // 00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, 00008 // the U.S. Government retains certain rights in this software. 00009 // 00010 // Redistribution and use in source and binary forms, with or without 00011 // modification, are permitted provided that the following conditions are 00012 // met: 00013 // 00014 // 1. Redistributions of source code must retain the above copyright 00015 // notice, this list of conditions and the following disclaimer. 00016 // 00017 // 2. Redistributions in binary form must reproduce the above copyright 00018 // notice, this list of conditions and the following disclaimer in the 00019 // documentation and/or other materials provided with the distribution. 00020 // 00021 // 3. Neither the name of the Corporation nor the names of the 00022 // contributors may be used to endorse or promote products derived from 00023 // this software without specific prior written permission. 00024 // 00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00036 // 00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 00038 // 00039 // ************************************************************************ 00040 //@HEADER 00041 00042 #ifndef BELOS_RCG_SOLMGR_HPP 00043 #define BELOS_RCG_SOLMGR_HPP 00044 00049 #include "BelosConfigDefs.hpp" 00050 #include "BelosTypes.hpp" 00051 00052 #include "BelosLinearProblem.hpp" 00053 #include "BelosSolverManager.hpp" 00054 00055 #include "BelosRCGIter.hpp" 00056 #include "BelosDGKSOrthoManager.hpp" 00057 #include "BelosICGSOrthoManager.hpp" 00058 #include "BelosIMGSOrthoManager.hpp" 00059 #include "BelosStatusTestMaxIters.hpp" 00060 #include "BelosStatusTestGenResNorm.hpp" 00061 #include "BelosStatusTestCombo.hpp" 00062 #include "BelosStatusTestOutputFactory.hpp" 00063 #include "BelosOutputManager.hpp" 00064 #include "Teuchos_BLAS.hpp" 00065 #include "Teuchos_LAPACK.hpp" 00066 #include "Teuchos_as.hpp" 00067 #ifdef BELOS_TEUCHOS_TIME_MONITOR 00068 #include "Teuchos_TimeMonitor.hpp" 00069 #endif 00070 00110 namespace Belos { 00111 00113 00114 00121 class RCGSolMgrLinearProblemFailure : public BelosError {public: 00122 RCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg) 00123 {}}; 00124 00131 class RCGSolMgrLAPACKFailure : public BelosError {public: 00132 RCGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg) 00133 {}}; 00134 00141 class RCGSolMgrRecyclingFailure : public BelosError {public: 00142 RCGSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg) 00143 {}}; 00144 00146 00147 00148 // Partial specialization for complex ScalarType. 00149 // This contains a trivial implementation. 00150 // See discussion in the class documentation above. 00151 template<class ScalarType, class MV, class OP, 00152 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex> 00153 class RCGSolMgr : 00154 public Details::RealSolverManager<ScalarType, MV, OP, 00155 Teuchos::ScalarTraits<ScalarType>::isComplex> 00156 { 00157 static const bool isComplex = Teuchos::ScalarTraits<ScalarType>::isComplex; 00158 typedef Details::RealSolverManager<ScalarType, MV, OP, isComplex> base_type; 00159 00160 public: 00161 RCGSolMgr () : 00162 base_type () 00163 {} 00164 RCGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 00165 const Teuchos::RCP<Teuchos::ParameterList> &pl) : 00166 base_type () 00167 {} 00168 virtual ~RCGSolMgr () {} 00169 }; 00170 00171 00172 // Partial specialization for real ScalarType. 00173 // This contains the actual working implementation of RCG. 00174 // See discussion in the class documentation above. 00175 template<class ScalarType, class MV, class OP> 00176 class RCGSolMgr<ScalarType, MV, OP, false> : 00177 public Details::RealSolverManager<ScalarType, MV, OP, false> { 00178 private: 00179 typedef MultiVecTraits<ScalarType,MV> MVT; 00180 typedef MultiVecTraitsExt<ScalarType,MV> MVText; 00181 typedef OperatorTraits<ScalarType,MV,OP> OPT; 00182 typedef Teuchos::ScalarTraits<ScalarType> SCT; 00183 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType; 00184 typedef Teuchos::ScalarTraits<MagnitudeType> MT; 00185 00186 public: 00187 00189 00190 00196 RCGSolMgr(); 00197 00219 RCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 00220 const Teuchos::RCP<Teuchos::ParameterList> &pl ); 00221 00223 virtual ~RCGSolMgr() {}; 00225 00227 00228 00229 const LinearProblem<ScalarType,MV,OP>& getProblem() const { 00230 return *problem_; 00231 } 00232 00234 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const; 00235 00237 Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; } 00238 00244 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const { 00245 return Teuchos::tuple(timerSolve_); 00246 } 00247 00252 MagnitudeType achievedTol() const { 00253 return achievedTol_; 00254 } 00255 00257 int getNumIters() const { 00258 return numIters_; 00259 } 00260 00262 bool isLOADetected() const { return false; } 00263 00265 00267 00268 00270 void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; } 00271 00273 void setParameters( const Teuchos::RCP<Teuchos::ParameterList> ¶ms ); 00274 00276 00278 00279 00285 void reset( const ResetType type ) { 00286 if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); 00287 else if (type & Belos::RecycleSubspace) existU_ = false; 00288 } 00290 00292 00293 00311 ReturnType solve(); 00312 00314 00317 00319 std::string description() const; 00320 00322 00323 private: 00324 00325 // Called by all constructors; Contains init instructions common to all constructors 00326 void init(); 00327 00328 // Computes harmonic eigenpairs of projected matrix created during one cycle. 00329 // Y contains the harmonic Ritz vectors corresponding to the recycleBlocks eigenvalues of smallest magnitude. 00330 void getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType> &F, 00331 const Teuchos::SerialDenseMatrix<int,ScalarType> &G, 00332 Teuchos::SerialDenseMatrix<int,ScalarType>& Y); 00333 00334 // Sort list of n floating-point numbers and return permutation vector 00335 void sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm); 00336 00337 // Initialize solver state storage 00338 void initializeStateStorage(); 00339 00340 // Linear problem. 00341 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_; 00342 00343 // Output manager. 00344 Teuchos::RCP<OutputManager<ScalarType> > printer_; 00345 Teuchos::RCP<std::ostream> outputStream_; 00346 00347 // Status test. 00348 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_; 00349 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_; 00350 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_; 00351 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_; 00352 00353 // Current parameter list. 00354 Teuchos::RCP<Teuchos::ParameterList> params_; 00355 00356 // Default solver values. 00357 static const MagnitudeType convtol_default_; 00358 static const int maxIters_default_; 00359 static const int blockSize_default_; 00360 static const int numBlocks_default_; 00361 static const int recycleBlocks_default_; 00362 static const bool showMaxResNormOnly_default_; 00363 static const int verbosity_default_; 00364 static const int outputStyle_default_; 00365 static const int outputFreq_default_; 00366 static const std::string label_default_; 00367 static const Teuchos::RCP<std::ostream> outputStream_default_; 00368 00369 // 00370 // Current solver values. 00371 // 00372 00374 MagnitudeType convtol_; 00375 00380 MagnitudeType achievedTol_; 00381 00383 int maxIters_; 00384 00386 int numIters_; 00387 00388 int numBlocks_, recycleBlocks_; 00389 bool showMaxResNormOnly_; 00390 int verbosity_, outputStyle_, outputFreq_; 00391 00393 // Solver State Storage 00395 // Search vectors 00396 Teuchos::RCP<MV> P_; 00397 // 00398 // A times current search direction 00399 Teuchos::RCP<MV> Ap_; 00400 // 00401 // Residual vector 00402 Teuchos::RCP<MV> r_; 00403 // 00404 // Preconditioned residual 00405 Teuchos::RCP<MV> z_; 00406 // 00407 // Flag indicating that the recycle space should be used 00408 bool existU_; 00409 // 00410 // Flag indicating that the updated recycle space has been created 00411 bool existU1_; 00412 // 00413 // Recycled subspace and its image 00414 Teuchos::RCP<MV> U_, AU_; 00415 // 00416 // Recycled subspace for next system and its image 00417 Teuchos::RCP<MV> U1_; 00418 // 00419 // Coefficients arising in RCG iteration 00420 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_; 00421 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Beta_; 00422 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_; 00423 // 00424 // Solutions to local least-squares problems 00425 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_; 00426 // 00427 // The matrix U^T A U 00428 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > UTAU_; 00429 // 00430 // LU factorization of U^T A U 00431 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_; 00432 // 00433 // Data from LU factorization of UTAU 00434 Teuchos::RCP<std::vector<int> > ipiv_; 00435 // 00436 // The matrix (AU)^T AU 00437 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAU_; 00438 // 00439 // The scalar r'*z 00440 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_; 00441 // 00442 // Matrices needed for calculation of harmonic Ritz eigenproblem 00443 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > F_,G_,Y_; 00444 // 00445 // Matrices needed for updating recycle space 00446 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > L2_,DeltaL2_,AU1TUDeltaL2_; 00447 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AU1TAU1_, AU1TU1_, AU1TAP_; 00448 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > FY_,GY_; 00449 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > APTAP_; 00450 Teuchos::RCP<MV> U1Y1_, PY2_; 00451 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAP_, AU1TU_; 00452 ScalarType dold; 00454 00455 // Timers. 00456 std::string label_; 00457 Teuchos::RCP<Teuchos::Time> timerSolve_; 00458 00459 // Internal state variables. 00460 bool params_Set_; 00461 }; 00462 00463 00464 // Default solver values. 00465 template<class ScalarType, class MV, class OP> 00466 const typename RCGSolMgr<ScalarType,MV,OP,false>::MagnitudeType 00467 RCGSolMgr<ScalarType,MV,OP,false>::convtol_default_ = 1e-8; 00468 00469 template<class ScalarType, class MV, class OP> 00470 const int RCGSolMgr<ScalarType,MV,OP,false>::maxIters_default_ = 1000; 00471 00472 template<class ScalarType, class MV, class OP> 00473 const int RCGSolMgr<ScalarType,MV,OP,false>::numBlocks_default_ = 25; 00474 00475 template<class ScalarType, class MV, class OP> 00476 const int RCGSolMgr<ScalarType,MV,OP,false>::blockSize_default_ = 1; 00477 00478 template<class ScalarType, class MV, class OP> 00479 const int RCGSolMgr<ScalarType,MV,OP,false>::recycleBlocks_default_ = 3; 00480 00481 template<class ScalarType, class MV, class OP> 00482 const bool RCGSolMgr<ScalarType,MV,OP,false>::showMaxResNormOnly_default_ = false; 00483 00484 template<class ScalarType, class MV, class OP> 00485 const int RCGSolMgr<ScalarType,MV,OP,false>::verbosity_default_ = Belos::Errors; 00486 00487 template<class ScalarType, class MV, class OP> 00488 const int RCGSolMgr<ScalarType,MV,OP,false>::outputStyle_default_ = Belos::General; 00489 00490 template<class ScalarType, class MV, class OP> 00491 const int RCGSolMgr<ScalarType,MV,OP,false>::outputFreq_default_ = -1; 00492 00493 template<class ScalarType, class MV, class OP> 00494 const std::string RCGSolMgr<ScalarType,MV,OP,false>::label_default_ = "Belos"; 00495 00496 template<class ScalarType, class MV, class OP> 00497 const Teuchos::RCP<std::ostream> RCGSolMgr<ScalarType,MV,OP,false>::outputStream_default_ = Teuchos::rcp(&std::cout,false); 00498 00499 // Empty Constructor 00500 template<class ScalarType, class MV, class OP> 00501 RCGSolMgr<ScalarType,MV,OP,false>::RCGSolMgr(): 00502 achievedTol_(0.0), 00503 numIters_(0) 00504 { 00505 init(); 00506 } 00507 00508 // Basic Constructor 00509 template<class ScalarType, class MV, class OP> 00510 RCGSolMgr<ScalarType,MV,OP,false>::RCGSolMgr( 00511 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 00512 const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 00513 problem_(problem), 00514 achievedTol_(0.0), 00515 numIters_(0) 00516 { 00517 init(); 00518 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager."); 00519 00520 // If the parameter list pointer is null, then set the current parameters to the default parameter list. 00521 if ( !is_null(pl) ) { 00522 setParameters( pl ); 00523 } 00524 } 00525 00526 // Common instructions executed in all constructors 00527 template<class ScalarType, class MV, class OP> 00528 void RCGSolMgr<ScalarType,MV,OP,false>::init() 00529 { 00530 outputStream_ = outputStream_default_; 00531 convtol_ = convtol_default_; 00532 maxIters_ = maxIters_default_; 00533 numBlocks_ = numBlocks_default_; 00534 recycleBlocks_ = recycleBlocks_default_; 00535 verbosity_ = verbosity_default_; 00536 outputStyle_= outputStyle_default_; 00537 outputFreq_= outputFreq_default_; 00538 showMaxResNormOnly_ = showMaxResNormOnly_default_; 00539 label_ = label_default_; 00540 params_Set_ = false; 00541 P_ = Teuchos::null; 00542 Ap_ = Teuchos::null; 00543 r_ = Teuchos::null; 00544 z_ = Teuchos::null; 00545 existU_ = false; 00546 existU1_ = false; 00547 U_ = Teuchos::null; 00548 AU_ = Teuchos::null; 00549 U1_ = Teuchos::null; 00550 Alpha_ = Teuchos::null; 00551 Beta_ = Teuchos::null; 00552 D_ = Teuchos::null; 00553 Delta_ = Teuchos::null; 00554 UTAU_ = Teuchos::null; 00555 LUUTAU_ = Teuchos::null; 00556 ipiv_ = Teuchos::null; 00557 AUTAU_ = Teuchos::null; 00558 rTz_old_ = Teuchos::null; 00559 F_ = Teuchos::null; 00560 G_ = Teuchos::null; 00561 Y_ = Teuchos::null; 00562 L2_ = Teuchos::null; 00563 DeltaL2_ = Teuchos::null; 00564 AU1TUDeltaL2_ = Teuchos::null; 00565 AU1TAU1_ = Teuchos::null; 00566 AU1TU1_ = Teuchos::null; 00567 AU1TAP_ = Teuchos::null; 00568 FY_ = Teuchos::null; 00569 GY_ = Teuchos::null; 00570 APTAP_ = Teuchos::null; 00571 U1Y1_ = Teuchos::null; 00572 PY2_ = Teuchos::null; 00573 AUTAP_ = Teuchos::null; 00574 AU1TU_ = Teuchos::null; 00575 dold = 0.; 00576 } 00577 00578 template<class ScalarType, class MV, class OP> 00579 void RCGSolMgr<ScalarType,MV,OP,false>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> ¶ms ) 00580 { 00581 // Create the internal parameter list if ones doesn't already exist. 00582 if (params_ == Teuchos::null) { 00583 params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) ); 00584 } 00585 else { 00586 params->validateParameters(*getValidParameters()); 00587 } 00588 00589 // Check for maximum number of iterations 00590 if (params->isParameter("Maximum Iterations")) { 00591 maxIters_ = params->get("Maximum Iterations",maxIters_default_); 00592 00593 // Update parameter in our list and in status test. 00594 params_->set("Maximum Iterations", maxIters_); 00595 if (maxIterTest_!=Teuchos::null) 00596 maxIterTest_->setMaxIters( maxIters_ ); 00597 } 00598 00599 // Check for the maximum number of blocks. 00600 if (params->isParameter("Num Blocks")) { 00601 numBlocks_ = params->get("Num Blocks",numBlocks_default_); 00602 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument, 00603 "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive."); 00604 00605 // Update parameter in our list. 00606 params_->set("Num Blocks", numBlocks_); 00607 } 00608 00609 // Check for the maximum number of blocks. 00610 if (params->isParameter("Num Recycled Blocks")) { 00611 recycleBlocks_ = params->get("Num Recycled Blocks",recycleBlocks_default_); 00612 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument, 00613 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive."); 00614 00615 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument, 00616 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\"."); 00617 00618 // Update parameter in our list. 00619 params_->set("Num Recycled Blocks", recycleBlocks_); 00620 } 00621 00622 // Check to see if the timer label changed. 00623 if (params->isParameter("Timer Label")) { 00624 std::string tempLabel = params->get("Timer Label", label_default_); 00625 00626 // Update parameter in our list and solver timer 00627 if (tempLabel != label_) { 00628 label_ = tempLabel; 00629 params_->set("Timer Label", label_); 00630 std::string solveLabel = label_ + ": RCGSolMgr total solve time"; 00631 #ifdef BELOS_TEUCHOS_TIME_MONITOR 00632 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel); 00633 #endif 00634 } 00635 } 00636 00637 // Check for a change in verbosity level 00638 if (params->isParameter("Verbosity")) { 00639 if (Teuchos::isParameterType<int>(*params,"Verbosity")) { 00640 verbosity_ = params->get("Verbosity", verbosity_default_); 00641 } else { 00642 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity"); 00643 } 00644 00645 // Update parameter in our list. 00646 params_->set("Verbosity", verbosity_); 00647 if (printer_ != Teuchos::null) 00648 printer_->setVerbosity(verbosity_); 00649 } 00650 00651 // Check for a change in output style 00652 if (params->isParameter("Output Style")) { 00653 if (Teuchos::isParameterType<int>(*params,"Output Style")) { 00654 outputStyle_ = params->get("Output Style", outputStyle_default_); 00655 } else { 00656 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style"); 00657 } 00658 00659 // Reconstruct the convergence test if the explicit residual test is not being used. 00660 params_->set("Output Style", outputStyle_); 00661 outputTest_ = Teuchos::null; 00662 } 00663 00664 // output stream 00665 if (params->isParameter("Output Stream")) { 00666 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream"); 00667 00668 // Update parameter in our list. 00669 params_->set("Output Stream", outputStream_); 00670 if (printer_ != Teuchos::null) 00671 printer_->setOStream( outputStream_ ); 00672 } 00673 00674 // frequency level 00675 if (verbosity_ & Belos::StatusTestDetails) { 00676 if (params->isParameter("Output Frequency")) { 00677 outputFreq_ = params->get("Output Frequency", outputFreq_default_); 00678 } 00679 00680 // Update parameter in out list and output status test. 00681 params_->set("Output Frequency", outputFreq_); 00682 if (outputTest_ != Teuchos::null) 00683 outputTest_->setOutputFrequency( outputFreq_ ); 00684 } 00685 00686 // Create output manager if we need to. 00687 if (printer_ == Teuchos::null) { 00688 printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) ); 00689 } 00690 00691 // Convergence 00692 typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t; 00693 typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t; 00694 00695 // Check for convergence tolerance 00696 if (params->isParameter("Convergence Tolerance")) { 00697 convtol_ = params->get("Convergence Tolerance",convtol_default_); 00698 00699 // Update parameter in our list and residual tests. 00700 params_->set("Convergence Tolerance", convtol_); 00701 if (convTest_ != Teuchos::null) 00702 convTest_->setTolerance( convtol_ ); 00703 } 00704 00705 if (params->isParameter("Show Maximum Residual Norm Only")) { 00706 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only"); 00707 00708 // Update parameter in our list and residual tests 00709 params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_); 00710 if (convTest_ != Teuchos::null) 00711 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ ); 00712 } 00713 00714 // Create status tests if we need to. 00715 00716 // Basic test checks maximum iterations and native residual. 00717 if (maxIterTest_ == Teuchos::null) 00718 maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) ); 00719 00720 // Implicit residual test, using the native residual to determine if convergence was achieved. 00721 if (convTest_ == Teuchos::null) 00722 convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) ); 00723 00724 if (sTest_ == Teuchos::null) 00725 sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) ); 00726 00727 if (outputTest_ == Teuchos::null) { 00728 00729 // Create the status test output class. 00730 // This class manages and formats the output from the status test. 00731 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ ); 00732 outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined ); 00733 00734 // Set the solver string for the output test 00735 std::string solverDesc = " Recycling CG "; 00736 outputTest_->setSolverDesc( solverDesc ); 00737 } 00738 00739 // Create the timer if we need to. 00740 if (timerSolve_ == Teuchos::null) { 00741 std::string solveLabel = label_ + ": RCGSolMgr total solve time"; 00742 #ifdef BELOS_TEUCHOS_TIME_MONITOR 00743 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel); 00744 #endif 00745 } 00746 00747 // Inform the solver manager that the current parameters were set. 00748 params_Set_ = true; 00749 } 00750 00751 00752 template<class ScalarType, class MV, class OP> 00753 Teuchos::RCP<const Teuchos::ParameterList> 00754 RCGSolMgr<ScalarType,MV,OP,false>::getValidParameters() const 00755 { 00756 static Teuchos::RCP<const Teuchos::ParameterList> validPL; 00757 00758 // Set all the valid parameters and their default values. 00759 if(is_null(validPL)) { 00760 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList(); 00761 pl->set("Convergence Tolerance", convtol_default_, 00762 "The relative residual tolerance that needs to be achieved by the\n" 00763 "iterative solver in order for the linear system to be declared converged."); 00764 pl->set("Maximum Iterations", maxIters_default_, 00765 "The maximum number of block iterations allowed for each\n" 00766 "set of RHS solved."); 00767 pl->set("Block Size", blockSize_default_, 00768 "Block Size Parameter -- currently must be 1 for RCG"); 00769 pl->set("Num Blocks", numBlocks_default_, 00770 "The length of a cycle (and this max number of search vectors kept)\n"); 00771 pl->set("Num Recycled Blocks", recycleBlocks_default_, 00772 "The number of vectors in the recycle subspace."); 00773 pl->set("Verbosity", verbosity_default_, 00774 "What type(s) of solver information should be outputted\n" 00775 "to the output stream."); 00776 pl->set("Output Style", outputStyle_default_, 00777 "What style is used for the solver information outputted\n" 00778 "to the output stream."); 00779 pl->set("Output Frequency", outputFreq_default_, 00780 "How often convergence information should be outputted\n" 00781 "to the output stream."); 00782 pl->set("Output Stream", outputStream_default_, 00783 "A reference-counted pointer to the output stream where all\n" 00784 "solver output is sent."); 00785 pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_, 00786 "When convergence information is printed, only show the maximum\n" 00787 "relative residual norm when the block size is greater than one."); 00788 pl->set("Timer Label", label_default_, 00789 "The string to use as a prefix for the timer labels."); 00790 // pl->set("Restart Timers", restartTimers_); 00791 validPL = pl; 00792 } 00793 return validPL; 00794 } 00795 00796 // initializeStateStorage 00797 template<class ScalarType, class MV, class OP> 00798 void RCGSolMgr<ScalarType,MV,OP,false>::initializeStateStorage() { 00799 00800 // Check if there is any multivector to clone from. 00801 Teuchos::RCP<const MV> rhsMV = problem_->getRHS(); 00802 if (rhsMV == Teuchos::null) { 00803 // Nothing to do 00804 return; 00805 } 00806 else { 00807 00808 // Initialize the state storage 00809 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVText::GetGlobalLength(*rhsMV),std::invalid_argument, 00810 "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!"); 00811 00812 // If the subspace has not been initialized before, generate it using the RHS from lp_. 00813 if (P_ == Teuchos::null) { 00814 P_ = MVT::Clone( *rhsMV, numBlocks_+2 ); 00815 } 00816 else { 00817 // Generate P_ by cloning itself ONLY if more space is needed. 00818 if (MVT::GetNumberVecs(*P_) < numBlocks_+2) { 00819 Teuchos::RCP<const MV> tmp = P_; 00820 P_ = MVT::Clone( *tmp, numBlocks_+2 ); 00821 } 00822 } 00823 00824 // Generate Ap_ only if it doesn't exist 00825 if (Ap_ == Teuchos::null) 00826 Ap_ = MVT::Clone( *rhsMV, 1 ); 00827 00828 // Generate r_ only if it doesn't exist 00829 if (r_ == Teuchos::null) 00830 r_ = MVT::Clone( *rhsMV, 1 ); 00831 00832 // Generate z_ only if it doesn't exist 00833 if (z_ == Teuchos::null) 00834 z_ = MVT::Clone( *rhsMV, 1 ); 00835 00836 // If the recycle space has not been initialized before, generate it using the RHS from lp_. 00837 if (U_ == Teuchos::null) { 00838 U_ = MVT::Clone( *rhsMV, recycleBlocks_ ); 00839 } 00840 else { 00841 // Generate U_ by cloning itself ONLY if more space is needed. 00842 if (MVT::GetNumberVecs(*U_) < recycleBlocks_) { 00843 Teuchos::RCP<const MV> tmp = U_; 00844 U_ = MVT::Clone( *tmp, recycleBlocks_ ); 00845 } 00846 } 00847 00848 // If the recycle space has not be initialized before, generate it using the RHS from lp_. 00849 if (AU_ == Teuchos::null) { 00850 AU_ = MVT::Clone( *rhsMV, recycleBlocks_ ); 00851 } 00852 else { 00853 // Generate AU_ by cloning itself ONLY if more space is needed. 00854 if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) { 00855 Teuchos::RCP<const MV> tmp = AU_; 00856 AU_ = MVT::Clone( *tmp, recycleBlocks_ ); 00857 } 00858 } 00859 00860 // If the recycle space has not been initialized before, generate it using the RHS from lp_. 00861 if (U1_ == Teuchos::null) { 00862 U1_ = MVT::Clone( *rhsMV, recycleBlocks_ ); 00863 } 00864 else { 00865 // Generate U1_ by cloning itself ONLY if more space is needed. 00866 if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) { 00867 Teuchos::RCP<const MV> tmp = U1_; 00868 U1_ = MVT::Clone( *tmp, recycleBlocks_ ); 00869 } 00870 } 00871 00872 // Generate Alpha_ only if it doesn't exist, otherwise resize it. 00873 if (Alpha_ == Teuchos::null) 00874 Alpha_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, 1 ) ); 00875 else { 00876 if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) ) 00877 Alpha_->reshape( numBlocks_, 1 ); 00878 } 00879 00880 // Generate Beta_ only if it doesn't exist, otherwise resize it. 00881 if (Beta_ == Teuchos::null) 00882 Beta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + 1, 1 ) ); 00883 else { 00884 if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) ) 00885 Beta_->reshape( numBlocks_ + 1, 1 ); 00886 } 00887 00888 // Generate D_ only if it doesn't exist, otherwise resize it. 00889 if (D_ == Teuchos::null) 00890 D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ , 1 ) ); 00891 else { 00892 if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) ) 00893 D_->reshape( numBlocks_, 1 ); 00894 } 00895 00896 // Generate Delta_ only if it doesn't exist, otherwise resize it. 00897 if (Delta_ == Teuchos::null) 00898 Delta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ + 1 ) ); 00899 else { 00900 if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) ) 00901 Delta_->reshape( recycleBlocks_, numBlocks_ + 1 ); 00902 } 00903 00904 // Generate UTAU_ only if it doesn't exist, otherwise resize it. 00905 if (UTAU_ == Teuchos::null) 00906 UTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 00907 else { 00908 if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) ) 00909 UTAU_->reshape( recycleBlocks_, recycleBlocks_ ); 00910 } 00911 00912 // Generate LUUTAU_ only if it doesn't exist, otherwise resize it. 00913 if (LUUTAU_ == Teuchos::null) 00914 LUUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 00915 else { 00916 if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) ) 00917 LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ ); 00918 } 00919 00920 // Generate ipiv_ only if it doesn't exist, otherwise resize it. 00921 if (ipiv_ == Teuchos::null) 00922 ipiv_ = Teuchos::rcp( new std::vector<int>(recycleBlocks_) ); 00923 else { 00924 if ( (int)ipiv_->size() != recycleBlocks_ ) // if ipiv not correct size, always resize it 00925 ipiv_->resize(recycleBlocks_); 00926 } 00927 00928 // Generate AUTAU_ only if it doesn't exist, otherwise resize it. 00929 if (AUTAU_ == Teuchos::null) 00930 AUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 00931 else { 00932 if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) ) 00933 AUTAU_->reshape( recycleBlocks_, recycleBlocks_ ); 00934 } 00935 00936 // Generate rTz_old_ only if it doesn't exist 00937 if (rTz_old_ == Teuchos::null) 00938 rTz_old_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 ) ); 00939 else { 00940 if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) ) 00941 rTz_old_->reshape( 1, 1 ); 00942 } 00943 00944 // Generate F_ only if it doesn't exist 00945 if (F_ == Teuchos::null) 00946 F_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) ); 00947 else { 00948 if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) ) 00949 F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ); 00950 } 00951 00952 // Generate G_ only if it doesn't exist 00953 if (G_ == Teuchos::null) 00954 G_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) ); 00955 else { 00956 if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) ) 00957 G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ); 00958 } 00959 00960 // Generate Y_ only if it doesn't exist 00961 if (Y_ == Teuchos::null) 00962 Y_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, recycleBlocks_ ) ); 00963 else { 00964 if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) ) 00965 Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ); 00966 } 00967 00968 // Generate L2_ only if it doesn't exist 00969 if (L2_ == Teuchos::null) 00970 L2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+1, numBlocks_ ) ); 00971 else { 00972 if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) ) 00973 L2_->reshape( numBlocks_+1, numBlocks_ ); 00974 } 00975 00976 // Generate DeltaL2_ only if it doesn't exist 00977 if (DeltaL2_ == Teuchos::null) 00978 DeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) ); 00979 else { 00980 if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) ) 00981 DeltaL2_->reshape( recycleBlocks_, numBlocks_ ); 00982 } 00983 00984 // Generate AU1TUDeltaL2_ only if it doesn't exist 00985 if (AU1TUDeltaL2_ == Teuchos::null) 00986 AU1TUDeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) ); 00987 else { 00988 if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) ) 00989 AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ ); 00990 } 00991 00992 // Generate AU1TAU1_ only if it doesn't exist 00993 if (AU1TAU1_ == Teuchos::null) 00994 AU1TAU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 00995 else { 00996 if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) ) 00997 AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ ); 00998 } 00999 01000 // Generate GY_ only if it doesn't exist 01001 if (GY_ == Teuchos::null) 01002 GY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) ); 01003 else { 01004 if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) ) 01005 GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ ); 01006 } 01007 01008 // Generate AU1TU1_ only if it doesn't exist 01009 if (AU1TU1_ == Teuchos::null) 01010 AU1TU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 01011 else { 01012 if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) ) 01013 AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ ); 01014 } 01015 01016 // Generate FY_ only if it doesn't exist 01017 if (FY_ == Teuchos::null) 01018 FY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) ); 01019 else { 01020 if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) ) 01021 FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ ); 01022 } 01023 01024 // Generate AU1TAP_ only if it doesn't exist 01025 if (AU1TAP_ == Teuchos::null) 01026 AU1TAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) ); 01027 else { 01028 if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) ) 01029 AU1TAP_->reshape( recycleBlocks_, numBlocks_ ); 01030 } 01031 01032 // Generate APTAP_ only if it doesn't exist 01033 if (APTAP_ == Teuchos::null) 01034 APTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, numBlocks_ ) ); 01035 else { 01036 if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) ) 01037 APTAP_->reshape( numBlocks_, numBlocks_ ); 01038 } 01039 01040 // If the subspace has not been initialized before, generate it using the RHS from lp_. 01041 if (U1Y1_ == Teuchos::null) { 01042 U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ ); 01043 } 01044 else { 01045 // Generate U1Y1_ by cloning itself ONLY if more space is needed. 01046 if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) { 01047 Teuchos::RCP<const MV> tmp = U1Y1_; 01048 U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ ); 01049 } 01050 } 01051 01052 // If the subspace has not been initialized before, generate it using the RHS from lp_. 01053 if (PY2_ == Teuchos::null) { 01054 PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ ); 01055 } 01056 else { 01057 // Generate PY2_ by cloning itself ONLY if more space is needed. 01058 if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) { 01059 Teuchos::RCP<const MV> tmp = PY2_; 01060 PY2_ = MVT::Clone( *tmp, recycleBlocks_ ); 01061 } 01062 } 01063 01064 // Generate AUTAP_ only if it doesn't exist 01065 if (AUTAP_ == Teuchos::null) 01066 AUTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) ); 01067 else { 01068 if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) ) 01069 AUTAP_->reshape( recycleBlocks_, numBlocks_ ); 01070 } 01071 01072 // Generate AU1TU_ only if it doesn't exist 01073 if (AU1TU_ == Teuchos::null) 01074 AU1TU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 01075 else { 01076 if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) ) 01077 AU1TU_->reshape( recycleBlocks_, recycleBlocks_ ); 01078 } 01079 01080 01081 } 01082 } 01083 01084 template<class ScalarType, class MV, class OP> 01085 ReturnType RCGSolMgr<ScalarType,MV,OP,false>::solve() { 01086 01087 Teuchos::LAPACK<int,ScalarType> lapack; 01088 std::vector<int> index(recycleBlocks_); 01089 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one(); 01090 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero(); 01091 01092 // Count of number of cycles performed on current rhs 01093 int cycle = 0; 01094 01095 // Set the current parameters if they were not set before. 01096 // NOTE: This may occur if the user generated the solver manager with the default constructor and 01097 // then didn't set any parameters using setParameters(). 01098 if (!params_Set_) { 01099 setParameters(Teuchos::parameterList(*getValidParameters())); 01100 } 01101 01102 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,RCGSolMgrLinearProblemFailure, 01103 "Belos::RCGSolMgr::solve(): Linear problem is not a valid object."); 01104 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),RCGSolMgrLinearProblemFailure, 01105 "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called."); 01106 TEUCHOS_TEST_FOR_EXCEPTION((problem_->getLeftPrec() != Teuchos::null)&&(problem_->getRightPrec() != Teuchos::null), 01107 RCGSolMgrLinearProblemFailure, 01108 "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner."); 01109 01110 // Grab the preconditioning object 01111 Teuchos::RCP<OP> precObj; 01112 if (problem_->getLeftPrec() != Teuchos::null) { 01113 precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec()); 01114 } 01115 else if (problem_->getRightPrec() != Teuchos::null) { 01116 precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec()); 01117 } 01118 01119 // Create indices for the linear systems to be solved. 01120 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) ); 01121 std::vector<int> currIdx(1); 01122 currIdx[0] = 0; 01123 01124 // Inform the linear problem of the current linear system to solve. 01125 problem_->setLSIndex( currIdx ); 01126 01127 // Check the number of blocks and change them if necessary. 01128 ptrdiff_t dim = MVText::GetGlobalLength( *(problem_->getRHS()) ); 01129 if (numBlocks_ > dim) { 01130 numBlocks_ = Teuchos::asSafe<int>(dim); 01131 params_->set("Num Blocks", numBlocks_); 01132 printer_->stream(Warnings) << 01133 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl << 01134 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl; 01135 } 01136 01137 // Initialize storage for all state variables 01138 initializeStateStorage(); 01139 01140 // Parameter list 01141 Teuchos::ParameterList plist; 01142 plist.set("Num Blocks",numBlocks_); 01143 plist.set("Recycled Blocks",recycleBlocks_); 01144 01145 // Reset the status test. 01146 outputTest_->reset(); 01147 01148 // Assume convergence is achieved, then let any failed convergence set this to false. 01149 bool isConverged = true; 01150 01151 // Compute AU = A*U, UTAU = U'*AU, AUTAU = (AU)'*(AU) 01152 if (existU_) { 01153 index.resize(recycleBlocks_); 01154 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01155 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index ); 01156 index.resize(recycleBlocks_); 01157 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01158 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index ); 01159 // Initialize AU 01160 problem_->applyOp( *Utmp, *AUtmp ); 01161 // Initialize UTAU 01162 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ ); 01163 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp ); 01164 // Initialize AUTAU ( AUTAU = AU'*(M\AU) ) 01165 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ ); 01166 if ( precObj != Teuchos::null ) { 01167 index.resize(recycleBlocks_); 01168 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01169 index.resize(recycleBlocks_); 01170 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01171 Teuchos::RCP<MV> PCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage 01172 OPT::Apply( *precObj, *AUtmp, *PCAU ); 01173 MVT::MvTransMv( one, *AUtmp, *PCAU, AUTAUtmp ); 01174 } else { 01175 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp ); 01176 } 01177 } 01178 01180 // RCG solver 01181 01182 Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter; 01183 rcg_iter = Teuchos::rcp( new RCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) ); 01184 01185 // Enter solve() iterations 01186 { 01187 #ifdef BELOS_TEUCHOS_TIME_MONITOR 01188 Teuchos::TimeMonitor slvtimer(*timerSolve_); 01189 #endif 01190 01191 while ( numRHS2Solve > 0 ) { 01192 01193 // Debugging output to tell use if recycle space exists and will be used 01194 if (printer_->isVerbosity( Debug ) ) { 01195 if (existU_) printer_->print( Debug, "Using recycle space generated from previous call to solve()." ); 01196 else printer_->print( Debug, "No recycle space exists." ); 01197 } 01198 01199 // Reset the number of iterations. 01200 rcg_iter->resetNumIters(); 01201 01202 // Set the current number of recycle blocks and subspace dimension with the RCG iteration. 01203 rcg_iter->setSize( recycleBlocks_, numBlocks_ ); 01204 01205 // Reset the number of calls that the status test output knows about. 01206 outputTest_->resetNumCalls(); 01207 01208 // indicate that updated recycle space has not yet been generated for this linear system 01209 existU1_ = false; 01210 01211 // reset cycle count 01212 cycle = 0; 01213 01214 // Get the current residual 01215 problem_->computeCurrResVec( &*r_ ); 01216 01217 // If U exists, find best soln over this space first 01218 if (existU_) { 01219 // Solve linear system UTAU * y = (U'*r) 01220 Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1); 01221 index.resize(recycleBlocks_); 01222 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01223 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index ); 01224 MVT::MvTransMv( one, *Utmp, *r_, Utr ); 01225 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ ); 01226 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ); 01227 LUUTAUtmp.assign(UTAUtmp); 01228 int info = 0; 01229 lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info); 01230 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure, 01231 "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution."); 01232 01233 // Update solution (x = x + U*y) 01234 MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() ); 01235 01236 // Update residual ( r = r - AU*y ) 01237 index.resize(recycleBlocks_); 01238 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01239 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index ); 01240 MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ ); 01241 } 01242 01243 if ( precObj != Teuchos::null ) { 01244 OPT::Apply( *precObj, *r_, *z_ ); 01245 } else { 01246 z_ = r_; 01247 } 01248 01249 // rTz_old = r'*z 01250 MVT::MvTransMv( one, *r_, *z_, *rTz_old_ ); 01251 01252 if ( existU_ ) { 01253 // mu = UTAU\(AU'*z); 01254 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1); 01255 index.resize(recycleBlocks_); 01256 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01257 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index ); 01258 MVT::MvTransMv( one, *AUtmp, *z_, mu ); 01259 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ); 01260 char TRANS = 'N'; 01261 int info; 01262 lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info ); 01263 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure, 01264 "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution."); 01265 // p = z - U*mu; 01266 index.resize( 1 ); 01267 index[0] = 0; 01268 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index ); 01269 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp); 01270 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp ); 01271 } else { 01272 // p = z; 01273 index.resize( 1 ); 01274 index[0] = 0; 01275 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index ); 01276 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp); 01277 } 01278 01279 // Set the new state and initialize the solver. 01280 RCGIterState<ScalarType,MV> newstate; 01281 01282 // Create RCP views here 01283 index.resize( numBlocks_+1 ); 01284 for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; } 01285 newstate.P = MVT::CloneViewNonConst( *P_, index ); 01286 index.resize( recycleBlocks_ ); 01287 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01288 newstate.U = MVT::CloneViewNonConst( *U_, index ); 01289 index.resize( recycleBlocks_ ); 01290 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01291 newstate.AU = MVT::CloneViewNonConst( *AU_, index ); 01292 newstate.Alpha = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) ); 01293 newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) ); 01294 newstate.D = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) ); 01295 newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) ); 01296 newstate.LUUTAU = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) ); 01297 // assign the rest of the values in the struct 01298 newstate.curDim = 1; // We have initialized the first search vector 01299 newstate.Ap = Ap_; 01300 newstate.r = r_; 01301 newstate.z = z_; 01302 newstate.existU = existU_; 01303 newstate.ipiv = ipiv_; 01304 newstate.rTz_old = rTz_old_; 01305 01306 rcg_iter->initialize(newstate); 01307 01308 while(1) { 01309 01310 // tell rcg_iter to iterate 01311 try { 01312 rcg_iter->iterate(); 01313 01315 // 01316 // check convergence first 01317 // 01319 if ( convTest_->getStatus() == Passed ) { 01320 // We have convergence 01321 break; // break from while(1){rcg_iter->iterate()} 01322 } 01324 // 01325 // check for maximum iterations 01326 // 01328 else if ( maxIterTest_->getStatus() == Passed ) { 01329 // we don't have convergence 01330 isConverged = false; 01331 break; // break from while(1){rcg_iter->iterate()} 01332 } 01334 // 01335 // check if cycle complete; update for next cycle 01336 // 01338 else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) { 01339 // index into P_ of last search vector generated this cycle 01340 int lastp = -1; 01341 // index into Beta_ of last entry generated this cycle 01342 int lastBeta = -1; 01343 if (recycleBlocks_ > 0) { 01344 if (!existU_) { 01345 if (cycle == 0) { // No U, no U1 01346 01347 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ ); 01348 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ ); 01349 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 ); 01350 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 ); 01351 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 ); 01352 Ftmp.putScalar(zero); 01353 Gtmp.putScalar(zero); 01354 for (int ii=0;ii<numBlocks_;ii++) { 01355 Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0)); 01356 if (ii > 0) { 01357 Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01358 Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01359 } 01360 Ftmp(ii,ii) = Dtmp(ii,0); 01361 } 01362 01363 // compute harmonic Ritz vectors 01364 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ ); 01365 getHarmonicVecs(Ftmp,Gtmp,Ytmp); 01366 01367 // U1 = [P(:,1:end-1)*Y]; 01368 index.resize( numBlocks_ ); 01369 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; } 01370 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index ); 01371 index.resize( recycleBlocks_ ); 01372 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01373 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index ); 01374 MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp ); 01375 01376 // Precompute some variables for next cycle 01377 01378 // AU1TAU1 = Y'*G*Y; 01379 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ ); 01380 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero); 01381 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ ); 01382 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero); 01383 01384 // AU1TU1 = Y'*F*Y; 01385 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ ); 01386 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero); 01387 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ ); 01388 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero); 01389 01390 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 ); 01391 // Must reinitialize AU1TAP; can become dense later 01392 AU1TAPtmp.putScalar(zero); 01393 // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end)); 01394 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0); 01395 for (int ii=0; ii<recycleBlocks_; ++ii) { 01396 AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp; 01397 } 01398 01399 // indicate that updated recycle space now defined 01400 existU1_ = true; 01401 01402 // Indicate the size of the P, Beta structures generated this cycle 01403 lastp = numBlocks_; 01404 lastBeta = numBlocks_-1; 01405 01406 } // if (cycle == 0) 01407 else { // No U, but U1 guaranteed to exist now 01408 01409 // Finish computation of subblocks 01410 // AU1TAP = AU1TAP * D(1); 01411 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 ); 01412 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ ); 01413 AU1TAPtmp.scale(Dtmp(0,0)); 01414 01415 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 ); 01416 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 ); 01417 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ ); 01418 APTAPtmp.putScalar(zero); 01419 for (int ii=0; ii<numBlocks_; ii++) { 01420 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0)); 01421 if (ii > 0) { 01422 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01423 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01424 } 01425 } 01426 01427 // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)]; 01428 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01429 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ ); 01430 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01431 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01432 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01433 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ ); 01434 F11.assign(AU1TU1tmp); 01435 F12.putScalar(zero); 01436 F21.putScalar(zero); 01437 F22.putScalar(zero); 01438 for(int ii=0;ii<numBlocks_;ii++) { 01439 F22(ii,ii) = Dtmp(ii,0); 01440 } 01441 01442 // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP]; 01443 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01444 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ ); 01445 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ ); 01446 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01447 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01448 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01449 G11.assign(AU1TAU1tmp); 01450 G12.assign(AU1TAPtmp); 01451 // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually) 01452 for (int ii=0;ii<recycleBlocks_;++ii) 01453 for (int jj=0;jj<numBlocks_;++jj) 01454 G21(jj,ii) = G12(ii,jj); 01455 G22.assign(APTAPtmp); 01456 01457 // compute harmonic Ritz vectors 01458 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ ); 01459 getHarmonicVecs(Ftmp,Gtmp,Ytmp); 01460 01461 // U1 = [U1 P(:,2:end-1)]*Y; 01462 index.resize( numBlocks_ ); 01463 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; } 01464 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index ); 01465 index.resize( recycleBlocks_ ); 01466 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01467 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index ); 01468 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01469 index.resize( recycleBlocks_ ); 01470 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01471 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index ); 01472 index.resize( recycleBlocks_ ); 01473 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01474 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index ); 01475 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ ); 01476 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp ); 01477 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp ); 01478 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp); 01479 01480 // Precompute some variables for next cycle 01481 01482 // AU1TAU1 = Y'*G*Y; 01483 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01484 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero); 01485 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero); 01486 01487 // AU1TAP = zeros(k,m); 01488 // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end)); 01489 AU1TAPtmp.putScalar(zero); 01490 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0); 01491 for (int ii=0; ii<recycleBlocks_; ++ii) { 01492 AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp; 01493 } 01494 01495 // AU1TU1 = Y'*F*Y; 01496 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01497 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero); 01498 //Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ ); 01499 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero); 01500 01501 // Indicate the size of the P, Beta structures generated this cycle 01502 lastp = numBlocks_+1; 01503 lastBeta = numBlocks_; 01504 01505 } // if (cycle != 1) 01506 } // if (!existU_) 01507 else { // U exists 01508 if (cycle == 0) { // No U1, but U exists 01509 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 ); 01510 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 ); 01511 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 ); 01512 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ ); 01513 APTAPtmp.putScalar(zero); 01514 for (int ii=0; ii<numBlocks_; ii++) { 01515 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0)); 01516 if (ii > 0) { 01517 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01518 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01519 } 01520 } 01521 01522 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ ); 01523 L2tmp.putScalar(zero); 01524 for(int ii=0;ii<numBlocks_;ii++) { 01525 L2tmp(ii,ii) = 1./Alphatmp(ii,0); 01526 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0); 01527 } 01528 01529 // AUTAP = UTAU*Delta*L2; 01530 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ ); 01531 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ ); 01532 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 ); 01533 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ ); 01534 DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero); 01535 AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero); 01536 01537 // F = [UTAU zeros(k,m); zeros(m,k) diag(D)]; 01538 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01539 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ ); 01540 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01541 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01542 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01543 F11.assign(UTAUtmp); 01544 F12.putScalar(zero); 01545 F21.putScalar(zero); 01546 F22.putScalar(zero); 01547 for(int ii=0;ii<numBlocks_;ii++) { 01548 F22(ii,ii) = Dtmp(ii,0); 01549 } 01550 01551 // G = [AUTAU AUTAP; AUTAP' APTAP]; 01552 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01553 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ ); 01554 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01555 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01556 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01557 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ ); 01558 G11.assign(AUTAUtmp); 01559 G12.assign(AUTAPtmp); 01560 // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually) 01561 for (int ii=0;ii<recycleBlocks_;++ii) 01562 for (int jj=0;jj<numBlocks_;++jj) 01563 G21(jj,ii) = G12(ii,jj); 01564 G22.assign(APTAPtmp); 01565 01566 // compute harmonic Ritz vectors 01567 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ ); 01568 getHarmonicVecs(Ftmp,Gtmp,Ytmp); 01569 01570 // U1 = [U P(:,1:end-1)]*Y; 01571 index.resize( recycleBlocks_ ); 01572 for (int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; } 01573 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index ); 01574 index.resize( numBlocks_ ); 01575 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; } 01576 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index ); 01577 index.resize( recycleBlocks_ ); 01578 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01579 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index ); 01580 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01581 index.resize( recycleBlocks_ ); 01582 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01583 Teuchos::RCP<MV> UY1tmp = MVT::CloneViewNonConst( *U1Y1_, index ); 01584 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ ); 01585 index.resize( recycleBlocks_ ); 01586 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01587 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index ); 01588 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp ); 01589 MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp ); 01590 MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp); 01591 01592 // Precompute some variables for next cycle 01593 01594 // AU1TAU1 = Y'*G*Y; 01595 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01596 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero); 01597 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ ); 01598 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero); 01599 01600 // AU1TU1 = Y'*F*Y; 01601 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01602 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero); 01603 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ ); 01604 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero); 01605 01606 // AU1TU = UTAU; 01607 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ ); 01608 AU1TUtmp.assign(UTAUtmp); 01609 01610 // dold = D(end); 01611 dold = Dtmp(numBlocks_-1,0); 01612 01613 // indicate that updated recycle space now defined 01614 existU1_ = true; 01615 01616 // Indicate the size of the P, Beta structures generated this cycle 01617 lastp = numBlocks_; 01618 lastBeta = numBlocks_-1; 01619 } 01620 else { // Have U and U1 01621 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 ); 01622 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 ); 01623 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 ); 01624 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ ); 01625 for (int ii=0; ii<numBlocks_; ii++) { 01626 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0)); 01627 if (ii > 0) { 01628 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01629 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01630 } 01631 } 01632 01633 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ ); 01634 for(int ii=0;ii<numBlocks_;ii++) { 01635 L2tmp(ii,ii) = 1./Alphatmp(ii,0); 01636 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0); 01637 } 01638 01639 // M(end,1) = dold*(-Beta(1)/Alpha(1)); 01640 // AU1TAP = Y'*[AU1TU*Delta*L2; M]; 01641 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ ); 01642 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 ); 01643 DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero); 01644 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ ); 01645 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ ); 01646 AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero); 01647 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ ); 01648 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ ); 01649 AU1TAPtmp.putScalar(zero); 01650 AU1TAPtmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero); 01651 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01652 ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0)); 01653 for(int ii=0;ii<recycleBlocks_;ii++) { 01654 AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val; 01655 } 01656 01657 // AU1TU = Y1'*AU1TU 01658 Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ ); 01659 Y1TAU1TU.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero); 01660 AU1TUtmp.assign(Y1TAU1TU); 01661 01662 // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)]; 01663 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01664 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ ); 01665 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01666 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01667 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01668 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ ); 01669 F11.assign(AU1TU1tmp); 01670 for(int ii=0;ii<numBlocks_;ii++) { 01671 F22(ii,ii) = Dtmp(ii,0); 01672 } 01673 01674 // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP]; 01675 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01676 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ ); 01677 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01678 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01679 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01680 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ ); 01681 G11.assign(AU1TAU1tmp); 01682 G12.assign(AU1TAPtmp); 01683 // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually) 01684 for (int ii=0;ii<recycleBlocks_;++ii) 01685 for (int jj=0;jj<numBlocks_;++jj) 01686 G21(jj,ii) = G12(ii,jj); 01687 G22.assign(APTAPtmp); 01688 01689 // compute harmonic Ritz vectors 01690 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ ); 01691 getHarmonicVecs(Ftmp,Gtmp,Ytmp); 01692 01693 // U1 = [U1 P(:,2:end-1)]*Y; 01694 index.resize( numBlocks_ ); 01695 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; } 01696 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index ); 01697 index.resize( recycleBlocks_ ); 01698 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01699 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index ); 01700 index.resize( recycleBlocks_ ); 01701 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01702 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index ); 01703 index.resize( recycleBlocks_ ); 01704 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01705 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index ); 01706 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp ); 01707 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp ); 01708 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp); 01709 01710 // Precompute some variables for next cycle 01711 01712 // AU1TAU1 = Y'*G*Y; 01713 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01714 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero); 01715 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero); 01716 01717 // AU1TU1 = Y'*F*Y; 01718 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01719 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero); 01720 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero); 01721 01722 // dold = D(end); 01723 dold = Dtmp(numBlocks_-1,0); 01724 01725 // Indicate the size of the P, Beta structures generated this cycle 01726 lastp = numBlocks_+1; 01727 lastBeta = numBlocks_; 01728 01729 } 01730 } 01731 } // if (recycleBlocks_ > 0) 01732 01733 // Cleanup after end of cycle 01734 01735 // P = P(:,end-1:end); 01736 index.resize( 2 ); 01737 index[0] = lastp-1; index[1] = lastp; 01738 Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_, index ); 01739 index[0] = 0; index[1] = 1; 01740 MVT::SetBlock(*Ptmp2,index,*P_); 01741 01742 // Beta = Beta(end); 01743 (*Beta_)(0,0) = (*Beta_)(lastBeta,0); 01744 01745 // Delta = Delta(:,end); 01746 if (existU_) { // Delta only initialized if U exists 01747 Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 ); 01748 Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ ); 01749 mu1.assign(mu2); 01750 } 01751 01752 // Now reinitialize state variables for next cycle 01753 newstate.P = Teuchos::null; 01754 index.resize( numBlocks_+1 ); 01755 for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; } 01756 newstate.P = MVT::CloneViewNonConst( *P_, index ); 01757 01758 newstate.Beta = Teuchos::null; 01759 newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) ); 01760 01761 newstate.Delta = Teuchos::null; 01762 newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) ); 01763 01764 newstate.curDim = 1; // We have initialized the first search vector 01765 01766 // Pass to iteration object 01767 rcg_iter->initialize(newstate); 01768 01769 // increment cycle count 01770 cycle = cycle + 1; 01771 01772 } 01774 // 01775 // we returned from iterate(), but none of our status tests Passed. 01776 // something is wrong, and it is probably our fault. 01777 // 01779 else { 01780 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error, 01781 "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate()."); 01782 } 01783 } 01784 catch (const std::exception &e) { 01785 printer_->stream(Errors) << "Error! Caught std::exception in RCGIter::iterate() at iteration " 01786 << rcg_iter->getNumIters() << std::endl 01787 << e.what() << std::endl; 01788 throw; 01789 } 01790 } 01791 01792 // Inform the linear problem that we are finished with this block linear system. 01793 problem_->setCurrLS(); 01794 01795 // Update indices for the linear systems to be solved. 01796 numRHS2Solve -= 1; 01797 if ( numRHS2Solve > 0 ) { 01798 currIdx[0]++; 01799 // Set the next indices. 01800 problem_->setLSIndex( currIdx ); 01801 } 01802 else { 01803 currIdx.resize( numRHS2Solve ); 01804 } 01805 01806 // Update the recycle space for the next linear system 01807 if (existU1_) { // be sure updated recycle space was created 01808 // U = U1 01809 index.resize(recycleBlocks_); 01810 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01811 MVT::SetBlock(*U1_,index,*U_); 01812 // Set flag indicating recycle space is now defined 01813 existU_ = true; 01814 if (numRHS2Solve > 0) { // also update AU, UTAU, and AUTAU 01815 // Free pointers in newstate 01816 newstate.P = Teuchos::null; 01817 newstate.Ap = Teuchos::null; 01818 newstate.r = Teuchos::null; 01819 newstate.z = Teuchos::null; 01820 newstate.U = Teuchos::null; 01821 newstate.AU = Teuchos::null; 01822 newstate.Alpha = Teuchos::null; 01823 newstate.Beta = Teuchos::null; 01824 newstate.D = Teuchos::null; 01825 newstate.Delta = Teuchos::null; 01826 newstate.LUUTAU = Teuchos::null; 01827 newstate.ipiv = Teuchos::null; 01828 newstate.rTz_old = Teuchos::null; 01829 01830 // Reinitialize AU, UTAU, AUTAU 01831 index.resize(recycleBlocks_); 01832 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01833 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index ); 01834 index.resize(recycleBlocks_); 01835 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01836 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index ); 01837 // Initialize AU 01838 problem_->applyOp( *Utmp, *AUtmp ); 01839 // Initialize UTAU 01840 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ ); 01841 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp ); 01842 // Initialize AUTAU ( AUTAU = AU'*(M\AU) ) 01843 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ ); 01844 if ( precObj != Teuchos::null ) { 01845 index.resize(recycleBlocks_); 01846 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01847 index.resize(recycleBlocks_); 01848 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01849 Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage 01850 OPT::Apply( *precObj, *AUtmp, *LeftPCAU ); 01851 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp ); 01852 } else { 01853 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp ); 01854 } 01855 } // if (numRHS2Solve > 0) 01856 01857 } // if (existU1) 01858 }// while ( numRHS2Solve > 0 ) 01859 01860 } 01861 01862 // print final summary 01863 sTest_->print( printer_->stream(FinalSummary) ); 01864 01865 // print timing information 01866 #ifdef BELOS_TEUCHOS_TIME_MONITOR 01867 // Calling summarize() can be expensive, so don't call unless the 01868 // user wants to print out timing details. summarize() will do all 01869 // the work even if it's passed a "black hole" output stream. 01870 if (verbosity_ & TimingDetails) 01871 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) ); 01872 #endif 01873 01874 // get iteration information for this solve 01875 numIters_ = maxIterTest_->getNumIters(); 01876 01877 // Save the convergence test value ("achieved tolerance") for this solve. 01878 { 01879 using Teuchos::rcp_dynamic_cast; 01880 typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type; 01881 // testValues is nonnull and not persistent. 01882 const std::vector<MagnitudeType>* pTestValues = 01883 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue(); 01884 01885 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error, 01886 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 01887 "method returned NULL. Please report this bug to the Belos developers."); 01888 01889 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error, 01890 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() " 01891 "method returned a vector of length zero. Please report this bug to the " 01892 "Belos developers."); 01893 01894 // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the 01895 // achieved tolerances for all vectors in the current solve(), or 01896 // just for the vectors from the last deflation? 01897 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end()); 01898 } 01899 01900 if (!isConverged) { 01901 return Unconverged; // return from RCGSolMgr::solve() 01902 } 01903 return Converged; // return from RCGSolMgr::solve() 01904 } 01905 01906 // Compute the harmonic eigenpairs of the projected, dense system. 01907 template<class ScalarType, class MV, class OP> 01908 void RCGSolMgr<ScalarType,MV,OP,false>::getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType>& F, 01909 const Teuchos::SerialDenseMatrix<int,ScalarType>& G, 01910 Teuchos::SerialDenseMatrix<int,ScalarType>& Y) { 01911 // order of F,G 01912 int n = F.numCols(); 01913 01914 // The LAPACK interface 01915 Teuchos::LAPACK<int,ScalarType> lapack; 01916 01917 // Magnitude of harmonic Ritz values 01918 std::vector<MagnitudeType> w(n); 01919 01920 // Sorted order of harmonic Ritz values 01921 std::vector<int> iperm(n); 01922 01923 // Compute k smallest harmonic Ritz pairs 01924 // SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, INFO ) 01925 int itype = 1; // solve A*x = (lambda)*B*x 01926 char jobz='V'; // compute eigenvalues and eigenvectors 01927 char uplo='U'; // since F,G symmetric, reference only their upper triangular data 01928 std::vector<ScalarType> work(1); 01929 int lwork = -1; 01930 int info = 0; 01931 // since SYGV destroys workspace, create copies of F,G 01932 Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_, F_->numRows(), F_->numCols() ); 01933 Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *G_, G_->numRows(), G_->numCols() ); 01934 01935 // query for optimal workspace size 01936 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info); 01937 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure, 01938 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size."); 01939 lwork = (int)work[0]; 01940 work.resize(lwork); 01941 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info); 01942 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure, 01943 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions."); 01944 01945 01946 // Construct magnitude of each harmonic Ritz value 01947 this->sort(w,n,iperm); 01948 01949 // Select recycledBlocks_ smallest eigenvectors 01950 for( int i=0; i<recycleBlocks_; i++ ) { 01951 for( int j=0; j<n; j++ ) { 01952 Y(j,i) = G2(j,iperm[i]); 01953 } 01954 } 01955 01956 } 01957 01958 // This method sorts list of n floating-point numbers and return permutation vector 01959 template<class ScalarType, class MV, class OP> 01960 void RCGSolMgr<ScalarType,MV,OP,false>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm) 01961 { 01962 int l, r, j, i, flag; 01963 int RR2; 01964 double dRR, dK; 01965 01966 // Initialize the permutation vector. 01967 for(j=0;j<n;j++) 01968 iperm[j] = j; 01969 01970 if (n <= 1) return; 01971 01972 l = n / 2 + 1; 01973 r = n - 1; 01974 l = l - 1; 01975 dRR = dlist[l - 1]; 01976 dK = dlist[l - 1]; 01977 01978 RR2 = iperm[l - 1]; 01979 while (r != 0) { 01980 j = l; 01981 flag = 1; 01982 01983 while (flag == 1) { 01984 i = j; 01985 j = j + j; 01986 01987 if (j > r + 1) 01988 flag = 0; 01989 else { 01990 if (j < r + 1) 01991 if (dlist[j] > dlist[j - 1]) j = j + 1; 01992 01993 if (dlist[j - 1] > dK) { 01994 dlist[i - 1] = dlist[j - 1]; 01995 iperm[i - 1] = iperm[j - 1]; 01996 } 01997 else { 01998 flag = 0; 01999 } 02000 } 02001 } 02002 dlist[i - 1] = dRR; 02003 iperm[i - 1] = RR2; 02004 if (l == 1) { 02005 dRR = dlist [r]; 02006 RR2 = iperm[r]; 02007 dK = dlist[r]; 02008 dlist[r] = dlist[0]; 02009 iperm[r] = iperm[0]; 02010 r = r - 1; 02011 } 02012 else { 02013 l = l - 1; 02014 dRR = dlist[l - 1]; 02015 RR2 = iperm[l - 1]; 02016 dK = dlist[l - 1]; 02017 } 02018 } 02019 dlist[0] = dRR; 02020 iperm[0] = RR2; 02021 } 02022 02023 // This method requires the solver manager to return a std::string that describes itself. 02024 template<class ScalarType, class MV, class OP> 02025 std::string RCGSolMgr<ScalarType,MV,OP,false>::description() const 02026 { 02027 std::ostringstream oss; 02028 oss << "Belos::RCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">"; 02029 return oss.str(); 02030 } 02031 02032 } // end Belos namespace 02033 02034 #endif /* BELOS_RCG_SOLMGR_HPP */
1.7.6.1
