Belos Version of the Day
BelosPCPGSolMgr.hpp
Go to the documentation of this file.
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_PCPG_SOLMGR_HPP
00043 #define BELOS_PCPG_SOLMGR_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 
00052 #include "BelosLinearProblem.hpp"
00053 #include "BelosSolverManager.hpp"
00054 
00055 #include "BelosPCPGIter.hpp"
00056 
00057 #include "BelosDGKSOrthoManager.hpp"
00058 #include "BelosICGSOrthoManager.hpp"
00059 #include "BelosIMGSOrthoManager.hpp"
00060 #include "BelosStatusTestMaxIters.hpp"
00061 #include "BelosStatusTestGenResNorm.hpp"
00062 #include "BelosStatusTestCombo.hpp"
00063 #include "BelosStatusTestOutputFactory.hpp"
00064 #include "BelosOutputManager.hpp"
00065 #include "Teuchos_BLAS.hpp"
00066 #include "Teuchos_LAPACK.hpp"
00067 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00068 #include "Teuchos_TimeMonitor.hpp"
00069 #endif
00070 
00071 //#include <vector>              //getPermThatSorts
00072 //#include <algorithm>           //getPermThatSorts
00073 
00074 
00091 namespace Belos {
00092   
00094 
00095   
00102   class PCPGSolMgrLinearProblemFailure : public BelosError {public:
00103     PCPGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00104     {}};
00105   
00111   class PCPGSolMgrOrthoFailure : public BelosError {public:
00112     PCPGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00113     {}};
00114   
00121   class PCPGSolMgrLAPACKFailure : public BelosError {public:
00122     PCPGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
00123     {}};
00124 
00131   class PCPGSolMgrRecyclingFailure : public BelosError {public:
00132     PCPGSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg)
00133     {}};
00134 
00136 
00137 
00138   template<class ScalarType, class MV, class OP>
00139   class PCPGSolMgr : public SolverManager<ScalarType,MV,OP> {
00140     
00141   private:
00142     typedef MultiVecTraits<ScalarType,MV> MVT;
00143     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00144     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00145     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00146     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00147     
00148   public:
00149     
00151 
00152    
00159      PCPGSolMgr();
00160  
00196     PCPGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00197           const Teuchos::RCP<Teuchos::ParameterList> &pl );
00198     
00200     virtual ~PCPGSolMgr() {};
00202     
00204 
00205     
00208     const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00209       return *problem_;
00210     }
00211 
00214     Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00215 
00218     Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00219  
00225     Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00226       return Teuchos::tuple(timerSolve_);
00227     }
00228 
00234     MagnitudeType achievedTol() const {
00235       return achievedTol_;
00236     }
00237 
00239     int getNumIters() const {
00240       return numIters_;
00241     }
00242  
00245     bool isLOADetected() const { return false; }
00246  
00248     
00250 
00251    
00253     void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00254    
00256     void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
00257     
00259    
00261 
00262 
00266     void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00268  
00270 
00271     
00289     ReturnType solve();
00290     
00292     
00295     
00297     std::string description() const;
00298     
00300     
00301   private:
00302 
00303     // In the A-inner product, perform an RRQR decomposition without using A unless absolutely necessary.  Given 
00304     // the seed space U and C = A U, find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU.
00305     int ARRQR(int numVecs, int numOrthVecs, const Teuchos::SerialDenseMatrix<int,ScalarType>& D);
00306 
00307     // Linear problem.
00308     Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00309     
00310     // Output manager.
00311     Teuchos::RCP<OutputManager<ScalarType> > printer_;
00312     Teuchos::RCP<std::ostream> outputStream_;
00313 
00314     // Status test.
00315     Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00316     Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00317     Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
00318     Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00319 
00320     // Orthogonalization manager.
00321     Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_; 
00322     
00323     // Current parameter list.
00324     Teuchos::RCP<Teuchos::ParameterList> params_;
00325 
00326     // Default solver values.
00327     static const MagnitudeType convtol_default_;
00328     static const MagnitudeType orthoKappa_default_;
00329     static const int maxIters_default_;
00330     static const int deflatedBlocks_default_;
00331     static const int savedBlocks_default_;
00332     static const int verbosity_default_;
00333     static const int outputStyle_default_;
00334     static const int outputFreq_default_;
00335     static const std::string label_default_;
00336     static const std::string orthoType_default_;
00337     static const Teuchos::RCP<std::ostream> outputStream_default_;
00338 
00339     //
00340     // Current solver values.
00341     //
00342 
00344     MagnitudeType convtol_;
00345 
00347     MagnitudeType orthoKappa_;
00348 
00350     MagnitudeType achievedTol_;
00351 
00353     int numIters_;
00354 
00356     int maxIters_;
00357 
00358     int deflatedBlocks_, savedBlocks_, verbosity_, outputStyle_, outputFreq_;
00359     std::string orthoType_; 
00360 
00361     // Recycled subspace, its image and the residual
00362     Teuchos::RCP<MV> U_, C_, R_;
00363 
00364     // Actual dimension of current recycling subspace (<= savedBlocks_ )
00365     int dimU_; 
00366 
00367     // Timers.
00368     std::string label_;
00369     Teuchos::RCP<Teuchos::Time> timerSolve_;
00370 
00371     // Internal state variables.
00372     bool isSet_;
00373   };
00374 
00375 
00376 // Default solver values.
00377 template<class ScalarType, class MV, class OP>
00378 const typename PCPGSolMgr<ScalarType,MV,OP>::MagnitudeType PCPGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00379 
00380 template<class ScalarType, class MV, class OP>
00381 const typename PCPGSolMgr<ScalarType,MV,OP>::MagnitudeType PCPGSolMgr<ScalarType,MV,OP>::orthoKappa_default_ = -1.0;
00382 
00383 template<class ScalarType, class MV, class OP>
00384 const int PCPGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00385 
00386 template<class ScalarType, class MV, class OP>
00387 const int PCPGSolMgr<ScalarType,MV,OP>::deflatedBlocks_default_ = 2;
00388 
00389 template<class ScalarType, class MV, class OP>
00390 const int PCPGSolMgr<ScalarType,MV,OP>::savedBlocks_default_ = 16;
00391 
00392 template<class ScalarType, class MV, class OP>
00393 const int PCPGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00394 
00395 template<class ScalarType, class MV, class OP>
00396 const int PCPGSolMgr<ScalarType,MV,OP>::outputStyle_default_ = Belos::General;
00397 
00398 template<class ScalarType, class MV, class OP>
00399 const int PCPGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00400 
00401 template<class ScalarType, class MV, class OP>
00402 const std::string PCPGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00403 
00404 template<class ScalarType, class MV, class OP>
00405 const std::string PCPGSolMgr<ScalarType,MV,OP>::orthoType_default_ = "DGKS";
00406 
00407 template<class ScalarType, class MV, class OP>
00408 const Teuchos::RCP<std::ostream> PCPGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00409 
00410 
00411 // Empty Constructor
00412 template<class ScalarType, class MV, class OP>
00413 PCPGSolMgr<ScalarType,MV,OP>::PCPGSolMgr() :
00414   outputStream_(outputStream_default_),
00415   convtol_(convtol_default_),
00416   orthoKappa_(orthoKappa_default_),
00417   achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
00418   numIters_(0),
00419   maxIters_(maxIters_default_),
00420   deflatedBlocks_(deflatedBlocks_default_),
00421   savedBlocks_(savedBlocks_default_),
00422   verbosity_(verbosity_default_),
00423   outputStyle_(outputStyle_default_),
00424   outputFreq_(outputFreq_default_),
00425   orthoType_(orthoType_default_),
00426   label_(label_default_),
00427   isSet_(false)
00428 {}
00429 
00430 
00431 // Basic Constructor
00432 template<class ScalarType, class MV, class OP>
00433 PCPGSolMgr<ScalarType,MV,OP>::PCPGSolMgr( 
00434                const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00435                const Teuchos::RCP<Teuchos::ParameterList> &pl ) : 
00436   problem_(problem),
00437   outputStream_(outputStream_default_),
00438   convtol_(convtol_default_),
00439   orthoKappa_(orthoKappa_default_),
00440   achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
00441   numIters_(0),
00442   maxIters_(maxIters_default_),
00443   deflatedBlocks_(deflatedBlocks_default_),
00444   savedBlocks_(savedBlocks_default_),
00445   verbosity_(verbosity_default_),
00446   outputStyle_(outputStyle_default_),
00447   outputFreq_(outputFreq_default_),
00448   orthoType_(orthoType_default_),
00449   label_(label_default_),
00450   isSet_(false)
00451 {
00452   TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00453 
00454   if (!is_null(pl)) {
00455     // Set the parameters using the list that was passed in.
00456     setParameters( pl );  
00457   }
00458 }
00459 
00460 
00461 template<class ScalarType, class MV, class OP>
00462 void PCPGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
00463 {
00464   // Create the internal parameter list if ones doesn't already exist.
00465   if (params_ == Teuchos::null) {
00466     params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00467   }
00468   else {
00469     params->validateParameters(*getValidParameters());
00470   }
00471 
00472   // Check for maximum number of iterations
00473   if (params->isParameter("Maximum Iterations")) {
00474     maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00475 
00476     // Update parameter in our list and in status test.
00477     params_->set("Maximum Iterations", maxIters_);
00478     if (maxIterTest_!=Teuchos::null)
00479       maxIterTest_->setMaxIters( maxIters_ );
00480   }
00481 
00482   // Check for the maximum numbers of saved and deflated blocks.
00483   if (params->isParameter("Num Saved Blocks")) {
00484     savedBlocks_ = params->get("Num Saved Blocks",savedBlocks_default_);
00485     TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ <= 0, std::invalid_argument,
00486            "Belos::PCPGSolMgr: \"Num Saved Blocks\" must be strictly positive.");
00487 
00488     // savedBlocks > number of matrix rows and columns, not known in parameters.
00489     //TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks_ >= maxIters_, std::invalid_argument,
00490     //"Belos::PCPGSolMgr: \"Num Saved Blocks\" must be less than \"Maximum Iterations\".");
00491 
00492     // Update parameter in our list.
00493     params_->set("Num Saved Blocks", savedBlocks_);
00494   }
00495   if (params->isParameter("Num Deflated Blocks")) {
00496     deflatedBlocks_ = params->get("Num Deflated Blocks",deflatedBlocks_default_);
00497     TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ < 0, std::invalid_argument,
00498            "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be positive.");
00499 
00500     TEUCHOS_TEST_FOR_EXCEPTION(deflatedBlocks_ > savedBlocks_, std::invalid_argument,
00501            "Belos::PCPGSolMgr: \"Num Deflated Blocks\" must be <= \"Num Saved Blocks\".");
00502 
00503     // Update parameter in our list.
00504     params_->set("Num Deflated Blocks", deflatedBlocks_);
00505   }
00506 
00507   // Check to see if the timer label changed.
00508   if (params->isParameter("Timer Label")) {
00509     std::string tempLabel = params->get("Timer Label", label_default_);
00510 
00511     // Update parameter in our list and solver timer
00512     if (tempLabel != label_) {
00513       label_ = tempLabel;
00514       params_->set("Timer Label", label_);
00515       std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
00516 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00517       timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00518 #endif
00519     }
00520   }
00521 
00522   // Check if the orthogonalization changed.
00523   if (params->isParameter("Orthogonalization")) {
00524     std::string tempOrthoType = params->get("Orthogonalization",orthoType_default_);
00525     TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType != "DGKS" && tempOrthoType != "ICGS" && tempOrthoType != "IMGS", 
00526       std::invalid_argument,
00527       "Belos::PCPGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
00528     if (tempOrthoType != orthoType_) {
00529       orthoType_ = tempOrthoType;
00530       // Create orthogonalization manager
00531       if (orthoType_=="DGKS") {
00532   if (orthoKappa_ <= 0) {
00533     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00534   }
00535   else {
00536     ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00537     Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00538   }
00539       }
00540       else if (orthoType_=="ICGS") {
00541   ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00542       } 
00543       else if (orthoType_=="IMGS") {
00544   ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00545       } 
00546     }  
00547   }
00548 
00549   // Check which orthogonalization constant to use.
00550   if (params->isParameter("Orthogonalization Constant")) {
00551     orthoKappa_ = params->get("Orthogonalization Constant",orthoKappa_default_);
00552 
00553     // Update parameter in our list.
00554     params_->set("Orthogonalization Constant",orthoKappa_);
00555     if (orthoType_=="DGKS") {
00556       if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
00557   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00558       }
00559     } 
00560   }
00561 
00562   // Check for a change in verbosity level
00563   if (params->isParameter("Verbosity")) {
00564     if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00565       verbosity_ = params->get("Verbosity", verbosity_default_);
00566     } else {
00567       verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00568     }
00569 
00570     // Update parameter in our list.
00571     params_->set("Verbosity", verbosity_);
00572     if (printer_ != Teuchos::null)
00573       printer_->setVerbosity(verbosity_);
00574   }
00575 
00576   // Check for a change in output style
00577   if (params->isParameter("Output Style")) {
00578     if (Teuchos::isParameterType<int>(*params,"Output Style")) {
00579       outputStyle_ = params->get("Output Style", outputStyle_default_);
00580     } else {
00581       outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
00582     }
00583 
00584     // Reconstruct the convergence test if the explicit residual test is not being used.
00585     params_->set("Output Style", outputStyle_);
00586     outputTest_ = Teuchos::null;
00587   }
00588 
00589   // output stream
00590   if (params->isParameter("Output Stream")) {
00591     outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00592 
00593     // Update parameter in our list.
00594     params_->set("Output Stream", outputStream_);
00595     if (printer_ != Teuchos::null)
00596       printer_->setOStream( outputStream_ );
00597   }
00598 
00599   // frequency level
00600   if (verbosity_ & Belos::StatusTestDetails) {
00601     if (params->isParameter("Output Frequency")) {
00602       outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00603     }
00604 
00605     // Update parameter in out list and output status test.
00606     params_->set("Output Frequency", outputFreq_);
00607     if (outputTest_ != Teuchos::null)
00608       outputTest_->setOutputFrequency( outputFreq_ );
00609   }
00610 
00611   // Create output manager if we need to.
00612   if (printer_ == Teuchos::null) {
00613     printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00614   }  
00615   
00616   // Convergence
00617   typedef Belos::StatusTestCombo<ScalarType,MV,OP>  StatusTestCombo_t;
00618   typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP>  StatusTestResNorm_t;
00619 
00620   // Check for convergence tolerance
00621   if (params->isParameter("Convergence Tolerance")) {
00622     convtol_ = params->get("Convergence Tolerance",convtol_default_);
00623 
00624     // Update parameter in our list and residual tests.
00625     params_->set("Convergence Tolerance", convtol_);
00626     if (convTest_ != Teuchos::null)
00627       convTest_->setTolerance( convtol_ );
00628   }
00629 
00630   // Create status tests if we need to.
00631 
00632   // Basic test checks maximum iterations and native residual.
00633   if (maxIterTest_ == Teuchos::null)
00634     maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00635 
00636   if (convTest_ == Teuchos::null)
00637     convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
00638 
00639   sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00640   
00641   // Create the status test output class.
00642   // This class manages and formats the output from the status test.
00643   StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
00644   outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
00645 
00646   // Set the solver string for the output test
00647   std::string solverDesc = " PCPG ";
00648   outputTest_->setSolverDesc( solverDesc );
00649 
00650 
00651   // Create orthogonalization manager if we need to.
00652   if (ortho_ == Teuchos::null) {
00653     if (orthoType_=="DGKS") {
00654       if (orthoKappa_ <= 0) {
00655   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00656       }
00657       else {
00658   ortho_ = Teuchos::rcp( new DGKSOrthoManager<ScalarType,MV,OP>( label_ ) );
00659   Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
00660       }
00661     }
00662     else if (orthoType_=="ICGS") {
00663       ortho_ = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00664     } 
00665     else if (orthoType_=="IMGS") {
00666       ortho_ = Teuchos::rcp( new IMGSOrthoManager<ScalarType,MV,OP>( label_ ) );
00667     } 
00668     else {
00669       TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!="ICGS"&&orthoType_!="DGKS"&&orthoType_!="IMGS",std::logic_error,
00670        "Belos::PCPGSolMgr(): Invalid orthogonalization type.");
00671     }  
00672   }
00673 
00674   // Create the timer if we need to.
00675   if (timerSolve_ == Teuchos::null) {
00676     std::string solveLabel = label_ + ": PCPGSolMgr total solve time";
00677 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00678     timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00679 #endif
00680   }
00681 
00682   // Inform the solver manager that the current parameters were set.
00683   isSet_ = true;
00684 }
00685 
00686     
00687 template<class ScalarType, class MV, class OP>
00688 Teuchos::RCP<const Teuchos::ParameterList>
00689 PCPGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00690 {
00691   static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00692   if (is_null(validPL)) {
00693     Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00694     // Set all the valid parameters and their default values.
00695     pl->set("Convergence Tolerance", convtol_default_,
00696       "The relative residual tolerance that needs to be achieved by the\n"
00697       "iterative solver in order for the linear system to be declared converged.");
00698     pl->set("Maximum Iterations", maxIters_default_,
00699       "The maximum number of iterations allowed for each\n"
00700       "set of RHS solved.");
00701     pl->set("Num Deflated Blocks", deflatedBlocks_default_,
00702       "The maximum number of vectors in the seed subspace." );
00703     pl->set("Num Saved Blocks", savedBlocks_default_,
00704       "The maximum number of vectors saved from old Krylov subspaces." );
00705     pl->set("Verbosity", verbosity_default_,
00706       "What type(s) of solver information should be outputted\n"
00707       "to the output stream.");
00708     pl->set("Output Style", outputStyle_default_,
00709       "What style is used for the solver information outputted\n"
00710       "to the output stream.");
00711     pl->set("Output Frequency", outputFreq_default_,
00712       "How often convergence information should be outputted\n"
00713       "to the output stream.");  
00714     pl->set("Output Stream", outputStream_default_,
00715       "A reference-counted pointer to the output stream where all\n"
00716       "solver output is sent.");
00717     pl->set("Timer Label", label_default_,
00718       "The string to use as a prefix for the timer labels.");
00719     //  pl->set("Restart Timers", restartTimers_);
00720     pl->set("Orthogonalization", orthoType_default_,
00721       "The type of orthogonalization to use: DGKS, ICGS, IMGS");
00722     pl->set("Orthogonalization Constant",orthoKappa_default_,
00723       "The constant used by DGKS orthogonalization to determine\n"
00724       "whether another step of classical Gram-Schmidt is necessary.");
00725     validPL = pl;
00726   }
00727   return validPL;
00728 }
00729 
00730   
00731 // solve()
00732 template<class ScalarType, class MV, class OP>
00733 ReturnType PCPGSolMgr<ScalarType,MV,OP>::solve() {
00734 
00735   // Set the current parameters if are not set already.
00736   if (!isSet_) { setParameters( params_ ); }
00737 
00738   Teuchos::BLAS<int,ScalarType> blas;
00739   Teuchos::LAPACK<int,ScalarType> lapack;
00740   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00741   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00742   
00743   TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,PCPGSolMgrLinearProblemFailure,
00744                      "Belos::PCPGSolMgr::solve(): Linear problem is not a valid object.");
00745 
00746   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PCPGSolMgrLinearProblemFailure,
00747                      "Belos::PCPGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00748 
00749   // Create indices for the linear systems to be solved.
00750   int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00751   std::vector<int> currIdx(1);
00752   currIdx[0] = 0;
00753 
00754    bool debug = false;
00755 
00756   // Inform the linear problem of the current linear system to solve.
00757   problem_->setLSIndex( currIdx ); // block size == 1
00758 
00759   // Assume convergence is achieved, then let any failed convergence set this to false.
00760   bool isConverged = true;  
00761 
00763   // PCPG iteration parameter list
00764   Teuchos::ParameterList plist;
00765   plist.set("Saved Blocks", savedBlocks_);
00766   plist.set("Block Size", 1);
00767   plist.set("Keep Diagonal", true);
00768   plist.set("Initialize Diagonal", true);
00769   
00771   // PCPG solver
00772   
00773   Teuchos::RCP<PCPGIter<ScalarType,MV,OP> > pcpg_iter;
00774   pcpg_iter = Teuchos::rcp( new PCPGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
00775   // Number of iterations required to generate initial recycle space (if needed)
00776 
00777   // Enter solve() iterations
00778   {
00779 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00780     Teuchos::TimeMonitor slvtimer(*timerSolve_);
00781 #endif
00782     while ( numRHS2Solve > 0 ) {  // test for quick return
00783       
00784       // Reset the status test.
00785       outputTest_->reset();
00786 
00787       // Create the first block in the current Krylov basis (residual).
00788       if (R_ == Teuchos::null)
00789         R_ = MVT::Clone( *(problem_->getRHS()), 1 ); 
00790 
00791       problem_->computeCurrResVec( &*R_ );
00792 
00793       
00794       // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
00795       // TODO: ensure hypothesis right here ... I have to think about use cases.
00796 
00797       if( U_ != Teuchos::null ){
00798         // Hypothesis: if U_ is not null, then neither is C_ and furthermore U'C= I.
00799 
00800         // possibly over solved equation ...  I want residual norms
00801         // relative to the initial residual, not what I am about to compute.
00802         Teuchos::RCP<MV> cur_soln_vec = problem_->getCurrLHSVec();
00803         std::vector<MagnitudeType> rnorm0(1);
00804         MVT::MvNorm( *R_, rnorm0 ); // rnorm0  = norm(R_);
00805 
00806         // Z := U_'*R_; xo += U_*Z ;R_ -= C_*Z
00807         std::cout  << "Solver Manager:  dimU_ = " << dimU_ << std::endl;
00808         Teuchos::SerialDenseMatrix<int,ScalarType> Z( dimU_, 1 );
00809 
00810         Teuchos::RCP<const MV> Uactive, Cactive;
00811         std::vector<int> active_columns( dimU_ );
00812         for (int i=0; i < dimU_; ++i) active_columns[i] = i;
00813         Uactive = MVT::CloneView(*U_, active_columns);
00814         Cactive = MVT::CloneView(*C_, active_columns);
00815 
00816         if( debug ){
00817           std::cout << " Solver Manager : check duality of seed basis " << std::endl;
00818           Teuchos::SerialDenseMatrix<int,ScalarType> H( dimU_, dimU_ );
00819           MVT::MvTransMv( one, *Uactive, *Cactive, H );
00820           H.print( std::cout );
00821         }
00822         
00823         MVT::MvTransMv( one, *Uactive, *R_, Z );
00824         Teuchos::RCP<MV> tempU = MVT::Clone( *R_, 1 );
00825         MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );  // UZ
00826         MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );  // xo += tmp;
00827         MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );  // CZ
00828         MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );  // R_ -= tmp;
00829         std::vector<MagnitudeType> rnorm(1);
00830         MVT::MvNorm( *R_, rnorm );
00831         if( rnorm[0] < rnorm0[0] * .001 ){  //reorthogonalize
00832           MVT::MvTransMv( one, *Uactive, *R_, Z );
00833           MVT::MvTimesMatAddMv( one, *Uactive, Z, zero, *tempU );
00834           MVT::MvAddMv( one, *tempU, one, *cur_soln_vec, *cur_soln_vec );  // xo += UZ;
00835           MVT::MvTimesMatAddMv( one, *Cactive, Z, zero, *tempU );
00836           MVT::MvAddMv( -one, *tempU, one, *R_, *R_ );  // R_ -= CZ;
00837         }
00838         Uactive = Teuchos::null;
00839         Cactive = Teuchos::null;
00840         tempU = Teuchos::null;
00841       }
00842       else { 
00843         dimU_ = 0;
00844       }
00845 
00846 
00847       // Set the new state and initialize the solver.
00848       PCPGIterState<ScalarType,MV> pcpgState; // fails if R == null.
00849 
00850       pcpgState.R = R_;
00851       if( U_ != Teuchos::null ) pcpgState.U = U_;
00852       if( C_ != Teuchos::null ) pcpgState.C = C_;
00853       if( dimU_ > 0 ) pcpgState.curDim = dimU_;
00854       pcpg_iter->initialize(pcpgState);
00855 
00856       // treat initialize() exceptions here?  how to use try-catch-throw? DMD
00857 
00858       // Get the current number of deflated blocks with the PCPG iteration
00859       dimU_ = pcpgState.curDim;
00860       if( !dimU_ ) printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
00861       pcpg_iter->resetNumIters();
00862 
00863       if( dimU_ > savedBlocks_ )
00864         std::cout << "Error: dimU_  = " << dimU_ << " > savedBlocks_ = " << savedBlocks_ << std::endl; 
00865 
00866       while(1) { // dummy loop for break
00867 
00868         // tell pcpg_iter to iterate
00869         try {
00870           if( debug ) printf("********** Calling iterate...\n");
00871           pcpg_iter->iterate();
00872 
00874           //
00875           // check convergence first
00876           //
00878           if ( convTest_->getStatus() == Passed ) {
00879             // we have convergence
00880             break;  // break from while(1){pcpg_iter->iterate()}
00881           }
00883           //
00884           // check for maximum iterations
00885           //
00887           else if ( maxIterTest_->getStatus() == Passed ) {
00888             // we don't have convergence
00889             isConverged = false;
00890             break;  // break from while(1){pcpg_iter->iterate()}
00891           }
00892           else {
00893 
00895           //
00896           // we returned from iterate(), but none of our status tests Passed.
00897           // Something is wrong, and it is probably the developers fault.
00898           //
00900 
00901             TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
00902                                "Belos::PCPGSolMgr::solve(): Invalid return from PCPGIter::iterate().");
00903           } // end if
00904         } // end try
00905         catch (const PCPGIterOrthoFailure &e) {
00906 
00907           // Check to see if the most recent solution yielded convergence.
00908           sTest_->checkStatus( &*pcpg_iter );
00909           if (convTest_->getStatus() != Passed)
00910             isConverged = false;
00911           break;
00912         }
00913         catch (const std::exception &e) {
00914           printer_->stream(Errors) << "Error! Caught exception in PCPGIter::iterate() at iteration "
00915                                    << pcpg_iter->getNumIters() << std::endl
00916                                    << e.what() << std::endl;
00917           throw;
00918         }
00919       } // end of while(1)
00920 
00921       // Update the linear problem.
00922       Teuchos::RCP<MV> update = pcpg_iter->getCurrentUpdate();
00923       problem_->updateSolution( update, true );
00924 
00925       // Inform the linear problem that we are finished with this block linear system.
00926       problem_->setCurrLS();
00927 
00928       // Get the state.   How did pcpgState die?
00929       PCPGIterState<ScalarType,MV> oldState = pcpg_iter->getState();
00930 
00931       dimU_ = oldState.curDim;
00932       int q = oldState.prevUdim;
00933 
00934       std::cout << "SolverManager: dimU_ " << dimU_ << "   prevUdim= " << q << std::endl;
00935 
00936       if( q > deflatedBlocks_ ) 
00937         std::cout << "SolverManager: Error deflatedBlocks = " << deflatedBlocks_ << std::endl;
00938 
00939       int rank; 
00940       if( dimU_ > q ){ // Orthogonalize [U;C](:,prevUdim:dimU_)
00941         //Given the seed space U and C = A U for some symmetric positive definite A,
00942         //find U1 and C1 with span(U1)=span(U) such that C1'U1 = I maintaining C=AU
00943 
00944         //oldState.D->print( std::cout ); D = diag( C'*U )
00945 
00946         U_ = oldState.U; //MVT::MvPrint( *U_, std::cout ); 
00947         C_ = oldState.C; //MVT::MvPrint( *C_, std::cout ); 
00948         rank = ARRQR(dimU_,q, *oldState.D );
00949         if( rank < dimU_ ) {
00950                 std::cout << " rank decreased in ARRQR, something to do? " << std::endl;
00951         }  
00952         dimU_ = rank;
00953 
00954       } // Now U_ and C_ = AU are dual bases.
00955 
00956       if( dimU_ > deflatedBlocks_ ){
00957 
00958         if( !deflatedBlocks_ ){
00959            U_ = Teuchos::null;
00960            C_ = Teuchos::null;
00961            dimU_ = deflatedBlocks_;
00962            break;
00963         }
00964 
00965         bool Harmonic = false; // (Harmonic) Ritz vectors
00966      
00967         Teuchos::RCP<MV> Uorth;
00968 
00969         std::vector<int> active_cols( dimU_ ); 
00970         for (int i=0; i < dimU_; ++i) active_cols[i] = i;
00971 
00972         if( Harmonic ){
00973           Uorth = MVT::CloneCopy(*C_, active_cols); 
00974         }
00975         else{
00976           Uorth = MVT::CloneCopy(*U_, active_cols); 
00977         }
00978 
00979         // Explicitly construct Q and R factors 
00980         Teuchos::SerialDenseMatrix<int,ScalarType> R(dimU_,dimU_);
00981         rank = ortho_->normalize(*Uorth, Teuchos::rcp(&R,false));
00982         Uorth = Teuchos::null;
00983         // TODO:  During the previous solve, the matrix that normalizes U(1:q) was computed and discarded.  
00984         // One might save it, reuse it here, and just normalize columns U(q+1:dimU_) here.
00985 
00986         // throw an error if U is both A-orthonormal and rank deficient
00987         TEUCHOS_TEST_FOR_EXCEPTION(rank < dimU_,PCPGSolMgrOrthoFailure,
00988                            "Belos::PCPGSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
00989 
00990 
00991         // R VT' = Ur S,   
00992         Teuchos::SerialDenseMatrix<int,ScalarType> VT; // Not referenced
00993         Teuchos::SerialDenseMatrix<int,ScalarType> Ur; // Not referenced
00994         int lwork = 5*dimU_;                           // minimal, extra computation < 67*dimU_
00995         int info = 0;  // Hermite
00996         int lrwork = 1;
00997         if( problem_->isHermitian() ) lrwork = dimU_;
00998         std::vector<ScalarType> work(lwork); // 
00999         std::vector<ScalarType> Svec(dimU_); // 
01000         std::vector<ScalarType> rwork(lrwork); 
01001         lapack.GESVD('N', 'O',
01002                    R.numRows(),R.numCols(),R.values(), R.numRows(),
01003                    &Svec[0],
01004                    Ur.values(),1,
01005                    VT.values(),1, // Output: VT stored in R
01006                    &work[0], lwork,
01007                    &rwork[0], &info);
01008 
01009         TEUCHOS_TEST_FOR_EXCEPTION(info != 0, PCPGSolMgrLAPACKFailure,
01010            "Belos::PCPGSolMgr::solve(): LAPACK _GESVD failed to compute singular values.");
01011 
01012         if( work[0] !=  67. * dimU_ )
01013            std::cout << " SVD " << dimU_ <<  " lwork " << work[0]  << std::endl;
01014         for( int i=0; i< dimU_; i++)
01015            std::cout << i << " " << Svec[i] << std::endl;
01016            
01017         Teuchos::SerialDenseMatrix<int,ScalarType> wholeV( R, Teuchos::TRANS);
01018 
01019         int startRow = 0, startCol = 0;
01020         if( Harmonic )
01021           startCol = dimU_ - deflatedBlocks_;
01022 
01023         Teuchos::SerialDenseMatrix<int,ScalarType> V(Teuchos::Copy,
01024                                                      wholeV,
01025                                                      wholeV.numRows(),
01026                                                      deflatedBlocks_,
01027                                          startRow,
01028                                          startCol);
01029         std::vector<int> active_columns( dimU_ );
01030         std::vector<int> def_cols( deflatedBlocks_ );
01031         for (int i=0; i < dimU_; ++i) active_columns[i] = i;
01032         for (int i=0; i < deflatedBlocks_; ++i) def_cols[i] = i;
01033 
01034         Teuchos::RCP<MV> Uactive = MVT::CloneViewNonConst(*U_, def_cols);
01035         Teuchos::RCP<MV> Ucopy = MVT::CloneCopy( *U_, active_columns );
01036         MVT::MvTimesMatAddMv( one, *Ucopy, V, zero, *Uactive ); //  U:= U*V
01037         Ucopy   = Teuchos::null;
01038         Uactive = Teuchos::null;
01039         Teuchos::RCP<MV> Cactive = MVT::CloneViewNonConst(*C_, def_cols);
01040         Teuchos::RCP<MV> Ccopy = MVT::CloneCopy( *C_, active_columns );
01041         MVT::MvTimesMatAddMv( one, *Ccopy, V, zero, *Cactive ); //  C:= C*V
01042         Ccopy  = Teuchos::null;
01043         Cactive = Teuchos::null;
01044         dimU_ = deflatedBlocks_;
01045       }
01046       printer_->stream(Debug) << " Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << dimU_ << std::endl << std::endl;
01047 
01048       // Inform the linear problem that we are finished with this block linear system.
01049       problem_->setCurrLS();
01050 
01051       // Update indices for the linear systems to be solved.
01052       numRHS2Solve -= 1;
01053       if ( numRHS2Solve > 0 ) {
01054   currIdx[0]++;
01055 
01056         // Set the next indices.
01057         problem_->setLSIndex( currIdx );
01058       }
01059       else {
01060         currIdx.resize( numRHS2Solve );
01061       }
01062     }// while ( numRHS2Solve > 0 )
01063   }
01064     
01065   // print final summary
01066   sTest_->print( printer_->stream(FinalSummary) );
01067   
01068   // print timing information
01069 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01070   // Calling summarize() can be expensive, so don't call unless the
01071   // user wants to print out timing details.  summarize() will do all
01072   // the work even if it's passed a "black hole" output stream.
01073   if (verbosity_ & TimingDetails)
01074     Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
01075 #endif
01076 
01077   // Save the convergence test value ("achieved tolerance") for this solve.
01078   {
01079     using Teuchos::rcp_dynamic_cast;
01080     typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
01081     // testValues is nonnull and not persistent.
01082     const std::vector<MagnitudeType>* pTestValues = 
01083       rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
01084     
01085     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
01086       "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
01087       "method returned NULL.  Please report this bug to the Belos developers.");
01088     
01089     TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
01090       "Belos::PCPGSolMgr::solve(): The convergence test's getTestValue() "
01091       "method returned a vector of length zero.  Please report this bug to the "
01092       "Belos developers.");
01093 
01094     // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
01095     // achieved tolerances for all vectors in the current solve(), or
01096     // just for the vectors from the last deflation?
01097     achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
01098   }
01099  
01100   // get iteration information for this solve
01101   numIters_ = maxIterTest_->getNumIters();
01102  
01103   if (!isConverged) {
01104     return Unconverged; // return from PCPGSolMgr::solve() 
01105   }
01106   return Converged; // return from PCPGSolMgr::solve() 
01107 }
01108 
01109 // A-orthogonalize the Seed Space
01110 // Note that Anasazi::GenOrthoManager provides simplified versions of the algorithm,
01111 // that are not rank revealing, and are not designed for PCPG in other ways too.
01112 template<class ScalarType, class MV, class OP>
01113 int PCPGSolMgr<ScalarType,MV,OP>::ARRQR(int p, int q, const Teuchos::SerialDenseMatrix<int,ScalarType>& D)
01114 {
01115   using Teuchos::RCP;
01116   ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
01117   ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
01118 
01119   // Allocate memory for scalars.
01120   Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
01121   Teuchos::SerialDenseMatrix<int,ScalarType> gamma( 1, 1 );
01122   Teuchos::SerialDenseMatrix<int,ScalarType> anorm( 1, 1 );
01123   std::vector<int> curind(1);
01124   std::vector<int> ipiv(p - q); // RRQR Pivot indices
01125   std::vector<ScalarType> Pivots(p); // RRQR Pivots
01126   int i, imax, j, k, l;
01127   ScalarType rteps = 1.5e-8;
01128 
01129   // Scale such that diag( U'C) = I
01130   for( int i = q ; i < p ; i++ ){
01131     ipiv[i-q] = i;
01132     curind[0] = i;
01133     RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
01134     RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
01135     anorm(0,0) = one / Teuchos::ScalarTraits<ScalarType>::squareroot( D(i-q,i-q) ) ;
01136     MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P );
01137     MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP );
01138     Pivots[i]  = one;
01139   }
01140 
01141   for( i = q ; i < p ; i++ ){  
01142     if( q < i && i < p-1 ){ // Find the largest pivot
01143       imax = i;
01144       l = ipiv[imax-q];
01145       for( j = i+1 ; j < p ; j++ ){  
01146          k = ipiv[j-q];  
01147          if( Pivots[k] > Pivots[l] ){
01148            imax = j;  
01149            l = k;
01150          }
01151       } // end for
01152       if( imax > i ){ 
01153           l = ipiv[imax-q]; // swap ipiv( imax ) and ipiv(i+1)
01154           ipiv[imax-q] = ipiv[i-q];
01155           ipiv[i-q] = l;
01156       }
01157     } // largest pivot found
01158     int k = ipiv[i-q];
01159  
01160     if( Pivots[k]  > 1.5625e-2 ){
01161       anorm(0,0) =  Pivots[k]; // A-norm of u
01162     }
01163     else{ // anorm(0,0) = sqrt( U(:,k)'*C(:,k) );
01164       curind[0] = k;
01165       RCP<const MV> P = MVT::CloneView(*U_,curind);
01166       RCP<const MV> AP = MVT::CloneView(*C_,curind);
01167       MVT::MvTransMv( one, *P, *AP, anorm );
01168       anorm(0,0) = Teuchos::ScalarTraits<ScalarType>::squareroot( anorm(0,0) ) ;
01169     }
01170     if( rteps <= anorm(0,0) && anorm(0,0) < 9.765625e-4){
01171        /*
01172        C(:,k) = A*U(:,k);  % Change C
01173        fixC = U(:, ipiv(1:i-1) )'*C(:,k);
01174        U(:,k) = U(:,k) - U(:, ipiv(1:i-1) )*fixC;
01175        C(:,k) = C(:,k) - C(:, ipiv(1:i-1) )*fixC;
01176        anorm = sqrt( U(:,k)'*C(:,k) );
01177        */
01178        std::cout << "ARRQR: Bad case not implemented" << std::endl;
01179     }
01180     if( anorm(0,0) < rteps ){ // rank [U;C] = i-1
01181        std::cout << "ARRQR : deficient case not implemented " << std::endl;
01182        //U = U(:, ipiv(1:i-1) );
01183        //C = C(:, ipiv(1:i-1) );
01184        p = q + i; 
01185        // update curDim_ in State
01186        break; 
01187     }
01188     curind[0] = k;
01189     RCP<MV> P = MVT::CloneViewNonConst(*U_,curind);
01190     RCP<MV> AP = MVT::CloneViewNonConst(*C_,curind);
01191     MVT::MvAddMv( anorm(0,0), *P, zero, *AP, *P ); // U(:,k) = U(:,k)/anorm;
01192     MVT::MvAddMv( zero, *P, anorm(0,0), *AP, *AP ); // C(:,k) = C(:,k)/anorm;
01193     P = Teuchos::null;
01194     AP = Teuchos::null;
01195     Pivots[k] = one;                 // delete,  for diagonostic purposes
01196     P = MVT::CloneViewNonConst(*U_,curind);  // U(:,k)
01197     AP = MVT::CloneViewNonConst(*C_,curind); // C(:,k)
01198     for( j = i+1 ; j < p ; j++ ){
01199       l = ipiv[j-q];   // ahhh
01200       curind[0] = l;
01201       RCP<MV> Q = MVT::CloneViewNonConst(*U_,curind); // segmentation fault,  j=i+1=5
01202       MVT::MvTransMv( one, *Q, *AP, alpha); // alpha(0,0) = U(:,l)'*C(:,k);
01203       MVT::MvAddMv( -alpha(0,0), *P, one, *Q, *Q ); // U(:,l) -= U(:,k) * alpha(0,0);
01204       Q = Teuchos::null;
01205       RCP<MV> AQ = MVT::CloneViewNonConst(*C_,curind);
01206       MVT::MvAddMv( -alpha(0,0), *AP, one, *AQ, *AQ ); // C(:,l) -= C(:,l) - C(:,k) * alpha(0,0);
01207       AQ = Teuchos::null;
01208       gamma(0,0) = ( Pivots[l] - alpha(0,0))*( Pivots[l] + alpha(0,0));
01209       if( gamma(0,0) > 0){
01210         Pivots[l] = Teuchos::ScalarTraits<ScalarType>::squareroot( gamma(0,0) );
01211       }
01212       else {
01213         Pivots[l] = zero; //rank deficiency revealed
01214       }
01215     }
01216   }
01217   return p;
01218 }
01219 
01220 //  The method returns a string describing the solver manager.
01221 template<class ScalarType, class MV, class OP>
01222 std::string PCPGSolMgr<ScalarType,MV,OP>::description() const
01223 {
01224   std::ostringstream oss;
01225   oss << "Belos::PCPGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
01226   oss << "{";
01227   oss << "Ortho Type='"<<orthoType_;
01228   oss << "}";
01229   return oss.str();
01230 }
01231   
01232 } // end Belos namespace
01233 
01234 #endif /* BELOS_PCPG_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines