Anasazi Version of the Day
AnasaziGeneralizedDavidsonSolMgr.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP
00030 #define ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP
00031 
00036 #include "Teuchos_ParameterList.hpp"
00037 #include "Teuchos_RCPDecl.hpp"
00038 
00039 #include "AnasaziConfigDefs.hpp"
00040 #include "AnasaziTypes.hpp"
00041 #include "AnasaziEigenproblem.hpp"
00042 #include "AnasaziSolverManager.hpp"
00043 #include "AnasaziBasicOrthoManager.hpp"
00044 #include "AnasaziSVQBOrthoManager.hpp"
00045 #include "AnasaziICGSOrthoManager.hpp"
00046 #include "AnasaziBasicOutputManager.hpp"
00047 #include "AnasaziBasicSort.hpp"
00048 #include "AnasaziGeneralizedDavidson.hpp"
00049 #include "AnasaziStatusTestResNorm.hpp"
00050 #include "AnasaziStatusTestWithOrdering.hpp"
00051 
00052 using Teuchos::RCP;
00053 
00057 namespace Anasazi {
00058 
00077 template <class ScalarType, class MV, class OP>
00078 class GeneralizedDavidsonSolMgr : public SolverManager<ScalarType,MV,OP>
00079 {
00080     public:
00081 
00110         GeneralizedDavidsonSolMgr( const RCP< Eigenproblem<ScalarType,MV,OP> > &problem,
00111                                    Teuchos::ParameterList &pl );
00112 
00116         const Eigenproblem<ScalarType,MV,OP> & getProblem() const { return *d_problem; }
00117 
00121         int getNumIters() const { return d_solver->getNumIters(); }
00122 
00127         ReturnType solve();
00128 
00129     private:
00130 
00131         void getRestartState( GeneralizedDavidsonState<ScalarType,MV> &state );
00132 
00133         typedef MultiVecTraits<ScalarType,MV>        MVT;
00134         typedef Teuchos::ScalarTraits<ScalarType>    ST;
00135         typedef typename ST::magnitudeType           MagnitudeType;
00136         typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00137 
00138         RCP< Eigenproblem<ScalarType,MV,OP> >           d_problem;
00139         RCP< GeneralizedDavidson<ScalarType,MV,OP> >    d_solver;
00140         RCP< OutputManager<ScalarType> >                d_outputMan;
00141         RCP< OrthoManager<ScalarType,MV> >              d_orthoMan;
00142         RCP< SortManager<MagnitudeType> >               d_sortMan;
00143         RCP< StatusTest<ScalarType,MV,OP> >             d_tester;
00144         int d_maxRestarts;
00145         int d_restartDim;
00146 
00147 }; // class GeneralizedDavidsonSolMgr
00148 
00149 //---------------------------------------------------------------------------//
00150 // Prevent instantiation on complex scalar type
00151 //---------------------------------------------------------------------------//
00152 template <class MagnitudeType, class MV, class OP>
00153 class GeneralizedDavidsonSolMgr<std::complex<MagnitudeType>,MV,OP>
00154 {
00155   public:
00156 
00157     typedef std::complex<MagnitudeType> ScalarType;
00158     GeneralizedDavidsonSolMgr(
00159             const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00160             Teuchos::ParameterList &pl )
00161     {
00162         // Provide a compile error when attempting to instantiate on complex type
00163         MagnitudeType::this_class_is_missing_a_specialization();
00164     }
00165 };
00166 
00167 //---------------------------------------------------------------------------//
00168 // Start member definitions
00169 //---------------------------------------------------------------------------//
00170 
00171 //---------------------------------------------------------------------------//
00172 // Constructor
00173 //---------------------------------------------------------------------------//
00174 template <class ScalarType, class MV, class OP>
00175 GeneralizedDavidsonSolMgr<ScalarType,MV,OP>::GeneralizedDavidsonSolMgr(
00176         const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00177         Teuchos::ParameterList &pl )
00178    : d_problem(problem)
00179 {
00180     TEUCHOS_TEST_FOR_EXCEPTION( d_problem == Teuchos::null,                std::invalid_argument, "Problem not given to solver manager." );
00181     TEUCHOS_TEST_FOR_EXCEPTION( !d_problem->isProblemSet(),                std::invalid_argument, "Problem not set." );
00182     TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getA() == Teuchos::null &&
00183                                 d_problem->getOperator() == Teuchos::null, std::invalid_argument, "A operator not supplied on Eigenproblem." );
00184     TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getInitVec() == Teuchos::null,  std::invalid_argument, "No vector to clone from on Eigenproblem." );
00185     TEUCHOS_TEST_FOR_EXCEPTION( d_problem->getNEV() <= 0,                  std::invalid_argument, "Number of requested eigenvalues must be positive.");
00186 
00187     if( !pl.isType<int>("Block Size") )
00188     {
00189         pl.set<int>("Block Size",1);
00190     }
00191 
00192     if( !pl.isType<int>("Maximum Subspace Dimension") )
00193     {
00194         pl.set<int>("Maximum Subspace Dimension",3*problem->getNEV()*pl.get<int>("Block Size"));
00195     }
00196 
00197     if( !pl.isType<int>("Print Number of Ritz Values") )
00198     {
00199         int numToPrint = std::max( pl.get<int>("Block Size"), d_problem->getNEV() );
00200         pl.set<int>("Print Number of Ritz Values",numToPrint);
00201     }
00202 
00203     // Get convergence info
00204     MagnitudeType tol = pl.get<MagnitudeType>("Convergence Tolerance", MT::eps() );
00205     TEUCHOS_TEST_FOR_EXCEPTION( pl.get<MagnitudeType>("Convergence Tolerance") <= MT::zero(),
00206                                 std::invalid_argument, "Convergence Tolerance must be greater than zero." );
00207 
00208     // Get maximum restarts
00209     if( pl.isType<int>("Maximum Restarts") )
00210     {
00211         d_maxRestarts = pl.get<int>("Maximum Restarts");
00212         TEUCHOS_TEST_FOR_EXCEPTION( d_maxRestarts < 0, std::invalid_argument, "Maximum Restarts must be non-negative" );
00213     }
00214     else
00215     {
00216         d_maxRestarts = 20;
00217     }
00218 
00219     // Get maximum restarts
00220     d_restartDim = pl.get<int>("Restart Dimension",d_problem->getNEV());
00221     TEUCHOS_TEST_FOR_EXCEPTION( d_restartDim < d_problem->getNEV(),
00222             std::invalid_argument, "Restart Dimension must be at least NEV" );
00223 
00224     // Get initial guess type
00225     std::string initType;
00226     if( pl.isType<std::string>("Initial Guess") )
00227     {
00228         initType = pl.get<std::string>("Initial Guess");
00229         TEUCHOS_TEST_FOR_EXCEPTION( initType!="User" && initType!="Random", std::invalid_argument,
00230                                     "Initial Guess type must be 'User' or 'Random'." );
00231     }
00232     else
00233     {
00234         initType = "User";
00235     }
00236 
00237     // Get sort type
00238     std::string which;
00239     if( pl.isType<std::string>("Which") )
00240     {
00241         which = pl.get<std::string>("Which");
00242         TEUCHOS_TEST_FOR_EXCEPTION( which!="LM" && which!="SM" && which!="LR" && which!="SR" && which!="LI" && which!="SI",
00243                                     std::invalid_argument,
00244                                     "Which must be one of LM,SM,LR,SR,LI,SI." );
00245     }
00246     else
00247     {
00248         which = "LM";
00249     }
00250 
00251     // Build sort manager (currently must be stored as pointer to derived class)
00252     d_sortMan = Teuchos::rcp( new BasicSort<MagnitudeType>(which) );
00253 
00254     // Build orthogonalization manager
00255     std::string ortho = pl.get<std::string>("Orthogonalization","SVQB");
00256     TEUCHOS_TEST_FOR_EXCEPTION( ortho!="DGKS" && ortho!= "SVQB" && ortho!="ICGS", std::invalid_argument,
00257                                 "Anasazi::GeneralizedDavidsonSolMgr::constructor: Invalid orthogonalization type" );
00258 
00259     if( ortho=="DGKS" )
00260     {
00261         d_orthoMan = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>() );
00262     }
00263     else if( ortho=="SVQB" )
00264     {
00265         d_orthoMan = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>() );
00266     }
00267     else if( ortho=="ICGS" )
00268     {
00269         d_orthoMan = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>() );
00270     }
00271 
00272     // Build StatusTest
00273     bool scaleRes  = false; // Always false, scaling the residual is handled by the solver
00274     bool failOnNaN = false;
00275     RCP<StatusTest<ScalarType,MV,OP> > resNormTest = Teuchos::rcp(
00276             new StatusTestResNorm<ScalarType,MV,OP>(tol,d_problem->getNEV(),
00277                                     StatusTestResNorm<ScalarType,MV,OP>::RES_2NORM,scaleRes,failOnNaN) );
00278     d_tester = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(resNormTest,d_sortMan,d_problem->getNEV()) );
00279 
00280     // Build output manager
00281     int verbosity = pl.get<int>("Verbosity",Errors);
00282     d_outputMan = Teuchos::rcp( new BasicOutputManager<ScalarType>() );
00283     d_outputMan->setVerbosity( verbosity );
00284 
00285     // Build solver
00286     d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidsonSolMgr: Building solver" << std::endl;
00287     d_solver = Teuchos::rcp( new GeneralizedDavidson<ScalarType,MV,OP>( problem, d_sortMan, d_outputMan, d_tester, d_orthoMan, pl ) );
00288 }
00289 
00290 //---------------------------------------------------------------------------//
00291 // Solve
00292 //---------------------------------------------------------------------------//
00293 template <class ScalarType, class MV, class OP>
00294 ReturnType GeneralizedDavidsonSolMgr<ScalarType,MV,OP>::solve()
00295 {
00296     Eigensolution<ScalarType,MV> sol;
00297     sol.numVecs = 0;
00298     d_problem->setSolution(sol);
00299 
00300     d_solver->initialize();
00301     int restarts = 0;
00302     while( 1 )
00303     {
00304         // Call iterate on the solver
00305         d_solver->iterate();
00306 
00307         // If the solver converged, we're done
00308         if( d_tester->getStatus() == Passed )
00309             break;
00310 
00311         // If we're already at maximum number of restarts, wrap it up
00312         if( restarts == d_maxRestarts )
00313             break;
00314 
00315         // We need to restart
00316         d_solver->sortProblem( d_restartDim );
00317         GeneralizedDavidsonState<ScalarType,MV> state = d_solver->getState();
00318         getRestartState( state );
00319         d_solver->initialize( state );
00320         restarts++;
00321     }
00322 
00323     // Output final state
00324     if( d_outputMan->isVerbosity(FinalSummary) )
00325         d_solver->currentStatus(d_outputMan->stream(FinalSummary));
00326 
00327     // Fill solution struct
00328     sol.numVecs = d_tester->howMany();
00329     if( sol.numVecs > 0 )
00330     {
00331         std::vector<int> whichVecs = d_tester->whichVecs();
00332         std::vector<int> origIndex = d_solver->getRitzIndex();
00333 
00334         // Make sure no conjugate pairs are split
00335         // Because these are not sorted we have to check all values
00336         for( int i=0; i<sol.numVecs; ++i )
00337         {
00338             if( origIndex[ whichVecs[i] ] == 1 )
00339             {
00340                 if( std::find( whichVecs.begin(), whichVecs.end(), whichVecs[i]+1 ) == whichVecs.end() )
00341                 {
00342                     whichVecs.push_back( whichVecs[i]+1 );
00343                     sol.numVecs++;
00344                 }
00345             }
00346             else if( origIndex[ whichVecs[i] ] == -1 )
00347             {
00348                 if( std::find( whichVecs.begin(), whichVecs.end(), whichVecs[i]-1 ) == whichVecs.end() )
00349                 {
00350                     whichVecs.push_back( whichVecs[i]-1 );
00351                     sol.numVecs++;
00352                 }
00353             }
00354         }
00355 
00356         if( d_outputMan->isVerbosity(Debug) )
00357         {
00358             d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidsonSolMgr: "
00359                 << sol.numVecs << " eigenpairs converged" << std::endl;
00360         }
00361 
00362         // Sort converged values
00363         std::vector< Value<ScalarType> > origVals = d_solver->getRitzValues();
00364         std::vector<MagnitudeType> realParts;
00365         std::vector<MagnitudeType> imagParts;
00366         for( int i=0; i<sol.numVecs; ++i )
00367         {
00368             realParts.push_back( origVals[whichVecs[i]].realpart );
00369             imagParts.push_back( origVals[whichVecs[i]].imagpart );
00370         }
00371 
00372         std::vector<int> permVec(sol.numVecs);
00373         d_sortMan->sort( realParts, imagParts, Teuchos::rcpFromRef(permVec), sol.numVecs );
00374 
00375         // Create new which vector
00376         std::vector<int> newWhich;
00377         for( int i=0; i<sol.numVecs; ++i )
00378             newWhich.push_back( whichVecs[permVec[i]] );
00379 
00380         // Check if converged vectors are ordered
00381         bool ordered = true;
00382         for( int i=0; i<sol.numVecs; ++i )
00383         {
00384             if( newWhich[i]!=i )
00385             {
00386                 ordered = false;
00387                 break;
00388             }
00389         }
00390 
00391         if( ordered  )
00392         {
00393             // Everything is ordered, pull directly from solver and resize
00394             sol.index = origIndex;
00395             sol.index.resize(sol.numVecs);
00396             sol.Evals = d_solver->getRitzValues();
00397             sol.Evals.resize(sol.numVecs);
00398         }
00399         else
00400         {
00401             // Manually copy values into sol
00402 
00403             sol.index.resize(sol.numVecs);
00404             sol.Evals.resize(sol.numVecs);
00405 
00406             for( int i=0; i<sol.numVecs; ++i )
00407             {
00408                 sol.index[i] = origIndex[ newWhich[i] ];
00409                 sol.Evals[i] = origVals[  newWhich[i] ];
00410             }
00411         }
00412         sol.Evecs = MVT::CloneCopy( *(d_solver->getRitzVectors()), newWhich );
00413     }
00414     d_problem->setSolution(sol);
00415 
00416     // Return convergence status
00417     if( sol.numVecs < d_problem->getNEV() )
00418         return Unconverged;
00419 
00420     return Converged;
00421 }
00422 
00423 //---------------------------------------------------------------------------//
00424 // Update GeneralizedDavidson state for restarting
00425 //---------------------------------------------------------------------------//
00426 template <class ScalarType, class MV, class OP>
00427 void GeneralizedDavidsonSolMgr<ScalarType,MV,OP>::getRestartState(
00428         GeneralizedDavidsonState<ScalarType,MV> &state )
00429 {
00430     TEUCHOS_TEST_FOR_EXCEPTION( state.curDim <= d_restartDim, std::runtime_error,
00431             "Anasazi::GeneralizedDavidsonSolMgr: State dimension at restart is smaller than Restart Dimension" );
00432 
00433     std::vector<int> ritzIndex = d_solver->getRitzIndex();
00434 
00435     // Don't split conjugate pair when restarting
00436     int restartDim = d_restartDim;
00437     if( ritzIndex[d_restartDim-1]==1 )
00438         restartDim++;
00439 
00440     d_outputMan->stream(Debug) << " >> Anasazi::GeneralizedDavidsonSolMgr: Restarting with "
00441         << restartDim << " vectors" << std::endl;
00442 
00443     // We have already sorted the problem with d_restartDim "best" values
00444     //  in the leading position.  If we partition the Schur vectors (Z)
00445     //  of the projected problem as Z = [Z_wanted Z_unwanted], then the
00446     //  search subspace after the restart is V_restart = V*Z_wanted
00447     //  (same for AV,BV)
00448 
00449     // Get view of wanted portion of Z
00450     const Teuchos::SerialDenseMatrix<int,ScalarType> Z_wanted =
00451         Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*state.Z,state.curDim,restartDim);
00452 
00453     // Get indices for restart
00454     std::vector<int> allIndices(state.curDim);
00455     for( int i=0; i<state.curDim; ++i )
00456         allIndices[i] = i;
00457 
00458     RCP<const MV>  V_orig = MVT::CloneView( *state.V,  allIndices );
00459 
00460     // Get indices for restart
00461     std::vector<int> restartIndices(restartDim);
00462     for( int i=0; i<restartDim; ++i )
00463         restartIndices[i] = i;
00464 
00465     // Views of subspace vectors to be updated
00466     RCP<MV>  V_restart  = MVT::CloneViewNonConst( *state.V,  restartIndices );
00467 
00468     // Temp storage
00469     RCP<MV> restartVecs = MVT::Clone(*state.V,restartDim);
00470 
00471     // Reset V
00472     MVT::MvTimesMatAddMv(ST::one(),*V_orig,Z_wanted,ST::zero(),*restartVecs);
00473     MVT::SetBlock(*restartVecs,restartIndices,*V_restart);
00474 
00475     // V, Z each have orthonormal columns, therefore V*Z should as well
00476     if( d_outputMan->isVerbosity(Debug) )
00477     {
00478         MagnitudeType orthErr = d_orthoMan->orthonormError(*V_restart);
00479         std::stringstream os;
00480         os << " >> Anasazi::GeneralizedDavidsonSolMgr: Error in V^T V == I after restart : " << orthErr << std::endl;
00481         d_outputMan->print(Debug,os.str());
00482     }
00483 
00484     // Reset AV
00485     RCP<MV> AV_restart  = MVT::CloneViewNonConst( *state.AV, restartIndices );
00486     RCP<const MV> AV_orig = MVT::CloneView( *state.AV, allIndices );
00487 
00488     MVT::MvTimesMatAddMv(ST::one(),*AV_orig,Z_wanted,ST::zero(),*restartVecs);
00489     MVT::SetBlock(*restartVecs,restartIndices,*AV_restart);
00490 
00491     int err;
00492 
00493     // Update matrix projection as Z^{*}(V^{*}AV)Z
00494     const Teuchos::SerialDenseMatrix<int,ScalarType> VAV_orig( Teuchos::View, *state.VAV, state.curDim, state.curDim );
00495     Teuchos::SerialDenseMatrix<int,ScalarType> tmpMat(state.curDim, restartDim);
00496     err = tmpMat.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), VAV_orig, Z_wanted, ST::zero() );
00497     TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, "GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
00498 
00499     Teuchos::SerialDenseMatrix<int,ScalarType> VAV_restart( Teuchos::View, *state.VAV, restartDim, restartDim );
00500     err = VAV_restart.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, ST::one(), Z_wanted, tmpMat, ST::zero() );
00501     TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, "GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
00502 
00503     if( d_problem->getM() != Teuchos::null )
00504     {
00505         // Reset BV
00506         RCP<const MV> BV_orig     = MVT::CloneView( *state.BV, allIndices );
00507         RCP<MV>       BV_restart  = MVT::CloneViewNonConst( *state.BV, restartIndices );
00508 
00509         MVT::MvTimesMatAddMv(ST::one(),*BV_orig,Z_wanted,ST::zero(),*restartVecs);
00510         MVT::SetBlock(*restartVecs,restartIndices,*BV_restart);
00511 
00512 
00513         // Update matrix projection as Z^{*}(V^{*}BV)Z
00514         const Teuchos::SerialDenseMatrix<int,ScalarType> VBV_orig( Teuchos::View, *state.VBV, state.curDim, state.curDim );
00515         err = tmpMat.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, ST::one(), VBV_orig, Z_wanted, ST::zero() );
00516         TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, "GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
00517 
00518         Teuchos::SerialDenseMatrix<int,ScalarType> VBV_restart( Teuchos::View, *state.VBV, restartDim, restartDim );
00519         VBV_restart.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, ST::one(), Z_wanted, tmpMat, ST::zero() );
00520         TEUCHOS_TEST_FOR_EXCEPTION( err!=0, std::runtime_error, "GeneralizedDavidsonSolMgr::getRestartState: multiply returned nonzero error code" );
00521     }
00522 
00523     // Set Q,Z to identity
00524     state.Q->putScalar( ST::zero() );
00525     state.Z->putScalar( ST::zero() );
00526     for( int ii=0; ii<restartDim; ii++ )
00527     {
00528        (*state.Q)(ii,ii)= ST::one();
00529        (*state.Z)(ii,ii)= ST::one();
00530     }
00531 
00532     // Update current dimension
00533     state.curDim = restartDim;
00534 }
00535 
00536 } // namespace Anasazi
00537 
00538 #endif // ANASAZI_GENERALIZED_DAVIDSON_SOLMGR_HPP
00539 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends