NLPThyraEpetraMultiPointModelEval4DOptMain.cpp

00001 #include "NLPInterfacePack_NLPFirstOrderThyraModelEvaluator.hpp"
00002 #ifdef HAVE_MPI
00003 #include "EpetraExt_MultiPointModelEvaluator.h"
00004 #include "EpetraExt_MultiMpiComm.h"
00005 #endif
00006 #include "EpetraMultiPointModelEval4DOpt.hpp"
00007 #include "MoochoPack_MoochoThyraSolver.hpp"
00008 #include "Stratimikos_DefaultLinearSolverBuilder.hpp"
00009 #include "Thyra_EpetraModelEvaluator.hpp"
00010 #include "Teuchos_RCP.hpp"
00011 #include "Teuchos_GlobalMPISession.hpp"
00012 #include "Teuchos_VerboseObject.hpp"
00013 #include "Teuchos_CommandLineProcessor.hpp"
00014 #include "Teuchos_StandardCatchMacros.hpp"
00015 
00016 int main( int argc, char* argv[] )
00017 #ifdef HAVE_MPI
00018 {
00019   using Teuchos::CommandLineProcessor;
00020   typedef AbstractLinAlgPack::value_type  Scalar;
00021   using MoochoPack::MoochoSolver;
00022   using MoochoPack::MoochoThyraSolver;
00023 
00024   bool dummySuccess = true;
00025 
00026   Teuchos::GlobalMPISession mpiSession(&argc,&argv);
00027 
00028   Teuchos::RCP<Teuchos::FancyOStream>
00029     out = Teuchos::VerboseObjectBase::getDefaultOStream();
00030 
00031   try {
00032 
00033     Stratimikos::DefaultLinearSolverBuilder lowsfCreator;
00034     MoochoThyraSolver                     solver;
00035   
00036     //
00037     // Get options from the command line
00038     //
00039     
00040     Scalar       xt0         = 1.0;
00041     Scalar       xt1         = 1.0;
00042     Scalar       pt0         = 2.0;
00043     Scalar       pt1         = 0.0;
00044     Scalar       d           = 10.0;
00045     Scalar       x00         = 1.0;
00046     Scalar       x01         = 1.0;
00047     Scalar       p00         = 2.0;
00048     Scalar       p01         = 0.0;
00049     Scalar       pL0         = -1e+50;
00050     Scalar       pL1         = -1e+50;
00051     Scalar       pU0         = +1e+50;
00052     Scalar       pU1         = +1e+50;
00053     Scalar       q0          = 0.0;
00054 
00055     Scalar       xL0         = -1e+50;
00056     Scalar       xL1         = -1e+50;
00057     Scalar       xU0         = +1e+50;
00058     Scalar       xU1         = +1e+50;
00059 
00060     CommandLineProcessor  clp(false); // Don't throw exceptions
00061 
00062     lowsfCreator.setupCLP(&clp);
00063     solver.setupCLP(&clp);
00064 
00065     clp.setOption( "xt0", &xt0 );
00066     clp.setOption( "xt1", &xt1 );
00067     clp.setOption( "pt0", &pt0 );
00068     clp.setOption( "pt1", &pt1 );
00069     clp.setOption( "d", &d );
00070     clp.setOption( "x00", &x00 );
00071     clp.setOption( "x01", &x01 );
00072     clp.setOption( "p00", &p00 );
00073     clp.setOption( "q0", &q0 );
00074     clp.setOption( "p01", &p01 );
00075     clp.setOption( "pL0", &pL0 );
00076     clp.setOption( "pL1", &pL1 );
00077     clp.setOption( "pU0", &pU0 );
00078     clp.setOption( "pU1", &pU1 );
00079     clp.setOption( "xL0", &xL0 );
00080     clp.setOption( "xL1", &xL1 );
00081     clp.setOption( "xU0", &xU0 );
00082     clp.setOption( "xU1", &xU1 );
00083  
00084     CommandLineProcessor::EParseCommandLineReturn
00085       parse_return = clp.parse(argc,argv,&std::cerr);
00086 
00087     if( parse_return != CommandLineProcessor::PARSE_SUCCESSFUL )
00088       return parse_return;
00089 
00090     lowsfCreator.readParameters(out.get());
00091     solver.readParameters(out.get());
00092 
00093     //
00094     // Create the NLP
00095     //
00096     
00097     // Construct global and split (individual point) communicators
00098     int numPoints = 20;
00099     Teuchos::RCP<EpetraExt::MultiMpiComm> globalComm = 
00100       Teuchos::rcp(new EpetraExt::MultiMpiComm(MPI_COMM_WORLD, 1,  numPoints));
00101     int myPoints = globalComm->NumTimeStepsOnDomain();
00102     int myFirstPoint = globalComm->FirstTimeStepOnDomain();
00103 
00104     Teuchos::RCP<Epetra_MpiComm> epetra_comm =
00105       Teuchos::rcp_dynamic_cast<Epetra_MpiComm>(
00106         Teuchos::rcp(&globalComm->SubDomainComm(),false)
00107         );
00108 
00109     // Create the single-point EpetraExt::ModelEvaluator object
00110     Teuchos::RCP<EpetraMultiPointModelEval4DOpt>
00111       epetraModel = Teuchos::rcp(new EpetraMultiPointModelEval4DOpt(epetra_comm,xt0,xt1,pt0,pt1,d,x00,x01,p00,p01, q0));
00112     epetraModel->set_p_bounds(pL0,pL1,pU0,pU1);
00113     epetraModel->set_x_bounds(xL0,xL1,xU0,xU1);
00114 
00115     // Fill a vector of pointers to initial guesses for each point
00116     std::vector<Epetra_Vector*> multi_x_init(myPoints);
00117     Epetra_Vector init_vec(*(epetraModel->get_x_init().get()));
00118     for (int i=0; i<myPoints; i++) { multi_x_init[i] = &init_vec;}
00119 
00120     // Fill a vector of pointers to parameter vectors that define the multiple points
00121     Teuchos::RCP<std::vector< Teuchos::RCP<Epetra_Vector> > >  q_vec
00122         = Teuchos::rcp(new std::vector< Teuchos::RCP<Epetra_Vector> >(myPoints));
00123     for (int i=0; i<myPoints; i++) {
00124        q_vec->operator[](i) = Teuchos::rcp(new Epetra_Vector( *(epetraModel->get_p_map(1))));
00125        q_vec->operator[](i)->operator[](0) = 0.0 + 0.1*(i+myFirstPoint);
00126     }
00127 
00128     // Create the EpetraExt::MultiPointModelEvaluator object
00129     Teuchos::RCP<EpetraExt::MultiPointModelEvaluator>
00130        multiPointModel = rcp(new EpetraExt::MultiPointModelEvaluator(
00131                              epetraModel, globalComm, multi_x_init, q_vec));
00132 
00133     // Create the Thyra::EpetraModelEvaluator object
00134 
00135     Teuchos::RCP<Thyra::LinearOpWithSolveFactoryBase<double> >
00136       lowsFactory = lowsfCreator.createLinearSolveStrategy("");
00137 
00138     Teuchos::RCP<Thyra::EpetraModelEvaluator>
00139       epetraThyraModel = rcp(new Thyra::EpetraModelEvaluator());
00140     
00141     epetraThyraModel->initialize(multiPointModel,lowsFactory);
00142     
00143     //
00144     // Solve the NLP
00145     //
00146     
00147     // Set the model
00148     solver.setModel(epetraThyraModel);
00149 
00150     // Read the initial guess if one exists
00151     solver.readInitialGuess(out.get());
00152 
00153     // Solve the NLP
00154     const MoochoSolver::ESolutionStatus solution_status = solver.solve();
00155 
00156     // Write the parameters that where read
00157     lowsfCreator.writeParamsFile(*lowsFactory);
00158     solver.writeParamsFile();
00159     
00160     //
00161     // Return the solution status (0 if sucessfull)
00162     //
00163 
00164     return solution_status;
00165 
00166   }
00167   TEUCHOS_STANDARD_CATCH_STATEMENTS(true,*out,dummySuccess)
00168 
00169   return MoochoSolver::SOLVE_RETURN_EXCEPTION;
00170 
00171 }
00172 #else
00173 {
00174   cout << "NLPThyraEpetraModelEval4DOpt does nothing when HAVE_MPI is not defined" << endl;
00175   return 0;
00176 }
00177 #endif

Generated on Wed May 12 21:57:50 2010 for MOOCHO by  doxygen 1.4.7