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_Hopf_MooreSpence_SalingerBordering.H"
00043 #include "LOCA_Hopf_MooreSpence_ExtendedGroup.H"
00044 #include "LOCA_Hopf_MooreSpence_AbstractGroup.H"
00045 #include "LOCA_GlobalData.H"
00046 #include "LOCA_ErrorCheck.H"
00047
00048 LOCA::Hopf::MooreSpence::SalingerBordering::SalingerBordering(
00049 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00050 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00051 const Teuchos::RCP<Teuchos::ParameterList>& slvrParams) :
00052 globalData(global_data),
00053 solverParams(slvrParams),
00054 group(),
00055 hopfGroup(),
00056 yVector(),
00057 zVector(),
00058 CeRealVector(),
00059 CeImagVector(),
00060 dfdp(),
00061 dCedpReal(),
00062 dCedpImag(),
00063 ByVector(),
00064 minusBzVector()
00065 {
00066 }
00067
00068 LOCA::Hopf::MooreSpence::SalingerBordering::~SalingerBordering()
00069 {
00070 }
00071
00072 void
00073 LOCA::Hopf::MooreSpence::SalingerBordering::setBlocks(
00074 const Teuchos::RCP<LOCA::Hopf::MooreSpence::AbstractGroup>& group_,
00075 const Teuchos::RCP<LOCA::Hopf::MooreSpence::ExtendedGroup>& hopfGroup_,
00076 const Teuchos::RCP<const NOX::Abstract::Vector>& yVector_,
00077 const Teuchos::RCP<const NOX::Abstract::Vector>& zVector_,
00078 const Teuchos::RCP<const NOX::Abstract::Vector>& CeRealVector_,
00079 const Teuchos::RCP<const NOX::Abstract::Vector>& CeImagVector_,
00080 const Teuchos::RCP<const NOX::Abstract::Vector>& dfdp_,
00081 const Teuchos::RCP<const NOX::Abstract::Vector>& dCedpReal_,
00082 const Teuchos::RCP<const NOX::Abstract::Vector>& dCedpImag_,
00083 const Teuchos::RCP<const NOX::Abstract::Vector>& ByVector_,
00084 const Teuchos::RCP<const NOX::Abstract::Vector>& mBzVector_,
00085 double w_)
00086 {
00087 group = group_;
00088 hopfGroup = hopfGroup_;
00089 yVector = yVector_;
00090 zVector = zVector_;
00091 CeRealVector = CeRealVector_;
00092 CeImagVector = CeImagVector_;
00093 dfdp = dfdp_;
00094 dCedpReal = dCedpReal_;
00095 dCedpImag = dCedpImag_;
00096 ByVector = ByVector_;
00097 minusBzVector = mBzVector_;
00098 w = w_;
00099 }
00100
00101 NOX::Abstract::Group::ReturnType
00102 LOCA::Hopf::MooreSpence::SalingerBordering::solve(
00103 Teuchos::ParameterList& params,
00104 const LOCA::Hopf::MooreSpence::ExtendedMultiVector& input,
00105 LOCA::Hopf::MooreSpence::ExtendedMultiVector& result) const
00106 {
00107 string callingFunction =
00108 "LOCA::Hopf::MooreSpence::SalingerBordering::solve()";
00109 NOX::Abstract::Group::ReturnType status;
00110
00111
00112 Teuchos::RCP<const NOX::Abstract::MultiVector> input_x =
00113 input.getXMultiVec();
00114 Teuchos::RCP<const NOX::Abstract::MultiVector> input_y =
00115 input.getRealEigenMultiVec();
00116 Teuchos::RCP<const NOX::Abstract::MultiVector> input_z =
00117 input.getImagEigenMultiVec();
00118 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_w
00119 = input.getFrequencies();
00120 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_p
00121 = input.getBifParams();
00122
00123
00124 Teuchos::RCP<NOX::Abstract::MultiVector> result_x =
00125 result.getXMultiVec();
00126 Teuchos::RCP<NOX::Abstract::MultiVector> result_y =
00127 result.getRealEigenMultiVec();
00128 Teuchos::RCP<NOX::Abstract::MultiVector> result_z =
00129 result.getImagEigenMultiVec();
00130 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_w =
00131 result.getFrequencies();
00132 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_p =
00133 result.getBifParams();
00134
00135 int m = input.numVectors();
00136
00137 vector<int> index_input(m);
00138 for (int i=0; i<m; i++)
00139 index_input[i] = i;
00140
00141
00142
00143
00144
00145
00146 Teuchos::RCP<NOX::Abstract::MultiVector> cont_input_x =
00147 input_x->clone(m+1);
00148 Teuchos::RCP<NOX::Abstract::MultiVector> cont_input_y =
00149 input_y->clone(m+2);
00150 Teuchos::RCP<NOX::Abstract::MultiVector> cont_input_z =
00151 input_z->clone(m+2);
00152
00153 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_x =
00154 result_x->clone(m+1);
00155 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_y =
00156 result_y->clone(m+2);
00157 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_z =
00158 result_y->clone(m+2);
00159
00160
00161 cont_input_x->setBlock(*input_x, index_input);
00162
00163
00164 (*cont_input_x)[m] = *dfdp;
00165
00166
00167 cont_input_y->setBlock(*input_y, index_input);
00168
00169
00170 (*cont_input_y)[m] = *dCedpReal;
00171
00172
00173 (*cont_input_y)[m+1] = *minusBzVector;
00174
00175
00176 cont_input_z->setBlock(*input_z, index_input);
00177
00178
00179 (*cont_input_z)[m] = *dCedpImag;
00180
00181
00182 (*cont_input_z)[m+1] = *ByVector;
00183
00184
00185 cont_result_x->init(0.0);
00186 cont_result_y->init(0.0);
00187 cont_result_z->init(0.0);
00188
00189
00190 status = solveContiguous(params, *cont_input_x, *cont_input_y, *cont_input_z,
00191 *input_w, *input_p, *cont_result_x, *cont_result_y,
00192 *cont_result_z, *result_w, *result_p);
00193
00194
00195 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_x_view =
00196 cont_result_x->subView(index_input);
00197 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_y_view =
00198 cont_result_y->subView(index_input);
00199 Teuchos::RCP<NOX::Abstract::MultiVector> cont_result_z_view =
00200 cont_result_z->subView(index_input);
00201
00202
00203 *result_x = *cont_result_x_view;
00204 *result_y = *cont_result_y_view;
00205 *result_z = *cont_result_z_view;
00206
00207 return status;
00208 }
00209
00210
00211
00212
00213
00214
00215
00216
00217 NOX::Abstract::Group::ReturnType
00218 LOCA::Hopf::MooreSpence::SalingerBordering::solveContiguous(
00219 Teuchos::ParameterList& params,
00220 const NOX::Abstract::MultiVector& input_x,
00221 const NOX::Abstract::MultiVector& input_y,
00222 const NOX::Abstract::MultiVector& input_z,
00223 const NOX::Abstract::MultiVector::DenseMatrix& input_w,
00224 const NOX::Abstract::MultiVector::DenseMatrix& input_p,
00225 NOX::Abstract::MultiVector& result_x,
00226 NOX::Abstract::MultiVector& result_y,
00227 NOX::Abstract::MultiVector& result_z,
00228 NOX::Abstract::MultiVector::DenseMatrix& result_w,
00229 NOX::Abstract::MultiVector::DenseMatrix& result_p) const
00230 {
00231 string callingFunction =
00232 "LOCA::Hopf::MooreSpence::SalingerBordering::solveContiguous()";
00233 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00234 NOX::Abstract::Group::ReturnType status;
00235
00236 int m = input_x.numVectors()-1;
00237 vector<int> index_input(m);
00238 vector<int> index_dp(1);
00239 vector<int> index_B(1);
00240 vector<int> index_ip(m+1);
00241 for (int i=0; i<m; i++) {
00242 index_input[i] = i;
00243 index_ip[i] = i;
00244 }
00245 index_ip[m] = m;
00246 index_dp[0] = m;
00247 index_B[0] = m+1;
00248
00249
00250 if (!group->isJacobian()) {
00251 status = group->computeJacobian();
00252 finalStatus =
00253 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00254 finalStatus,
00255 callingFunction);
00256 }
00257
00258
00259 status = group->applyJacobianInverseMultiVector(params, input_x, result_x);
00260 finalStatus =
00261 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00262 callingFunction);
00263 Teuchos::RCP<NOX::Abstract::MultiVector> A =
00264 result_x.subView(index_input);
00265 Teuchos::RCP<NOX::Abstract::MultiVector> b =
00266 result_x.subView(index_dp);
00267
00268
00269 if (!group->isComplex()) {
00270 status = group->computeComplex(w);
00271 finalStatus =
00272 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00273 finalStatus,
00274 callingFunction);
00275 }
00276
00277
00278 Teuchos::RCP<NOX::Abstract::MultiVector> tmp_real =
00279 result_y.clone(NOX::ShapeCopy);
00280 Teuchos::RCP<NOX::Abstract::MultiVector> tmp_real_sub =
00281 tmp_real->subView(index_ip);
00282 Teuchos::RCP<NOX::Abstract::MultiVector> tmp_imag =
00283 result_y.clone(NOX::ShapeCopy);
00284 Teuchos::RCP<NOX::Abstract::MultiVector> tmp_imag_sub =
00285 tmp_imag->subView(index_ip);
00286 tmp_real->init(0.0);
00287 tmp_imag->init(0.0);
00288 status = group->computeDCeDxa(*yVector, *zVector, w, result_x,
00289 *CeRealVector, *CeImagVector, *tmp_real_sub,
00290 *tmp_imag_sub);
00291 finalStatus =
00292 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00293 callingFunction);
00294
00295
00296 tmp_real->update(1.0, input_y, -1.0);
00297 tmp_imag->update(1.0, input_z, -1.0);
00298
00299
00300 if (!group->isComplex()) {
00301 status = group->computeComplex(w);
00302 finalStatus =
00303 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00304 finalStatus,
00305 callingFunction);
00306 }
00307
00308
00309 status = group->applyComplexInverseMultiVector(params, *tmp_real, *tmp_imag,
00310 result_y, result_z);
00311 finalStatus =
00312 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00313 callingFunction);
00314 Teuchos::RCP<NOX::Abstract::MultiVector> C =
00315 result_y.subView(index_input);
00316 Teuchos::RCP<NOX::Abstract::MultiVector> D =
00317 result_z.subView(index_input);
00318 Teuchos::RCP<NOX::Abstract::MultiVector> e =
00319 result_y.subView(index_dp);
00320 Teuchos::RCP<NOX::Abstract::MultiVector> f =
00321 result_z.subView(index_dp);
00322 Teuchos::RCP<NOX::Abstract::MultiVector> g =
00323 result_y.subView(index_B);
00324 Teuchos::RCP<NOX::Abstract::MultiVector> h =
00325 result_z.subView(index_B);
00326
00327
00328
00329 NOX::Abstract::MultiVector::DenseMatrix ltC(1,m);
00330 NOX::Abstract::MultiVector::DenseMatrix ltD(1,m);
00331 double lte = hopfGroup->lTransNorm((*e)[0]);
00332 double ltf = hopfGroup->lTransNorm((*f)[0]);
00333 double ltg = hopfGroup->lTransNorm((*g)[0]);
00334 double lth = hopfGroup->lTransNorm((*h)[0]);
00335 double denom = lth*lte - ltg*ltf;
00336 hopfGroup->lTransNorm(*C, ltC);
00337 ltC -= input_w;
00338 ltC.scale(lth);
00339 hopfGroup->lTransNorm(*D, ltD);
00340 ltD -= input_p;
00341 result_p.assign(ltD);
00342 result_p.scale(-ltg);
00343 result_p += ltC;
00344 result_p.scale(1.0/denom);
00345
00346
00347 result_w.assign(result_p);
00348 result_w.scale(-ltf);
00349 result_w += ltD;
00350 result_w.scale(1.0/lth);
00351
00352
00353 A->update(Teuchos::NO_TRANS, -1.0, *b, result_p, 1.0);
00354
00355
00356 C->update(Teuchos::NO_TRANS, -1.0, *e, result_p, 1.0);
00357 C->update(Teuchos::NO_TRANS, -1.0, *g, result_w, 1.0);
00358
00359
00360 D->update(Teuchos::NO_TRANS, -1.0, *f, result_p, 1.0);
00361 D->update(Teuchos::NO_TRANS, -1.0, *h, result_w, 1.0);
00362
00363 return finalStatus;
00364 }