00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "NOX_Abstract_MultiVector.H"
00043 #include "Teuchos_ParameterList.hpp"
00044 #include "LOCA_Parameter_SublistParser.H"
00045 #include "LOCA_GlobalData.H"
00046 #include "LOCA_ErrorCheck.H"
00047 #include "NOX_Utils.H"
00048 #include "LOCA_Factory.H"
00049 #include "LOCA_Eigensolver_AnasaziStrategy.H"
00050 #include "LOCA_EigenvalueSort_Strategies.H"
00051
00052 #ifdef HAVE_LOCA_ANASAZI
00053 #include "Anasazi_LOCA_MultiVecTraits.H"
00054 #include "Anasazi_LOCA_OperatorTraits.H"
00055 #include "AnasaziBlockKrylovSchurSolMgr.hpp"
00056 #include "AnasaziBasicEigenproblem.hpp"
00057 #endif
00058
00059 LOCA::Eigensolver::AnasaziStrategy::AnasaziStrategy(
00060 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00061 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams_,
00062 const Teuchos::RCP<Teuchos::ParameterList>& eigParams) :
00063 globalData(global_data),
00064 topParams(topParams_),
00065 eigenParams(eigParams),
00066 solverParams(),
00067 blksz(1),
00068 nev(4),
00069 isSymmetric(false),
00070 locaSort()
00071 {
00072
00073 solverParams = Teuchos::rcp(new Teuchos::ParameterList);
00074 *solverParams = *(topParams->getSublist("Linear Solver"));
00075 if (Teuchos::isParameterType<double>(*eigenParams,"Linear Solve Tolerance"))
00076 solverParams->set("Tolerance", Teuchos::get<double>
00077 (*eigenParams, "Linear Solve Tolerance"));
00078
00079
00080 blksz = eigenParams->get("Block Size", 1);
00081 nev = eigenParams->get("Num Eigenvalues", 4);
00082 isSymmetric = eigenParams->get("Symmetric", false);
00083
00084
00085 eigenParams->get("Convergence Tolerance", 1.0e-7);
00086 eigenParams->get("Maximum Restarts", 1);
00087 eigenParams->get("Num Blocks", 30);
00088 eigenParams->get("Step Size", 1);
00089 if (!eigenParams->isParameter("Verbosity"))
00090 eigenParams->set("Verbosity",
00091 Anasazi::Errors +
00092 Anasazi::Warnings +
00093 Anasazi::FinalSummary);
00094
00095
00096 Teuchos::RCP<LOCA::EigenvalueSort::AbstractStrategy> sortingStrategy
00097 = globalData->locaFactory->createEigenvalueSortStrategy(topParams,
00098 eigenParams);
00099 locaSort =
00100 Teuchos::rcp(new Anasazi::LOCASort(globalData, sortingStrategy));
00101 eigenParams->set( "Sort Manager", locaSort );
00102 }
00103
00104 LOCA::Eigensolver::AnasaziStrategy::~AnasaziStrategy()
00105 {
00106 }
00107
00108 NOX::Abstract::Group::ReturnType
00109 LOCA::Eigensolver::AnasaziStrategy::computeEigenvalues(
00110 NOX::Abstract::Group& group,
00111 Teuchos::RCP< std::vector<double> >& evals_r,
00112 Teuchos::RCP< std::vector<double> >& evals_i,
00113 Teuchos::RCP< NOX::Abstract::MultiVector >& evecs_r,
00114 Teuchos::RCP< NOX::Abstract::MultiVector >& evecs_i)
00115 {
00116 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperIteration)) {
00117 globalData->locaUtils->out() << "\n" <<
00118 globalData->locaUtils->fill(64,'=') <<
00119 "\nAnasazi Eigensolver starting with block size " << blksz <<
00120 "\n" << std::endl;
00121 }
00122
00123
00124 const NOX::Abstract::Vector& xVector = group.getX();
00125
00126
00127 Teuchos::RCP<LOCA::AnasaziOperator::AbstractStrategy> anasaziOp
00128 = globalData->locaFactory->createAnasaziOperatorStrategy(
00129 topParams,
00130 eigenParams,
00131 solverParams,
00132 Teuchos::rcp(&group,false));
00133 Teuchos::RCP<MV> ivec = xVector.createMultiVector(blksz);
00134 ivec->random();
00135
00136
00137 anasaziOp->preProcessSeedVector(*ivec);
00138
00139
00140 Teuchos::RCP<Anasazi::BasicEigenproblem<double, MV, OP> >
00141 LOCAProblem =
00142 Teuchos::rcp( new Anasazi::BasicEigenproblem<double, MV, OP>(anasaziOp,
00143 ivec) );
00144
00145
00146 LOCAProblem->setNEV( nev );
00147
00148
00149 LOCAProblem->setHermitian(isSymmetric);
00150
00151
00152
00153 LOCAProblem->setProblem();
00154
00155
00156 Anasazi::BlockKrylovSchurSolMgr<double, MV, OP>
00157 LOCABlockKrylovSchur(LOCAProblem, *eigenParams);
00158
00159
00160 Anasazi::ReturnType returnCode = LOCABlockKrylovSchur.solve();
00161
00162
00163 const Anasazi::Eigensolution<double,MV>& anasaziSolution =
00164 LOCAProblem->getSolution();
00165 int numVecs = anasaziSolution.numVecs;
00166 evals_r =
00167 Teuchos::rcp(new std::vector<double>(numVecs));
00168 evals_i =
00169 Teuchos::rcp(new std::vector<double>(numVecs));
00170 for (int i=0; i<numVecs; i++) {
00171 (*evals_r)[i] = anasaziSolution.Evals[i].realpart;
00172 (*evals_i)[i] = anasaziSolution.Evals[i].imagpart;
00173 }
00174
00175 if (returnCode != Anasazi::Converged)
00176 globalData->locaUtils->out() <<
00177 globalData->locaUtils->fill(72, '*') << std::endl <<
00178 "WARNING: Anasazi eigensolver did not converge." << std::endl <<
00179 " Only " << numVecs << " of " << nev <<
00180 " eigenvalues were computed!" << std::endl <<
00181 globalData->locaUtils->fill(72, '*') << std::endl << std::endl;
00182
00183 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperIteration))
00184 globalData->locaUtils->out() <<
00185 "Untransformed eigenvalues (since the operator was " <<
00186 anasaziOp->label() << ")" << std::endl;
00187
00188
00189 Teuchos::RCP<MV> evecs = anasaziSolution.Evecs;
00190
00191 evecs_r = evecs->clone(numVecs);
00192 evecs_i = evecs->clone(numVecs);
00193 for (int i=0; i<numVecs; i++) {
00194
00195
00196 if (anasaziSolution.index[i] == 0) {
00197 (*evecs_r)[i] = (*evecs)[i];
00198 (*evecs_i)[i].init(0.0);
00199 }
00200
00201
00202 else if (anasaziSolution.index[i] == 1) {
00203 (*evecs_r)[i] = (*evecs)[i];
00204 (*evecs_i)[i] = (*evecs)[i+1];
00205 }
00206
00207
00208
00209 else if (anasaziSolution.index[i] == -1) {
00210 (*evecs_r)[i] = (*evecs)[i-1];
00211 (*evecs_i)[i].update(-1.0, (*evecs)[i], 0.0);
00212 }
00213 else {
00214 string func = "LOCA::Eigensolver::AnasaziStrategy::computeEigenvalues()";
00215 stringstream ss;
00216 ss << "Unknown anasazi index " << anasaziSolution.index[i];
00217 globalData->locaErrorCheck->throwError(func, ss.str());
00218 }
00219 }
00220
00221
00222 double rq_r, rq_i;
00223
00224 for (int i=0; i<numVecs; i++) {
00225
00226
00227 anasaziOp->transformEigenvalue((*evals_r)[i], (*evals_i)[i]);
00228
00229
00230 anasaziOp->rayleighQuotient((*evecs_r)[i], (*evecs_i)[i], rq_r, rq_i);
00231
00232
00233
00234 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperIteration)) {
00235 globalData->locaUtils->out() << "Eigenvalue " << i << " : " <<
00236 globalData->locaUtils->sciformat((*evals_r)[i]) << " " <<
00237 globalData->locaUtils->sciformat((*evals_i)[i]) <<
00238 " i : RQresid " <<
00239 globalData->locaUtils->sciformat(fabs((*evals_r)[i] - rq_r)) <<
00240 " " <<
00241 globalData->locaUtils->sciformat(fabs((*evals_i)[i] - rq_i)) <<
00242 " i" << std::endl;
00243 }
00244
00245 }
00246
00247
00248 std::vector<Anasazi::Value<double> > ritzValues =
00249 LOCABlockKrylovSchur.getRitzValues();
00250 int numRitz = ritzValues.size();
00251 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperIteration) &&
00252 numRitz>numVecs) {
00253 globalData->locaUtils->out() <<
00254 "~~~~~~~ remaining eigenvalue approximations ~~~~~~~~~~~~" << std::endl;
00255 }
00256 for (int i=numVecs; i<numRitz; i++) {
00257
00258
00259 anasaziOp->transformEigenvalue(ritzValues[i].realpart,
00260 ritzValues[i].imagpart);
00261
00262 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperIteration)) {
00263 globalData->locaUtils->out() <<
00264 "Eigenvalue " << i << " : " <<
00265 globalData->locaUtils->sciformat(ritzValues[i].realpart) << " " <<
00266 globalData->locaUtils->sciformat(ritzValues[i].imagpart) << " i" <<
00267 std::endl;
00268 }
00269
00270 }
00271
00272 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperIteration)) {
00273 globalData->locaUtils->out() <<
00274 "\nAnasazi Eigensolver finished.\n" <<
00275 globalData->locaUtils->fill(64,'=') << "\n" << std::endl;
00276 }
00277
00278 if (returnCode == Anasazi::Converged)
00279 return NOX::Abstract::Group::Ok;
00280 else
00281 return NOX::Abstract::Group::NotConverged;
00282 }
00283