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
00043 #include "LOCA_DerivUtils.H"
00044
00045 #include "LOCA_MultiContinuation_AbstractGroup.H"
00046 #include "LOCA_Hopf_MooreSpence_AbstractGroup.H"
00047 #include "LOCA_Hopf_MinimallyAugmented_AbstractGroup.H"
00048 #include "NOX_Abstract_Vector.H"
00049 #include "NOX_Abstract_MultiVector.H"
00050 #include "LOCA_Parameter_Vector.H"
00051 #include "NOX_Common.H"
00052 #include "LOCA_GlobalData.H"
00053 #include "LOCA_ErrorCheck.H"
00054
00055 LOCA::DerivUtils::DerivUtils(
00056 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00057 double perturb) :
00058 globalData(global_data),
00059 perturb(perturb)
00060 {
00061
00062 }
00063
00064 LOCA::DerivUtils::DerivUtils(const DerivUtils& source) :
00065 globalData(source.globalData),
00066 perturb(source.perturb)
00067 {
00068
00069 }
00070
00071 LOCA::DerivUtils::~DerivUtils()
00072 {
00073
00074 }
00075
00076 Teuchos::RCP<LOCA::DerivUtils>
00077 LOCA::DerivUtils::clone(NOX::CopyType type) const
00078 {
00079 return Teuchos::rcp(new DerivUtils(*this));
00080 }
00081
00082 NOX::Abstract::Group::ReturnType
00083 LOCA::DerivUtils::computeDfDp(LOCA::MultiContinuation::AbstractGroup& grp,
00084 const vector<int>& param_ids,
00085 NOX::Abstract::MultiVector& result,
00086 bool isValidF) const
00087 {
00088 string callingFunction =
00089 "LOCA::DerivUtils::computeDfDp()";
00090 NOX::Abstract::Group::ReturnType status, finalStatus;
00091
00092
00093 NOX::Abstract::Vector *f = &result[0];
00094 NOX::Abstract::Vector *dfdp = NULL;
00095
00096
00097 if (!isValidF) {
00098 finalStatus = grp.computeF();
00099 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00100 *f = grp.getF();
00101 }
00102 else
00103 finalStatus = NOX::Abstract::Group::Ok;
00104
00105 double param;
00106 double eps;
00107
00108
00109 for (unsigned int i=0; i<param_ids.size(); i++) {
00110
00111
00112 eps = perturbParam(grp, param, param_ids[i]);
00113
00114
00115 status = grp.computeF();
00116 finalStatus =
00117 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00118 finalStatus,
00119 callingFunction);
00120
00121
00122 dfdp = &result[i+1];
00123 dfdp->update(1.0, grp.getF(), -1.0, *f, 0.0);
00124 dfdp->scale(1.0/eps);
00125
00126
00127 grp.setParam(param_ids[i], param);
00128
00129 }
00130
00131 return finalStatus;
00132 }
00133
00134 NOX::Abstract::Group::ReturnType
00135 LOCA::DerivUtils::computeDJnDp(LOCA::MultiContinuation::AbstractGroup& grp,
00136 const vector<int>& paramIDs,
00137 const NOX::Abstract::Vector& nullVector,
00138 NOX::Abstract::MultiVector& result,
00139 bool isValid) const
00140 {
00141 string callingFunction =
00142 "LOCA::DerivUtils::computeDJnDp()";
00143 NOX::Abstract::Group::ReturnType status, finalStatus;
00144
00145
00146 NOX::Abstract::Vector *Jn = &result[0];
00147 NOX::Abstract::Vector *dJndp = NULL;
00148
00149
00150 if (!isValid) {
00151 finalStatus = grp.computeJacobian();
00152 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00153
00154 status = grp.applyJacobian(nullVector, *Jn);
00155 finalStatus =
00156 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00157 finalStatus,
00158 callingFunction);
00159 }
00160 else
00161 finalStatus = NOX::Abstract::Group::Ok;
00162
00163 double param;
00164 double eps;
00165
00166
00167 for (unsigned int i=0; i<paramIDs.size(); i++) {
00168
00169
00170 eps = perturbParam(grp, param, paramIDs[i]);
00171
00172
00173 status = grp.computeJacobian();
00174 finalStatus =
00175 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00176 callingFunction);
00177
00178 dJndp = &result[i+1];
00179 status = grp.applyJacobian(nullVector, *dJndp);
00180 finalStatus =
00181 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00182 finalStatus,
00183 callingFunction);
00184
00185
00186 dJndp->update(-1.0, *Jn, 1.0);
00187 dJndp->scale(1.0/eps);
00188
00189
00190 grp.setParam(paramIDs[i], param);
00191
00192 }
00193
00194 return finalStatus;
00195 }
00196
00197 NOX::Abstract::Group::ReturnType
00198 LOCA::DerivUtils::computeDJnDxa(LOCA::MultiContinuation::AbstractGroup& grp,
00199 const NOX::Abstract::Vector& nullVector,
00200 const NOX::Abstract::MultiVector& aVector,
00201 NOX::Abstract::MultiVector& result) const
00202 {
00203 string callingFunction =
00204 "LOCA::DerivUtils::computeDJnDxa()";
00205 NOX::Abstract::Group::ReturnType status, finalStatus;
00206
00207
00208 Teuchos::RCP<NOX::Abstract::Vector> baseJnVectorPtr =
00209 nullVector.clone(NOX::ShapeCopy);
00210
00211 if (!grp.isJacobian()) {
00212 finalStatus = grp.computeJacobian();
00213 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00214 }
00215 else
00216 finalStatus = NOX::Abstract::Group::Ok;
00217
00218 status = grp.applyJacobian(nullVector, *baseJnVectorPtr);
00219 finalStatus =
00220 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00221 callingFunction);
00222
00223
00224 status = computeDJnDxa(grp, nullVector, aVector, *baseJnVectorPtr, result);
00225 finalStatus =
00226 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00227 callingFunction);
00228
00229 return finalStatus;
00230 }
00231
00232 NOX::Abstract::Group::ReturnType
00233 LOCA::DerivUtils::computeDJnDxa(LOCA::MultiContinuation::AbstractGroup& grp,
00234 const NOX::Abstract::Vector& nullVector,
00235 const NOX::Abstract::MultiVector& aVector,
00236 const NOX::Abstract::Vector& JnVector,
00237 NOX::Abstract::MultiVector& result) const
00238 {
00239 string callingFunction =
00240 "LOCA::DerivUtils::computeDJnDxa()";
00241 NOX::Abstract::Group::ReturnType status, finalStatus;
00242
00243
00244 Teuchos::RCP<NOX::Abstract::Vector> Xvec =
00245 grp.getX().clone(NOX::DeepCopy);
00246
00247
00248 for (int i=0; i<aVector.numVectors(); i++) {
00249
00250
00251 double eps = perturbXVec(grp, *Xvec, aVector[i]);
00252
00253
00254 finalStatus = grp.computeJacobian();
00255 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00256
00257 status = grp.applyJacobian(nullVector, result[i]);
00258 finalStatus =
00259 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00260 finalStatus,
00261 callingFunction);
00262
00263
00264 result[i].update(-1.0, JnVector, 1.0);
00265 result[i].scale(1.0/eps);
00266
00267 }
00268
00269
00270 grp.setX(*Xvec);
00271
00272 return finalStatus;
00273 }
00274
00275 NOX::Abstract::Group::ReturnType
00276 LOCA::DerivUtils::computeDwtJnDp(
00277 LOCA::MultiContinuation::AbstractGroup& grp,
00278 const vector<int>& paramIDs,
00279 const NOX::Abstract::Vector& w,
00280 const NOX::Abstract::Vector& nullVector,
00281 NOX::Abstract::MultiVector::DenseMatrix& result,
00282 bool isValid) const
00283 {
00284 string callingFunction =
00285 "LOCA::DerivUtils::computeDwtJnDp()";
00286 NOX::Abstract::Group::ReturnType status, finalStatus;
00287
00288
00289 Teuchos::RCP<NOX::Abstract::Vector> Jn =
00290 w.clone(NOX::ShapeCopy);
00291 double base_wtJn;
00292
00293
00294 if (!isValid) {
00295
00296
00297 finalStatus = grp.computeJacobian();
00298 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00299
00300
00301 status = grp.applyJacobian(nullVector, *Jn);
00302 finalStatus =
00303 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00304 finalStatus,
00305 callingFunction);
00306
00307
00308 base_wtJn = w.innerProduct(*Jn);
00309 result(0,0) = base_wtJn;
00310 }
00311 else {
00312 base_wtJn = result(0,0);
00313 finalStatus = NOX::Abstract::Group::Ok;
00314 }
00315
00316 double param;
00317 double eps;
00318 double perturb_wtJn;
00319
00320
00321 for (unsigned int i=0; i<paramIDs.size(); i++) {
00322
00323
00324 eps = perturbParam(grp, param, paramIDs[i]);
00325
00326
00327 status = grp.computeJacobian();
00328 finalStatus =
00329 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00330 callingFunction);
00331
00332 status = grp.applyJacobian(nullVector, *Jn);
00333 finalStatus =
00334 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00335 finalStatus,
00336 callingFunction);
00337
00338
00339 perturb_wtJn = w.innerProduct(*Jn);
00340
00341
00342 result(0,i+1) = (perturb_wtJn - base_wtJn) / eps;
00343
00344
00345 grp.setParam(paramIDs[i], param);
00346
00347 }
00348
00349 return finalStatus;
00350 }
00351
00352 NOX::Abstract::Group::ReturnType
00353 LOCA::DerivUtils::computeDwtJDp(LOCA::MultiContinuation::AbstractGroup& grp,
00354 const vector<int>& paramIDs,
00355 const NOX::Abstract::Vector& w,
00356 NOX::Abstract::MultiVector& result,
00357 bool isValid) const
00358 {
00359 string callingFunction =
00360 "LOCA::DerivUtils::computeDwtJDp()";
00361 NOX::Abstract::Group::ReturnType status, finalStatus;
00362
00363
00364 NOX::Abstract::Vector *wtJ = &result[0];
00365 NOX::Abstract::Vector *dwtJdp = NULL;
00366
00367
00368 if (!isValid) {
00369 finalStatus = grp.computeJacobian();
00370 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00371
00372 status = grp.applyJacobianTranspose(w, *wtJ);
00373 finalStatus =
00374 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00375 finalStatus,
00376 callingFunction);
00377 }
00378 else
00379 finalStatus = NOX::Abstract::Group::Ok;
00380
00381 double param;
00382 double eps;
00383
00384
00385 for (unsigned int i=0; i<paramIDs.size(); i++) {
00386
00387
00388 eps = perturbParam(grp, param, paramIDs[i]);
00389
00390
00391 status = grp.computeJacobian();
00392 finalStatus =
00393 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00394 callingFunction);
00395
00396 dwtJdp = &result[i+1];
00397 status = grp.applyJacobianTranspose(w, *dwtJdp);
00398 finalStatus =
00399 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00400 finalStatus,
00401 callingFunction);
00402
00403
00404 dwtJdp->update(-1.0, *wtJ, 1.0);
00405 dwtJdp->scale(1.0/eps);
00406
00407
00408 grp.setParam(paramIDs[i], param);
00409
00410 }
00411
00412 return finalStatus;
00413 }
00414
00415 NOX::Abstract::Group::ReturnType
00416 LOCA::DerivUtils::computeDwtJnDx(LOCA::MultiContinuation::AbstractGroup& grp,
00417 const NOX::Abstract::Vector& w,
00418 const NOX::Abstract::Vector& nullVector,
00419 NOX::Abstract::Vector& result) const
00420 {
00421 string callingFunction =
00422 "LOCA::DerivUtils::computeDwtJnDx()";
00423 NOX::Abstract::Group::ReturnType status, finalStatus;
00424
00425
00426 Teuchos::RCP<NOX::Abstract::Vector> wtJ =
00427 w.clone(NOX::ShapeCopy);
00428
00429
00430 finalStatus = grp.computeJacobian();
00431 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00432
00433 status = grp.applyJacobianTranspose(w, *wtJ);
00434 finalStatus =
00435 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00436 callingFunction);
00437
00438
00439 Teuchos::RCP<NOX::Abstract::Vector> Xvec =
00440 grp.getX().clone(NOX::DeepCopy);
00441
00442
00443 double eps = perturbXVec(grp, *Xvec, nullVector);
00444
00445
00446 finalStatus = grp.computeJacobian();
00447 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00448
00449 status = grp.applyJacobianTranspose(w, result);
00450 finalStatus =
00451 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00452 finalStatus,
00453 callingFunction);
00454
00455
00456 result.update(-1.0, *wtJ, 1.0);
00457 result.scale(1.0/eps);
00458
00459
00460 grp.setX(*Xvec);
00461
00462 return finalStatus;
00463 }
00464
00465 NOX::Abstract::Group::ReturnType
00466 LOCA::DerivUtils::computeDwtJnDx(LOCA::MultiContinuation::AbstractGroup& grp,
00467 const NOX::Abstract::MultiVector& w,
00468 const NOX::Abstract::Vector& nullVector,
00469 NOX::Abstract::MultiVector& result) const
00470 {
00471 string callingFunction =
00472 "LOCA::DerivUtils::computeDwtJnDx()";
00473 NOX::Abstract::Group::ReturnType status, finalStatus;
00474
00475
00476 Teuchos::RCP<NOX::Abstract::MultiVector> wtJ =
00477 w.clone(NOX::ShapeCopy);
00478
00479
00480 finalStatus = grp.computeJacobian();
00481 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00482
00483 status = grp.applyJacobianTransposeMultiVector(w, *wtJ);
00484 finalStatus =
00485 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00486 callingFunction);
00487
00488
00489 Teuchos::RCP<NOX::Abstract::Vector> Xvec =
00490 grp.getX().clone(NOX::DeepCopy);
00491
00492
00493 double eps = perturbXVec(grp, *Xvec, nullVector);
00494
00495
00496 finalStatus = grp.computeJacobian();
00497 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00498
00499 status = grp.applyJacobianTransposeMultiVector(w, result);
00500 finalStatus =
00501 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00502 finalStatus,
00503 callingFunction);
00504
00505
00506 result.update(-1.0, *wtJ, 1.0);
00507 result.scale(1.0/eps);
00508
00509
00510 grp.setX(*Xvec);
00511
00512 return finalStatus;
00513 }
00514
00515 NOX::Abstract::Group::ReturnType
00516 LOCA::DerivUtils::computeDCeDp(LOCA::Hopf::MooreSpence::AbstractGroup& grp,
00517 const vector<int>& paramIDs,
00518 const NOX::Abstract::Vector& yVector,
00519 const NOX::Abstract::Vector& zVector,
00520 double w,
00521 NOX::Abstract::MultiVector& result_real,
00522 NOX::Abstract::MultiVector& result_imag,
00523 bool isValid) const
00524 {
00525 string callingFunction =
00526 "LOCA::DerivUtils::computeDCeDp()";
00527 NOX::Abstract::Group::ReturnType status, finalStatus;
00528
00529
00530 NOX::Abstract::Vector& CeReal = result_real[0];
00531 NOX::Abstract::Vector& CeImag = result_imag[0];
00532 NOX::Abstract::Vector* dCedpReal;
00533 NOX::Abstract::Vector* dCedpImag;
00534
00535
00536 if (!isValid) {
00537 finalStatus = grp.computeComplex(w);
00538 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00539
00540 status = grp.applyComplex(yVector, zVector, CeReal, CeImag);
00541 finalStatus =
00542 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00543 finalStatus,
00544 callingFunction);
00545 }
00546 else
00547 finalStatus = NOX::Abstract::Group::Ok;
00548
00549 double param;
00550 double eps;
00551
00552
00553 for (unsigned int i=0; i<paramIDs.size(); i++) {
00554
00555
00556 eps = perturbParam(grp, param, paramIDs[i]);
00557
00558
00559 status = grp.computeComplex(w);
00560 finalStatus =
00561 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00562 callingFunction);
00563
00564 dCedpReal = &(result_real[i+1]);
00565 dCedpImag = &(result_imag[i+1]);
00566 status = grp.applyComplex(yVector, zVector, *dCedpReal, *dCedpImag);
00567 finalStatus =
00568 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00569 finalStatus,
00570 callingFunction);
00571
00572
00573 dCedpReal->update(-1.0, CeReal, 1.0);
00574 dCedpReal->scale(1.0/eps);
00575 dCedpImag->update(-1.0, CeImag, 1.0);
00576 dCedpImag->scale(1.0/eps);
00577
00578
00579 grp.setParam(paramIDs[i], param);
00580
00581 }
00582
00583 return finalStatus;
00584 }
00585
00586 NOX::Abstract::Group::ReturnType
00587 LOCA::DerivUtils::computeDCeDxa(
00588 LOCA::Hopf::MooreSpence::AbstractGroup& grp,
00589 const NOX::Abstract::Vector& yVector,
00590 const NOX::Abstract::Vector& zVector,
00591 double w,
00592 const NOX::Abstract::MultiVector& aVector,
00593 NOX::Abstract::MultiVector& result_real,
00594 NOX::Abstract::MultiVector& result_imag) const
00595 {
00596 string callingFunction =
00597 "LOCA::DerivUtils::computeDCeDxa()";
00598 NOX::Abstract::Group::ReturnType status, finalStatus;
00599
00600
00601 Teuchos::RCP<NOX::Abstract::Vector> CeReal =
00602 yVector.clone(NOX::ShapeCopy);
00603 Teuchos::RCP<NOX::Abstract::Vector> CeImag =
00604 zVector.clone(NOX::ShapeCopy);
00605
00606
00607 finalStatus = grp.computeComplex(w);
00608 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00609
00610 status = grp.applyComplex(yVector, zVector, *CeReal, *CeImag);
00611 finalStatus =
00612 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00613 finalStatus,
00614 callingFunction);
00615
00616
00617 status =
00618 computeDCeDxa(grp, yVector, zVector, w, aVector, *CeReal, *CeImag,
00619 result_real, result_imag);
00620 finalStatus =
00621 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00622 callingFunction);
00623
00624 return finalStatus;
00625 }
00626
00627 NOX::Abstract::Group::ReturnType
00628 LOCA::DerivUtils::computeDCeDxa(
00629 LOCA::Hopf::MooreSpence::AbstractGroup& grp,
00630 const NOX::Abstract::Vector& yVector,
00631 const NOX::Abstract::Vector& zVector,
00632 double w,
00633 const NOX::Abstract::MultiVector& aVector,
00634 const NOX::Abstract::Vector& Ce_real,
00635 const NOX::Abstract::Vector& Ce_imag,
00636 NOX::Abstract::MultiVector& result_real,
00637 NOX::Abstract::MultiVector& result_imag) const
00638 {
00639 string callingFunction =
00640 "LOCA::DerivUtils::computeDCeDxa()";
00641 NOX::Abstract::Group::ReturnType status, finalStatus =
00642 NOX::Abstract::Group::Ok;
00643
00644
00645 Teuchos::RCP<NOX::Abstract::Vector> Xvec =
00646 grp.getX().clone(NOX::DeepCopy);
00647
00648
00649 for (int i=0; i<aVector.numVectors(); i++) {
00650
00651
00652 double eps = perturbXVec(grp, *Xvec, aVector[i]);
00653
00654
00655 status = grp.computeComplex(w);
00656 finalStatus =
00657 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00658 finalStatus,
00659 callingFunction);
00660
00661 status =
00662 grp.applyComplex(yVector, zVector, result_real[i], result_imag[i]);
00663 finalStatus =
00664 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00665 finalStatus,
00666 callingFunction);
00667
00668
00669 result_real[i].update(-1.0, Ce_real, 1.0); result_real[i].scale(1.0/eps);
00670 result_imag[i].update(-1.0, Ce_imag, 1.0); result_imag[i].scale(1.0/eps);
00671
00672 }
00673
00674
00675 grp.setX(*Xvec);
00676
00677 return finalStatus;
00678 }
00679
00680 NOX::Abstract::Group::ReturnType
00681 LOCA::DerivUtils::computeDwtCeDp(
00682 LOCA::Hopf::MinimallyAugmented::AbstractGroup& grp,
00683 const vector<int>& paramIDs,
00684 const NOX::Abstract::Vector& w1,
00685 const NOX::Abstract::Vector& w2,
00686 const NOX::Abstract::Vector& yVector,
00687 const NOX::Abstract::Vector& zVector,
00688 double omega,
00689 NOX::Abstract::MultiVector::DenseMatrix& result_real,
00690 NOX::Abstract::MultiVector::DenseMatrix& result_imag,
00691 bool isValid) const
00692 {
00693 string callingFunction =
00694 "LOCA::DerivUtils::computeDwtCeDp()";
00695 NOX::Abstract::Group::ReturnType status, finalStatus;
00696
00697
00698 Teuchos::RCP<NOX::Abstract::Vector> CeReal =
00699 w1.clone(NOX::ShapeCopy);
00700 Teuchos::RCP<NOX::Abstract::Vector> CeImag =
00701 w2.clone(NOX::ShapeCopy);
00702
00703
00704 if (!isValid) {
00705 finalStatus = grp.computeComplex(omega);
00706 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00707
00708 status = grp.applyComplex(yVector, zVector, *CeReal, *CeImag);
00709 finalStatus =
00710 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00711 finalStatus,
00712 callingFunction);
00713
00714
00715 result_real(0,0) = w1.innerProduct(*CeReal) + w2.innerProduct(*CeImag);
00716 result_imag(0,0) = w1.innerProduct(*CeImag) - w2.innerProduct(*CeReal);
00717 }
00718 else
00719 finalStatus = NOX::Abstract::Group::Ok;
00720
00721 double param;
00722 double eps;
00723
00724
00725 for (unsigned int i=0; i<paramIDs.size(); i++) {
00726
00727
00728 eps = perturbParam(grp, param, paramIDs[i]);
00729
00730
00731 status = grp.computeComplex(omega);
00732 finalStatus =
00733 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00734 callingFunction);
00735
00736 status = grp.applyComplex(yVector, zVector, *CeReal, *CeImag);
00737 finalStatus =
00738 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00739 finalStatus,
00740 callingFunction);
00741
00742
00743
00744 result_real(0,i+1) = (w1.innerProduct(*CeReal) + w2.innerProduct(*CeImag) -
00745 result_real(0,0)) / eps;
00746 result_imag(0,i+1) = (w1.innerProduct(*CeImag) - w2.innerProduct(*CeReal) -
00747 result_imag(0,0)) /eps;
00748
00749
00750 grp.setParam(paramIDs[i], param);
00751
00752 }
00753
00754 return finalStatus;
00755 }
00756
00757 NOX::Abstract::Group::ReturnType
00758 LOCA::DerivUtils::computeDwtCeDx(
00759 LOCA::Hopf::MinimallyAugmented::AbstractGroup& grp,
00760 const NOX::Abstract::Vector& w1,
00761 const NOX::Abstract::Vector& w2,
00762 const NOX::Abstract::Vector& yVector,
00763 const NOX::Abstract::Vector& zVector,
00764 double omega,
00765 NOX::Abstract::Vector& result_real,
00766 NOX::Abstract::Vector& result_imag) const
00767 {
00768 string callingFunction =
00769 "LOCA::DerivUtils::computeDwtCeDxa()";
00770 NOX::Abstract::Group::ReturnType status, finalStatus =
00771 NOX::Abstract::Group::Ok;
00772
00773
00774 Teuchos::RCP<NOX::Abstract::Vector> wtC_real =
00775 w1.clone(NOX::ShapeCopy);
00776 Teuchos::RCP<NOX::Abstract::Vector> wtC_imag =
00777 w2.clone(NOX::ShapeCopy);
00778
00779
00780 finalStatus = grp.computeComplex(omega);
00781 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00782
00783 status = grp.applyComplexTranspose(w1, w2, *wtC_real, *wtC_imag);
00784 finalStatus =
00785 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00786 callingFunction);
00787
00788
00789 Teuchos::RCP<NOX::Abstract::Vector> Xvec =
00790 grp.getX().clone(NOX::DeepCopy);
00791
00792
00793 double eps = perturbXVec(grp, *Xvec, yVector);
00794
00795
00796 status = grp.computeComplex(omega);
00797 finalStatus =
00798 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00799 finalStatus,
00800 callingFunction);
00801
00802 status =
00803 grp.applyComplexTranspose(w1, w2, result_real, result_imag);
00804 finalStatus =
00805 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00806 finalStatus,
00807 callingFunction);
00808
00809
00810 result_real.update(-1.0, *wtC_real, 1.0); result_real.scale(1.0/eps);
00811 result_imag.update(-1.0, *wtC_imag, 1.0); result_imag.scale(1.0/eps);
00812
00813
00814 result_imag.scale(-1.0);
00815
00816
00817 grp.setX(*Xvec);
00818
00819
00820 eps = perturbXVec(grp, *Xvec, zVector);
00821
00822
00823 status = grp.computeComplex(omega);
00824 finalStatus =
00825 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00826 finalStatus,
00827 callingFunction);
00828
00829 Teuchos::RCP<NOX::Abstract::Vector> tmp_r =
00830 result_real.clone(NOX::ShapeCopy);
00831 Teuchos::RCP<NOX::Abstract::Vector> tmp_i =
00832 result_imag.clone(NOX::ShapeCopy);
00833 status =
00834 grp.applyComplexTranspose(w1, w2, *tmp_r, *tmp_i);
00835 finalStatus =
00836 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00837 finalStatus,
00838 callingFunction);
00839
00840
00841 tmp_r->update(-1.0, *wtC_real, 1.0); tmp_r->scale(1.0/eps);
00842 tmp_i->update(-1.0, *wtC_imag, 1.0); tmp_i->scale(1.0/eps);
00843
00844
00845 tmp_i->scale(-1.0);
00846
00847 result_real.update(-1.0, *tmp_i, 1.0);
00848 result_imag.update( 1.0, *tmp_r, 1.0);
00849
00850
00851 grp.setX(*Xvec);
00852
00853 return finalStatus;
00854 }
00855
00856
00857
00858
00859
00860 double
00861 LOCA::DerivUtils::perturbParam(LOCA::MultiContinuation::AbstractGroup& grp,
00862 double& paramOrig,
00863 int param_id) const
00864 {
00865 paramOrig = grp.getParam(param_id);
00866
00867
00868 double eps = epsScalar(paramOrig);
00869 double param = paramOrig + eps;
00870
00871
00872 grp.setParam(param_id, param);
00873
00874
00875 return eps;
00876 }
00877
00878 double
00879 LOCA::DerivUtils::perturbXVec(LOCA::MultiContinuation::AbstractGroup& grp,
00880 const NOX::Abstract::Vector& xVector,
00881 const NOX::Abstract::Vector& aVector) const
00882 {
00883
00884 Teuchos::RCP<NOX::Abstract::Vector> tmpXVecPtr =
00885 xVector.clone(NOX::DeepCopy);
00886
00887
00888 double eps = epsVector(*tmpXVecPtr, aVector);
00889
00890
00891 grp.setX(tmpXVecPtr->update(eps, aVector, 1.0));
00892
00893
00894 return eps;
00895 }
00896
00897 double
00898 LOCA::DerivUtils::epsScalar(double p) const
00899 {
00900 return perturb * (perturb + fabs(p));
00901 }
00902
00903 double
00904 LOCA::DerivUtils::epsVector(const NOX::Abstract::Vector& xVector,
00905 const NOX::Abstract::Vector& aVector) const
00906 {
00907 return perturb * (perturb + xVector.norm(NOX::Abstract::Vector::TwoNorm)
00908 / (aVector.norm(NOX::Abstract::Vector::TwoNorm) + perturb));
00909 }