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 "LOCA_Thyra_Group.H"
00043 #include "NOX_Thyra_MultiVector.H"
00044 #include "Teuchos_TestForException.hpp"
00045 #include "Thyra_ModelEvaluator.hpp"
00046 #include "Thyra_SolveSupportTypes.hpp"
00047 #include "Thyra_DetachedMultiVectorView.hpp"
00048 #include "Thyra_VectorStdOps.hpp"
00049 #include "Teuchos_ParameterList.hpp"
00050 #include "LOCA_GlobalData.H"
00051 #include "LOCA_ErrorCheck.H"
00052
00053 LOCA::Thyra::Group::Group(
00054 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00055 const NOX::Thyra::Vector& initial_guess,
00056 const Teuchos::RCP< ::Thyra::ModelEvaluator<double> >& model,
00057 const LOCA::ParameterVector& p,
00058 int p_index,
00059 bool impl_dfdp) :
00060 NOX::Thyra::Group(initial_guess, model),
00061 LOCA::Abstract::Group(global_data),
00062 globalData(global_data),
00063 params(p),
00064 param_index(p_index),
00065 saveDataStrategy(),
00066 implement_dfdp(impl_dfdp)
00067 {
00068
00069
00070 RTOpPack::SubVectorView<double> pv(0,params.length(),
00071 params.getDoubleArrayPointer(),1);
00072 Teuchos::RCP<const ::Thyra::VectorSpaceBase<double> > ps =
00073 model_->get_p_space(param_index);
00074 param_thyra_vec = ::Thyra::createMemberView(ps,pv);
00075
00076
00077 Teuchos::RCP<const ::Thyra::VectorSpaceBase<double> > xs =
00078 model_->get_x_space();
00079 x_dot_vec = ::Thyra::createMember(xs);
00080 ::Thyra::put_scalar(0.0, x_dot_vec.get());
00081 }
00082
00083 LOCA::Thyra::Group::Group(const LOCA::Thyra::Group& source,
00084 NOX::CopyType type) :
00085 NOX::Thyra::Group(source, type),
00086 LOCA::Abstract::Group(source, type),
00087 globalData(source.globalData),
00088 params(source.params),
00089 param_index(source.param_index),
00090 saveDataStrategy(source.saveDataStrategy),
00091 implement_dfdp(source.implement_dfdp)
00092 {
00093
00094
00095 RTOpPack::SubVectorView<double> pv(0,params.length(),
00096 params.getDoubleArrayPointer(),1);
00097 Teuchos::RCP<const ::Thyra::VectorSpaceBase<double> > ps =
00098 model_->get_p_space(param_index);
00099 param_thyra_vec = ::Thyra::createMemberView(ps,pv);
00100
00101
00102 Teuchos::RCP<const ::Thyra::VectorSpaceBase<double> > xs =
00103 model_->get_x_space();
00104 x_dot_vec = ::Thyra::createMember(xs);
00105 ::Thyra::put_scalar(0.0, x_dot_vec.get());
00106 }
00107
00108 LOCA::Thyra::Group::~Group()
00109 {
00110 }
00111
00112 LOCA::Thyra::Group&
00113 LOCA::Thyra::Group::operator=(const LOCA::Thyra::Group& source)
00114 {
00115 if (this != &source) {
00116 NOX::Thyra::Group::operator=(source);
00117 LOCA::Abstract::Group::copy(source);
00118 params = source.params;
00119 param_index = source.param_index;
00120 saveDataStrategy = source.saveDataStrategy;
00121 implement_dfdp = source.implement_dfdp;
00122
00123
00124 }
00125 return *this;
00126 }
00127
00128 NOX::Abstract::Group&
00129 LOCA::Thyra::Group::operator=(const NOX::Abstract::Group& source)
00130 {
00131 operator=(dynamic_cast<const Group&> (source));
00132 return *this;
00133 }
00134
00135 NOX::Abstract::Group&
00136 LOCA::Thyra::Group::operator=(const NOX::Thyra::Group& source)
00137 {
00138 operator=(dynamic_cast<const Group&> (source));
00139 return *this;
00140 }
00141
00142 Teuchos::RCP<NOX::Abstract::Group>
00143 LOCA::Thyra::Group::clone(NOX::CopyType type) const
00144 {
00145 return Teuchos::rcp(new LOCA::Thyra::Group(*this, type));
00146 }
00147
00148 NOX::Abstract::Group::ReturnType
00149 LOCA::Thyra::Group::computeF()
00150 {
00151 if (this->isF())
00152 return NOX::Abstract::Group::Ok;
00153
00154 in_args_.set_x(x_vec_->getThyraRCPVector().assert_not_null());
00155 if (in_args_.supports(::Thyra::ModelEvaluatorBase::IN_ARG_x_dot))
00156 in_args_.set_x_dot(x_dot_vec);
00157 in_args_.set_p(param_index, param_thyra_vec);
00158 out_args_.set_f(f_vec_->getThyraRCPVector().assert_not_null());
00159
00160 model_->evalModel(in_args_, out_args_);
00161
00162 in_args_.set_x(Teuchos::null);
00163 in_args_.set_p(param_index, Teuchos::null);
00164 out_args_.set_f(Teuchos::null);
00165
00166 is_valid_f_ = true;
00167
00168 if (out_args_.isFailed())
00169 return NOX::Abstract::Group::Failed;
00170
00171 return NOX::Abstract::Group::Ok;
00172 }
00173
00174 NOX::Abstract::Group::ReturnType
00175 LOCA::Thyra::Group::computeJacobian()
00176 {
00177 if (this->isJacobian())
00178 return NOX::Abstract::Group::Ok;
00179
00180 in_args_.set_x(x_vec_->getThyraRCPVector());
00181 if (in_args_.supports(::Thyra::ModelEvaluatorBase::IN_ARG_x_dot))
00182 in_args_.set_x_dot(x_dot_vec);
00183 if (in_args_.supports(::Thyra::ModelEvaluatorBase::IN_ARG_alpha))
00184 in_args_.set_alpha(0.0);
00185 if (in_args_.supports(::Thyra::ModelEvaluatorBase::IN_ARG_beta))
00186 in_args_.set_beta(1.0);
00187 in_args_.set_p(param_index, param_thyra_vec);
00188 out_args_.set_W(shared_jacobian_->getObject(this));
00189
00190 model_->evalModel(in_args_, out_args_);
00191
00192 in_args_.set_x(Teuchos::null);
00193 in_args_.set_p(param_index, Teuchos::null);
00194 out_args_.set_W(Teuchos::null);
00195
00196 is_valid_jacobian_ = true;
00197
00198 if (out_args_.isFailed())
00199 return NOX::Abstract::Group::Failed;
00200
00201 return NOX::Abstract::Group::Ok;
00202 }
00203
00204 void
00205 LOCA::Thyra::Group::copy(const NOX::Abstract::Group& source)
00206 {
00207 *this = source;
00208 }
00209
00210 void
00211 LOCA::Thyra::Group::setParams(const LOCA::ParameterVector& p)
00212 {
00213 this->resetIsValidFlags();
00214 params = p;
00215 }
00216
00217 void
00218 LOCA::Thyra::Group::setParam(int paramID, double val)
00219 {
00220 this->resetIsValidFlags();
00221 params.setValue(paramID, val);
00222 }
00223
00224 void
00225 LOCA::Thyra::Group::setParam(string paramID, double val)
00226 {
00227 this->resetIsValidFlags();
00228 params.setValue(paramID, val);
00229 }
00230
00231 const LOCA::ParameterVector&
00232 LOCA::Thyra::Group::getParams() const
00233 {
00234 return params;
00235 }
00236
00237 double
00238 LOCA::Thyra::Group::getParam(int paramID) const
00239 {
00240 return params.getValue(paramID);
00241 }
00242
00243 double
00244 LOCA::Thyra::Group::getParam(string paramID) const
00245 {
00246 return params.getValue(paramID);
00247 }
00248
00249 NOX::Abstract::Group::ReturnType
00250 LOCA::Thyra::Group::computeDfDpMulti(const vector<int>& paramIDs,
00251 NOX::Abstract::MultiVector& fdfdp,
00252 bool isValidF)
00253 {
00254
00255
00256
00257 implement_dfdp = false;
00258
00259
00260
00261 if (!implement_dfdp ||
00262 !out_args_.supports(::Thyra::ModelEvaluatorBase::OUT_ARG_DfDp,
00263 param_index).supports(::Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL)) {
00264 NOX::Abstract::Group::ReturnType res =
00265 LOCA::Abstract::Group::computeDfDpMulti(paramIDs, fdfdp, isValidF);
00266 return res;
00267 }
00268
00269
00270 int num_vecs = fdfdp.numVectors()-1;
00271 std::vector<int> index_dfdp(num_vecs);
00272 for (int i=0; i<num_vecs; i++)
00273 index_dfdp[i] = i+1;
00274 NOX::Thyra::Vector& f = dynamic_cast<NOX::Thyra::Vector&>(fdfdp[0]);
00275 Teuchos::RCP<NOX::Abstract::MultiVector> dfdp =
00276 fdfdp.subView(index_dfdp);
00277
00278
00279
00280
00281
00282 int np = params.length();
00283 Teuchos::RCP<NOX::Thyra::MultiVector> dfdp_full =
00284 Teuchos::rcp_dynamic_cast<NOX::Thyra::MultiVector>(dfdp->clone(np));
00285
00286 ::Thyra::ModelEvaluatorBase::DerivativeMultiVector<double> dmv(dfdp_full->getThyraMultiVector(), ::Thyra::ModelEvaluatorBase::DERIV_MV_BY_COL);
00287 ::Thyra::ModelEvaluatorBase::Derivative<double> deriv(dmv);
00288
00289 in_args_.set_x(x_vec_->getThyraRCPVector().assert_not_null());
00290 if (in_args_.supports(::Thyra::ModelEvaluatorBase::IN_ARG_x_dot))
00291 in_args_.set_x_dot(x_dot_vec);
00292 in_args_.set_p(param_index, param_thyra_vec);
00293 if (!isValidF)
00294 out_args_.set_f(f.getThyraRCPVector().assert_not_null());
00295 out_args_.set_DfDp(param_index, deriv);
00296
00297
00298 model_->evalModel(in_args_, out_args_);
00299
00300
00301 for (int i=0; i<num_vecs; i++)
00302 (*dfdp)[i] = (*dfdp_full)[paramIDs[i]];
00303
00304
00305 in_args_.set_x(Teuchos::null);
00306 in_args_.set_p(param_index, Teuchos::null);
00307 out_args_.set_f(Teuchos::null);
00308 out_args_.set_DfDp(param_index,
00309 ::Thyra::ModelEvaluatorBase::Derivative<double>());
00310
00311 if (out_args_.isFailed())
00312 return NOX::Abstract::Group::Failed;
00313
00314 return NOX::Abstract::Group::Ok;
00315 }
00316
00317 void
00318 LOCA::Thyra::Group::preProcessContinuationStep(
00319 LOCA::Abstract::Iterator::StepStatus stepStatus)
00320 {
00321 if (saveDataStrategy != Teuchos::null)
00322 saveDataStrategy->preProcessContinuationStep(stepStatus);
00323 }
00324
00325 void
00326 LOCA::Thyra::Group::postProcessContinuationStep(
00327 LOCA::Abstract::Iterator::StepStatus stepStatus)
00328 {
00329 if (saveDataStrategy != Teuchos::null)
00330 saveDataStrategy->postProcessContinuationStep(stepStatus);
00331 }
00332
00333 void
00334 LOCA::Thyra::Group::projectToDraw(const NOX::Abstract::Vector& x,
00335 double *px) const
00336 {
00337 if (saveDataStrategy != Teuchos::null)
00338 saveDataStrategy->projectToDraw(x, px);
00339 }
00340
00341 int
00342 LOCA::Thyra::Group::projectToDrawDimension() const
00343 {
00344 if (saveDataStrategy != Teuchos::null)
00345 return saveDataStrategy->projectToDrawDimension();
00346 return 0;
00347 }
00348
00349 double
00350 LOCA::Thyra::Group::computeScaledDotProduct(
00351 const NOX::Abstract::Vector& a,
00352 const NOX::Abstract::Vector& b) const
00353 {
00354 return a.innerProduct(b) / a.length();
00355 }
00356
00357 void
00358 LOCA::Thyra::Group::printSolution(const double conParam) const
00359 {
00360 printSolution(*x_vec_, conParam);
00361 }
00362
00363 void
00364 LOCA::Thyra::Group::printSolution(const NOX::Abstract::Vector& x_,
00365 const double conParam) const
00366 {
00367 if (saveDataStrategy != Teuchos::null)
00368 saveDataStrategy->saveSolution(x_, conParam);
00369 }
00370
00371 void
00372 LOCA::Thyra::Group::scaleVector(NOX::Abstract::Vector& x) const
00373 {
00374 x.scale(1.0 / sqrt(static_cast<double>(x.length())));
00375 }
00376
00377 NOX::Abstract::Group::ReturnType
00378 LOCA::Thyra::Group::computeShiftedMatrix(double alpha, double beta)
00379 {
00380 in_args_.set_x(x_vec_->getThyraRCPVector());
00381 if (in_args_.supports(::Thyra::ModelEvaluatorBase::IN_ARG_x_dot))
00382 in_args_.set_x_dot(x_dot_vec);
00383 in_args_.set_p(param_index, param_thyra_vec);
00384 in_args_.set_alpha(-beta);
00385 in_args_.set_beta(alpha);
00386 out_args_.set_W(shared_jacobian_->getObject(this));
00387
00388 model_->evalModel(in_args_, out_args_);
00389
00390 in_args_.set_x(Teuchos::null);
00391 in_args_.set_p(param_index, Teuchos::null);
00392 in_args_.set_alpha(0.0);
00393 in_args_.set_beta(1.0);
00394 out_args_.set_W(Teuchos::null);
00395
00396 is_valid_jacobian_ = false;
00397
00398 if (out_args_.isFailed())
00399 return NOX::Abstract::Group::Failed;
00400
00401 return NOX::Abstract::Group::Ok;
00402 }
00403
00404 NOX::Abstract::Group::ReturnType
00405 LOCA::Thyra::Group::applyShiftedMatrix(const NOX::Abstract::Vector& input,
00406 NOX::Abstract::Vector& result) const
00407 {
00408 const NOX::Thyra::Vector& thyra_input =
00409 dynamic_cast<const NOX::Thyra::Vector&>(input);
00410 NOX::Thyra::Vector& thyra_result =
00411 dynamic_cast<NOX::Thyra::Vector&>(result);
00412
00413 ::Thyra::apply(*shared_jacobian_->getObject(), ::Thyra::NOTRANS,
00414 thyra_input.getThyraVector(), &thyra_result.getThyraVector());
00415
00416 return NOX::Abstract::Group::Ok;
00417 }
00418
00419 NOX::Abstract::Group::ReturnType
00420 LOCA::Thyra::Group::applyShiftedMatrixMultiVector(
00421 const NOX::Abstract::MultiVector& input,
00422 NOX::Abstract::MultiVector& result) const
00423 {
00424 const NOX::Thyra::MultiVector& nt_input =
00425 Teuchos::dyn_cast<const NOX::Thyra::MultiVector>(input);
00426 NOX::Thyra::MultiVector& nt_result =
00427 Teuchos::dyn_cast<NOX::Thyra::MultiVector>(result);
00428
00429 ::Thyra::apply(*shared_jacobian_->getObject(),
00430 ::Thyra::NOTRANS,
00431 *nt_input.getThyraMultiVector(),
00432 nt_result.getThyraMultiVector().get());
00433
00434 return NOX::Abstract::Group::Ok;
00435 }
00436
00437 NOX::Abstract::Group::ReturnType
00438 LOCA::Thyra::Group::applyShiftedMatrixInverseMultiVector(
00439 Teuchos::ParameterList& lsParams,
00440 const NOX::Abstract::MultiVector& input,
00441 NOX::Abstract::MultiVector& result) const
00442 {
00443 return this->applyJacobianInverseMultiVector(lsParams, input, result);
00444 }
00445
00446 void
00447 LOCA::Thyra::Group::setSaveDataStrategy(
00448 const Teuchos::RCP<LOCA::Thyra::SaveDataStrategy>& s)
00449 {
00450 saveDataStrategy = s;
00451 }