NOX_Playa_Group.cpp

00001 // $Id$ 
00002 // $Source$ 
00003 
00004 //@HEADER
00005 //   
00006 //@HEADER
00007 
00008 #include "NOX_Common.H"
00009 #include "NOX_Playa_Group.hpp"  // class definition
00010 #include "PlayaMPIComm.hpp"
00011 #include "PlayaOut.hpp"
00012 #include "PlayaTabs.hpp"
00013 
00014 using Playa::Out;
00015 using Playa::Tabs;
00016 using Teuchos::RCP;
00017 
00018 
00019 NOX::NOXPlaya::Group::Group(const Playa::Vector<double>& initcond, 
00020   const Playa::NonlinearOperator<double>& nonlinOp,
00021   const Playa::LinearSolver<double>& linsolver) 
00022   :
00023   precision(3),  // 3 digits of accuracy is default
00024   xVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, DeepCopy))),
00025   fVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00026   newtonVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00027   gradientVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00028   solver(linsolver),
00029   jacobian(),
00030   nonlinearOp(nonlinOp),
00031   normF(0.0)
00032 {  
00033   nonlinearOp.setEvalPt(xVector->getPlayaVector());
00034   resetIsValid();
00035 }
00036 
00037 NOX::NOXPlaya::Group::Group(const Playa::NonlinearOperator<double>& nonlinOp,
00038   const Playa::LinearSolver<double>& linsolver) 
00039   :
00040   precision(3), // 3 digits of accuracy is default
00041   xVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, DeepCopy))),
00042   fVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00043   newtonVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00044   gradientVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00045   solver(linsolver),
00046   jacobian(),
00047   nonlinearOp(nonlinOp),
00048   normF(0.0)
00049 {  
00050   nonlinearOp.setEvalPt(xVector->getPlayaVector());
00051   resetIsValid();
00052 }
00053 
00054 NOX::NOXPlaya::Group::Group(const Playa::Vector<double>& initcond, 
00055   const Playa::NonlinearOperator<double>& nonlinOp,
00056   const Playa::LinearSolver<double>& linsolver,
00057   int numdigits) 
00058   :
00059   precision(numdigits),
00060   xVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, DeepCopy))),
00061   fVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00062   newtonVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00063   gradientVector(rcp(new NOX::NOXPlaya::Vector(initcond, precision, ShapeCopy))),
00064   solver(linsolver),
00065   jacobian(),
00066   nonlinearOp(nonlinOp),
00067   normF(0.0)
00068 {  
00069   nonlinearOp.setEvalPt(xVector->getPlayaVector());
00070   resetIsValid();
00071 }
00072 
00073 NOX::NOXPlaya::Group::Group(const Playa::NonlinearOperator<double>& nonlinOp,
00074   const Playa::LinearSolver<double>& linsolver,
00075   int numdigits) 
00076   :
00077   precision(numdigits),
00078   xVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, DeepCopy))),
00079   fVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00080   newtonVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00081   gradientVector(rcp(new NOX::NOXPlaya::Vector(nonlinOp.getInitialGuess(), precision, ShapeCopy))),
00082   solver(linsolver),
00083   jacobian(),
00084   nonlinearOp(nonlinOp),
00085   normF(0.0)
00086 {  
00087   nonlinearOp.setEvalPt(xVector->getPlayaVector());
00088   resetIsValid();
00089 }
00090 
00091 
00092 NOX::NOXPlaya::Group::Group(const NOX::NOXPlaya::Group& source, NOX::CopyType type) :
00093   precision(source.precision),
00094   xVector(rcp(new NOX::NOXPlaya::Vector(*(source.xVector), precision, type))), 
00095   fVector(rcp(new NOX::NOXPlaya::Vector(*(source.fVector), precision, type))),  
00096   newtonVector(rcp(new NOX::NOXPlaya::Vector(*(source.newtonVector), precision, type))),
00097   gradientVector(rcp(new NOX::NOXPlaya::Vector(*(source.gradientVector), precision, type))),
00098   solver(source.solver),
00099   jacobian(source.jacobian),
00100   nonlinearOp(source.nonlinearOp),
00101   isValidF(false),
00102   isValidJacobian(false),
00103   isValidGradient(false),
00104   isValidNewton(false),
00105   normF(0.0)
00106 {
00107   switch (type) 
00108   {
00109     
00110     case NOX::DeepCopy:
00111     
00112       isValidF = source.isValidF;
00113       isValidGradient = source.isValidGradient;
00114       isValidNewton = source.isValidNewton;
00115       isValidJacobian = source.isValidJacobian;
00116       normF = source.normF;
00117       break;
00118 
00119     case NOX::ShapeCopy:
00120       resetIsValid();
00121       normF = 0.0;
00122       break;
00123 
00124     default:
00125       Out::os() << "NOX:Playa::Group - invalid CopyType for copy constructor." << std::endl;
00126       throw "NOX Playa Error";
00127   }
00128 
00129 }
00130 
00131 NOX::NOXPlaya::Group::~Group() 
00132 {
00133 }
00134 
00135 void NOX::NOXPlaya::Group::resetIsValid() //private
00136 {
00137   isValidF = false;
00138   isValidJacobian = false;
00139   isValidGradient = false;
00140   isValidNewton = false;
00141 }
00142 
00143 
00144 RCP<NOX::Abstract::Group> NOX::NOXPlaya::Group::clone(NOX::CopyType type) const 
00145 {
00146   return rcp(new NOX::NOXPlaya::Group(*this, type));
00147 }
00148 
00149 NOX::Abstract::Group& NOX::NOXPlaya::Group::operator=(const NOX::Abstract::Group& source)
00150 {
00151   return operator=(dynamic_cast<const NOX::NOXPlaya::Group&> (source));
00152 }
00153 
00154 NOX::Abstract::Group& NOX::NOXPlaya::Group::operator=(const NOX::NOXPlaya::Group& source)
00155 {
00156   if (this != &source) 
00157   {
00158 
00159     // Deep Copy of the xVector
00160     *xVector = *(source.xVector);
00161     nonlinearOp = source.nonlinearOp;
00162     solver = source.solver;
00163     jacobian = source.jacobian;
00164     precision = source.precision;
00165 
00166     // Update the isValidVectors
00167     isValidF = source.isValidF;
00168     isValidGradient = source.isValidGradient;
00169     isValidNewton = source.isValidNewton;
00170     isValidJacobian = source.isValidJacobian;
00171     
00172     // Only copy vectors that are valid
00173     if (isValidF) 
00174     {
00175       *fVector = *(source.fVector);
00176       normF = source.normF;
00177     }
00178 
00179     if (isValidGradient)
00180       *gradientVector = *(source.gradientVector);
00181     
00182     if (isValidNewton)
00183       *newtonVector = *(source.newtonVector);
00184     
00185     if (isValidJacobian)
00186       jacobian = source.jacobian;
00187   }
00188   return *this;
00189 }
00190 
00191 void NOX::NOXPlaya::Group::setX(const NOX::Abstract::Vector& y) 
00192 {
00193   setX(dynamic_cast<const NOX::NOXPlaya::Vector&> (y));
00194 }
00195 
00196 void NOX::NOXPlaya::Group::setX(const NOX::NOXPlaya::Vector& y) 
00197 {
00198   resetIsValid();
00199   nonlinearOp.setEvalPt(y.getPlayaVector());
00200   *xVector = y;
00201 }
00202 
00203 void NOX::NOXPlaya::Group::computeX(const NOX::Abstract::Group& grp, 
00204   const NOX::Abstract::Vector& d, 
00205   double step) 
00206 {
00207   // Cast to appropriate type, then call the "native" computeX
00208   const Group& tsfgrp = dynamic_cast<const NOX::NOXPlaya::Group&> (grp);
00209   const Vector& tsfd = dynamic_cast<const NOX::NOXPlaya::Vector&> (d);
00210   computeX(tsfgrp, tsfd, step); 
00211 }
00212 
00213 void NOX::NOXPlaya::Group::computeX(const Group& grp, const Vector& d, double step) 
00214 {
00215   Tabs tab;
00216   resetIsValid();
00217   xVector->update(1.0, *(grp.xVector), step, d);
00218 }
00219 
00220 NOX::Abstract::Group::ReturnType NOX::NOXPlaya::Group::computeF() 
00221 {
00222  
00223   if (nonlinearOp.verb() > 2)
00224   {
00225     Out::os() << "calling computeF()" << std::endl;
00226   }
00227 
00228   if (isValidF)
00229   {
00230     if (nonlinearOp.verb() > 2)
00231     {
00232       Out::os() << "reusing F" << std::endl;
00233     }
00234     return NOX::Abstract::Group::Ok;
00235   }
00236   else
00237   {
00238     if (nonlinearOp.verb() > 2)
00239     {
00240       Out::os() << "computing new F" << std::endl;
00241     }
00242     nonlinearOp.setEvalPt(xVector->getPlayaVector());
00243     *fVector = nonlinearOp.getFunctionValue();
00244     isValidF = true;
00245     normF = fVector->norm();
00246   }
00247 
00248   return (isValidF) ? (NOX::Abstract::Group::Ok) : (NOX::Abstract::Group::Failed);
00249 }
00250 
00251 NOX::Abstract::Group::ReturnType NOX::NOXPlaya::Group::computeJacobian() 
00252 {
00253 
00254   if (nonlinearOp.verb() > 2)
00255   {
00256     Out::os() << "calling computeJ()" << std::endl;
00257   }
00258 
00259   // Skip if the Jacobian is already valid
00260   if (isValidJacobian)
00261   {
00262     return NOX::Abstract::Group::Ok;
00263   }
00264   else
00265   {
00266     nonlinearOp.setEvalPt(xVector->getPlayaVector());
00267     jacobian = nonlinearOp.getJacobian();
00268 
00269     isValidJacobian = true;
00270   }
00271   return (isValidJacobian) ? (NOX::Abstract::Group::Ok) : (NOX::Abstract::Group::Failed);
00272 }
00273 
00274 NOX::Abstract::Group::ReturnType NOX::NOXPlaya::Group::computeGradient() 
00275 {
00276   if (nonlinearOp.verb() > 2)
00277   {
00278     Out::os() << "calling computeGrad()" << std::endl;
00279   }
00280   if (isValidGradient)
00281   {
00282     return NOX::Abstract::Group::Ok;
00283   }
00284   else
00285   {
00286     if (!isF()) 
00287     {
00288       Out::os() << "ERROR: NOX::NOXPlaya::Group::computeGrad() - F is out of date wrt X!" << std::endl;
00289       return NOX::Abstract::Group::BadDependency;
00290     }
00291 
00292     if (!isJacobian()) 
00293     {
00294       Out::os() << "ERROR: NOX::NOXPlaya::Group::computeGrad() - Jacobian is out of date wrt X!" << std::endl;
00295       return NOX::Abstract::Group::BadDependency;
00296     }
00297   
00298     // Compute Gradient = J' * F
00299 
00300     NOX::Abstract::Group::ReturnType status 
00301       = applyJacobianTranspose(*fVector,*gradientVector);
00302     isValidGradient = (status == NOX::Abstract::Group::Ok);
00303 
00304     // Return result
00305     return status;
00306   }
00307 }
00308 
00309 NOX::Abstract::Group::ReturnType 
00310 NOX::NOXPlaya::Group::computeNewton(Teuchos::ParameterList& p) 
00311 {
00312   if (isNewton())
00313   {
00314     return NOX::Abstract::Group::Ok;
00315   }
00316   else
00317   {
00318     if (!isF()) 
00319     {
00320       Out::os() << "ERROR: NOX::Example::Group::computeNewton() - invalid F" 
00321                 << std::endl;
00322       throw "NOX Error";
00323     }
00324 
00325     if (!isJacobian()) 
00326     {
00327       Out::os() << "ERROR: NOX::Example::Group::computeNewton() - invalid Jacobian" << std::endl;
00328       throw "NOX Error";
00329     }
00330 
00331 /*
00332   if (p.isParameter("Tolerance"))
00333   {
00334   double tol = p.get("Tolerance", tol);
00335   solver.updateTolerance(tol);
00336   }
00337 */
00338 
00339     NOX::Abstract::Group::ReturnType status 
00340       = applyJacobianInverse(p, *fVector, *newtonVector);
00341     isValidNewton = (status == NOX::Abstract::Group::Ok);
00342 
00343     
00344     // Scale soln by -1
00345     newtonVector->scale(-1.0);
00346 
00347     if (nonlinearOp.verb() > 0)
00348     {
00349       Out::os() << "newton step" << std::endl;
00350       newtonVector->getPlayaVector().print(Out::os());
00351     }
00352       
00353     // Return solution
00354     return status;
00355   }
00356 }
00357 
00358 
00359 NOX::Abstract::Group::ReturnType 
00360 NOX::NOXPlaya::Group::applyJacobian(const Abstract::Vector& input, 
00361   NOX::Abstract::Vector& result) const
00362 {
00363   const NOX::NOXPlaya::Vector& tsfinput = dynamic_cast<const NOX::NOXPlaya::Vector&> (input);
00364   NOX::NOXPlaya::Vector& tsfresult = dynamic_cast<NOX::NOXPlaya::Vector&> (result);
00365   return applyJacobian(tsfinput, tsfresult);
00366 }
00367 
00368 NOX::Abstract::Group::ReturnType 
00369 NOX::NOXPlaya::Group::applyJacobian(const NOX::NOXPlaya::Vector& input, 
00370   NOX::NOXPlaya::Vector& result) const
00371 {
00372   // Check validity of the Jacobian
00373   if (!isJacobian()) 
00374   {
00375     return NOX::Abstract::Group::BadDependency;
00376   }
00377   else
00378   {
00379     // Compute result = J * input
00380     jacobian.apply(input.getPlayaVector(),result.getPlayaVector());
00381     return NOX::Abstract::Group::Ok;
00382   }
00383 }
00384 
00385 NOX::Abstract::Group::ReturnType 
00386 NOX::NOXPlaya::Group::applyJacobianTranspose(const NOX::Abstract::Vector& input, 
00387   NOX::Abstract::Vector& result) const
00388 {
00389   const NOX::NOXPlaya::Vector& tsfinput = dynamic_cast<const NOX::NOXPlaya::Vector&> (input);
00390   NOX::NOXPlaya::Vector& tsfresult = dynamic_cast<NOX::NOXPlaya::Vector&> (result);
00391   return applyJacobianTranspose(tsfinput, tsfresult);
00392 }
00393 
00394 NOX::Abstract::Group::ReturnType 
00395 NOX::NOXPlaya::Group::applyJacobianTranspose(const NOX::NOXPlaya::Vector& input, NOX::NOXPlaya::Vector& result) const
00396 {
00397   // Check validity of the Jacobian
00398   if (!isJacobian()) 
00399   {
00400     return NOX::Abstract::Group::BadDependency;
00401   }
00402   else
00403   {
00404     // Compute result = J^T * input
00405     jacobian.applyTranspose(input.getPlayaVector(), result.getPlayaVector());
00406       
00407     return NOX::Abstract::Group::Ok;
00408   }
00409 }
00410 
00411 NOX::Abstract::Group::ReturnType 
00412 NOX::NOXPlaya::Group::applyJacobianInverse(Teuchos::ParameterList& p, 
00413   const Abstract::Vector& input, 
00414   NOX::Abstract::Vector& result) const 
00415 {
00416   const NOX::NOXPlaya::Vector& tsfinput = dynamic_cast<const NOX::NOXPlaya::Vector&> (input);
00417   NOX::NOXPlaya::Vector& tsfresult = dynamic_cast<NOX::NOXPlaya::Vector&> (result); 
00418   return applyJacobianInverse(p, tsfinput, tsfresult);
00419 }
00420 
00421 NOX::Abstract::Group::ReturnType 
00422 NOX::NOXPlaya::Group::applyJacobianInverse(Teuchos::ParameterList& p, 
00423   const NOX::NOXPlaya::Vector& input, 
00424   NOX::NOXPlaya::Vector& result) const 
00425 {
00426 
00427   if (!isJacobian()) 
00428   {
00429     Out::os() << "ERROR: NOX::NOXPlaya::Group::applyJacobianInverse() - invalid Jacobian" << std::endl;
00430     throw "NOX Error";
00431 
00432   }
00433 /*
00434   if (p.isParameter("Tolerance"))
00435   {
00436   double tol = p.get("Tolerance", tol);
00437   solver.updateTolerance(tol);
00438   }
00439 */
00440   if (nonlinearOp.verb() > 4)
00441   {
00442     Out::os() << "---------------- applying J^-1 ------------------" << std::endl;
00443     Out::os() << "J=" << std::endl;
00444     jacobian.print(Out::os());
00445     Out::os() << "F=" << std::endl;
00446     input.getPlayaVector().print(Out::os());
00447   }
00448 
00449   Playa::SolverState<double> status 
00450     = solver.solve(jacobian, input.getPlayaVector(),
00451       result.getPlayaVector());
00452   
00453   if (status.finalState() != Playa::SolveConverged)
00454   {
00455     return NOX::Abstract::Group::Failed;
00456   }
00457   else
00458   {
00459     if (nonlinearOp.verb() > 2)
00460     {
00461       Out::os() << "soln=" << std::endl;
00462       result.getPlayaVector().print(Out::os());
00463     }
00464     return NOX::Abstract::Group::Ok;
00465   }
00466 
00467 }
00468 
00469 bool NOX::NOXPlaya::Group::isF() const 
00470 {   
00471   return isValidF;
00472 }
00473 
00474 bool NOX::NOXPlaya::Group::isJacobian() const 
00475 {  
00476   return isValidJacobian;
00477 }
00478 
00479 bool NOX::NOXPlaya::Group::isGradient() const 
00480 {   
00481   return isValidGradient;
00482 }
00483 
00484 bool NOX::NOXPlaya::Group::isNewton() const 
00485 {   
00486   return isValidNewton;
00487 }
00488 
00489 const NOX::Abstract::Vector& NOX::NOXPlaya::Group::getX() const 
00490 {
00491   return *xVector;
00492 }
00493 
00494 const NOX::Abstract::Vector& NOX::NOXPlaya::Group::getF() const 
00495 {  
00496   if (nonlinearOp.verb() > 2)
00497   {
00498     Out::os() << "calling getF()" << std::endl;
00499   }
00500   TEUCHOS_TEST_FOR_EXCEPTION(!isF(), runtime_error, 
00501     "calling getF() with invalid function value");
00502   return *fVector;
00503 }
00504 
00505 double NOX::NOXPlaya::Group::getNormF() const
00506 {
00507   if (nonlinearOp.verb() > 2)
00508   {
00509     Out::os() << "normF = " << normF << std::endl;
00510   }
00511   TEUCHOS_TEST_FOR_EXCEPTION(!isF(), runtime_error, 
00512     "calling normF() with invalid function value");
00513   return normF;
00514 }
00515 
00516 const NOX::Abstract::Vector& NOX::NOXPlaya::Group::getGradient() const 
00517 { 
00518   return *gradientVector;
00519 }
00520 
00521 const NOX::Abstract::Vector& NOX::NOXPlaya::Group::getNewton() const 
00522 {
00523   return *newtonVector;
00524 }
00525 
00526 void NOX::NOXPlaya::Group::print() const
00527 {
00528   cout << "x = " << *xVector << "\n";
00529 
00530   if (isValidF) {
00531     cout << "F(x) = " << *fVector << "\n";
00532     cout << "|| F(x) || = " << normF << "\n";
00533   }
00534   else
00535     cout << "F(x) has not been computed" << "\n";
00536   
00537   cout << std::endl;
00538 }
00539 
00540 
00541 

doxygen