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_TurningPoint_MooreSpence_SalingerBordering.H"
00043 #include "LOCA_TurningPoint_MooreSpence_ExtendedGroup.H"
00044 #include "LOCA_TurningPoint_MooreSpence_AbstractGroup.H"
00045 #include "LOCA_GlobalData.H"
00046 #include "LOCA_ErrorCheck.H"
00047 #include "LOCA_Abstract_TransposeSolveGroup.H"
00048
00049 LOCA::TurningPoint::MooreSpence::SalingerBordering::SalingerBordering(
00050 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00051 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00052 const Teuchos::RCP<Teuchos::ParameterList>& slvrParams) :
00053 globalData(global_data),
00054 solverParams(slvrParams),
00055 group(),
00056 tpGroup(),
00057 nullVector(),
00058 JnVector(),
00059 dfdp(),
00060 dJndp()
00061 {
00062 }
00063
00064 LOCA::TurningPoint::MooreSpence::SalingerBordering::~SalingerBordering()
00065 {
00066 }
00067
00068 void
00069 LOCA::TurningPoint::MooreSpence::SalingerBordering::setBlocks(
00070 const Teuchos::RCP<LOCA::TurningPoint::MooreSpence::AbstractGroup>& group_,
00071 const Teuchos::RCP<LOCA::TurningPoint::MooreSpence::ExtendedGroup>& tpGroup_,
00072 const Teuchos::RCP<const NOX::Abstract::Vector>& nullVector_,
00073 const Teuchos::RCP<const NOX::Abstract::Vector>& JnVector_,
00074 const Teuchos::RCP<const NOX::Abstract::MultiVector>& dfdp_,
00075 const Teuchos::RCP<const NOX::Abstract::MultiVector>& dJndp_)
00076 {
00077 group = group_;
00078 tpGroup = tpGroup_;
00079 nullVector = nullVector_;
00080 JnVector = JnVector_;
00081 dfdp = dfdp_;
00082 dJndp = dJndp_;
00083 }
00084
00085 NOX::Abstract::Group::ReturnType
00086 LOCA::TurningPoint::MooreSpence::SalingerBordering::solve(
00087 Teuchos::ParameterList& params,
00088 const LOCA::TurningPoint::MooreSpence::ExtendedMultiVector& input,
00089 LOCA::TurningPoint::MooreSpence::ExtendedMultiVector& result) const
00090 {
00091 string callingFunction =
00092 "LOCA::TurningPoint::MooreSpence::SalingerBordering::solve()";
00093 NOX::Abstract::Group::ReturnType status;
00094
00095
00096 Teuchos::RCP<const NOX::Abstract::MultiVector> input_x =
00097 input.getXMultiVec();
00098 Teuchos::RCP<const NOX::Abstract::MultiVector> input_null =
00099 input.getNullMultiVec();
00100 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_param = input.getScalars();
00101
00102
00103 Teuchos::RCP<NOX::Abstract::MultiVector> result_x =
00104 result.getXMultiVec();
00105 Teuchos::RCP<NOX::Abstract::MultiVector> result_null =
00106 result.getNullMultiVec();
00107 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_param =
00108 result.getScalars();
00109
00110 int m = input.numVectors();
00111
00112 vector<int> index_input(m);
00113 for (int i=0; i<m; i++)
00114 index_input[i] = i;
00115
00116
00117
00118
00119
00120 Teuchos::RCP<NOX::Abstract::MultiVector> cont_input_x =
00121 input_x->clone(m+1);
00122 Teuchos::RCP<NOX::Abstract::MultiVector> cont_input_null =
00123 input_null->clone(m+1);
00124
00125 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_x =
00126 result_x->clone(m+1);
00127 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_null =
00128 result_null->clone(m+1);
00129
00130
00131 cont_input_x->setBlock(*input_x, index_input);
00132
00133
00134 (*cont_input_x)[m] = (*dfdp)[0];
00135
00136
00137 cont_input_null->setBlock(*input_null, index_input);
00138
00139
00140 (*cont_input_null)[m] = (*dJndp)[0];
00141
00142
00143 cont_result_x->init(0.0);
00144 cont_result_null->init(0.0);
00145
00146
00147 status = solveContiguous(params, *cont_input_x, *cont_input_null,
00148 *input_param, *cont_result_x, *cont_result_null,
00149 *result_param);
00150
00151
00152 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_x_view =
00153 cont_result_x->subView(index_input);
00154 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_null_view =
00155 cont_result_null->subView(index_input);
00156
00157
00158 *result_x = *cont_result_x_view;
00159 *result_null = *cont_result_null_view;
00160
00161 return status;
00162 }
00163
00164 NOX::Abstract::Group::ReturnType
00165 LOCA::TurningPoint::MooreSpence::SalingerBordering::solveTranspose(
00166 Teuchos::ParameterList& params,
00167 const LOCA::TurningPoint::MooreSpence::ExtendedMultiVector& input,
00168 LOCA::TurningPoint::MooreSpence::ExtendedMultiVector& result) const
00169 {
00170 string callingFunction =
00171 "LOCA::TurningPoint::MooreSpence::SalingerBordering::solveTranspose()";
00172 NOX::Abstract::Group::ReturnType status;
00173
00174
00175 Teuchos::RCP<const NOX::Abstract::MultiVector> input_x =
00176 input.getXMultiVec();
00177 Teuchos::RCP<const NOX::Abstract::MultiVector> input_null =
00178 input.getNullMultiVec();
00179 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_param = input.getScalars();
00180
00181
00182 Teuchos::RCP<NOX::Abstract::MultiVector> result_x =
00183 result.getXMultiVec();
00184 Teuchos::RCP<NOX::Abstract::MultiVector> result_null =
00185 result.getNullMultiVec();
00186 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_param =
00187 result.getScalars();
00188
00189 int m = input.numVectors();
00190
00191 vector<int> index_input(m);
00192 for (int i=0; i<m; i++)
00193 index_input[i] = i;
00194
00195
00196
00197
00198
00199 Teuchos::RCP<NOX::Abstract::MultiVector> cont_input_x =
00200 input_x->clone(m+1);
00201 Teuchos::RCP<NOX::Abstract::MultiVector> cont_input_null =
00202 input_null->clone(m+1);
00203
00204 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_x =
00205 result_x->clone(m+1);
00206 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_null =
00207 result_null->clone(m+1);
00208
00209
00210 cont_input_x->setBlock(*input_x, index_input);
00211
00212
00213 (*cont_input_x)[m].init(0);
00214
00215
00216 cont_input_null->setBlock(*input_null, index_input);
00217
00218
00219 Teuchos::RCP<NOX::Abstract::Vector> phi = tpGroup->getLengthVector();
00220 (*cont_input_null)[m].update(-1.0, *phi, 0.0);
00221
00222
00223 cont_result_x->init(0.0);
00224 cont_result_null->init(0.0);
00225
00226
00227 status = solveTransposeContiguous(params, *cont_input_x, *cont_input_null,
00228 *input_param, *cont_result_x,
00229 *cont_result_null,
00230 *result_param);
00231
00232
00233 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_x_view =
00234 cont_result_x->subView(index_input);
00235 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_null_view =
00236 cont_result_null->subView(index_input);
00237
00238
00239 *result_x = *cont_result_x_view;
00240 *result_null = *cont_result_null_view;
00241
00242 return status;
00243 }
00244
00245
00246
00247
00248
00249
00250 NOX::Abstract::Group::ReturnType
00251 LOCA::TurningPoint::MooreSpence::SalingerBordering::solveContiguous(
00252 Teuchos::ParameterList& params,
00253 const NOX::Abstract::MultiVector& input_x,
00254 const NOX::Abstract::MultiVector& input_null,
00255 const NOX::Abstract::MultiVector::DenseMatrix& input_param,
00256 NOX::Abstract::MultiVector& result_x,
00257 NOX::Abstract::MultiVector& result_null,
00258 NOX::Abstract::MultiVector::DenseMatrix& result_param) const
00259 {
00260 string callingFunction =
00261 "LOCA::TurningPoint::MooreSpence::SalingerBordering::solveContiguous()";
00262 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00263 NOX::Abstract::Group::ReturnType status;
00264
00265 int m = input_x.numVectors()-1;
00266 vector<int> index_input(m);
00267 vector<int> index_dp(1);
00268 for (int i=0; i<m; i++)
00269 index_input[i] = i;
00270 index_dp[0] = m;
00271
00272
00273 if (!group->isJacobian()) {
00274 status = group->computeJacobian();
00275 finalStatus =
00276 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00277 finalStatus,
00278 callingFunction);
00279 }
00280
00281
00282 status = group->applyJacobianInverseMultiVector(params, input_x, result_x);
00283 finalStatus =
00284 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00285 callingFunction);
00286 Teuchos::RCP<NOX::Abstract::MultiVector> A =
00287 result_x.subView(index_input);
00288 Teuchos::RCP<NOX::Abstract::MultiVector> b =
00289 result_x.subView(index_dp);
00290
00291
00292 Teuchos::RCP<NOX::Abstract::MultiVector> tmp =
00293 result_x.clone(NOX::ShapeCopy);
00294 status = group->computeDJnDxaMulti(*nullVector, *JnVector, result_x,
00295 *tmp);
00296 finalStatus =
00297 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00298 callingFunction);
00299
00300
00301 tmp->update(-1.0, input_null, 1.0);
00302
00303
00304 if (!group->isJacobian()) {
00305 status = group->computeJacobian();
00306 finalStatus =
00307 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00308 finalStatus,
00309 callingFunction);
00310 }
00311
00312
00313 status = group->applyJacobianInverseMultiVector(params, *tmp, result_null);
00314 finalStatus =
00315 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00316 callingFunction);
00317 Teuchos::RCP<NOX::Abstract::MultiVector> C =
00318 result_null.subView(index_input);
00319 Teuchos::RCP<NOX::Abstract::MultiVector> d =
00320 result_null.subView(index_dp);
00321
00322
00323 tpGroup->lTransNorm(*C, result_param);
00324 result_param += input_param;
00325 double denom = tpGroup->lTransNorm((*d)[0]);
00326 result_param.scale(1.0/denom);
00327
00328
00329 A->update(Teuchos::NO_TRANS, -1.0, *b, result_param, 1.0);
00330
00331
00332 C->update(Teuchos::NO_TRANS, 1.0, *d, result_param, -1.0);
00333
00334 return finalStatus;
00335 }
00336
00337
00338
00339
00340
00341
00342 NOX::Abstract::Group::ReturnType
00343 LOCA::TurningPoint::MooreSpence::SalingerBordering::solveTransposeContiguous(
00344 Teuchos::ParameterList& params,
00345 const NOX::Abstract::MultiVector& input_x,
00346 const NOX::Abstract::MultiVector& input_null,
00347 const NOX::Abstract::MultiVector::DenseMatrix& input_param,
00348 NOX::Abstract::MultiVector& result_x,
00349 NOX::Abstract::MultiVector& result_null,
00350 NOX::Abstract::MultiVector::DenseMatrix& result_param) const
00351 {
00352 string callingFunction =
00353 "LOCA::TurningPoint::MooreSpence::SalingerBordering::solveTransposeContiguous()";
00354 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00355 NOX::Abstract::Group::ReturnType status;
00356
00357
00358 Teuchos::RCP<LOCA::Abstract::TransposeSolveGroup> tsGroup =
00359 Teuchos::rcp_dynamic_cast<LOCA::Abstract::TransposeSolveGroup>(group);
00360 if (tsGroup == Teuchos::null)
00361 globalData->locaErrorCheck->throwError(callingFunction,
00362 "Underlying group must be derived from NOX::Abstract::TransposeSolveGroup for transpose solve");
00363
00364 int m = input_x.numVectors()-1;
00365 vector<int> index_input(m);
00366 vector<int> index_dp(1);
00367 for (int i=0; i<m; i++)
00368 index_input[i] = i;
00369 index_dp[0] = m;
00370
00371
00372 if (!group->isJacobian()) {
00373 status = group->computeJacobian();
00374 finalStatus =
00375 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00376 finalStatus,
00377 callingFunction);
00378 }
00379
00380
00381 status = tsGroup->applyJacobianTransposeInverseMultiVector(params,
00382 input_null,
00383 result_null);
00384 finalStatus =
00385 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00386 callingFunction);
00387 Teuchos::RCP<NOX::Abstract::MultiVector> A =
00388 result_null.subView(index_input);
00389 Teuchos::RCP<NOX::Abstract::MultiVector> b =
00390 result_null.subView(index_dp);
00391
00392
00393 Teuchos::RCP<NOX::Abstract::MultiVector> tmp =
00394 result_null.clone(NOX::ShapeCopy);
00395 status = group->computeDwtJnDxMulti(result_null, *nullVector, *tmp);
00396 finalStatus =
00397 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00398 finalStatus,
00399 callingFunction);
00400
00401
00402 tmp->update(1.0, input_x, -1.0);
00403
00404
00405 if (!group->isJacobian()) {
00406 status = group->computeJacobian();
00407 finalStatus =
00408 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00409 finalStatus,
00410 callingFunction);
00411 }
00412
00413
00414 status = tsGroup->applyJacobianTransposeInverseMultiVector(params, *tmp,
00415 result_x);
00416 finalStatus =
00417 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00418 callingFunction);
00419 Teuchos::RCP<NOX::Abstract::MultiVector> C =
00420 result_x.subView(index_input);
00421 Teuchos::RCP<NOX::Abstract::MultiVector> d =
00422 result_x.subView(index_dp);
00423
00424
00425 NOX::Abstract::MultiVector::DenseMatrix t1(1,m+1);
00426 result_null.multiply(1.0, *dJndp, t1);
00427
00428
00429 NOX::Abstract::MultiVector::DenseMatrix t2(1,m+1);
00430 result_x.multiply(1.0, *dfdp, t2);
00431
00432
00433 double denom = t2(0,m) + t1(0,m);
00434 for (int i=0; i<m; i++)
00435 result_param(0,i) = (input_param(0,i) - t2(0,i) - t1(0,i))/denom;
00436
00437
00438 A->update(Teuchos::NO_TRANS, 1.0, *b, result_param, 1.0);
00439
00440
00441 C->update(Teuchos::NO_TRANS, 1.0, *d, result_param, 1.0);
00442
00443 return finalStatus;
00444 }