00001
00002
00003
00004
00005 #ifndef PLAYA_ANASAZIEIGENSOLVER_IMPL_HPP
00006 #define PLAYA_ANASAZIEIGENSOLVER_IMPL_HPP
00007
00008 #include "PlayaDefs.hpp"
00009 #include "PlayaAnasaziEigensolverDecl.hpp"
00010 #include "PlayaParameterListPreconditionerFactory.hpp"
00011 #include "PlayaPreconditionerFactory.hpp"
00012 #include "AnasaziBlockKrylovSchurSolMgr.hpp"
00013 #include "AnasaziBlockDavidsonSolMgr.hpp"
00014 #include "AnasaziSimpleLOBPCGSolMgr.hpp"
00015 #include "AnasaziLOBPCGSolMgr.hpp"
00016 #include "AnasaziBasicEigenproblem.hpp"
00017 #include "PlayaAnasaziAdapter.hpp"
00018
00019
00020 namespace Playa
00021 {
00022 using Teuchos::ParameterList;
00023 using Anasazi::SimpleMV;
00024
00025 template <class MV, class OP>
00026 class InitTraits
00027 {
00028 public:
00029
00030 static RCP<OP> opPtr(const LinearOperator<double>& A);
00031
00032
00033 static RCP<MV> makeMV(int numVecs, const VectorSpace<double>& space);
00034
00035
00036 static Vector<double> vec(const RCP<MV>& mv, int i);
00037 };
00038
00039
00040
00041 template <> class InitTraits<SimpleMV, LinearOperator<double> >
00042 {
00043 public:
00044 typedef SimpleMV MV;
00045 typedef LinearOperator<double> OP;
00046
00047
00048 static RCP<OP> opPtr(const LinearOperator<double>& A)
00049 {
00050 if (A.ptr().get() != 0)
00051 return rcp(new LinearOperator<double>(A));
00052 else
00053 {
00054 RCP<LinearOperator<double> > rtn;
00055 return rtn;
00056 }
00057 }
00058
00059
00060 static RCP<MV> makeMV(int blockSize, const VectorSpace<double>& space)
00061 {
00062 RCP<MV> mv = rcp(new MV(blockSize));
00063 for (int i=0; i<blockSize; i++) (*mv)[i] = space.createMember();
00064 return mv;
00065 }
00066
00067
00068 static Vector<double> vec(const RCP<MV>& mv, int i)
00069 {
00070 return (*mv)[i];
00071 }
00072
00073
00074 };
00075
00076
00077
00078
00079 template <class Scalar>
00080 inline void AnasaziEigensolver<Scalar>::solve(
00081 const LinearOperator<Scalar>& K,
00082 const LinearOperator<Scalar>& M,
00083 Array<Vector<Scalar> >& evecs,
00084 Array<std::complex<Scalar> >& ew) const
00085 {
00086 typedef SimpleMV MV;
00087 typedef LinearOperator<Scalar> OP;
00088
00089 typedef Anasazi::MultiVecTraits<Scalar,MV> MVT;
00090 typedef Anasazi::OperatorTraits<Scalar,MV,OP> OPT;
00091
00092 TimeMonitor timer(solveTimer());
00093 VectorSpace<Scalar> KDomain = K.domain();
00094
00095 RCP<OP> KPtr = InitTraits<MV, OP>::opPtr(K);
00096 RCP<OP> MPtr = InitTraits<MV, OP>::opPtr(M);
00097
00098
00099
00100
00101
00102 std::string method = this->params().get<string>("Method");
00103 int numEigs = this->params().get<int>("Number of Eigenvalues");
00104 int blockSize = this->params().get<int>("Block Size");
00105 bool usePrec = this->params().get<bool>("Use Preconditioner");
00106 bool hermitian = this->params().get<bool>("Is Hermitian");
00107
00108
00109
00110
00111
00112 RCP<MV> mv = InitTraits<MV, OP>::makeMV(blockSize, KDomain);
00113
00114
00115 MVT::MvRandom( *mv );
00116
00117
00118 ParameterList precParams = this->params().sublist("Preconditioner");
00119 PreconditionerFactory<double> precFactory
00120 = new ParameterListPreconditionerFactory(precParams);
00121
00122 LinearOperator<Scalar> P;
00123 if (usePrec)
00124 {
00125 TimeMonitor pTimer(precondBuildTimer());
00126 P = precFactory.createPreconditioner(K).right();
00127 }
00128
00129
00130 RCP<Anasazi::Eigenproblem<Scalar,MV,OP> > problem;
00131
00132 if (MPtr.get() != 0)
00133 {
00134 problem =
00135 rcp( new Anasazi::BasicEigenproblem<Scalar,MV,OP>(KPtr, MPtr, mv) );
00136 }
00137 else
00138 {
00139 problem =
00140 rcp( new Anasazi::BasicEigenproblem<Scalar,MV,OP>(KPtr, mv) );
00141 }
00142
00143 ParameterList eigParams = this->params();
00144 problem->setHermitian(hermitian);
00145 problem->setNEV(numEigs);
00146 if (usePrec) problem->setPrec(InitTraits<MV, OP>::opPtr(P));
00147
00148 bool ret = problem->setProblem();
00149 TEUCHOS_TEST_FOR_EXCEPTION(!ret, std::runtime_error,
00150 "Eigenproblem not setup correctly");
00151
00152
00153
00154 RCP<Anasazi::SolverManager<Scalar,MV,OP> > MySolverMan;
00155 if (method=="Block Davidson")
00156 {
00157 MySolverMan = rcp(new Anasazi::BlockDavidsonSolMgr<Scalar,MV,OP>(problem, eigParams));
00158 }
00159 else if (method=="Block Krylov Schur")
00160 {
00161 MySolverMan = rcp(new Anasazi::BlockKrylovSchurSolMgr<Scalar,MV,OP>(problem, eigParams));
00162 }
00163 else if (method=="Simple LOBPCG")
00164 {
00165 MySolverMan = rcp(new Anasazi::SimpleLOBPCGSolMgr<Scalar,MV,OP>(problem, eigParams));
00166 }
00167 else if (method=="LOBPCG")
00168 {
00169 MySolverMan = rcp(new Anasazi::LOBPCGSolMgr<Scalar,MV,OP>(problem, eigParams));
00170 }
00171 else
00172 {
00173 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00174 "solver method [" << method << "] not recognized");
00175 }
00176
00177
00178 Anasazi::ReturnType returnCode = MySolverMan->solve();
00179 Out::os() << "return code = " << returnCode << endl;
00180 TEUCHOS_TEST_FOR_EXCEPTION(returnCode != Anasazi::Converged,
00181 std::runtime_error, "Anasazi did not converge!");
00182
00183
00184 Anasazi::Eigensolution<Scalar,MV> sol = problem->getSolution();
00185 RCP<MV> evecs_mv = sol.Evecs;
00186 int numev = sol.numVecs;
00187
00188
00189 ew.resize(numev);
00190 evecs.resize(numev);
00191
00192 for (int i=0; i<numev; i++)
00193 {
00194 Vector<Scalar> tmp = InitTraits<MV, OP>::vec(evecs_mv, i);
00195
00196 evecs[i] = KDomain.createMember();
00197 evecs[i].acceptCopyOf(tmp);
00198
00199
00200
00201
00202 ew[i].real() = sol.Evals[i].realpart;
00203 ew[i].imag() = sol.Evals[i].imagpart;
00204 }
00205 }
00206
00207
00208 }
00209
00210
00211 #endif