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 #include "Epetra_Comm.h"
00040 #include "Epetra_Map.h"
00041 #include "Epetra_Import.h"
00042 #include "Epetra_Vector.h"
00043 #include "Epetra_CrsGraph.h"
00044 #include "Epetra_CrsMatrix.h"
00045
00046 #include "NOX_Abstract_Group.H"
00047 #include "NOX_Epetra_Vector.H"
00048 #include "Teuchos_ParameterList.hpp"
00049 #include "NOX_Utils.H"
00050
00051 #include "NOX_Epetra_FiniteDifference.H"
00052
00053 using namespace NOX;
00054 using namespace NOX::Epetra;
00055
00056 FiniteDifference::FiniteDifference(
00057 Teuchos::ParameterList& printingParams,
00058 const Teuchos::RCP<NOX::Epetra::Interface::Required>& i,
00059 const NOX::Epetra::Vector& x,
00060 double beta_,
00061 double alpha_) :
00062 utils(printingParams),
00063 interface(i),
00064 x_perturb(x.getEpetraVector()),
00065 fo(x.getEpetraVector()),
00066 fp(x.getEpetraVector()),
00067 Jc(x.getEpetraVector()),
00068 alpha(alpha_),
00069 beta(beta_),
00070 betaType(Scalar),
00071 diffType(Forward),
00072 label("NOX::FiniteDifference Jacobian"),
00073 useGroupForComputeF(false)
00074 {
00075
00076 jacobian = createGraphAndJacobian(*i, x.getEpetraVector());
00077 }
00078
00079 FiniteDifference::
00080 FiniteDifference(
00081 Teuchos::ParameterList& printingParams,
00082 const Teuchos::RCP<NOX::Epetra::Interface::Required>& i,
00083 const NOX::Epetra::Vector& x,
00084 const Teuchos::RCP<const Epetra_Vector>& beta_,
00085 double alpha_) :
00086 utils(printingParams),
00087 interface(i),
00088 x_perturb(x.getEpetraVector()),
00089 fo(x.getEpetraVector()),
00090 fp(x.getEpetraVector()),
00091 Jc(x.getEpetraVector()),
00092 alpha(alpha_),
00093 beta(0),
00094 betaVector(beta_),
00095 betaType(Vector),
00096 diffType(Forward),
00097 label("NOX::FiniteDifference Jacobian"),
00098 useGroupForComputeF(false)
00099 {
00100
00101 jacobian = createGraphAndJacobian(*i, x.getEpetraVector());
00102 }
00103
00104 FiniteDifference::FiniteDifference(
00105 Teuchos::ParameterList& printingParams,
00106 const Teuchos::RCP<NOX::Epetra::Interface::Required>& i,
00107 const NOX::Epetra::Vector& x,
00108 const Teuchos::RCP<Epetra_CrsGraph>& userGraph,
00109 double beta_,
00110 double alpha_) :
00111 utils(printingParams),
00112 graph(userGraph),
00113 interface(i),
00114 x_perturb(x.getEpetraVector()),
00115 fo(x.getEpetraVector()),
00116 fp(x.getEpetraVector()),
00117 Jc(x.getEpetraVector()),
00118 alpha(alpha_),
00119 beta(beta_),
00120 betaType(Scalar),
00121 diffType(Forward),
00122 label("NOX::FiniteDifference Jacobian"),
00123 useGroupForComputeF(false)
00124 {
00125
00126
00127 jacobian = Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
00128 jacobian->FillComplete();
00129 }
00130
00131 FiniteDifference::FiniteDifference(
00132 Teuchos::ParameterList& printingParams,
00133 const Teuchos::RCP<NOX::Epetra::Interface::Required>& i,
00134 const NOX::Epetra::Vector& x,
00135 const Teuchos::RCP<Epetra_CrsGraph>& userGraph,
00136 const Teuchos::RCP<const Epetra_Vector>& beta_,
00137 double alpha_) :
00138 utils(printingParams),
00139 graph(userGraph),
00140 interface(i),
00141 x_perturb(x.getEpetraVector()),
00142 fo(x.getEpetraVector()),
00143 fp(x.getEpetraVector()),
00144 Jc(x.getEpetraVector()),
00145 alpha(alpha_),
00146 beta(0),
00147 betaVector(beta_),
00148 betaType(Vector),
00149 diffType(Forward),
00150 label("NOX::FiniteDifference Jacobian"),
00151 useGroupForComputeF(false)
00152 {
00153
00154
00155 jacobian = Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
00156 jacobian->FillComplete();
00157 }
00158
00159 FiniteDifference::~FiniteDifference()
00160 {
00161
00162 }
00163
00164 const char* FiniteDifference::Label () const
00165 {
00166 return label.c_str();
00167 }
00168
00169 int FiniteDifference::SetUseTranspose(bool UseTranspose)
00170 {
00171 return jacobian->SetUseTranspose(UseTranspose);
00172 }
00173
00174 int FiniteDifference::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00175 {
00176 return jacobian->Apply(X, Y);
00177 }
00178
00179 int FiniteDifference::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00180 {
00181 return jacobian->ApplyInverse(X, Y);
00182 }
00183
00184 bool FiniteDifference::UseTranspose() const
00185 {
00186 return jacobian->UseTranspose();
00187 }
00188
00189 bool FiniteDifference::HasNormInf() const
00190 {
00191 return jacobian->HasNormInf();
00192 }
00193
00194 const Epetra_Map& FiniteDifference::OperatorDomainMap() const
00195 {
00196 return jacobian->OperatorDomainMap();
00197 }
00198
00199 const Epetra_Map& FiniteDifference::OperatorRangeMap() const
00200 {
00201 return jacobian->OperatorRangeMap();
00202 }
00203
00204 bool FiniteDifference::Filled() const
00205 {
00206 return jacobian->Filled();
00207 }
00208
00209 int FiniteDifference::NumMyRowEntries(int MyRow, int & NumEntries) const
00210 {
00211 return jacobian->NumMyRowEntries(MyRow, NumEntries);
00212 }
00213
00214 int FiniteDifference::MaxNumEntries() const
00215 {
00216 return jacobian->MaxNumEntries();
00217 }
00218
00219 int FiniteDifference::ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values, int * Indices) const
00220 {
00221 return jacobian->ExtractMyRowCopy(MyRow, Length, NumEntries, Values, Indices);
00222 }
00223
00224 int FiniteDifference::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00225 {
00226 return jacobian->ExtractDiagonalCopy(Diagonal);
00227 }
00228
00229 int FiniteDifference::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00230 {
00231 return jacobian->Multiply(TransA, X, Y);
00232 }
00233
00234 int FiniteDifference::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00235 {
00236 return jacobian->Solve(Upper, Trans, UnitDiagonal, X, Y);
00237 }
00238
00239 int FiniteDifference::InvRowSums(Epetra_Vector& x) const
00240 {
00241 return jacobian->InvRowSums(x);
00242 }
00243
00244 int FiniteDifference::LeftScale(const Epetra_Vector& x)
00245 {
00246 return jacobian->LeftScale(x);
00247 }
00248
00249 int FiniteDifference::InvColSums(Epetra_Vector& x) const
00250 {
00251 return jacobian->InvColSums(x);
00252 }
00253
00254 int FiniteDifference::RightScale(const Epetra_Vector& x)
00255 {
00256 return jacobian->RightScale(x);
00257 }
00258
00259 double FiniteDifference::NormInf() const
00260 {
00261 return jacobian->NormInf();
00262 }
00263
00264 double FiniteDifference::NormOne() const
00265 {
00266 return jacobian->NormOne();
00267 }
00268
00269 int FiniteDifference::NumGlobalNonzeros() const
00270 {
00271 return jacobian->NumGlobalNonzeros();
00272 }
00273
00274 int FiniteDifference::NumGlobalRows() const
00275 {
00276 return jacobian->NumGlobalRows();
00277 }
00278
00279 int FiniteDifference::NumGlobalCols() const
00280 {
00281 return jacobian->NumGlobalCols();
00282 }
00283
00284 int FiniteDifference::NumGlobalDiagonals() const
00285 {
00286 return jacobian->NumGlobalDiagonals();
00287 }
00288
00289 int FiniteDifference::NumMyNonzeros() const
00290 {
00291 return jacobian->NumMyNonzeros();
00292 }
00293
00294 int FiniteDifference::NumMyRows() const
00295 {
00296 return jacobian->NumMyRows();
00297 }
00298
00299 int FiniteDifference::NumMyCols() const
00300 {
00301 return jacobian->NumMyCols();
00302 }
00303
00304 int FiniteDifference::NumMyDiagonals() const
00305 {
00306 return jacobian->NumMyDiagonals();
00307 }
00308
00309 bool FiniteDifference::LowerTriangular() const
00310 {
00311 return jacobian->LowerTriangular();
00312 }
00313
00314 bool FiniteDifference::UpperTriangular() const
00315 {
00316 return jacobian->UpperTriangular();
00317 }
00318
00319 const Epetra_Comm& FiniteDifference::Comm() const
00320 {
00321 return jacobian->Comm();
00322 }
00323
00324 const Epetra_Map& FiniteDifference::RowMatrixRowMap() const
00325 {
00326 return jacobian->RowMatrixRowMap();
00327 }
00328
00329 const Epetra_Map& FiniteDifference::RowMatrixColMap() const
00330 {
00331 return jacobian->RowMatrixColMap();
00332 }
00333
00334 const Epetra_Import* FiniteDifference::RowMatrixImporter() const
00335 {
00336 return jacobian->RowMatrixImporter();
00337 }
00338
00339 const Epetra_BlockMap& FiniteDifference::Map() const
00340 {
00341 return jacobian->Map();
00342 }
00343
00344 bool FiniteDifference::computeJacobian(const Epetra_Vector& x)
00345 {
00346 return( computeJacobian(x, *this));
00347 }
00348
00349 bool FiniteDifference::computeJacobian(const Epetra_Vector& x, Epetra_Operator& Jac)
00350 {
00351
00352 FiniteDifference* testMatrix = dynamic_cast<FiniteDifference*>(&Jac);
00353 if (testMatrix == 0) {
00354 utils.out() << "ERROR: NOX::Epetra::FiniteDifference::computeJacobian() - "
00355 << "Jacobian to evaluate is not a FiniteDifference object!"
00356 << endl;
00357 throw "NOX Error";
00358 }
00359
00360 const Epetra_BlockMap& map = fo.Map();
00361
00362
00363
00364 Epetra_CrsMatrix& jac = *testMatrix->jacobian;
00365
00366
00367 jac.PutScalar(0.0);
00368
00369
00370 if ( diffType == Centered )
00371 if ( Teuchos::is_null(fmPtr) )
00372 fmPtr = Teuchos::rcp(new Epetra_Vector(x));
00373
00374 double scaleFactor = 1.0;
00375 if ( diffType == Backward )
00376 scaleFactor = -1.0;
00377
00378 double eta = 0.0;
00379
00380 int min = map.MinAllGID();
00381 int max = map.MaxAllGID();
00382 int myMin = map.MinMyGID();
00383 int myMax = map.MaxMyGID();
00384
00385
00386 computeF(x, fo, NOX::Epetra::Interface::Required::FD_Res);
00387
00388 x_perturb = x;
00389
00390
00391 for (int k = min; k < max+1; k++) {
00392
00393
00394
00395 int proc = 0;
00396 if (map.MyGID(k)) {
00397
00398 if (betaType == Scalar)
00399 eta = alpha*x[map.LID(k)] + beta;
00400 else
00401 eta = alpha*x[map.LID(k)] + (*betaVector)[map.LID(k)];
00402
00403 x_perturb[map.LID(k)] += scaleFactor * eta;
00404 proc = map.Comm().MyPID();
00405 }
00406
00407
00408 int broadcastProc = 0;
00409 map.Comm().SumAll(&proc ,&broadcastProc , 1);
00410
00411 map.Comm().Broadcast(&eta, 1, broadcastProc);
00412
00413
00414 computeF(x_perturb,fp, NOX::Epetra::Interface::Required::FD_Res);
00415
00416 if ( diffType == Centered ) {
00417 if (map.MyGID(k))
00418 x_perturb[map.LID(k)] -= 2.0 * eta;
00419 computeF(x_perturb,*fmPtr, NOX::Epetra::Interface::Required::FD_Res);
00420 }
00421
00422
00423 if ( diffType != Centered ) {
00424 Jc.Update(1.0, fp, -1.0, fo, 0.0);
00425 Jc.Scale( 1.0/(scaleFactor * eta) );
00426 }
00427 else {
00428 Jc.Update(1.0, fp, -1.0, *fmPtr, 0.0);
00429 Jc.Scale( 1.0/(2.0 * eta) );
00430 }
00431
00432
00433 for (int j = myMin; j < myMax+1; j++) {
00434 if (!map.MyGID(j))
00435 continue;
00436 if (Jc[map.LID(j)] != 0.0) {
00437 jac.ReplaceGlobalValues(j,1,&Jc[map.LID(j)],&k);
00438 }
00439 }
00440
00441
00442 x_perturb = x;
00443
00444 }
00445
00446 jac.FillComplete();
00447
00448 return true;
00449 }
00450
00451 bool FiniteDifference::computePreconditioner(const Epetra_Vector& x,
00452 Epetra_Operator& Prec,
00453 Teuchos::ParameterList* precParams)
00454 {
00455 return computeJacobian(x, *this);
00456 }
00457
00458 Teuchos::RCP<Epetra_CrsMatrix> FiniteDifference::
00459 createGraphAndJacobian(Interface::Required& i, const Epetra_Vector& x)
00460 {
00461
00462 const Epetra_BlockMap& map = fo.Map();
00463
00464 double eta = 0.0;
00465
00466 int min = map.MinAllGID();
00467 int max = map.MaxAllGID();
00468 int myMin = map.MinMyGID();
00469 int myMax = map.MaxMyGID();
00470
00471
00472 graph = Teuchos::rcp(new Epetra_CrsGraph(Copy,map,10));
00473
00474
00475 computeF(x, fo, NOX::Epetra::Interface::Required::FD_Res);
00476
00477
00478 for (int k = min; k < max+1; k++) {
00479
00480
00481
00482 if (map.MyGID(k)) {
00483
00484 if (betaType == Scalar)
00485 eta = alpha*x[map.LID(k)] + beta;
00486 else
00487 eta = alpha*x[map.LID(k)] + (*betaVector)[map.LID(k)];
00488
00489 x_perturb[map.LID(k)] += eta;
00490 }
00491
00492
00493 computeF(x_perturb, fp, NOX::Epetra::Interface::Required::FD_Res);
00494
00495
00496 Jc.Update(1.0, fp, -1.0, fo, 0.0);
00497
00498
00499
00500 for (int j = myMin; j < myMax+1; j++) {
00501
00502 if (!map.MyGID(j))
00503 continue;
00504 if (Jc[map.LID(j)] != 0.0) {
00505 graph->InsertGlobalIndices(j,1,&k);
00506 }
00507 }
00508
00509
00510 x_perturb = x;
00511
00512 }
00513
00514 graph->FillComplete();
00515 jacobian = Teuchos::rcp(new Epetra_CrsMatrix(Copy, *graph));
00516 jacobian->FillComplete();
00517
00518 return jacobian;
00519
00520 }
00521
00522 void FiniteDifference::setDifferenceMethod(DifferenceType diffType_)
00523 {
00524 diffType = diffType_;
00525 }
00526
00527 Epetra_CrsMatrix& FiniteDifference::getUnderlyingMatrix() const
00528 {
00529 return *jacobian;
00530 }
00531
00532 void FiniteDifference::Print(ostream& strm) const
00533 {
00534 jacobian->Print(strm);
00535 }
00536
00537 void FiniteDifference::setGroupForComputeF(NOX::Abstract::Group& group)
00538 {
00539 useGroupForComputeF = true;
00540 groupPtr = group.clone();
00541 return;
00542 }
00543
00544
00545 bool FiniteDifference::computeF(const Epetra_Vector& input,
00546 Epetra_Vector& result,
00547 NOX::Epetra::Interface::Required::FillType)
00548 {
00549 bool ok = false;
00550
00551 if (!useGroupForComputeF)
00552 ok = interface->computeF(input, result,
00553 NOX::Epetra::Interface::Required::FD_Res);
00554 else {
00555
00556
00557 Epetra_Vector& nonconstInput = const_cast<Epetra_Vector&>(input);
00558
00559 Teuchos::RCP<Epetra_Vector> tmpEpetraInputVec =
00560
00561 Teuchos::rcp(new Epetra_Vector(nonconstInput));
00562 NOX::Epetra::Vector noxX(tmpEpetraInputVec,
00563 NOX::Epetra::Vector::CreateCopy);
00564 groupPtr->setX(noxX);
00565 groupPtr->computeF();
00566 result = dynamic_cast<const NOX::Epetra::Vector&>
00567 (groupPtr->getF()).getEpetraVector();
00568 ok = true;
00569
00570 }
00571
00572 return ok;
00573 }