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 "Teuchos_BLAS.hpp"
00043 #include "LOCA_Pitchfork_MooreSpence_SalingerBordering.H"
00044 #include "LOCA_Pitchfork_MooreSpence_ExtendedGroup.H"
00045 #include "LOCA_Pitchfork_MooreSpence_AbstractGroup.H"
00046 #include "LOCA_GlobalData.H"
00047 #include "LOCA_ErrorCheck.H"
00048
00049 LOCA::Pitchfork::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 pfGroup(),
00057 asymMultiVector(),
00058 asymVector(),
00059 nullVector(),
00060 JnVector(),
00061 dfdp(),
00062 dJndp()
00063 {
00064 }
00065
00066 LOCA::Pitchfork::MooreSpence::SalingerBordering::~SalingerBordering()
00067 {
00068 }
00069
00070 void
00071 LOCA::Pitchfork::MooreSpence::SalingerBordering::setBlocks(
00072 const Teuchos::RCP<LOCA::Pitchfork::MooreSpence::AbstractGroup>& group_,
00073 const Teuchos::RCP<LOCA::Pitchfork::MooreSpence::ExtendedGroup>& pfGroup_,
00074 const Teuchos::RCP<const NOX::Abstract::MultiVector>& asymMultiVector_,
00075 const Teuchos::RCP<const NOX::Abstract::Vector>& nullVector_,
00076 const Teuchos::RCP<const NOX::Abstract::Vector>& JnVector_,
00077 const Teuchos::RCP<const NOX::Abstract::Vector>& dfdp_,
00078 const Teuchos::RCP<const NOX::Abstract::Vector>& dJndp_)
00079 {
00080 group = group_;
00081 pfGroup = pfGroup_;
00082 asymMultiVector = asymMultiVector_;
00083 asymVector = Teuchos::rcp(&(*asymMultiVector)[0], false);
00084 nullVector = nullVector_;
00085 JnVector = JnVector_;
00086 dfdp = dfdp_;
00087 dJndp = dJndp_;
00088 }
00089
00090 NOX::Abstract::Group::ReturnType
00091 LOCA::Pitchfork::MooreSpence::SalingerBordering::solve(
00092 Teuchos::ParameterList& params,
00093 const LOCA::Pitchfork::MooreSpence::ExtendedMultiVector& input,
00094 LOCA::Pitchfork::MooreSpence::ExtendedMultiVector& result) const
00095 {
00096 string callingFunction =
00097 "LOCA::Pitchfork::MooreSpence::SalingerBordering::solve()";
00098 NOX::Abstract::Group::ReturnType status;
00099
00100
00101 Teuchos::RCP<const NOX::Abstract::MultiVector> input_x =
00102 input.getXMultiVec();
00103 Teuchos::RCP<const NOX::Abstract::MultiVector> input_null =
00104 input.getNullMultiVec();
00105 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_slack = input.getSlacks();
00106 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_param = input.getBifParams();
00107
00108
00109 Teuchos::RCP<NOX::Abstract::MultiVector> result_x =
00110 result.getXMultiVec();
00111 Teuchos::RCP<NOX::Abstract::MultiVector> result_null =
00112 result.getNullMultiVec();
00113 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_slack =
00114 result.getSlacks();
00115 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_param =
00116 result.getBifParams();
00117
00118 int m = input.numVectors();
00119 vector<int> index_input(m);
00120 for (int i=0; i<m; i++)
00121 index_input[i] = i;
00122
00123
00124
00125
00126
00127 Teuchos::RCP<NOX::Abstract::MultiVector> cont_input_x =
00128 input_x->clone(m+2);
00129 Teuchos::RCP<NOX::Abstract::MultiVector> cont_input_null =
00130 input_null->clone(m+2);
00131
00132 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_x =
00133 result_x->clone(m+2);
00134 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_null =
00135 result_null->clone(m+2);
00136
00137
00138 cont_input_x->setBlock(*input_x, index_input);
00139
00140
00141 (*cont_input_x)[m] = *dfdp;
00142
00143
00144 (*cont_input_x)[m+1] = *asymVector;
00145
00146
00147 cont_input_null->setBlock(*input_null, index_input);
00148
00149
00150 (*cont_input_null)[m] = *dJndp;
00151
00152
00153 (*cont_input_null)[m+1].init(0.0);
00154
00155
00156 cont_result_x->init(0.0);
00157 cont_result_null->init(0.0);
00158
00159
00160 status = solveContiguous(params, *cont_input_x, *cont_input_null,
00161 *input_slack, *input_param,
00162 *cont_result_x, *cont_result_null,
00163 *result_slack, *result_param);
00164
00165
00166 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_x_view =
00167 cont_result_x->subView(index_input);
00168 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_null_view =
00169 cont_result_null->subView(index_input);
00170
00171
00172 *result_x = *cont_result_x_view;
00173 *result_null = *cont_result_null_view;
00174
00175 return status;
00176 }
00177
00178
00179
00180
00181
00182
00183
00184 NOX::Abstract::Group::ReturnType
00185 LOCA::Pitchfork::MooreSpence::SalingerBordering::solveContiguous(
00186 Teuchos::ParameterList& params,
00187 const NOX::Abstract::MultiVector& input_x,
00188 const NOX::Abstract::MultiVector& input_null,
00189 const NOX::Abstract::MultiVector::DenseMatrix& input_slack,
00190 const NOX::Abstract::MultiVector::DenseMatrix& input_param,
00191 NOX::Abstract::MultiVector& result_x,
00192 NOX::Abstract::MultiVector& result_null,
00193 NOX::Abstract::MultiVector::DenseMatrix& result_slack,
00194 NOX::Abstract::MultiVector::DenseMatrix& result_param) const
00195 {
00196 string callingFunction =
00197 "LOCA::Pitchfork::MooreSpence::SalingerBordering::solveContiguous()";
00198 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00199 NOX::Abstract::Group::ReturnType status;
00200
00201 int m = input_x.numVectors()-2;
00202 vector<int> index_input(m);
00203 vector<int> index_dp(1);
00204 vector<int> index_s(1);
00205 for (int i=0; i<m; i++)
00206 index_input[i] = i;
00207 index_dp[0] = m;
00208 index_s[0] = m+1;
00209
00210
00211 if (!group->isJacobian()) {
00212 status = group->computeJacobian();
00213 finalStatus =
00214 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00215 finalStatus,
00216 callingFunction);
00217 }
00218
00219
00220 status = group->applyJacobianInverseMultiVector(params, input_x, result_x);
00221 finalStatus =
00222 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00223 callingFunction);
00224 Teuchos::RCP<NOX::Abstract::MultiVector> A =
00225 result_x.subView(index_input);
00226 Teuchos::RCP<NOX::Abstract::MultiVector> b =
00227 result_x.subView(index_dp);
00228 Teuchos::RCP<NOX::Abstract::MultiVector> c =
00229 result_x.subView(index_s);
00230
00231
00232 Teuchos::RCP<NOX::Abstract::MultiVector> tmp =
00233 result_x.clone(NOX::ShapeCopy);
00234 status = group->computeDJnDxaMulti(*nullVector, *JnVector, result_x,
00235 *tmp);
00236 finalStatus =
00237 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00238 callingFunction);
00239
00240
00241 tmp->update(1.0, input_null, -1.0);
00242
00243
00244 if (!group->isJacobian()) {
00245 status = group->computeJacobian();
00246 finalStatus =
00247 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00248 finalStatus,
00249 callingFunction);
00250 }
00251
00252
00253 status = group->applyJacobianInverseMultiVector(params, *tmp, result_null);
00254 finalStatus =
00255 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00256 callingFunction);
00257 Teuchos::RCP<NOX::Abstract::MultiVector> D =
00258 result_null.subView(index_input);
00259 Teuchos::RCP<NOX::Abstract::MultiVector> e =
00260 result_null.subView(index_dp);
00261 Teuchos::RCP<NOX::Abstract::MultiVector> g =
00262 result_null.subView(index_s);
00263
00264
00265
00266 Teuchos::BLAS<int,double> dblas;
00267 double lte = pfGroup->lTransNorm((*e)[0]);
00268 double ltf = pfGroup->lTransNorm((*g)[0]);
00269 double ipb = group->innerProduct((*b)[0], *asymVector);
00270 double ipc = group->innerProduct((*c)[0], *asymVector);
00271 double denom = lte*ipc - ltf*ipb;
00272 group->innerProduct(*asymMultiVector, *A, result_slack);
00273 pfGroup->lTransNorm(*D, result_param);
00274
00275 for (int i=0; i<m; i++) {
00276 result_slack(0,i) = (lte*(result_slack(0,i) - input_slack(0,i)) -
00277 ipb*(result_param(0,i) - input_param(0,i))) / denom;
00278 result_param(0,i) = (result_param(0,i) - input_param(0,i) -
00279 ltf*result_slack(0,i)) / lte;
00280 }
00281
00282
00283 A->update(Teuchos::NO_TRANS, -1.0, *b, result_param, 1.0);
00284 A->update(Teuchos::NO_TRANS, -1.0, *c, result_slack, 1.0);
00285
00286
00287 D->update(Teuchos::NO_TRANS, -1.0, *e, result_param, 1.0);
00288 D->update(Teuchos::NO_TRANS, -1.0, *g, result_slack, 1.0);
00289
00290 return finalStatus;
00291 }