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 "NOX_Common.H"
00040 #include "NOX_Epetra_BroydenOperator.H"
00041
00042
00043 #ifdef HAVE_NOX_DEBUG
00044 #ifdef HAVE_NOX_EPETRAEXT
00045 #include "EpetraExt_BlockMapOut.h"
00046 #include "EpetraExt_MultiVectorOut.h"
00047 #include "EpetraExt_RowMatrixOut.h"
00048 #endif
00049 #endif
00050
00051 #include "Epetra_Map.h"
00052
00053 using namespace NOX;
00054 using namespace NOX::Epetra;
00055
00056 BroydenOperator::BroydenOperator(
00057 Teuchos::ParameterList & nlParams_ ,
00058 const Teuchos::RCP<NOX::Utils>& utils_ ,
00059 Epetra_Vector & solnVec ,
00060 const Teuchos::RCP<Epetra_CrsMatrix>& mat ,
00061 bool verbose_ ) :
00062 verbose ( verbose_ ) ,
00063 crsMatrix ( Teuchos::rcp( new Epetra_CrsMatrix(*mat)) ) ,
00064 nlParams ( nlParams_ ) ,
00065 utils ( utils_ ) ,
00066 prePostOperator( utils, nlParams.sublist("Solver Options") ) ,
00067 label ( "NOX::Epetra::BroydenOperator" ) ,
00068 isValidStep ( false ) ,
00069 isValidYield ( false ) ,
00070 isValidBroyden ( false ) ,
00071 entriesRemoved ( mat->NumMyRows(), false )
00072 {
00073 initialize( nlParams, solnVec );
00074 }
00075
00076
00077
00078 BroydenOperator::BroydenOperator(const BroydenOperator & bOp) :
00079 verbose ( bOp.verbose ) ,
00080 stepVec ( bOp.stepVec ) ,
00081 yieldVec ( bOp.yieldVec ) ,
00082 workVec ( bOp.workVec ) ,
00083 oldX ( bOp.oldX ) ,
00084 oldF ( bOp.oldF ) ,
00085 crsMatrix ( bOp.crsMatrix ) ,
00086 jacIntPtr ( bOp.jacIntPtr ) ,
00087 jacMatrixPtr ( bOp.jacMatrixPtr ) ,
00088 precIntPtr ( bOp.precIntPtr ) ,
00089 precMatrixPtr ( bOp.precMatrixPtr ) ,
00090 nlParams ( bOp.nlParams ) ,
00091 utils ( bOp.utils ) ,
00092 prePostOperator( utils, nlParams.sublist("Solver Options") ),
00093 label ( "NOX::Epetra::BroydenOperator" ),
00094 isValidStep ( bOp.isValidStep ) ,
00095 isValidYield ( bOp.isValidYield ) ,
00096 isValidBroyden ( bOp.isValidBroyden ) ,
00097 entriesRemoved ( bOp.entriesRemoved )
00098 {
00099 }
00100
00101
00102
00103 bool
00104 BroydenOperator::initialize( Teuchos::ParameterList & nlParams, const Epetra_Vector & vec )
00105 {
00106 stepVec = Teuchos::rcp( new NOX::Epetra::Vector(vec) );
00107 yieldVec = Teuchos::rcp( new NOX::Epetra::Vector(vec) );
00108 workVec = Teuchos::rcp( new NOX::Epetra::Vector(vec) );
00109 oldX = Teuchos::rcp( new NOX::Epetra::Vector(vec) );
00110 oldF = Teuchos::rcp( new NOX::Epetra::Vector(vec) );
00111
00112
00113
00114
00115
00116
00117
00118
00119 Teuchos::RCP<NOX::Abstract::PrePostOperator> me = Teuchos::rcp(this, false);
00120 nlParams.sublist("Solver Options").set("User Defined Pre/Post Operator", me);
00121
00122 return true;
00123 }
00124
00125
00126
00127 BroydenOperator::~BroydenOperator()
00128 {
00129 }
00130
00131 const char* BroydenOperator::Label () const
00132 {
00133 return label.c_str();
00134 }
00135
00136 int BroydenOperator::SetUseTranspose(bool UseTranspose)
00137 {
00138 return crsMatrix->SetUseTranspose(UseTranspose);
00139 }
00140
00141 int BroydenOperator::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00142 {
00143 return crsMatrix->Apply(X, Y);
00144 }
00145
00146 int BroydenOperator::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00147 {
00148 return crsMatrix->ApplyInverse(X, Y);
00149 }
00150
00151 bool BroydenOperator::UseTranspose() const
00152 {
00153 return crsMatrix->UseTranspose();
00154 }
00155
00156 bool BroydenOperator::HasNormInf() const
00157 {
00158 return crsMatrix->HasNormInf();
00159 }
00160
00161 const Epetra_Map& BroydenOperator::OperatorDomainMap() const
00162 {
00163 return crsMatrix->OperatorDomainMap();
00164 }
00165
00166 const Epetra_Map& BroydenOperator::OperatorRangeMap() const
00167 {
00168 return crsMatrix->OperatorRangeMap();
00169 }
00170
00171 bool BroydenOperator::Filled() const
00172 {
00173 return crsMatrix->Filled();
00174 }
00175
00176 int BroydenOperator::NumMyRowEntries(int MyRow, int & NumEntries) const
00177 {
00178 return crsMatrix->NumMyRowEntries(MyRow, NumEntries);
00179 }
00180
00181 int BroydenOperator::MaxNumEntries() const
00182 {
00183 return crsMatrix->MaxNumEntries();
00184 }
00185
00186 inline int BroydenOperator::ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values, int * Indices) const
00187 {
00188 return crsMatrix->ExtractMyRowCopy(MyRow, Length, NumEntries, Values, Indices);
00189 }
00190
00191 int BroydenOperator::ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
00192 {
00193 return crsMatrix->ExtractDiagonalCopy(Diagonal);
00194 }
00195
00196 int BroydenOperator::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00197 {
00198 return crsMatrix->Multiply(TransA, X, Y);
00199 }
00200
00201 int BroydenOperator::Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00202 {
00203 return crsMatrix->Solve(Upper, Trans, UnitDiagonal, X, Y);
00204 }
00205
00206 int BroydenOperator::InvRowSums(Epetra_Vector& x) const
00207 {
00208 return crsMatrix->InvRowSums(x);
00209 }
00210
00211 int BroydenOperator::LeftScale(const Epetra_Vector& x)
00212 {
00213 return crsMatrix->LeftScale(x);
00214 }
00215
00216 int BroydenOperator::InvColSums(Epetra_Vector& x) const
00217 {
00218 return crsMatrix->InvColSums(x);
00219 }
00220
00221 int BroydenOperator::RightScale(const Epetra_Vector& x)
00222 {
00223 return crsMatrix->RightScale(x);
00224 }
00225
00226 double BroydenOperator::NormInf() const
00227 {
00228 return crsMatrix->NormInf();
00229 }
00230
00231 double BroydenOperator::NormOne() const
00232 {
00233 return crsMatrix->NormOne();
00234 }
00235
00236 int BroydenOperator::NumGlobalNonzeros() const
00237 {
00238 return crsMatrix->NumGlobalNonzeros();
00239 }
00240
00241 int BroydenOperator::NumGlobalRows() const
00242 {
00243 return crsMatrix->NumGlobalRows();
00244 }
00245
00246 int BroydenOperator::NumGlobalCols() const
00247 {
00248 return crsMatrix->NumGlobalCols();
00249 }
00250
00251 int BroydenOperator::NumGlobalDiagonals() const
00252 {
00253 return crsMatrix->NumGlobalDiagonals();
00254 }
00255
00256 int BroydenOperator::NumMyNonzeros() const
00257 {
00258 return crsMatrix->NumMyNonzeros();
00259 }
00260
00261 int BroydenOperator::NumMyRows() const
00262 {
00263 return crsMatrix->NumMyRows();
00264 }
00265
00266 int BroydenOperator::NumMyCols() const
00267 {
00268 return crsMatrix->NumMyCols();
00269 }
00270
00271 int BroydenOperator::NumMyDiagonals() const
00272 {
00273 return crsMatrix->NumMyDiagonals();
00274 }
00275
00276 bool BroydenOperator::LowerTriangular() const
00277 {
00278 return crsMatrix->LowerTriangular();
00279 }
00280
00281 bool BroydenOperator::UpperTriangular() const
00282 {
00283 return crsMatrix->UpperTriangular();
00284 }
00285
00286 const Epetra_Comm& BroydenOperator::Comm() const
00287 {
00288 return crsMatrix->Comm();
00289 }
00290
00291 const Epetra_Map& BroydenOperator::RowMatrixRowMap() const
00292 {
00293 return crsMatrix->RowMatrixRowMap();
00294 }
00295
00296 const Epetra_Map& BroydenOperator::RowMatrixColMap() const
00297 {
00298 return crsMatrix->RowMatrixColMap();
00299 }
00300
00301 const Epetra_Import* BroydenOperator::RowMatrixImporter() const
00302 {
00303 return crsMatrix->RowMatrixImporter();
00304 }
00305
00306 const Epetra_BlockMap& BroydenOperator::Map() const
00307 {
00308 return crsMatrix->Map();
00309 }
00310
00311
00312
00313 bool
00314 BroydenOperator::computeJacobian( const Epetra_Vector & x, Epetra_Operator& Jac )
00315 {
00316 bool ok = false;
00317
00318 ok = computeSparseBroydenUpdate();
00319
00320 vector<ReplacementInterface *>::iterator iter = replacementInterfaces.begin() ,
00321 iter_end = replacementInterfaces.end() ;
00322
00323 for( ; iter_end != iter; ++iter )
00324 {
00325 Teuchos::RCP<const Epetra_CrsMatrix> pMat = (*iter)->getReplacementValuesMatrix(x, ReplacementInterface::JACOBIAN );
00326 replaceBroydenMatrixValues( *pMat );
00327 }
00328
00329 return ok;
00330
00331 }
00332
00333
00334
00335 bool
00336 BroydenOperator::computePreconditioner( const Epetra_Vector & x,
00337 Epetra_Operator& Prec,
00338 Teuchos::ParameterList * pList )
00339 {
00340 bool ok = false;
00341
00342 ok = computeSparseBroydenUpdate();
00343
00344 vector<ReplacementInterface *>::iterator iter = replacementInterfaces.begin() ,
00345 iter_end = replacementInterfaces.end() ;
00346
00347 for( ; iter_end != iter; ++iter )
00348 {
00349 Teuchos::RCP<const Epetra_CrsMatrix> pMat = (*iter)->getReplacementValuesMatrix(x, ReplacementInterface::PRECONDITIONER );
00350 replaceBroydenMatrixValues( *pMat );
00351 }
00352
00353 return ok;
00354
00355 }
00356
00357
00358
00359 void
00360 BroydenOperator::runPreSolve( const NOX::Solver::Generic & solver)
00361 {
00362
00363
00364 isValidStep = false;
00365 isValidYield = false;
00366
00367 prePostOperator.runPreSolve( solver );
00368
00369 return;
00370 }
00371
00372
00373
00374 void
00375 BroydenOperator::runPreIterate( const NOX::Solver::Generic & solver)
00376 {
00377 prePostOperator.runPreIterate( solver );
00378
00379 return;
00380 }
00381
00382
00383
00384 void
00385 BroydenOperator::runPostIterate( const NOX::Solver::Generic & solver)
00386 {
00387
00388
00389 if( solver.getNumIterations() > 0 )
00390 {
00391
00392
00393
00394
00395 const Abstract::Group & oldSolnGrp = solver.getPreviousSolutionGroup();
00396 const Abstract::Group & solnGrp = solver.getSolutionGroup();
00397
00398
00399 *workVec = solnGrp.getX();
00400 workVec->update(-1.0, oldSolnGrp.getX(), 1.0);
00401
00402 setStepVector( *workVec );
00403
00404
00405 *workVec = solnGrp.getF();
00406 workVec->update(-1.0, oldSolnGrp.getF(), 1.0);
00407
00408 setYieldVector( *workVec );
00409
00410
00411
00412 }
00413
00414 prePostOperator.runPostIterate( solver );
00415
00416 return;
00417 }
00418
00419
00420
00421 void
00422 BroydenOperator::runPostSolve( const NOX::Solver::Generic & solver)
00423 {
00424 prePostOperator.runPostSolve( solver );
00425
00426 return;
00427 }
00428
00429
00430
00431 void
00432 BroydenOperator::setStepVector( Epetra_Vector & vec )
00433 {
00434 stepVec->getEpetraVector() = vec;
00435 isValidStep = true ;
00436 isValidBroyden = false ;
00437
00438 return;
00439 }
00440
00441
00442
00443 void
00444 BroydenOperator::setStepVector( NOX::Epetra::Vector & vec )
00445 {
00446 *stepVec = vec;
00447 isValidStep = true ;
00448 isValidBroyden = false ;
00449
00450 return;
00451 }
00452
00453
00454
00455 void
00456 BroydenOperator::setYieldVector( NOX::Epetra::Vector & vec )
00457 {
00458 *yieldVec = vec;
00459 isValidYield = true ;
00460 isValidBroyden = false ;
00461
00462 return;
00463 }
00464
00465
00466
00467 void
00468 BroydenOperator::setYieldVector( Epetra_Vector & vec )
00469 {
00470 yieldVec->getEpetraVector() = vec;
00471 isValidYield = true ;
00472 isValidBroyden = false ;
00473
00474 return;
00475 }
00476
00477
00478
00479 bool
00480 BroydenOperator::computeSparseBroydenUpdate()
00481 {
00482
00483 if( isValidBroyden )
00484 return true;
00485 else if( isValidStep && isValidYield )
00486 {
00487
00488 int ierr = crsMatrix->Multiply( false, stepVec->getEpetraVector(), workVec->getEpetraVector() );
00489 if( ierr )
00490 {
00491 cout << "ERROR: NOX::Epetra::BroydenOperator::computeSparseBroydenUpdate(...) "
00492 << "- crsMatrix->Multiply() failed!!!" << endl;
00493 throw "NOX Error: Broyden Update Failed";
00494 }
00495
00496 int numEntries = crsMatrix->NumMyCols();
00497 int * indices ;
00498 double * values ;
00499
00500 for( int row = 0; row < crsMatrix->NumMyRows(); ++row )
00501 {
00502 ierr = crsMatrix->ExtractMyRowView( row, numEntries, values, indices );
00503 if( ierr )
00504 {
00505 cout << "ERROR (" << ierr << ") : NOX::Epetra::BroydenOperator::computeSparseBroydenUpdate() "
00506 << "- crsMatrix->ExtractGlobalRowView(...) failed for row --> " << row << endl;
00507 throw "NOX Error: Broyden Update Failed";
00508 }
00509
00510 double diffVal = yieldVec->getEpetraVector()[row] - workVec->getEpetraVector()[row];
00511
00512 double rowStepInnerProduct = 0.0;
00513
00514 if( entriesRemoved[row] )
00515 {
00516 list<int>::iterator iter = retainedEntries[row].begin() ,
00517 iter_end = retainedEntries[row].end() ;
00518
00519 for( ; iter_end != iter; ++iter )
00520 rowStepInnerProduct +=
00521 stepVec->getEpetraVector()[indices[*iter]] * stepVec->getEpetraVector()[indices[*iter]];
00522
00523 for( iter = retainedEntries[row].begin(); iter_end != iter; ++iter )
00524 values[*iter] += diffVal * stepVec->getEpetraVector()[indices[*iter]] / rowStepInnerProduct;
00525 }
00526 else
00527 {
00528 for( int col = 0; col < numEntries; ++col )
00529 rowStepInnerProduct +=
00530 stepVec->getEpetraVector()[indices[col]] * stepVec->getEpetraVector()[indices[col]];
00531
00532 for( int col = 0; col < numEntries; ++col )
00533 (*values++) += diffVal * stepVec->getEpetraVector()[(*indices++)] / rowStepInnerProduct;
00534 }
00535 }
00536
00537
00538
00539
00540 if( verbose )
00541 cout << "\t... BroydenOperator::computeSparseBroydenUpdate()..." << endl;
00542
00543 }
00544 else
00545 {
00546 cout << "Warning: NOX::Epetra::BroydenOperator::computeSparseBroydenUpdate(...) "
00547 << "- either the step vector or the yield vector or both is not valid." << endl;
00548 cout << "Leaving existing matrix unchanged." << endl;
00549 }
00550
00551
00552 #ifdef HAVE_NOX_DEBUG
00553 #ifdef HAVE_NOX_EPETRAEXT
00554
00555 static int broydenOutputCount;
00556
00557 Teuchos::ParameterList & broydenParamas =
00558 nlParams.sublist("Direction").sublist("Newton").sublist("Broyden Op");
00559
00560 if( broydenParamas.get("Write Broyden Info", false) )
00561 {
00562 std::ostringstream outputNumber;
00563 outputNumber << broydenOutputCount++ ;
00564
00565 std::string prefixName = broydenParamas.get("Write Broyden Info File Prefix", "BroydenOp");
00566 std::string postfixName = outputNumber.str();
00567 postfixName += ".mm";
00568
00569 std::string mapFileName = prefixName + "_Map_" + postfixName;
00570 std::string matFileName = prefixName + "_Matrix_" + postfixName;
00571 std::string dxFileName = prefixName + "_dX_" + postfixName;
00572 std::string dfFileName = prefixName + "_dF_" + postfixName;
00573
00574 Epetra_RowMatrix * printMatrix = NULL;
00575 printMatrix = dynamic_cast<Epetra_RowMatrix*>(crsMatrix.get());
00576
00577 if( NULL == printMatrix )
00578 {
00579 cout << "Error: NOX::Epetra::BroydenOperator::computeSparseBroydenUpdate() - "
00580 << "Could not get a valid crsMatrix!\n"
00581 << "Please set the \"Write Linear System\" parameter to false." << endl;
00582 throw "NOX Error";
00583 }
00584
00585 EpetraExt::BlockMapToMatrixMarketFile(mapFileName.c_str(),
00586 printMatrix->RowMatrixRowMap());
00587 EpetraExt::RowMatrixToMatrixMarketFile(matFileName.c_str(), *printMatrix,
00588 "test matrix", "Broyden Matrix XXX");
00589 EpetraExt::MultiVectorToMatrixMarketFile(dxFileName.c_str(),
00590 stepVec->getEpetraVector());
00591 EpetraExt::MultiVectorToMatrixMarketFile(dfFileName.c_str(),
00592 yieldVec->getEpetraVector());
00593
00594 }
00595 #endif
00596 #endif
00597
00598 return true;
00599 }
00600
00601
00602
00603 void
00604 BroydenOperator::resetBroydenMatrix( const Epetra_CrsMatrix & mat )
00605 {
00606 *crsMatrix = mat;
00607
00608 return;
00609 }
00610
00611
00612
00613 bool
00614 BroydenOperator::isStep() const
00615 {
00616 return isValidStep;
00617 }
00618
00619
00620
00621 bool
00622 BroydenOperator::isYield() const
00623 {
00624 return isValidYield;
00625 }
00626
00627
00628
00629 bool
00630 BroydenOperator::isBroyden() const
00631 {
00632 return isValidBroyden;
00633 }
00634
00635
00636
00637 void
00638 BroydenOperator::removeEntriesFromBroydenUpdate( const Epetra_CrsGraph & graph )
00639 {
00640
00641 int numRemoveIndices ;
00642 int * removeIndPtr ;
00643 int ierr ;
00644
00645 cout << graph << endl;
00646
00647 for( int row = 0; row < graph.NumMyRows(); ++row)
00648 {
00649 ierr = graph.ExtractMyRowView( row, numRemoveIndices, removeIndPtr );
00650 if( ierr )
00651 {
00652 cout << "ERROR (" << ierr << ") : "
00653 << "NOX::Epetra::BroydenOperator::removeEntriesFromBroydenUpdate(...)"
00654 << " - Extract indices error for row --> " << row << endl;
00655 throw "NOX Broyden Operator Error";
00656 }
00657
00658 if( 0 != numRemoveIndices )
00659 {
00660
00661 map<int, bool> removeIndTable;
00662 for( int k = 0; k < numRemoveIndices; ++k )
00663 removeIndTable[ graph.ColMap().GID(removeIndPtr[k]) ] = true;
00664
00665
00666 int numOrigIndices = 0;
00667 int * indPtr;
00668
00669 ierr = crsMatrix->Graph().ExtractMyRowView( row, numOrigIndices, indPtr );
00670 if( ierr )
00671 {
00672 cout << "ERROR (" << ierr << ") : "
00673 << "NOX::Epetra::BroydenOperator::removeEntriesFromBroydenUpdate(...)"
00674 << " - Extract indices error for row --> " << row << endl;
00675 throw "NOX Broyden Operator Error";
00676 }
00677
00678
00679 if( retainedEntries.end() == retainedEntries.find(row) )
00680 {
00681 list<int> inds;
00682
00683 for( int k = 0; k < numOrigIndices; ++k )
00684 {
00685 if( removeIndTable.end() == removeIndTable.find( crsMatrix->Graph().ColMap().GID(indPtr[k]) ) )
00686 inds.push_back(k);
00687 }
00688
00689 retainedEntries[row] = inds;
00690 }
00691 else
00692 {
00693 list<int> & inds = retainedEntries[row];
00694
00695 list<int>::iterator iter = inds.begin() ,
00696 iter_end = inds.end() ;
00697
00698 for( ; iter_end != iter; ++iter )
00699 {
00700 if( !removeIndTable[ *iter ] )
00701 inds.remove( *iter );
00702 }
00703 }
00704
00705 entriesRemoved[row] = true;
00706 }
00707 }
00708
00709 return;
00710 }
00711
00712
00713
00714 void
00715 BroydenOperator::replaceBroydenMatrixValues( const Epetra_CrsMatrix & mat)
00716 {
00717 double * values ;
00718 int * indices ;
00719 int numEntries ;
00720 int ierr ;
00721
00722 for( int row = 0; row < mat.NumMyRows(); ++row)
00723 {
00724 ierr = mat.ExtractMyRowView(row, numEntries, values, indices);
00725 ierr += crsMatrix->ReplaceGlobalValues(row, numEntries, values, indices);
00726 if( ierr )
00727 {
00728 cout << "ERROR (" << ierr << ") : "
00729 << "NOX::Epetra::BroydenOperator::replaceBroydenMatrixValues(...)"
00730 << " - Extract or Replace values error for row --> "
00731 << row << endl;
00732 throw "NOX Broyden Operator Error";
00733 }
00734 }
00735 }
00736
00737
00738
00739 #ifdef HAVE_NOX_DEBUG
00740
00741 void
00742 BroydenOperator::outputActiveEntries()
00743 {
00744
00745 const Epetra_Comm & comm = stepVec->getEpetraVector().Comm();
00746 int numProcs = comm.NumProc();
00747 int myPID = nlParams.sublist("Printing").get("MyPID", 0);
00748
00749 map<int, list<int> >::iterator iter = retainedEntries.begin(),
00750 iter_end = retainedEntries.end() ;
00751
00752 for( int pid = 0; pid < numProcs; ++pid )
00753 {
00754 if( pid == myPID )
00755 {
00756 for( ; iter_end != iter; ++iter )
00757 {
00758 int row = (*iter).first;
00759
00760
00761 int numIndices ;
00762 int * indPtr ;
00763
00764 crsMatrix->Graph().ExtractMyRowView( row, numIndices, indPtr );
00765
00766 list<int> & entries = (*iter).second;
00767
00768 list<int>::iterator colIter = entries.begin() ,
00769 colIter_end = entries.end() ;
00770
00771 for( ; colIter_end != colIter; ++colIter )
00772 cout << "[" << pid << "] \t " << row << " \t " << *colIter << " ("
00773 << indPtr[*colIter] << ")" << endl;
00774 }
00775 }
00776 comm.Barrier();
00777 }
00778
00779 return;
00780 }
00781
00782
00783
00784 #endif