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_PhippsBordering.H"
00043 #include "LOCA_TurningPoint_MooreSpence_ExtendedGroup.H"
00044 #include "LOCA_TurningPoint_MooreSpence_AbstractGroup.H"
00045 #include "LOCA_BorderedSolver_AbstractStrategy.H"
00046 #include "LOCA_GlobalData.H"
00047 #include "LOCA_Factory.H"
00048 #include "LOCA_ErrorCheck.H"
00049 #include "Teuchos_LAPACK.hpp"
00050 #include "LOCA_BorderedSolver_JacobianOperator.H"
00051
00052 LOCA::TurningPoint::MooreSpence::PhippsBordering::PhippsBordering(
00053 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00054 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00055 const Teuchos::RCP<Teuchos::ParameterList>& slvrParams) :
00056 globalData(global_data),
00057 solverParams(slvrParams),
00058 group(),
00059 tpGroup(),
00060 nullVector(),
00061 JnVector(),
00062 dfdp(),
00063 dJndp(),
00064 borderedSolver(),
00065 nullMultiVector(),
00066 JnMultiVector(),
00067 s(0.0)
00068 {
00069 borderedSolver =
00070 globalData->locaFactory->createBorderedSolverStrategy(topParams,
00071 solverParams);
00072 }
00073
00074 LOCA::TurningPoint::MooreSpence::PhippsBordering::~PhippsBordering()
00075 {
00076 }
00077
00078 void
00079 LOCA::TurningPoint::MooreSpence::PhippsBordering::setBlocks(
00080 const Teuchos::RCP<LOCA::TurningPoint::MooreSpence::AbstractGroup>& group_,
00081 const Teuchos::RCP<LOCA::TurningPoint::MooreSpence::ExtendedGroup>& tpGroup_,
00082 const Teuchos::RCP<const NOX::Abstract::Vector>& nullVector_,
00083 const Teuchos::RCP<const NOX::Abstract::Vector>& JnVector_,
00084 const Teuchos::RCP<const NOX::Abstract::MultiVector>& dfdp_,
00085 const Teuchos::RCP<const NOX::Abstract::MultiVector>& dJndp_)
00086 {
00087 group = group_;
00088 tpGroup = tpGroup_;
00089 nullVector = nullVector_;
00090 JnVector = JnVector_;
00091 dfdp = dfdp_;
00092 dJndp = dJndp_;
00093
00094
00095 nullMultiVector = nullVector->createMultiVector(1, NOX::DeepCopy);
00096 JnMultiVector = JnVector->createMultiVector(1, NOX::DeepCopy);
00097 s = JnVector->norm(NOX::Abstract::Vector::TwoNorm);
00098 JnMultiVector->scale(1.0/s);
00099
00100
00101 Teuchos::RCP<const LOCA::BorderedSolver::JacobianOperator> op =
00102 Teuchos::rcp(new LOCA::BorderedSolver::JacobianOperator(group));
00103 borderedSolver->setMatrixBlocksMultiVecConstraint(op,
00104 JnMultiVector,
00105 nullMultiVector,
00106 Teuchos::null);
00107 NOX::Abstract::Group::ReturnType status = borderedSolver->initForSolve();
00108 globalData->locaErrorCheck->checkReturnType(status,
00109 "LOCA::Pitchfork::MooreSpence::PhippsBordering::setBlocks()");
00110 }
00111
00112 NOX::Abstract::Group::ReturnType
00113 LOCA::TurningPoint::MooreSpence::PhippsBordering::solve(
00114 Teuchos::ParameterList& params,
00115 const LOCA::TurningPoint::MooreSpence::ExtendedMultiVector& input,
00116 LOCA::TurningPoint::MooreSpence::ExtendedMultiVector& result) const
00117 {
00118 string callingFunction =
00119 "LOCA::TurningPoint::MooreSpence::PhippsBordering::solve()";
00120 NOX::Abstract::Group::ReturnType status;
00121
00122
00123 Teuchos::RCP<const NOX::Abstract::MultiVector> input_x =
00124 input.getXMultiVec();
00125 Teuchos::RCP<const NOX::Abstract::MultiVector> input_null =
00126 input.getNullMultiVec();
00127 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_param = input.getScalars();
00128
00129
00130 Teuchos::RCP<NOX::Abstract::MultiVector> result_x =
00131 result.getXMultiVec();
00132 Teuchos::RCP<NOX::Abstract::MultiVector> result_null =
00133 result.getNullMultiVec();
00134 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_param =
00135 result.getScalars();
00136
00137 int m = input.numVectors();
00138 vector<int> index_input(m);
00139 for (int i=0; i<m; i++)
00140 index_input[i] = i;
00141
00142
00143
00144
00145
00146 Teuchos::RCP<NOX::Abstract::MultiVector> cont_input_x =
00147 input_x->clone(m+2);
00148 Teuchos::RCP<NOX::Abstract::MultiVector> cont_input_null =
00149 input_null->clone(m+2);
00150
00151 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_x =
00152 result_x->clone(m+2);
00153 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_null =
00154 result_null->clone(m+2);
00155
00156
00157 cont_input_x->setBlock(*input_x, index_input);
00158
00159
00160 (*cont_input_x)[m] = (*dfdp)[0];
00161
00162
00163 (*cont_input_x)[m+1].init(0.0);
00164
00165
00166 cont_input_null->setBlock(*input_null, index_input);
00167
00168
00169 (*cont_input_null)[m] = (*dJndp)[0];
00170
00171
00172 (*cont_input_null)[m+1].init(0.0);
00173
00174
00175 cont_result_x->init(0.0);
00176 cont_result_null->init(0.0);
00177
00178
00179 status = solveContiguous(params, *cont_input_x, *cont_input_null,
00180 *input_param, *cont_result_x, *cont_result_null,
00181 *result_param);
00182
00183
00184 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_x_view =
00185 cont_result_x->subView(index_input);
00186 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_null_view =
00187 cont_result_null->subView(index_input);
00188
00189
00190 *result_x = *cont_result_x_view;
00191 *result_null = *cont_result_null_view;
00192
00193 return status;
00194 }
00195
00196
00197
00198
00199
00200
00201 NOX::Abstract::Group::ReturnType
00202 LOCA::TurningPoint::MooreSpence::PhippsBordering::solveContiguous(
00203 Teuchos::ParameterList& params,
00204 const NOX::Abstract::MultiVector& input_x,
00205 const NOX::Abstract::MultiVector& input_null,
00206 const NOX::Abstract::MultiVector::DenseMatrix& input_param,
00207 NOX::Abstract::MultiVector& result_x,
00208 NOX::Abstract::MultiVector& result_null,
00209 NOX::Abstract::MultiVector::DenseMatrix& result_param) const
00210 {
00211 string callingFunction =
00212 "LOCA::TurningPoint::MooreSpence::PhippsBordering::solveContiguous()";
00213 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00214 NOX::Abstract::Group::ReturnType status;
00215
00216 int m = input_x.numVectors()-2;
00217 vector<int> index_input(m);
00218 vector<int> index_input_dp(m+1);
00219 vector<int> index_null(1);
00220 vector<int> index_dp(1);
00221 for (int i=0; i<m; i++) {
00222 index_input[i] = i;
00223 index_input_dp[i] = i;
00224 }
00225 index_input_dp[m] = m;
00226 index_dp[0] = m;
00227 index_null[0] = m+1;
00228
00229 NOX::Abstract::MultiVector::DenseMatrix tmp_mat_1(1, m+1);
00230 NOX::Abstract::MultiVector::DenseMatrix tmp_mat_2(1, m+2);
00231
00232
00233 Teuchos::RCP<NOX::Abstract::MultiVector> input_x_view =
00234 input_x.subView(index_input_dp);
00235 Teuchos::RCP<NOX::Abstract::MultiVector> result_x_view =
00236 result_x.subView(index_input_dp);
00237
00238
00239 if (!group->isJacobian()) {
00240 status = group->computeJacobian();
00241 finalStatus =
00242 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00243 finalStatus,
00244 callingFunction);
00245 }
00246
00247
00248
00249 status = borderedSolver->applyInverse(params, input_x_view.get(), NULL,
00250 *result_x_view, tmp_mat_1);
00251 finalStatus =
00252 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00253 callingFunction);
00254 Teuchos::RCP<NOX::Abstract::MultiVector> A =
00255 result_x.subView(index_input);
00256 Teuchos::RCP<NOX::Abstract::MultiVector> B =
00257 result_x.subView(index_dp);
00258 double b = tmp_mat_1(0,m);
00259
00260
00261 result_x[m+1] = *nullVector;
00262 Teuchos::RCP<NOX::Abstract::MultiVector> tmp =
00263 result_x.clone(NOX::ShapeCopy);
00264 status = group->computeDJnDxaMulti(*nullVector, *JnVector, result_x,
00265 *tmp);
00266 finalStatus =
00267 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00268 callingFunction);
00269
00270
00271 tmp->update(-1.0, input_null, 1.0);
00272
00273
00274 if (!group->isJacobian()) {
00275 status = group->computeJacobian();
00276 finalStatus =
00277 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00278 finalStatus,
00279 callingFunction);
00280 }
00281
00282
00283
00284 status = borderedSolver->applyInverse(params, tmp.get(), NULL, result_null,
00285 tmp_mat_2);
00286 finalStatus =
00287 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00288 callingFunction);
00289 Teuchos::RCP<NOX::Abstract::MultiVector> C =
00290 result_null.subView(index_input);
00291 Teuchos::RCP<NOX::Abstract::MultiVector> D =
00292 result_null.subView(index_dp);
00293 Teuchos::RCP<NOX::Abstract::MultiVector> E =
00294 result_null.subView(index_null);
00295 double d = tmp_mat_2(0, m);
00296 double e = tmp_mat_2(0, m+1);
00297
00298
00299 double M[9];
00300 M[0] = s; M[1] = e; M[2] = -tpGroup->lTransNorm((*E)[0]);
00301 M[3] = 0.0; M[4] = s; M[5] = tpGroup->lTransNorm(*nullVector);
00302 M[6] = b; M[7] = -d; M[8] = tpGroup->lTransNorm((*D)[0]);
00303
00304
00305 tpGroup->lTransNorm(*C, result_param);
00306 result_param += input_param;
00307
00308 double *R = new double[3*m];
00309 for (int i=0; i<m; i++) {
00310 R[3*i] = tmp_mat_1(0,i);
00311 R[3*i+1] = -tmp_mat_2(0,i);
00312 R[3*i+2] = result_param(0,i);
00313 }
00314
00315
00316 int three = 3;
00317 int piv[3];
00318 int info;
00319 Teuchos::LAPACK<int,double> L;
00320 L.GESV(three, m, M, three, piv, R, three, &info);
00321 if (info != 0) {
00322 globalData->locaErrorCheck->throwError(
00323 callingFunction,
00324 "Solve of 3x3 coefficient matrix failed!");
00325 return NOX::Abstract::Group::Failed;
00326 }
00327
00328 NOX::Abstract::MultiVector::DenseMatrix alpha(1,m);
00329 NOX::Abstract::MultiVector::DenseMatrix beta(1,m);
00330 for (int i=0; i<m; i++) {
00331 alpha(0,i) = R[3*i];
00332 beta(0,i) = R[3*i+1];
00333 result_param(0,i) = R[3*i+2];
00334 }
00335
00336
00337 A->update(Teuchos::NO_TRANS, -1.0, *B, result_param, 1.0);
00338 A->update(Teuchos::NO_TRANS, 1.0, *nullMultiVector, alpha, 1.0);
00339
00340
00341
00342 C->update(Teuchos::NO_TRANS, 1.0, *D, result_param, -1.0);
00343 C->update(Teuchos::NO_TRANS, -1.0, *E, alpha, 1.0);
00344 C->update(Teuchos::NO_TRANS, 1.0, *nullMultiVector, beta, 1.0);
00345
00346 delete [] R;
00347
00348 return finalStatus;
00349 }