|
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_TimeMonitor.hpp" 00067 00077 namespace Belos { 00078 00080 00081 00088 class RCGSolMgrLinearProblemFailure : public BelosError {public: 00089 RCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg) 00090 {}}; 00091 00098 class RCGSolMgrLAPACKFailure : public BelosError {public: 00099 RCGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg) 00100 {}}; 00101 00108 class RCGSolMgrRecyclingFailure : public BelosError {public: 00109 RCGSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg) 00110 {}}; 00111 00113 00114 template<class ScalarType, class MV, class OP> 00115 class RCGSolMgr : public SolverManager<ScalarType,MV,OP> { 00116 00117 private: 00118 typedef MultiVecTraits<ScalarType,MV> MVT; 00119 typedef OperatorTraits<ScalarType,MV,OP> OPT; 00120 typedef Teuchos::ScalarTraits<ScalarType> SCT; 00121 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType; 00122 typedef Teuchos::ScalarTraits<MagnitudeType> MT; 00123 00124 public: 00125 00127 00128 00134 RCGSolMgr(); 00135 00157 RCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 00158 const Teuchos::RCP<Teuchos::ParameterList> &pl ); 00159 00161 virtual ~RCGSolMgr() {}; 00163 00165 00166 00167 const LinearProblem<ScalarType,MV,OP>& getProblem() const { 00168 return *problem_; 00169 } 00170 00172 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const; 00173 00175 Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; } 00176 00182 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const { 00183 return tuple(timerSolve_); 00184 } 00185 00187 int getNumIters() const { 00188 return numIters_; 00189 } 00190 00192 bool isLOADetected() const { return false; } 00193 00195 00197 00198 00200 void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; } 00201 00203 void setParameters( const Teuchos::RCP<Teuchos::ParameterList> ¶ms ); 00204 00206 00208 00209 00215 void reset( const ResetType type ) { 00216 if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); 00217 else if (type & Belos::RecycleSubspace) existU_ = false; 00218 } 00220 00222 00223 00241 ReturnType solve(); 00242 00244 00247 00249 std::string description() const; 00250 00252 00253 private: 00254 00255 // Called by all constructors; Contains init instructions common to all constructors 00256 void init(); 00257 00258 // Computes harmonic eigenpairs of projected matrix created during one cycle. 00259 // Y contains the harmonic Ritz vectors corresponding to the recycleBlocks eigenvalues of smallest magnitude. 00260 void getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType> &F, 00261 const Teuchos::SerialDenseMatrix<int,ScalarType> &G, 00262 Teuchos::SerialDenseMatrix<int,ScalarType>& Y); 00263 00264 // Sort list of n floating-point numbers and return permutation vector 00265 void sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm); 00266 00267 // Initialize solver state storage 00268 void initializeStateStorage(); 00269 00270 // Linear problem. 00271 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_; 00272 00273 // Output manager. 00274 Teuchos::RCP<OutputManager<ScalarType> > printer_; 00275 Teuchos::RCP<std::ostream> outputStream_; 00276 00277 // Status test. 00278 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_; 00279 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_; 00280 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_; 00281 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_; 00282 00283 // Current parameter list. 00284 Teuchos::RCP<ParameterList> params_; 00285 00286 // Default solver values. 00287 static const MagnitudeType convtol_default_; 00288 static const int maxIters_default_; 00289 static const int numBlocks_default_; 00290 static const int recycleBlocks_default_; 00291 static const bool showMaxResNormOnly_default_; 00292 static const int verbosity_default_; 00293 static const int outputStyle_default_; 00294 static const int outputFreq_default_; 00295 static const std::string label_default_; 00296 static const Teuchos::RCP<std::ostream> outputStream_default_; 00297 00298 // Current solver values. 00299 MagnitudeType convtol_; 00300 int maxIters_, numIters_; 00301 int numBlocks_, recycleBlocks_; 00302 bool showMaxResNormOnly_; 00303 int verbosity_, outputStyle_, outputFreq_; 00304 00306 // Solver State Storage 00308 // Search vectors 00309 Teuchos::RCP<MV> P_; 00310 // 00311 // A times current search direction 00312 Teuchos::RCP<MV> Ap_; 00313 // 00314 // Residual vector 00315 Teuchos::RCP<MV> r_; 00316 // 00317 // Preconditioned residual 00318 Teuchos::RCP<MV> z_; 00319 // 00320 // Flag indicating that the recycle space should be used 00321 bool existU_; 00322 // 00323 // Flag indicating that the updated recycle space has been created 00324 bool existU1_; 00325 // 00326 // Recycled subspace and its image 00327 Teuchos::RCP<MV> U_, AU_; 00328 // 00329 // Recycled subspace for next system and its image 00330 Teuchos::RCP<MV> U1_; 00331 // 00332 // Coefficients arising in RCG iteration 00333 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Alpha_; 00334 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Beta_; 00335 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_; 00336 // 00337 // Solutions to local least-squares problems 00338 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > Delta_; 00339 // 00340 // The matrix U^T A U 00341 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > UTAU_; 00342 // 00343 // LU factorization of U^T A U 00344 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > LUUTAU_; 00345 // 00346 // Data from LU factorization of UTAU 00347 Teuchos::RCP<std::vector<int> > ipiv_; 00348 // 00349 // The matrix (AU)^T AU 00350 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAU_; 00351 // 00352 // The scalar r'*z 00353 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > rTz_old_; 00354 // 00355 // Matrices needed for calculation of harmonic Ritz eigenproblem 00356 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > F_,G_,Y_; 00357 // 00358 // Matrices needed for updating recycle space 00359 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > L2_,DeltaL2_,AU1TUDeltaL2_; 00360 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AU1TAU1_, AU1TU1_, AU1TAP_; 00361 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > FY_,GY_; 00362 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > APTAP_; 00363 Teuchos::RCP<MV> U1Y1_, PY2_; 00364 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > AUTAP_, AU1TU_; 00365 ScalarType dold; 00367 00368 // Timers. 00369 std::string label_; 00370 Teuchos::RCP<Teuchos::Time> timerSolve_; 00371 00372 // Internal state variables. 00373 bool params_Set_; 00374 }; 00375 00376 00377 // Default solver values. 00378 template<class ScalarType, class MV, class OP> 00379 const typename RCGSolMgr<ScalarType,MV,OP>::MagnitudeType RCGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8; 00380 00381 template<class ScalarType, class MV, class OP> 00382 const int RCGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000; 00383 00384 template<class ScalarType, class MV, class OP> 00385 const int RCGSolMgr<ScalarType,MV,OP>::numBlocks_default_ = 25; 00386 00387 template<class ScalarType, class MV, class OP> 00388 const int RCGSolMgr<ScalarType,MV,OP>::recycleBlocks_default_ = 3; 00389 00390 template<class ScalarType, class MV, class OP> 00391 const bool RCGSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false; 00392 00393 template<class ScalarType, class MV, class OP> 00394 const int RCGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors; 00395 00396 template<class ScalarType, class MV, class OP> 00397 const int RCGSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General; 00398 00399 template<class ScalarType, class MV, class OP> 00400 const int RCGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1; 00401 00402 template<class ScalarType, class MV, class OP> 00403 const std::string RCGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos"; 00404 00405 template<class ScalarType, class MV, class OP> 00406 const Teuchos::RCP<std::ostream> RCGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false); 00407 00408 // Empty Constructor 00409 template<class ScalarType, class MV, class OP> 00410 RCGSolMgr<ScalarType,MV,OP>::RCGSolMgr() { 00411 init(); 00412 } 00413 00414 // Basic Constructor 00415 template<class ScalarType, class MV, class OP> 00416 RCGSolMgr<ScalarType,MV,OP>::RCGSolMgr( 00417 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 00418 const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 00419 problem_(problem) 00420 { 00421 init(); 00422 TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager."); 00423 00424 // If the parameter list pointer is null, then set the current parameters to the default parameter list. 00425 if ( !is_null(pl) ) { 00426 setParameters( pl ); 00427 } 00428 } 00429 00430 // Common instructions executed in all constructors 00431 template<class ScalarType, class MV, class OP> 00432 void RCGSolMgr<ScalarType,MV,OP>::init() 00433 { 00434 outputStream_ = outputStream_default_; 00435 convtol_ = convtol_default_; 00436 maxIters_ = maxIters_default_; 00437 numBlocks_ = numBlocks_default_; 00438 recycleBlocks_ = recycleBlocks_default_; 00439 verbosity_ = verbosity_default_; 00440 outputStyle_= outputStyle_default_; 00441 outputFreq_= outputFreq_default_; 00442 showMaxResNormOnly_ = showMaxResNormOnly_default_; 00443 label_ = label_default_; 00444 params_Set_ = false; 00445 P_ = Teuchos::null; 00446 Ap_ = Teuchos::null; 00447 r_ = Teuchos::null; 00448 z_ = Teuchos::null; 00449 existU_ = false; 00450 existU1_ = false; 00451 U_ = Teuchos::null; 00452 AU_ = Teuchos::null; 00453 U1_ = Teuchos::null; 00454 Alpha_ = Teuchos::null; 00455 Beta_ = Teuchos::null; 00456 D_ = Teuchos::null; 00457 Delta_ = Teuchos::null; 00458 UTAU_ = Teuchos::null; 00459 LUUTAU_ = Teuchos::null; 00460 ipiv_ = Teuchos::null; 00461 AUTAU_ = Teuchos::null; 00462 rTz_old_ = Teuchos::null; 00463 F_ = Teuchos::null; 00464 G_ = Teuchos::null; 00465 Y_ = Teuchos::null; 00466 L2_ = Teuchos::null; 00467 DeltaL2_ = Teuchos::null; 00468 AU1TUDeltaL2_ = Teuchos::null; 00469 AU1TAU1_ = Teuchos::null; 00470 AU1TU1_ = Teuchos::null; 00471 AU1TAP_ = Teuchos::null; 00472 FY_ = Teuchos::null; 00473 GY_ = Teuchos::null; 00474 APTAP_ = Teuchos::null; 00475 U1Y1_ = Teuchos::null; 00476 PY2_ = Teuchos::null; 00477 AUTAP_ = Teuchos::null; 00478 AU1TU_ = Teuchos::null; 00479 dold = 0.; 00480 } 00481 00482 template<class ScalarType, class MV, class OP> 00483 void RCGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> ¶ms ) 00484 { 00485 // Create the internal parameter list if ones doesn't already exist. 00486 if (params_ == Teuchos::null) { 00487 params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) ); 00488 } 00489 else { 00490 params->validateParameters(*getValidParameters()); 00491 } 00492 00493 // Check for maximum number of iterations 00494 if (params->isParameter("Maximum Iterations")) { 00495 maxIters_ = params->get("Maximum Iterations",maxIters_default_); 00496 00497 // Update parameter in our list and in status test. 00498 params_->set("Maximum Iterations", maxIters_); 00499 if (maxIterTest_!=Teuchos::null) 00500 maxIterTest_->setMaxIters( maxIters_ ); 00501 } 00502 00503 // Check for the maximum number of blocks. 00504 if (params->isParameter("Num Blocks")) { 00505 numBlocks_ = params->get("Num Blocks",numBlocks_default_); 00506 TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument, 00507 "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive."); 00508 00509 // Update parameter in our list. 00510 params_->set("Num Blocks", numBlocks_); 00511 } 00512 00513 // Check for the maximum number of blocks. 00514 if (params->isParameter("Num Recycled Blocks")) { 00515 recycleBlocks_ = params->get("Num Recycled Blocks",recycleBlocks_default_); 00516 TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument, 00517 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive."); 00518 00519 TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument, 00520 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\"."); 00521 00522 // Update parameter in our list. 00523 params_->set("Num Recycled Blocks", recycleBlocks_); 00524 } 00525 00526 // Check to see if the timer label changed. 00527 if (params->isParameter("Timer Label")) { 00528 std::string tempLabel = params->get("Timer Label", label_default_); 00529 00530 // Update parameter in our list and solver timer 00531 if (tempLabel != label_) { 00532 label_ = tempLabel; 00533 params_->set("Timer Label", label_); 00534 std::string solveLabel = label_ + ": RCGSolMgr total solve time"; 00535 timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel); 00536 } 00537 } 00538 00539 // Check for a change in verbosity level 00540 if (params->isParameter("Verbosity")) { 00541 if (Teuchos::isParameterType<int>(*params,"Verbosity")) { 00542 verbosity_ = params->get("Verbosity", verbosity_default_); 00543 } else { 00544 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity"); 00545 } 00546 00547 // Update parameter in our list. 00548 params_->set("Verbosity", verbosity_); 00549 if (printer_ != Teuchos::null) 00550 printer_->setVerbosity(verbosity_); 00551 } 00552 00553 // Check for a change in output style 00554 if (params->isParameter("Output Style")) { 00555 if (Teuchos::isParameterType<int>(*params,"Output Style")) { 00556 outputStyle_ = params->get("Output Style", outputStyle_default_); 00557 } else { 00558 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style"); 00559 } 00560 00561 // Reconstruct the convergence test if the explicit residual test is not being used. 00562 params_->set("Output Style", outputStyle_); 00563 outputTest_ = Teuchos::null; 00564 } 00565 00566 // output stream 00567 if (params->isParameter("Output Stream")) { 00568 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream"); 00569 00570 // Update parameter in our list. 00571 params_->set("Output Stream", outputStream_); 00572 if (printer_ != Teuchos::null) 00573 printer_->setOStream( outputStream_ ); 00574 } 00575 00576 // frequency level 00577 if (verbosity_ & Belos::StatusTestDetails) { 00578 if (params->isParameter("Output Frequency")) { 00579 outputFreq_ = params->get("Output Frequency", outputFreq_default_); 00580 } 00581 00582 // Update parameter in out list and output status test. 00583 params_->set("Output Frequency", outputFreq_); 00584 if (outputTest_ != Teuchos::null) 00585 outputTest_->setOutputFrequency( outputFreq_ ); 00586 } 00587 00588 // Create output manager if we need to. 00589 if (printer_ == Teuchos::null) { 00590 printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) ); 00591 } 00592 00593 // Convergence 00594 typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t; 00595 typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t; 00596 00597 // Check for convergence tolerance 00598 if (params->isParameter("Convergence Tolerance")) { 00599 convtol_ = params->get("Convergence Tolerance",convtol_default_); 00600 00601 // Update parameter in our list and residual tests. 00602 params_->set("Convergence Tolerance", convtol_); 00603 if (convTest_ != Teuchos::null) 00604 convTest_->setTolerance( convtol_ ); 00605 } 00606 00607 if (params->isParameter("Show Maximum Residual Norm Only")) { 00608 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only"); 00609 00610 // Update parameter in our list and residual tests 00611 params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_); 00612 if (convTest_ != Teuchos::null) 00613 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ ); 00614 } 00615 00616 // Create status tests if we need to. 00617 00618 // Basic test checks maximum iterations and native residual. 00619 if (maxIterTest_ == Teuchos::null) 00620 maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) ); 00621 00622 // Implicit residual test, using the native residual to determine if convergence was achieved. 00623 if (convTest_ == Teuchos::null) 00624 convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) ); 00625 00626 if (sTest_ == Teuchos::null) 00627 sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) ); 00628 00629 if (outputTest_ == Teuchos::null) { 00630 00631 // Create the status test output class. 00632 // This class manages and formats the output from the status test. 00633 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ ); 00634 outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined ); 00635 00636 // Set the solver string for the output test 00637 std::string solverDesc = " Recycling CG "; 00638 outputTest_->setSolverDesc( solverDesc ); 00639 } 00640 00641 // Create the timer if we need to. 00642 if (timerSolve_ == Teuchos::null) { 00643 std::string solveLabel = label_ + ": RCGSolMgr total solve time"; 00644 timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel); 00645 } 00646 00647 // Inform the solver manager that the current parameters were set. 00648 params_Set_ = true; 00649 } 00650 00651 00652 template<class ScalarType, class MV, class OP> 00653 Teuchos::RCP<const Teuchos::ParameterList> 00654 RCGSolMgr<ScalarType,MV,OP>::getValidParameters() const 00655 { 00656 static Teuchos::RCP<const Teuchos::ParameterList> validPL; 00657 00658 // Set all the valid parameters and their default values. 00659 if(is_null(validPL)) { 00660 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList(); 00661 pl->set("Convergence Tolerance", convtol_default_, 00662 "The relative residual tolerance that needs to be achieved by the\n" 00663 "iterative solver in order for the linear system to be declared converged."); 00664 pl->set("Maximum Iterations", maxIters_default_, 00665 "The maximum number of block iterations allowed for each\n" 00666 "set of RHS solved."); 00667 pl->set("Num Blocks", numBlocks_default_, 00668 "The length of a cycle (and this max number of search vectors kept)\n"); 00669 pl->set("Num Recycled Blocks", recycleBlocks_default_, 00670 "The number of vectors in the recycle subspace."); 00671 pl->set("Verbosity", verbosity_default_, 00672 "What type(s) of solver information should be outputted\n" 00673 "to the output stream."); 00674 pl->set("Output Style", outputStyle_default_, 00675 "What style is used for the solver information outputted\n" 00676 "to the output stream."); 00677 pl->set("Output Frequency", outputFreq_default_, 00678 "How often convergence information should be outputted\n" 00679 "to the output stream."); 00680 pl->set("Output Stream", outputStream_default_, 00681 "A reference-counted pointer to the output stream where all\n" 00682 "solver output is sent."); 00683 pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_, 00684 "When convergence information is printed, only show the maximum\n" 00685 "relative residual norm when the block size is greater than one."); 00686 pl->set("Timer Label", label_default_, 00687 "The string to use as a prefix for the timer labels."); 00688 // pl->set("Restart Timers", restartTimers_); 00689 validPL = pl; 00690 } 00691 return validPL; 00692 } 00693 00694 // initializeStateStorage 00695 template<class ScalarType, class MV, class OP> 00696 void RCGSolMgr<ScalarType,MV,OP>::initializeStateStorage() { 00697 00698 // Check if there is any multivector to clone from. 00699 Teuchos::RCP<const MV> rhsMV = problem_->getRHS(); 00700 if (rhsMV == Teuchos::null) { 00701 // Nothing to do 00702 return; 00703 } 00704 else { 00705 00706 // Initialize the state storage 00707 TEST_FOR_EXCEPTION(numBlocks_ > MVT::GetVecLength(*rhsMV),std::invalid_argument, 00708 "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!"); 00709 00710 // If the subspace has not been initialized before, generate it using the RHS from lp_. 00711 if (P_ == Teuchos::null) { 00712 P_ = MVT::Clone( *rhsMV, numBlocks_+2 ); 00713 } 00714 else { 00715 // Generate P_ by cloning itself ONLY if more space is needed. 00716 if (MVT::GetNumberVecs(*P_) < numBlocks_+2) { 00717 Teuchos::RCP<const MV> tmp = P_; 00718 P_ = MVT::Clone( *tmp, numBlocks_+2 ); 00719 } 00720 } 00721 00722 // Generate Ap_ only if it doesn't exist 00723 if (Ap_ == Teuchos::null) 00724 Ap_ = MVT::Clone( *rhsMV, 1 ); 00725 00726 // Generate r_ only if it doesn't exist 00727 if (r_ == Teuchos::null) 00728 r_ = MVT::Clone( *rhsMV, 1 ); 00729 00730 // Generate z_ only if it doesn't exist 00731 if (z_ == Teuchos::null) 00732 z_ = MVT::Clone( *rhsMV, 1 ); 00733 00734 // If the recycle space has not been initialized before, generate it using the RHS from lp_. 00735 if (U_ == Teuchos::null) { 00736 U_ = MVT::Clone( *rhsMV, recycleBlocks_ ); 00737 } 00738 else { 00739 // Generate U_ by cloning itself ONLY if more space is needed. 00740 if (MVT::GetNumberVecs(*U_) < recycleBlocks_) { 00741 Teuchos::RCP<const MV> tmp = U_; 00742 U_ = MVT::Clone( *tmp, recycleBlocks_ ); 00743 } 00744 } 00745 00746 // If the recycle space has not be initialized before, generate it using the RHS from lp_. 00747 if (AU_ == Teuchos::null) { 00748 AU_ = MVT::Clone( *rhsMV, recycleBlocks_ ); 00749 } 00750 else { 00751 // Generate AU_ by cloning itself ONLY if more space is needed. 00752 if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) { 00753 Teuchos::RCP<const MV> tmp = AU_; 00754 AU_ = MVT::Clone( *tmp, recycleBlocks_ ); 00755 } 00756 } 00757 00758 // If the recycle space has not been initialized before, generate it using the RHS from lp_. 00759 if (U1_ == Teuchos::null) { 00760 U1_ = MVT::Clone( *rhsMV, recycleBlocks_ ); 00761 } 00762 else { 00763 // Generate U1_ by cloning itself ONLY if more space is needed. 00764 if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) { 00765 Teuchos::RCP<const MV> tmp = U1_; 00766 U1_ = MVT::Clone( *tmp, recycleBlocks_ ); 00767 } 00768 } 00769 00770 // Generate Alpha_ only if it doesn't exist, otherwise resize it. 00771 if (Alpha_ == Teuchos::null) 00772 Alpha_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, 1 ) ); 00773 else { 00774 if ( (Alpha_->numRows() != numBlocks_) || (Alpha_->numCols() != 1) ) 00775 Alpha_->reshape( numBlocks_, 1 ); 00776 } 00777 00778 // Generate Beta_ only if it doesn't exist, otherwise resize it. 00779 if (Beta_ == Teuchos::null) 00780 Beta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + 1, 1 ) ); 00781 else { 00782 if ( (Beta_->numRows() != (numBlocks_+1)) || (Beta_->numCols() != 1) ) 00783 Beta_->reshape( numBlocks_ + 1, 1 ); 00784 } 00785 00786 // Generate D_ only if it doesn't exist, otherwise resize it. 00787 if (D_ == Teuchos::null) 00788 D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ , 1 ) ); 00789 else { 00790 if ( (D_->numRows() != numBlocks_) || (D_->numCols() != 1) ) 00791 D_->reshape( numBlocks_, 1 ); 00792 } 00793 00794 // Generate Delta_ only if it doesn't exist, otherwise resize it. 00795 if (Delta_ == Teuchos::null) 00796 Delta_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ + 1 ) ); 00797 else { 00798 if ( (Delta_->numRows() != recycleBlocks_) || (Delta_->numCols() != (numBlocks_ + 1)) ) 00799 Delta_->reshape( recycleBlocks_, numBlocks_ + 1 ); 00800 } 00801 00802 // Generate UTAU_ only if it doesn't exist, otherwise resize it. 00803 if (UTAU_ == Teuchos::null) 00804 UTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 00805 else { 00806 if ( (UTAU_->numRows() != recycleBlocks_) || (UTAU_->numCols() != recycleBlocks_) ) 00807 UTAU_->reshape( recycleBlocks_, recycleBlocks_ ); 00808 } 00809 00810 // Generate LUUTAU_ only if it doesn't exist, otherwise resize it. 00811 if (LUUTAU_ == Teuchos::null) 00812 LUUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 00813 else { 00814 if ( (LUUTAU_->numRows() != recycleBlocks_) || (LUUTAU_->numCols() != recycleBlocks_) ) 00815 LUUTAU_->reshape( recycleBlocks_, recycleBlocks_ ); 00816 } 00817 00818 // Generate ipiv_ only if it doesn't exist, otherwise resize it. 00819 if (ipiv_ == Teuchos::null) 00820 ipiv_ = Teuchos::rcp( new std::vector<int>(recycleBlocks_) ); 00821 else { 00822 if ( (int)ipiv_->size() != recycleBlocks_ ) // if ipiv not correct size, always resize it 00823 ipiv_->resize(recycleBlocks_); 00824 } 00825 00826 // Generate AUTAU_ only if it doesn't exist, otherwise resize it. 00827 if (AUTAU_ == Teuchos::null) 00828 AUTAU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 00829 else { 00830 if ( (AUTAU_->numRows() != recycleBlocks_) || (AUTAU_->numCols() != recycleBlocks_) ) 00831 AUTAU_->reshape( recycleBlocks_, recycleBlocks_ ); 00832 } 00833 00834 // Generate rTz_old_ only if it doesn't exist 00835 if (rTz_old_ == Teuchos::null) 00836 rTz_old_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( 1, 1 ) ); 00837 else { 00838 if ( (rTz_old_->numRows() != 1) || (rTz_old_->numCols() != 1) ) 00839 rTz_old_->reshape( 1, 1 ); 00840 } 00841 00842 // Generate F_ only if it doesn't exist 00843 if (F_ == Teuchos::null) 00844 F_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) ); 00845 else { 00846 if ( (F_->numRows() != (numBlocks_+recycleBlocks_)) || (F_->numCols() != numBlocks_+recycleBlocks_) ) 00847 F_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ); 00848 } 00849 00850 // Generate G_ only if it doesn't exist 00851 if (G_ == Teuchos::null) 00852 G_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ) ); 00853 else { 00854 if ( (G_->numRows() != (numBlocks_+recycleBlocks_)) || (G_->numCols() != numBlocks_+recycleBlocks_) ) 00855 G_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ); 00856 } 00857 00858 // Generate Y_ only if it doesn't exist 00859 if (Y_ == Teuchos::null) 00860 Y_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycleBlocks_, recycleBlocks_ ) ); 00861 else { 00862 if ( (Y_->numRows() != (numBlocks_+recycleBlocks_)) || (Y_->numCols() != numBlocks_+recycleBlocks_) ) 00863 Y_->reshape( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ ); 00864 } 00865 00866 // Generate L2_ only if it doesn't exist 00867 if (L2_ == Teuchos::null) 00868 L2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+1, numBlocks_ ) ); 00869 else { 00870 if ( (L2_->numRows() != (numBlocks_+1)) || (L2_->numCols() != numBlocks_) ) 00871 L2_->reshape( numBlocks_+1, numBlocks_ ); 00872 } 00873 00874 // Generate DeltaL2_ only if it doesn't exist 00875 if (DeltaL2_ == Teuchos::null) 00876 DeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) ); 00877 else { 00878 if ( (DeltaL2_->numRows() != (recycleBlocks_)) || (DeltaL2_->numCols() != (numBlocks_) ) ) 00879 DeltaL2_->reshape( recycleBlocks_, numBlocks_ ); 00880 } 00881 00882 // Generate AU1TUDeltaL2_ only if it doesn't exist 00883 if (AU1TUDeltaL2_ == Teuchos::null) 00884 AU1TUDeltaL2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) ); 00885 else { 00886 if ( (AU1TUDeltaL2_->numRows() != (recycleBlocks_)) || (AU1TUDeltaL2_->numCols() != (numBlocks_) ) ) 00887 AU1TUDeltaL2_->reshape( recycleBlocks_, numBlocks_ ); 00888 } 00889 00890 // Generate AU1TAU1_ only if it doesn't exist 00891 if (AU1TAU1_ == Teuchos::null) 00892 AU1TAU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 00893 else { 00894 if ( (AU1TAU1_->numRows() != (recycleBlocks_)) || (AU1TAU1_->numCols() != (recycleBlocks_) ) ) 00895 AU1TAU1_->reshape( recycleBlocks_, recycleBlocks_ ); 00896 } 00897 00898 // Generate GY_ only if it doesn't exist 00899 if (GY_ == Teuchos::null) 00900 GY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) ); 00901 else { 00902 if ( (GY_->numRows() != (numBlocks_ + recycleBlocks_)) || (GY_->numCols() != (recycleBlocks_) ) ) 00903 GY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ ); 00904 } 00905 00906 // Generate AU1TU1_ only if it doesn't exist 00907 if (AU1TU1_ == Teuchos::null) 00908 AU1TU1_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 00909 else { 00910 if ( (AU1TU1_->numRows() != (recycleBlocks_)) || (AU1TU1_->numCols() != (recycleBlocks_) ) ) 00911 AU1TU1_->reshape( recycleBlocks_, recycleBlocks_ ); 00912 } 00913 00914 // Generate FY_ only if it doesn't exist 00915 if (FY_ == Teuchos::null) 00916 FY_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_ + recycleBlocks_, recycleBlocks_ ) ); 00917 else { 00918 if ( (FY_->numRows() != (numBlocks_ + recycleBlocks_)) || (FY_->numCols() != (recycleBlocks_) ) ) 00919 FY_->reshape( numBlocks_+recycleBlocks_, recycleBlocks_ ); 00920 } 00921 00922 // Generate AU1TAP_ only if it doesn't exist 00923 if (AU1TAP_ == Teuchos::null) 00924 AU1TAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) ); 00925 else { 00926 if ( (AU1TAP_->numRows() != (recycleBlocks_)) || (AU1TAP_->numCols() != (numBlocks_) ) ) 00927 AU1TAP_->reshape( recycleBlocks_, numBlocks_ ); 00928 } 00929 00930 // Generate APTAP_ only if it doesn't exist 00931 if (APTAP_ == Teuchos::null) 00932 APTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_, numBlocks_ ) ); 00933 else { 00934 if ( (APTAP_->numRows() != (numBlocks_)) || (APTAP_->numCols() != (numBlocks_) ) ) 00935 APTAP_->reshape( numBlocks_, numBlocks_ ); 00936 } 00937 00938 // If the subspace has not been initialized before, generate it using the RHS from lp_. 00939 if (U1Y1_ == Teuchos::null) { 00940 U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ ); 00941 } 00942 else { 00943 // Generate U1Y1_ by cloning itself ONLY if more space is needed. 00944 if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) { 00945 Teuchos::RCP<const MV> tmp = U1Y1_; 00946 U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ ); 00947 } 00948 } 00949 00950 // If the subspace has not been initialized before, generate it using the RHS from lp_. 00951 if (PY2_ == Teuchos::null) { 00952 PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ ); 00953 } 00954 else { 00955 // Generate PY2_ by cloning itself ONLY if more space is needed. 00956 if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) { 00957 Teuchos::RCP<const MV> tmp = PY2_; 00958 PY2_ = MVT::Clone( *tmp, recycleBlocks_ ); 00959 } 00960 } 00961 00962 // Generate AUTAP_ only if it doesn't exist 00963 if (AUTAP_ == Teuchos::null) 00964 AUTAP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, numBlocks_ ) ); 00965 else { 00966 if ( (AUTAP_->numRows() != (recycleBlocks_)) || (AUTAP_->numCols() != (numBlocks_) ) ) 00967 AUTAP_->reshape( recycleBlocks_, numBlocks_ ); 00968 } 00969 00970 // Generate AU1TU_ only if it doesn't exist 00971 if (AU1TU_ == Teuchos::null) 00972 AU1TU_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycleBlocks_, recycleBlocks_ ) ); 00973 else { 00974 if ( (AU1TU_->numRows() != (recycleBlocks_)) || (AU1TU_->numCols() != (recycleBlocks_) ) ) 00975 AU1TU_->reshape( recycleBlocks_, recycleBlocks_ ); 00976 } 00977 00978 00979 } 00980 } 00981 00982 template<class ScalarType, class MV, class OP> 00983 ReturnType RCGSolMgr<ScalarType,MV,OP>::solve() { 00984 00985 Teuchos::LAPACK<int,ScalarType> lapack; 00986 std::vector<int> index(recycleBlocks_); 00987 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one(); 00988 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero(); 00989 00990 // Count of number of cycles performed on current rhs 00991 int cycle = 0; 00992 00993 // Set the current parameters if they were not set before. 00994 // NOTE: This may occur if the user generated the solver manager with the default constructor and 00995 // then didn't set any parameters using setParameters(). 00996 if (!params_Set_) { 00997 setParameters(Teuchos::parameterList(*getValidParameters())); 00998 } 00999 01000 TEST_FOR_EXCEPTION(problem_ == Teuchos::null,RCGSolMgrLinearProblemFailure, 01001 "Belos::RCGSolMgr::solve(): Linear problem is not a valid object."); 01002 TEST_FOR_EXCEPTION(!problem_->isProblemSet(),RCGSolMgrLinearProblemFailure, 01003 "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called."); 01004 01005 // Create indices for the linear systems to be solved. 01006 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) ); 01007 std::vector<int> currIdx(1); 01008 currIdx[0] = 0; 01009 01010 // Inform the linear problem of the current linear system to solve. 01011 problem_->setLSIndex( currIdx ); 01012 01013 // Check the number of blocks and change them if necessary. 01014 int dim = MVT::GetVecLength( *(problem_->getRHS()) ); 01015 if (numBlocks_ > dim) { 01016 numBlocks_ = dim; 01017 params_->set("Num Blocks", numBlocks_); 01018 printer_->stream(Warnings) << 01019 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl << 01020 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl; 01021 } 01022 01023 // Initialize storage for all state variables 01024 initializeStateStorage(); 01025 01026 // Parameter list 01027 Teuchos::ParameterList plist; 01028 plist.set("Num Blocks",numBlocks_); 01029 plist.set("Recycled Blocks",recycleBlocks_); 01030 01031 // Reset the status test. 01032 outputTest_->reset(); 01033 01034 // Assume convergence is achieved, then let any failed convergence set this to false. 01035 bool isConverged = true; 01036 01037 // Compute AU = A*U, UTAU = U'*AU, AUTAU = (AU)'*(AU) 01038 if (existU_) { 01039 index.resize(recycleBlocks_); 01040 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01041 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index ); 01042 index.resize(recycleBlocks_); 01043 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01044 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index ); 01045 // Initialize AU 01046 problem_->applyOp( *Utmp, *AUtmp ); 01047 // Initialize UTAU 01048 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ ); 01049 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp ); 01050 // Initialize AUTAU ( AUTAU = AU'*(M\AU) ) 01051 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ ); 01052 if ( problem_->getLeftPrec() != Teuchos::null ) { 01053 index.resize(recycleBlocks_); 01054 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01055 index.resize(recycleBlocks_); 01056 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01057 Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage 01058 problem_->applyLeftPrec( *AUtmp, *LeftPCAU ); 01059 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp ); 01060 } else { 01061 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp ); 01062 } 01063 } 01064 01066 // RCG solver 01067 01068 Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter; 01069 rcg_iter = Teuchos::rcp( new RCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) ); 01070 01071 // Enter solve() iterations 01072 { 01073 Teuchos::TimeMonitor slvtimer(*timerSolve_); 01074 01075 while ( numRHS2Solve > 0 ) { 01076 01077 // Debugging output to tell use if recycle space exists and will be used 01078 if (printer_->isVerbosity( Debug ) ) { 01079 if (existU_) printer_->print( Debug, "Using recycle space generated from previous call to solve()." ); 01080 else printer_->print( Debug, "No recycle space exists." ); 01081 } 01082 01083 // Reset the number of iterations. 01084 rcg_iter->resetNumIters(); 01085 01086 // Set the current number of recycle blocks and subspace dimension with the RCG iteration. 01087 rcg_iter->setSize( recycleBlocks_, numBlocks_ ); 01088 01089 // Reset the number of calls that the status test output knows about. 01090 outputTest_->resetNumCalls(); 01091 01092 // indicate that updated recycle space has not yet been generated for this linear system 01093 existU1_ = false; 01094 01095 // reset cycle count 01096 cycle = 0; 01097 01098 // Get the current residual 01099 problem_->computeCurrResVec( &*r_ ); 01100 01101 // If U exists, find best soln over this space first 01102 if (existU_) { 01103 // Solve linear system UTAU * y = (U'*r) 01104 Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1); 01105 index.resize(recycleBlocks_); 01106 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01107 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index ); 01108 MVT::MvTransMv( one, *Utmp, *r_, Utr ); 01109 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ ); 01110 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ); 01111 LUUTAUtmp.assign(UTAUtmp); 01112 int info = 0; 01113 lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info); 01114 TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure, 01115 "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution."); 01116 01117 // Update solution (x = x + U*y) 01118 MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() ); 01119 01120 // Update residual ( r = r - AU*y ) 01121 index.resize(recycleBlocks_); 01122 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01123 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index ); 01124 MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ ); 01125 } 01126 01127 if ( problem_->getLeftPrec() != Teuchos::null ) { 01128 problem_->applyLeftPrec( *r_, *z_ ); 01129 } else { 01130 z_ = r_; 01131 } 01132 01133 // rTz_old = r'*z 01134 MVT::MvTransMv( one, *r_, *z_, *rTz_old_ ); 01135 01136 if ( existU_ ) { 01137 // mu = UTAU\(AU'*z); 01138 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1); 01139 index.resize(recycleBlocks_); 01140 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01141 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index ); 01142 MVT::MvTransMv( one, *AUtmp, *z_, mu ); 01143 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ); 01144 char TRANS = 'N'; 01145 int info; 01146 lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info ); 01147 TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure, 01148 "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution."); 01149 // p = z - U*mu; 01150 index.resize( 1 ); 01151 index[0] = 0; 01152 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index ); 01153 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp); 01154 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp ); 01155 } else { 01156 // p = z; 01157 index.resize( 1 ); 01158 index[0] = 0; 01159 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index ); 01160 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp); 01161 } 01162 01163 // Set the new state and initialize the solver. 01164 RCGIterState<ScalarType,MV> newstate; 01165 01166 // Create RCP views here 01167 index.resize( numBlocks_+1 ); 01168 for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; } 01169 newstate.P = MVT::CloneViewNonConst( *P_, index ); 01170 index.resize( recycleBlocks_ ); 01171 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01172 newstate.U = MVT::CloneViewNonConst( *U_, index ); 01173 index.resize( recycleBlocks_ ); 01174 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01175 newstate.AU = MVT::CloneViewNonConst( *AU_, index ); 01176 newstate.Alpha = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) ); 01177 newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) ); 01178 newstate.D = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) ); 01179 newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) ); 01180 newstate.LUUTAU = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) ); 01181 // assign the rest of the values in the struct 01182 newstate.curDim = 1; // We have initialized the first search vector 01183 newstate.Ap = Ap_; 01184 newstate.r = r_; 01185 newstate.z = z_; 01186 newstate.existU = existU_; 01187 newstate.ipiv = ipiv_; 01188 newstate.rTz_old = rTz_old_; 01189 01190 rcg_iter->initialize(newstate); 01191 01192 while(1) { 01193 01194 // tell rcg_iter to iterate 01195 try { 01196 rcg_iter->iterate(); 01197 01199 // 01200 // check convergence first 01201 // 01203 if ( convTest_->getStatus() == Passed ) { 01204 // We have convergence 01205 break; // break from while(1){rcg_iter->iterate()} 01206 } 01208 // 01209 // check for maximum iterations 01210 // 01212 else if ( maxIterTest_->getStatus() == Passed ) { 01213 // we don't have convergence 01214 isConverged = false; 01215 break; // break from while(1){rcg_iter->iterate()} 01216 } 01218 // 01219 // check if cycle complete; update for next cycle 01220 // 01222 else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) { 01223 // index into P_ of last search vector generated this cycle 01224 int lastp = -1; 01225 // index into Beta_ of last entry generated this cycle 01226 int lastBeta = -1; 01227 if (recycleBlocks_ > 0) { 01228 if (!existU_) { 01229 if (cycle == 0) { // No U, no U1 01230 01231 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ ); 01232 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ ); 01233 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 ); 01234 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 ); 01235 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 ); 01236 Ftmp.putScalar(zero); 01237 Gtmp.putScalar(zero); 01238 for (int ii=0;ii<numBlocks_;ii++) { 01239 Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0)); 01240 if (ii > 0) { 01241 Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01242 Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01243 } 01244 Ftmp(ii,ii) = Dtmp(ii,0); 01245 } 01246 01247 // compute harmonic Ritz vectors 01248 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ ); 01249 getHarmonicVecs(Ftmp,Gtmp,Ytmp); 01250 01251 // U1 = [P(:,1:end-1)*Y]; 01252 index.resize( numBlocks_ ); 01253 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; } 01254 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index ); 01255 index.resize( recycleBlocks_ ); 01256 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01257 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index ); 01258 MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp ); 01259 01260 // Precompute some variables for next cycle 01261 01262 // AU1TAU1 = Y'*G*Y; 01263 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ ); 01264 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero); 01265 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ ); 01266 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero); 01267 01268 // AU1TU1 = Y'*F*Y; 01269 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ ); 01270 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero); 01271 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ ); 01272 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero); 01273 01274 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 ); 01275 // Must reinitialize AU1TAP; can become dense later 01276 AU1TAPtmp.putScalar(zero); 01277 // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end)); 01278 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0); 01279 for (int ii=0; ii<recycleBlocks_; ++ii) { 01280 AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp; 01281 } 01282 01283 // indicate that updated recycle space now defined 01284 existU1_ = true; 01285 01286 // Indicate the size of the P, Beta structures generated this cycle 01287 lastp = numBlocks_; 01288 lastBeta = numBlocks_-1; 01289 01290 } // if (cycle == 0) 01291 else { // No U, but U1 guaranteed to exist now 01292 01293 // Finish computation of subblocks 01294 // AU1TAP = AU1TAP * D(1); 01295 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 ); 01296 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ ); 01297 AU1TAPtmp.scale(Dtmp(0,0)); 01298 01299 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 ); 01300 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 ); 01301 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ ); 01302 APTAPtmp.putScalar(zero); 01303 for (int ii=0; ii<numBlocks_; ii++) { 01304 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0)); 01305 if (ii > 0) { 01306 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01307 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01308 } 01309 } 01310 01311 // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)]; 01312 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01313 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ ); 01314 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01315 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01316 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01317 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ ); 01318 F11.assign(AU1TU1tmp); 01319 F12.putScalar(zero); 01320 F21.putScalar(zero); 01321 F22.putScalar(zero); 01322 for(int ii=0;ii<numBlocks_;ii++) { 01323 F22(ii,ii) = Dtmp(ii,0); 01324 } 01325 01326 // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP]; 01327 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01328 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ ); 01329 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ ); 01330 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01331 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01332 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01333 G11.assign(AU1TAU1tmp); 01334 G12.assign(AU1TAPtmp); 01335 // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually) 01336 for (int ii=0;ii<recycleBlocks_;++ii) 01337 for (int jj=0;jj<numBlocks_;++jj) 01338 G21(jj,ii) = G12(ii,jj); 01339 G22.assign(APTAPtmp); 01340 01341 // compute harmonic Ritz vectors 01342 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ ); 01343 getHarmonicVecs(Ftmp,Gtmp,Ytmp); 01344 01345 // U1 = [U1 P(:,2:end-1)]*Y; 01346 index.resize( numBlocks_ ); 01347 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; } 01348 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index ); 01349 index.resize( recycleBlocks_ ); 01350 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01351 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index ); 01352 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01353 index.resize( recycleBlocks_ ); 01354 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01355 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index ); 01356 index.resize( recycleBlocks_ ); 01357 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01358 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index ); 01359 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ ); 01360 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp ); 01361 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp ); 01362 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp); 01363 01364 // Precompute some variables for next cycle 01365 01366 // AU1TAU1 = Y'*G*Y; 01367 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01368 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero); 01369 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero); 01370 01371 // AU1TAP = zeros(k,m); 01372 // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end)); 01373 AU1TAPtmp.putScalar(zero); 01374 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0); 01375 for (int ii=0; ii<recycleBlocks_; ++ii) { 01376 AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp; 01377 } 01378 01379 // AU1TU1 = Y'*F*Y; 01380 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01381 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero); 01382 //Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ ); 01383 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero); 01384 01385 // Indicate the size of the P, Beta structures generated this cycle 01386 lastp = numBlocks_+1; 01387 lastBeta = numBlocks_; 01388 01389 } // if (cycle != 1) 01390 } // if (!existU_) 01391 else { // U exists 01392 if (cycle == 0) { // No U1, but U exists 01393 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 ); 01394 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 ); 01395 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 ); 01396 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ ); 01397 APTAPtmp.putScalar(zero); 01398 for (int ii=0; ii<numBlocks_; ii++) { 01399 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0)); 01400 if (ii > 0) { 01401 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01402 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01403 } 01404 } 01405 01406 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ ); 01407 L2tmp.putScalar(zero); 01408 for(int ii=0;ii<numBlocks_;ii++) { 01409 L2tmp(ii,ii) = 1./Alphatmp(ii,0); 01410 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0); 01411 } 01412 01413 // AUTAP = UTAU*Delta*L2; 01414 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ ); 01415 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ ); 01416 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 ); 01417 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ ); 01418 DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero); 01419 AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero); 01420 01421 // F = [UTAU zeros(k,m); zeros(m,k) diag(D)]; 01422 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01423 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ ); 01424 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01425 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01426 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01427 F11.assign(UTAUtmp); 01428 F12.putScalar(zero); 01429 F21.putScalar(zero); 01430 F22.putScalar(zero); 01431 for(int ii=0;ii<numBlocks_;ii++) { 01432 F22(ii,ii) = Dtmp(ii,0); 01433 } 01434 01435 // G = [AUTAU AUTAP; AUTAP' APTAP]; 01436 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01437 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ ); 01438 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01439 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01440 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01441 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ ); 01442 G11.assign(AUTAUtmp); 01443 G12.assign(AUTAPtmp); 01444 // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually) 01445 for (int ii=0;ii<recycleBlocks_;++ii) 01446 for (int jj=0;jj<numBlocks_;++jj) 01447 G21(jj,ii) = G12(ii,jj); 01448 G22.assign(APTAPtmp); 01449 01450 // compute harmonic Ritz vectors 01451 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ ); 01452 getHarmonicVecs(Ftmp,Gtmp,Ytmp); 01453 01454 // U1 = [U P(:,1:end-1)]*Y; 01455 index.resize( recycleBlocks_ ); 01456 for (int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; } 01457 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index ); 01458 index.resize( numBlocks_ ); 01459 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; } 01460 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index ); 01461 index.resize( recycleBlocks_ ); 01462 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01463 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index ); 01464 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01465 index.resize( recycleBlocks_ ); 01466 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01467 Teuchos::RCP<MV> UY1tmp = MVT::CloneViewNonConst( *U1Y1_, index ); 01468 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ ); 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 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp ); 01473 MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp ); 01474 MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp); 01475 01476 // Precompute some variables for next cycle 01477 01478 // AU1TAU1 = Y'*G*Y; 01479 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01480 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero); 01481 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ ); 01482 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero); 01483 01484 // AU1TU1 = Y'*F*Y; 01485 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01486 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero); 01487 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ ); 01488 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero); 01489 01490 // AU1TU = UTAU; 01491 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ ); 01492 AU1TUtmp.assign(UTAUtmp); 01493 01494 // dold = D(end); 01495 dold = Dtmp(numBlocks_-1,0); 01496 01497 // indicate that updated recycle space now defined 01498 existU1_ = true; 01499 01500 // Indicate the size of the P, Beta structures generated this cycle 01501 lastp = numBlocks_; 01502 lastBeta = numBlocks_-1; 01503 } 01504 else { // Have U and U1 01505 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 ); 01506 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 ); 01507 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 ); 01508 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ ); 01509 for (int ii=0; ii<numBlocks_; ii++) { 01510 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0)); 01511 if (ii > 0) { 01512 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01513 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0); 01514 } 01515 } 01516 01517 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ ); 01518 for(int ii=0;ii<numBlocks_;ii++) { 01519 L2tmp(ii,ii) = 1./Alphatmp(ii,0); 01520 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0); 01521 } 01522 01523 // M(end,1) = dold*(-Beta(1)/Alpha(1)); 01524 // AU1TAP = Y'*[AU1TU*Delta*L2; M]; 01525 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ ); 01526 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 ); 01527 DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero); 01528 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ ); 01529 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ ); 01530 AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero); 01531 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ ); 01532 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ ); 01533 AU1TAPtmp.putScalar(zero); 01534 AU1TAPtmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero); 01535 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01536 ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0)); 01537 for(int ii=0;ii<recycleBlocks_;ii++) { 01538 AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val; 01539 } 01540 01541 // AU1TU = Y1'*AU1TU 01542 Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ ); 01543 Y1TAU1TU.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero); 01544 AU1TUtmp.assign(Y1TAU1TU); 01545 01546 // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)]; 01547 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01548 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ ); 01549 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01550 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01551 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01552 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ ); 01553 F11.assign(AU1TU1tmp); 01554 for(int ii=0;ii<numBlocks_;ii++) { 01555 F22(ii,ii) = Dtmp(ii,0); 01556 } 01557 01558 // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP]; 01559 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) ); 01560 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ ); 01561 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ ); 01562 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 ); 01563 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ ); 01564 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ ); 01565 G11.assign(AU1TAU1tmp); 01566 G12.assign(AU1TAPtmp); 01567 // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually) 01568 for (int ii=0;ii<recycleBlocks_;++ii) 01569 for (int jj=0;jj<numBlocks_;++jj) 01570 G21(jj,ii) = G12(ii,jj); 01571 G22.assign(APTAPtmp); 01572 01573 // compute harmonic Ritz vectors 01574 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ ); 01575 getHarmonicVecs(Ftmp,Gtmp,Ytmp); 01576 01577 // U1 = [U1 P(:,2:end-1)]*Y; 01578 index.resize( numBlocks_ ); 01579 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; } 01580 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index ); 01581 index.resize( recycleBlocks_ ); 01582 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01583 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index ); 01584 index.resize( recycleBlocks_ ); 01585 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01586 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index ); 01587 index.resize( recycleBlocks_ ); 01588 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01589 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index ); 01590 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp ); 01591 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp ); 01592 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp); 01593 01594 // Precompute some variables for next cycle 01595 01596 // AU1TAU1 = Y'*G*Y; 01597 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01598 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero); 01599 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero); 01600 01601 // AU1TU1 = Y'*F*Y; 01602 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ ); 01603 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero); 01604 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero); 01605 01606 // dold = D(end); 01607 dold = Dtmp(numBlocks_-1,0); 01608 01609 // Indicate the size of the P, Beta structures generated this cycle 01610 lastp = numBlocks_+1; 01611 lastBeta = numBlocks_; 01612 01613 } 01614 } 01615 } // if (recycleBlocks_ > 0) 01616 01617 // Cleanup after end of cycle 01618 01619 // P = P(:,end-1:end); 01620 index.resize( 2 ); 01621 index[0] = lastp-1; index[1] = lastp; 01622 Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_, index ); 01623 index[0] = 0; index[1] = 1; 01624 MVT::SetBlock(*Ptmp2,index,*P_); 01625 01626 // Beta = Beta(end); 01627 (*Beta_)(0,0) = (*Beta_)(lastBeta,0); 01628 01629 // Delta = Delta(:,end); 01630 if (existU_) { // Delta only initialized if U exists 01631 Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 ); 01632 Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ ); 01633 mu1.assign(mu2); 01634 } 01635 01636 // Now reinitialize state variables for next cycle 01637 newstate.P = Teuchos::null; 01638 index.resize( numBlocks_+1 ); 01639 for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; } 01640 newstate.P = MVT::CloneViewNonConst( *P_, index ); 01641 01642 newstate.Beta = Teuchos::null; 01643 newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) ); 01644 01645 newstate.Delta = Teuchos::null; 01646 newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) ); 01647 01648 newstate.curDim = 1; // We have initialized the first search vector 01649 01650 // Pass to iteration object 01651 rcg_iter->initialize(newstate); 01652 01653 // increment cycle count 01654 cycle = cycle + 1; 01655 01656 } 01658 // 01659 // we returned from iterate(), but none of our status tests Passed. 01660 // something is wrong, and it is probably our fault. 01661 // 01663 else { 01664 TEST_FOR_EXCEPTION(true,std::logic_error, 01665 "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate()."); 01666 } 01667 } 01668 catch (const std::exception &e) { 01669 printer_->stream(Errors) << "Error! Caught std::exception in RCGIter::iterate() at iteration " 01670 << rcg_iter->getNumIters() << std::endl 01671 << e.what() << std::endl; 01672 throw; 01673 } 01674 } 01675 01676 // Inform the linear problem that we are finished with this block linear system. 01677 problem_->setCurrLS(); 01678 01679 // Update indices for the linear systems to be solved. 01680 numRHS2Solve -= 1; 01681 if ( numRHS2Solve > 0 ) { 01682 currIdx[0]++; 01683 // Set the next indices. 01684 problem_->setLSIndex( currIdx ); 01685 } 01686 else { 01687 currIdx.resize( numRHS2Solve ); 01688 } 01689 01690 // Update the recycle space for the next linear system 01691 if (existU1_) { // be sure updated recycle space was created 01692 // U = U1 01693 index.resize(recycleBlocks_); 01694 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01695 MVT::SetBlock(*U1_,index,*U_); 01696 // Set flag indicating recycle space is now defined 01697 existU_ = true; 01698 if (numRHS2Solve > 0) { // also update AU, UTAU, and AUTAU 01699 // Free pointers in newstate 01700 newstate.P = Teuchos::null; 01701 newstate.Ap = Teuchos::null; 01702 newstate.r = Teuchos::null; 01703 newstate.z = Teuchos::null; 01704 newstate.U = Teuchos::null; 01705 newstate.AU = Teuchos::null; 01706 newstate.Alpha = Teuchos::null; 01707 newstate.Beta = Teuchos::null; 01708 newstate.D = Teuchos::null; 01709 newstate.Delta = Teuchos::null; 01710 newstate.LUUTAU = Teuchos::null; 01711 newstate.ipiv = Teuchos::null; 01712 newstate.rTz_old = Teuchos::null; 01713 01714 // Reinitialize AU, UTAU, AUTAU 01715 index.resize(recycleBlocks_); 01716 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01717 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index ); 01718 index.resize(recycleBlocks_); 01719 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01720 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index ); 01721 // Initialize AU 01722 problem_->applyOp( *Utmp, *AUtmp ); 01723 // Initialize UTAU 01724 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ ); 01725 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp ); 01726 // Initialize AUTAU ( AUTAU = AU'*(M\AU) ) 01727 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ ); 01728 if ( problem_->getLeftPrec() != Teuchos::null ) { 01729 index.resize(recycleBlocks_); 01730 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; } 01731 index.resize(recycleBlocks_); 01732 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; } 01733 Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage 01734 problem_->applyLeftPrec( *AUtmp, *LeftPCAU ); 01735 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp ); 01736 } else { 01737 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp ); 01738 } 01739 } // if (numRHS2Solve > 0) 01740 01741 } // if (existU1) 01742 }// while ( numRHS2Solve > 0 ) 01743 01744 } 01745 01746 // print final summary 01747 sTest_->print( printer_->stream(FinalSummary) ); 01748 01749 // print timing information 01750 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) ); 01751 01752 // get iteration information for this solve 01753 numIters_ = maxIterTest_->getNumIters(); 01754 01755 if (!isConverged) { 01756 return Unconverged; // return from RCGSolMgr::solve() 01757 } 01758 return Converged; // return from RCGSolMgr::solve() 01759 } 01760 01761 // Compute the harmonic eigenpairs of the projected, dense system. 01762 template<class ScalarType, class MV, class OP> 01763 void RCGSolMgr<ScalarType,MV,OP>::getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType>& F, 01764 const Teuchos::SerialDenseMatrix<int,ScalarType>& G, 01765 Teuchos::SerialDenseMatrix<int,ScalarType>& Y) { 01766 // order of F,G 01767 int n = F.numCols(); 01768 01769 // The LAPACK interface 01770 Teuchos::LAPACK<int,ScalarType> lapack; 01771 01772 // Magnitude of harmonic Ritz values 01773 std::vector<MagnitudeType> w(n); 01774 01775 // Sorted order of harmonic Ritz values 01776 std::vector<int> iperm(n); 01777 01778 // Compute k smallest harmonic Ritz pairs 01779 // SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, INFO ) 01780 int itype = 1; // solve A*x = (lambda)*B*x 01781 char jobz='V'; // compute eigenvalues and eigenvectors 01782 char uplo='U'; // since F,G symmetric, reference only their upper triangular data 01783 std::vector<ScalarType> work(1); 01784 int lwork = -1; 01785 int info = 0; 01786 // since SYGV destroys workspace, create copies of F,G 01787 Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_, F_->numRows(), F_->numCols() ); 01788 Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *G_, G_->numRows(), G_->numCols() ); 01789 01790 // query for optimal workspace size 01791 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info); 01792 TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure, 01793 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size."); 01794 lwork = (int)work[0]; 01795 work.resize(lwork); 01796 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info); 01797 TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure, 01798 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions."); 01799 01800 01801 // Construct magnitude of each harmonic Ritz value 01802 this->sort(w,n,iperm); 01803 01804 // Select recycledBlocks_ smallest eigenvectors 01805 for( int i=0; i<recycleBlocks_; i++ ) { 01806 for( int j=0; j<n; j++ ) { 01807 Y(j,i) = G2(j,iperm[i]); 01808 } 01809 } 01810 01811 } 01812 01813 // This method sorts list of n floating-point numbers and return permutation vector 01814 template<class ScalarType, class MV, class OP> 01815 void RCGSolMgr<ScalarType,MV,OP>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm) 01816 { 01817 int l, r, j, i, flag; 01818 int RR2; 01819 double dRR, dK; 01820 01821 // Initialize the permutation vector. 01822 for(j=0;j<n;j++) 01823 iperm[j] = j; 01824 01825 if (n <= 1) return; 01826 01827 l = n / 2 + 1; 01828 r = n - 1; 01829 l = l - 1; 01830 dRR = dlist[l - 1]; 01831 dK = dlist[l - 1]; 01832 01833 RR2 = iperm[l - 1]; 01834 while (r != 0) { 01835 j = l; 01836 flag = 1; 01837 01838 while (flag == 1) { 01839 i = j; 01840 j = j + j; 01841 01842 if (j > r + 1) 01843 flag = 0; 01844 else { 01845 if (j < r + 1) 01846 if (dlist[j] > dlist[j - 1]) j = j + 1; 01847 01848 if (dlist[j - 1] > dK) { 01849 dlist[i - 1] = dlist[j - 1]; 01850 iperm[i - 1] = iperm[j - 1]; 01851 } 01852 else { 01853 flag = 0; 01854 } 01855 } 01856 } 01857 dlist[i - 1] = dRR; 01858 iperm[i - 1] = RR2; 01859 if (l == 1) { 01860 dRR = dlist [r]; 01861 RR2 = iperm[r]; 01862 dK = dlist[r]; 01863 dlist[r] = dlist[0]; 01864 iperm[r] = iperm[0]; 01865 r = r - 1; 01866 } 01867 else { 01868 l = l - 1; 01869 dRR = dlist[l - 1]; 01870 RR2 = iperm[l - 1]; 01871 dK = dlist[l - 1]; 01872 } 01873 } 01874 dlist[0] = dRR; 01875 iperm[0] = RR2; 01876 } 01877 01878 // This method requires the solver manager to return a std::string that describes itself. 01879 template<class ScalarType, class MV, class OP> 01880 std::string RCGSolMgr<ScalarType,MV,OP>::description() const 01881 { 01882 std::ostringstream oss; 01883 oss << "Belos::RCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">"; 01884 return oss.str(); 01885 } 01886 01887 } // end Belos namespace 01888 01889 #endif /* BELOS_RCG_SOLMGR_HPP */
1.7.4
