diff --git a/CMakeLists.txt b/CMakeLists.txt index 8b2db6541..f356cc3c0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,7 +27,8 @@ option(HIOP_USE_EIGEN "Build with Eigen support" ON) option(HIOP_USE_MPI "Build with MPI support" ON) option(HIOP_USE_GPU "Build with support for GPUs - CUDA or HIP libraries" OFF) option(HIOP_TEST_WITH_BSUB "Use `jsrun` instead of `mpirun` commands when running tests" OFF) -option(HIOP_USE_RAJA "Build with portability abstraction library RAJA" OFF) +option(HIOP_USE_RAJA "Build with portability abstraction library RAJA" OFF) +option(HIOP_USE_AXOM "Build with AXOM to use Sidre for scalable checkpointing" OFF) option(HIOP_DEEPCHECKS "Extra checks and asserts in the code with a high penalty on performance" OFF) option(HIOP_WITH_KRON_REDUCTION "Build Kron Reduction code (requires UMFPACK)" OFF) option(HIOP_DEVELOPER_MODE "Build with extended warnings and options" OFF) @@ -289,6 +290,19 @@ if(HIOP_USE_RAJA) message(STATUS "Found umpire pkg-config: ${umpire_CONFIG}") endif() +if(HIOP_USE_AXOM) + if(HIOP_USE_MPI) + find_package(AXOM CONFIG + PATHS ${AXOM_DIR} ${AXOM_DIR}/lib/cmake/ + REQUIRED) + target_link_libraries(hiop_tpl INTERFACE axom) + message(STATUS "Found AXOM pkg-config: ${AXOM_CONFIG}") + elseif(HIOP_USE_MPI) + message(FATAL_ERROR "Error: HIOP_USE_MPI is required when HIOP_USE_AXOM is ON") + endif() +endif() + + if(HIOP_WITH_KRON_REDUCTION) set(HIOP_UMFPACK_DIR CACHE PATH "Path to UMFPACK directory") include(FindUMFPACK) diff --git a/doc/hiop_usermanual.pdf b/doc/hiop_usermanual.pdf index f69a6e790..ded095461 100644 Binary files a/doc/hiop_usermanual.pdf and b/doc/hiop_usermanual.pdf differ diff --git a/doc/src/sections/solver_options.tex b/doc/src/sections/solver_options.tex index 64057a221..a2a73be42 100755 --- a/doc/src/sections/solver_options.tex +++ b/doc/src/sections/solver_options.tex @@ -423,6 +423,16 @@ \subsubsection{Problem preprocessing} \medskip +\subsubsection{Checkpointing of the solver state and restarting}\label{sec:checkpoint} +As detailed in Section~\ref{sec:checkpoint_API}, \Hi can save/load its internal state to/from disk. All the options in this section require an Axom-enabled build (use ``-DHIOP\_USE\_AXOM=ON'' with cmake) and are supported only by the quasi-Newton IPM solver (\texttt{hiopAlgFilterIPMQuasiNewton} class) for the \texttt{hiopInterfaceDenseConstraints} NLP formulation/interface. + +\noindent \textbf{checkpoint\_save}: Save state of NLP solver to file indicated by value of option ``checkpoint\_file''. String values ``yes'' or ``no'', default ``no''. + +\noindent \textbf{checkpoint\_load\_on\_start} On (re)start the NLP solver will load checkpoint file specified by ``checkpoint\_file`` option. String values ``yes'' or ``no'', default ``no''. + +\noindent \textbf{checkpoint\_file} Path to checkpoint file to load from or save to. If present, the character ``\#'' is replaced with the iteration number at which the checkpointing is saved (but \textit{not} when loaded). \Hi adds a ``.root'' extension internally if the value of the option is a directory. If this option is not specified and loading or saving checkpoints is enabled, \Hi will use a file named ``hiop\_state\_chk''. + +\noindent \textbf{checkpoint\_save\_every\_N\_iter} Iteration frequency of saving checkpoints to disk if ``checkpoint\_save'' is ``yes''. Takes positive integer values with a default value $10$. \subsubsection{Miscellaneous options} diff --git a/doc/src/techrep_main.tex b/doc/src/techrep_main.tex index 2edb3f11b..f16f36fe7 100755 --- a/doc/src/techrep_main.tex +++ b/doc/src/techrep_main.tex @@ -133,7 +133,7 @@ \vspace{3cm} {\huge\bfseries \Hi\ -- User Guide} \\[14pt] - {\large\bfseries version 1.03} + {\large\bfseries version 1.1.0} \vspace{3cm} @@ -155,7 +155,7 @@ \vspace{4.75cm} \textcolor{violet}{{\large\bfseries Oct 15, 2017} \\ -{\large\bfseries Updated Feb 5, 2024}} +{\large\bfseries Updated Sept 22, 2024}} \vspace{0.75cm} @@ -474,6 +474,29 @@ \subsubsection{Calling \Hi for a \texttt{hiopInterfaceDenseConstraints} formulat \end{lstlisting} The standalone drivers \texttt{NlpDenseConsEx1}, \texttt{NlpDenseConsEx2}, and \texttt{NlpDenseConsEx3} inside directory \texttt{src/Drivers/} under the \Hi's root directory contain more detailed examples of the use of \Hi. +\subsubsection{Checkpointing}\label{sec:checkpoint_API} +File checkpointing is available for \Hi's quasi-Newton IPM solver, which is used exclusively to solve \texttt{hiopInterfaceDenseConstraints} formulation. This can be helpful when running a job on +a cluster that enforces limits on the job’s running time. +Later, this feature will also be provided for other solvers, such as the Newton IPM (used exclusively with sparse NLP) and HiOp-PriDec. + +The checkpointing I/O is based on Axom's scalable Sidre data manager (see \url{https://axom.readthedocs.io/en/develop/axom/sidre/docs/sphinx/index.html} for more information) and, thus, requires an Axom-enabled build (use ``-DHIOP\_USE\_AXOM=ON'' with cmake). + +There are two ways to use \Hi's checkpointing. The first is via the quasi-Newton solver's API, namely, the methods +\begin{lstlisting} +void load_state_from_sidre_group(const ::axom::sidre::Group& group); +void save_state_to_sidre_group(::axom::sidre::Group& group); +\end{lstlisting} +of \texttt{hiopAlgFilterIPMQuasiNewton} solver class. New Sidre views will be created (or reused) within the group passed as argument to load / save state variables of the quasi-Newton solver. Alternatively, \texttt{hiopAlgFilterIPMQuasiNewton} solver class offers similar methods to work directly with a file, namely, +\begin{lstlisting} +bool load_state_from_file(const ::std::string& path) noexcept; +bool save_state_to_file(const ::std::string& path) noexcept; +\end{lstlisting} +These two methods will create the Sidre group internally and checkpoint to/from it using the first two methods. + +A second avenue to checkpoint is via user options. This is detailed in Section~\ref{sec:checkpoint}. + + \warningcp{Note:} A couple of particularities stemming from the use of Sidre must be acknowledged. First, a checkpoint file should be loaded using HiOp with the same number of MPI ranks as when it was saved. Second, checkpointing is not available for non-MPI builds due to Axom having MPI as a dependency. Finally, when loading from or saving to a checkpoint file, the sizes of the file's variables (Sidre views) must match the sizes of the HiOp variables to which the data is loaded or saved, meaning \Hi will throw an exception if an existing file is (re)used to load or save a algorithm state for a problem that changed sizes since the file was created. + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% NLP Sparse diff --git a/src/Drivers/Dense/NlpDenseConsEx1.cpp b/src/Drivers/Dense/NlpDenseConsEx1.cpp index df97ee21e..cd0a31676 100644 --- a/src/Drivers/Dense/NlpDenseConsEx1.cpp +++ b/src/Drivers/Dense/NlpDenseConsEx1.cpp @@ -4,6 +4,14 @@ #include #include +#ifdef HIOP_USE_AXOM +#include +#include +#include +#include +using namespace axom; +#endif + using namespace hiop; Ex1Meshing1D::Ex1Meshing1D(double a, double b, size_type glob_n, double r, MPI_Comm comm_) @@ -178,10 +186,59 @@ void DiscretizedFunction::setFunctionValue(index_type i_global, const double& va this->data_[i_local]=value; } - - /* DenseConsEx1 class implementation */ +bool DenseConsEx1::iterate_callback(int iter, + double obj_value, + double logbar_obj_value, + int n, + const double* x, + const double* z_L, + const double* z_U, + int m_ineq, + const double* s, + int m, + const double* g, + const double* lambda, + double inf_pr, + double inf_du, + double onenorm_pr, + double mu, + double alpha_du, + double alpha_pr, + int ls_trials) +{ +#ifdef HIOP_USE_AXOM + //save state to sidre::Group every 5 iterations if a solver/algorithm object was provided + if(iter > 0 && (iter % 5 == 0) && nullptr!=solver_) { + // + //Example of how to save HiOp state to axom::sidre::Group + // + + //We first manufacture a Group. User code supposedly already has one. + sidre::DataStore ds; + sidre::Group* group = ds.getRoot()->createGroup("HiOp quasi-Newton alg state"); + + //the actual saving of state to group + try { + solver_->save_state_to_sidre_group(*group); + } catch(std::runtime_error& e) { + //user chooses action when an error occured in saving the state... + //we choose to stop HiOp + return false; + } + + //User code can further inspect the Group or add addtl info to DataStore, with the end goal + //of saving it to file before HiOp starts next iteration. Here we just save it. + sidre::IOManager writer(comm); + int n_files; + MPI_Comm_size(comm, &n_files); + writer.write(ds.getRoot(), n_files, "hiop_state_ex1", sidre::Group::getDefaultIOProtocol()); + } +#endif + return true; +} + /*set c to * c(t) = 1-10*t, for 0<=t<=1/10, * 0, for 1/10<=t<=1. diff --git a/src/Drivers/Dense/NlpDenseConsEx1.hpp b/src/Drivers/Dense/NlpDenseConsEx1.hpp index d9b0235f7..967255f04 100644 --- a/src/Drivers/Dense/NlpDenseConsEx1.hpp +++ b/src/Drivers/Dense/NlpDenseConsEx1.hpp @@ -78,7 +78,7 @@ class Ex1Meshing1D MPI_Comm comm; int my_rank, comm_size; index_type* col_partition; - + friend class DiscretizedFunction; private: @@ -112,7 +112,9 @@ class DenseConsEx1 : public hiop::hiopInterfaceDenseConstraints { public: DenseConsEx1(int n_mesh_elem=100, double mesh_ratio=1.0) - : n_vars(n_mesh_elem), comm(MPI_COMM_WORLD) + : n_vars(n_mesh_elem), + comm(MPI_COMM_WORLD), + solver_(nullptr) { //create the members _mesh = new Ex1Meshing1D(0.0,1.0, n_vars, mesh_ratio, comm); @@ -218,6 +220,31 @@ class DenseConsEx1 : public hiop::hiopInterfaceDenseConstraints } return true; } + + inline void set_solver(hiop::hiopAlgFilterIPM* alg_obj) + { + solver_ = alg_obj; + } + + bool iterate_callback(int iter, + double obj_value, + double logbar_obj_value, + int n, + const double* x, + const double* z_L, + const double* z_U, + int m_ineq, + const double* s, + int m, + const double* g, + const double* lambda, + double inf_pr, + double inf_du, + double onenorm_pr, + double mu, + double alpha_du, + double alpha_pr, + int ls_trials); private: int n_vars; MPI_Comm comm; @@ -228,6 +255,10 @@ class DenseConsEx1 : public hiop::hiopInterfaceDenseConstraints DiscretizedFunction* c; DiscretizedFunction* x; //proxy for taking hiop's variable in and working with it as a function + /// Pointer to the solver, to be used to checkpoint + hiop::hiopAlgFilterIPM* solver_; + +private: //populates the linear term c void set_c(); }; diff --git a/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp b/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp index 4fbebc4b0..4bd5c5c1b 100644 --- a/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp +++ b/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp @@ -7,10 +7,23 @@ #include #include +#ifdef HIOP_USE_AXOM +#include +#include +#include +#include +using namespace axom; +#endif + + using namespace hiop; static bool self_check(size_type n, double obj_value); - +#ifdef HIOP_USE_AXOM +static bool do_load_checkpoint_test(const size_type& mesh_size, + const double& ratio, + const double& obj_val_expected); +#endif static bool parse_arguments(int argc, char **argv, size_type& n, double& distortion_ratio, bool& self_check) { n = 20000; distortion_ratio=1.; self_check=false; //default options @@ -67,24 +80,27 @@ int main(int argc, char **argv) err = MPI_Init(&argc, &argv); assert(MPI_SUCCESS==err); err = MPI_Comm_rank(MPI_COMM_WORLD,&rank); assert(MPI_SUCCESS==err); err = MPI_Comm_size(MPI_COMM_WORLD,&numRanks); assert(MPI_SUCCESS==err); - if(0==rank) printf("Support for MPI is enabled\n"); + if(0==rank) { + printf("Support for MPI is enabled\n"); + } #endif - bool selfCheck; size_type mesh_size; double ratio; - if(!parse_arguments(argc, argv, mesh_size, ratio, selfCheck)) { usage(argv[0]); return 1;} - + bool selfCheck; + size_type mesh_size; + double ratio; + double objective = 0.; + if(!parse_arguments(argc, argv, mesh_size, ratio, selfCheck)) { + usage(argv[0]); + return 1; + } + DenseConsEx1 problem(mesh_size, ratio); - //if(rank==0) printf("interface created\n"); hiop::hiopNlpDenseConstraints nlp(problem); - //if(rank==0) printf("nlp formulation created\n"); - - //nlp.options->SetIntegerValue("verbosity_level", 4); - //nlp.options->SetNumericValue("tolerance", 1e-4); - //nlp.options->SetStringValue("duals_init", "zero"); - //nlp.options->SetIntegerValue("max_iter", 2); hiop::hiopAlgFilterIPM solver(&nlp); + problem.set_solver(&solver); + hiop::hiopSolveStatus status = solver.run(); - double objective = solver.getObjective(); + objective = solver.getObjective(); //this is used for testing when the driver is called with -selfcheck if(selfCheck) { @@ -97,7 +113,19 @@ int main(int argc, char **argv) } } - if(0==rank) printf("Objective: %18.12e\n", objective); + if(0==rank) { + printf("Objective: %18.12e\n", objective); + } + +#ifdef HIOP_USE_AXOM + // example/test for HiOp's load checkpoint API. + if(!do_load_checkpoint_test(mesh_size, ratio, objective)) { + if(rank==0) { + printf("Load checkpoint and restart test failed."); + } + return -1; + } +#endif #ifdef HIOP_USE_MPI MPI_Finalize(); #endif @@ -134,3 +162,56 @@ static bool self_check(size_type n, double objval) return true; } + +#ifdef HIOP_USE_AXOM +/** + * An illustration on how to use load_state_from_sidre_group API method of HiOp's algorithm class. + * + * + */ +static bool do_load_checkpoint_test(const size_type& mesh_size, + const double& ratio, + const double& obj_val_expected) +{ + //Pretend this is new job and recreate the HiOp objects. + DenseConsEx1 problem(mesh_size, ratio); + hiop::hiopNlpDenseConstraints nlp(problem); + + hiop::hiopAlgFilterIPM solver(&nlp); + + // + // example of how to use load_state_sidre_group to warm-start + // + + //Supposedly, the user code should have the group in hand before asking HiOp to load from it. + //We will manufacture it by loading a sidre checkpoint file. Here the checkpoint file + //"hiop_state_ex1.root" was created from the interface class' iterate_callback method + //(saved every 5 iterations) + sidre::DataStore ds; + + try { + sidre::IOManager reader(MPI_COMM_WORLD); + reader.read(ds.getRoot(), "hiop_state_ex1.root", false); + } catch(std::exception& e) { + printf("Failed to read checkpoint file. Error: [%s]", e.what()); + return false; + } + + + //the actual API call + try { + const sidre::Group* group = ds.getRoot()->getGroup("HiOp quasi-Newton alg state"); + solver.load_state_from_sidre_group(*group); + } catch(std::runtime_error& e) { + printf("Failed to load from sidre::group. Error: [%s]", e.what()); + return false; + } + + hiop::hiopSolveStatus status = solver.run(); + double obj_val = solver.getObjective(); + if(obj_val != obj_val_expected) { + return false; + } + return true; +} +#endif // HIOP_USE_AXOM diff --git a/src/Interface/hiopInterface.hpp b/src/Interface/hiopInterface.hpp index f0b6ceced..47c54d3ca 100644 --- a/src/Interface/hiopInterface.hpp +++ b/src/Interface/hiopInterface.hpp @@ -467,8 +467,8 @@ class hiopInterfaceBase } /** - * This method is used to provide an user all the hiop iterate - * procedure. @see solution_callback() for an explanation of the parameters. + * This method is used to provide user all the internal hiop iterates. @see solution_callback() + * for an explanation of the parameters. * * @param[in] x array of (local) entries of the primal variables (managed by Umpire, see note below) * @param[in] z_L array of (local) entries of the dual variables for lower bounds (managed by Umpire, see note below) @@ -496,7 +496,7 @@ class hiopInterfaceBase { return true; } - + /** * A wildcard function used to change the primal variables. * diff --git a/src/Interface/hiop_defs.hpp.in b/src/Interface/hiop_defs.hpp.in index 0af21826d..68570c8b2 100644 --- a/src/Interface/hiop_defs.hpp.in +++ b/src/Interface/hiop_defs.hpp.in @@ -12,6 +12,7 @@ #cmakedefine HIOP_USE_PARDISO #cmakedefine HIOP_USE_RESOLVE #cmakedefine HIOP_USE_GINKGO +#cmakedefine HIOP_USE_AXOM #define HIOP_VERSION "@PROJECT_VERSION@" #define HIOP_VERSION_MAJOR @PROJECT_VERSION_MAJOR@ #define HIOP_VERSION_MINOR @PROJECT_VERSION_MINOR@ diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index 2ac5fcc83..973c6d475 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -66,6 +66,11 @@ #include "hiopCppStdUtils.hpp" +#ifdef HIOP_USE_AXOM +#include "SidreHelper.hpp" +using namespace axom; +#endif + #include #include #include @@ -76,15 +81,23 @@ namespace hiop { hiopAlgFilterIPMBase::hiopAlgFilterIPMBase(hiopNlpFormulation* nlp_in, const bool within_FR) - : soc_dir(nullptr), - onenorm_pr_curr_{0.0}, - c_soc(nullptr), - d_soc(nullptr), - within_FR_{within_FR}, - pd_perturb_{nullptr}, - fact_acceptor_{nullptr} + : nlp(nlp_in), + logbar(nullptr), + it_curr(nullptr), + it_trial(nullptr), + dir(nullptr), + soc_dir(nullptr), + resid(nullptr), + resid_trial(nullptr), + iter_num_(0), + iter_num_total_(0), + onenorm_pr_curr_(0.0), + c_soc(nullptr), + d_soc(nullptr), + within_FR_(within_FR), + pd_perturb_(nullptr), + fact_acceptor_(nullptr) { - nlp = nlp_in; //force completion of the nlp's initialization nlp->finalizeInitialization(); } @@ -915,7 +928,8 @@ void hiopAlgFilterIPMBase::displayTerminationMsg() /////////////////////////////////////////////////////////////////////////////////////////////////// hiopAlgFilterIPMQuasiNewton::hiopAlgFilterIPMQuasiNewton(hiopNlpDenseConstraints* nlp_in, const bool within_FR) - : hiopAlgFilterIPMBase(nlp_in, within_FR) + : hiopAlgFilterIPMBase(nlp_in, within_FR), + load_state_api_called_(false) { nlpdc = nlp_in; reload_options(); @@ -976,8 +990,48 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() nlp->runStats.tmOptimizTotal.start(); - startingProcedure(*it_curr, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); //this also evaluates the nlp - _mu=mu0; + iter_num_ = 0; + + // + // starting point: + // - user provided (with slack adjustments and lsq eq. duals initialization) + // - load checkpoint API (method load_state_from_sidre_group) called before calling this method + // - checkpoint from file (option "checkpoint_load_on_start") + // + if(nlp->options->GetString("checkpoint_load_on_start") != "yes" && !load_state_api_called_) { + //this also evaluates the nlp + startingProcedure(*it_curr, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); + _mu=mu0; + iter_num_total_ = 0; + } else { + if(!load_state_api_called_) { + // + //checkpoint load from file + // +#ifdef HIOP_USE_AXOM + //load from file: will populate it_curr, _Hess_lagr, and algorithmic parameters + auto chkpnt_ok = load_state_from_file(nlp->options->GetString("checkpoint_file")); + if(!chkpnt_ok) { + nlp->log->printf(hovWarning, "Using default starting procedure (no checkpoint load!).\n"); + iter_num_total_ = 0; + //fall back on the default starting procedure (it also evaluates the nlp) + startingProcedure(*it_curr, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); + _mu=mu0; + iter_num_total_ = 0; + } +#else + nlp->log->printf(hovWarning, + "Unexpected checkpoint misconfiguration. " + "Will use user-provided starting point.\n"); +#endif + } + //additionally: need to evaluate the nlp + if(!this->evalNlp_noHess(*it_curr, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d)) { + nlp->log->printf(hovError, "Failure in evaluating user NLP functions at loaded checkpoint."); + return Error_In_User_Function; + } + solver_status_ = NlpSolve_SolveNotCalled; + } //update log bar logbar->updateWithNlpInfo(*it_curr, _mu, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); @@ -987,7 +1041,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() nlp->log->write("First residual-------------", *resid, hovIteration); - iter_num=0; nlp->runStats.nIter=iter_num; + nlp->runStats.nIter = iter_num_; bool disableLS = nlp->options->GetString("accept_every_trial_step")=="yes"; theta_max = theta_max_fact_*fmax(1.0,resid->get_theta()); @@ -1074,7 +1128,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() } //user callback - if(!nlp->user_callback_iterate(iter_num, + if(!nlp->user_callback_iterate(iter_num_, _f_nlp, logbar->f_logbar, *it_curr->get_x(), @@ -1095,10 +1149,15 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() solver_status_ = User_Stopped; break; } +#ifdef HIOP_USE_AXOM + //checkpointing - based on options provided by the user + checkpointing_stuff(); +#endif + /************************************************* * Termination check ************************************************/ - if(checkTermination(_err_nlp, iter_num, solver_status_)) { + if(checkTermination(_err_nlp, iter_num_, solver_status_)) { break; } if(NlpSolve_Pending!=solver_status_) break; //failure of the line search or user stopped. @@ -1112,7 +1171,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() if(!mu_updated) { break; } - nlp->log->printf(hovScalars, "Iter[%d] barrier params reduced: mu=%g tau=%g\n", iter_num, _mu, _tau); + nlp->log->printf(hovScalars, "Iter[%d] barrier params reduced: mu=%g tau=%g\n", iter_num_, _mu, _tau); //update only logbar problem and residual (the NLP didn't change) logbar->updateWithNlpInfo(*it_curr, _mu, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); @@ -1143,7 +1202,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() break; } } - nlp->log->printf(hovScalars, "Iter[%d] logbarObj=%23.17e (mu=%12.5e)\n", iter_num, logbar->f_logbar,_mu); + nlp->log->printf(hovScalars, "Iter[%d] logbarObj=%23.17e (mu=%12.5e)\n", iter_num_, logbar->f_logbar,_mu); /**************************************************** * Search direction calculation ***************************************************/ @@ -1152,7 +1211,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() nlp->runStats.kkt.start_optimiz_iteration(); - //first update the Hessian and kkt system + //update the Hessian and kkt system Hess->update(*it_curr,*_grad_f,*_Jac_c,*_Jac_d); if(!kkt->update(it_curr, _grad_f, _Jac_c, _Jac_d, _Hess_Lagr)) { nlp->log->write("Unrecoverable error in step computation (factorization) [1]. Will exit here.", @@ -1164,7 +1223,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() auto* fact_acceptor_ic = dynamic_cast (fact_acceptor_); if(fact_acceptor_ic) { //compute_search_direction call below updates linsol safe mode flag and linsol_safe_mode_lastiter - if(!compute_search_direction(kkt, linsol_safe_mode_on, linsol_forcequick, iter_num)) { + if(!compute_search_direction(kkt, linsol_safe_mode_on, linsol_forcequick, iter_num_)) { delete kkt; return solver_status_ = Err_Step_Computation; } @@ -1172,7 +1231,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() auto* fact_acceptor_dwd = dynamic_cast (fact_acceptor_); assert(fact_acceptor_dwd); //compute_search_direction call below updates linsol safe mode flag and linsol_safe_mode_lastiter - if(!compute_search_direction_inertia_free(kkt, linsol_safe_mode_on, linsol_forcequick, iter_num)) { + if(!compute_search_direction_inertia_free(kkt, linsol_safe_mode_on, linsol_forcequick, iter_num_)) { //it failed under safe mode delete kkt; return solver_status_ = Err_Step_Computation; @@ -1184,7 +1243,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() nlp->log->printf(hovSummary, "%s", nlp->runStats.kkt.get_summary_last_iter().c_str()); } - nlp->log->printf(hovIteration, "Iter[%d] full search direction -------------\n", iter_num); + nlp->log->printf(hovIteration, "Iter[%d] full search direction -------------\n", iter_num_); nlp->log->write("", *dir, hovIteration); /*************************************************************** * backtracking line search @@ -1325,7 +1384,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() //post line-search stuff //filter is augmented whenever the switching condition or Armijo rule do not hold for the trial point that was just accepted - if(nlp->options->GetString("force_resto")=="yes" && !within_FR_ && iter_num == 1) { + if(nlp->options->GetString("force_resto")=="yes" && !within_FR_ && iter_num_ == 1) { use_fr = apply_feasibility_restoration(kkt); if(use_fr) { // continue iterations if FR is accepted @@ -1385,10 +1444,14 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() } nlp->log->printf(hovScalars, - "Iter[%d] -> accepted step primal=[%17.11e] dual=[%17.11e]\n", iter_num, _alpha_primal, + "Iter[%d] -> accepted step primal=[%17.11e] dual=[%17.11e]\n", + iter_num_, + _alpha_primal, _alpha_dual); - iter_num++; - nlp->runStats.nIter=iter_num; + nlp->log->printf(hovScalars, "Finished iter[%d] global iter[%d]\n.", iter_num_, iter_num_total_); + iter_num_++; + iter_num_total_++; + nlp->runStats.nIter = iter_num_; // fr problem has already updated dual, slacks and NLP functions if(!use_fr) { @@ -1420,7 +1483,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() //update current iterate (do a fast swap of the pointers) hiopIterate* pit=it_curr; it_curr=it_trial; it_trial=pit; - nlp->log->printf(hovIteration, "Iter[%d] -> full iterate:", iter_num); + nlp->log->printf(hovIteration, "Iter[%d] -> full iterate:", iter_num_); nlp->log->write("", *it_curr, hovIteration); nlp->runStats.tmSolverInternal.stop(); //----- @@ -1431,7 +1494,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() //update residual resid->update(*it_curr,_f_nlp, *_c, *_d,*_grad_f,*_Jac_c,*_Jac_d, *logbar); - nlp->log->printf(hovIteration, "Iter[%d] full residual:-------------\n", iter_num); + nlp->log->printf(hovIteration, "Iter[%d] full residual:-------------\n", iter_num_); nlp->log->write("", *resid, hovIteration); } @@ -1457,12 +1520,12 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() void hiopAlgFilterIPMQuasiNewton::outputIteration(int lsStatus, int lsNum, int use_soc, int use_fr) { - if(iter_num/10*10==iter_num) + if(iter_num_/10*10==iter_num_) nlp->log->printf(hovSummary, "iter objective inf_pr inf_du lg(mu) alpha_du alpha_pr linesrch\n"); if(lsStatus==-1) nlp->log->printf(hovSummary, "%4d %14.7e %7.3e %7.3e %6.2f %7.3e %7.3e -(-)\n", - iter_num, _f_nlp/nlp->get_obj_scale(), _err_nlp_feas, _err_nlp_optim, + iter_num_total_, _f_nlp/nlp->get_obj_scale(), _err_nlp_feas, _err_nlp_optim, log10(_mu), _alpha_dual, _alpha_primal); else { char stepType[2]; @@ -1480,11 +1543,248 @@ void hiopAlgFilterIPMQuasiNewton::outputIteration(int lsStatus, int lsNum, int u } nlp->log->printf(hovSummary, "%4d %14.7e %7.3e %7.3e %6.2f %7.3e %7.3e %d(%s)\n", - iter_num, _f_nlp/nlp->get_obj_scale(), _err_nlp_feas, _err_nlp_optim, + iter_num_total_, _f_nlp/nlp->get_obj_scale(), _err_nlp_feas, _err_nlp_optim, log10(_mu), _alpha_dual, _alpha_primal, lsNum, stepType); } } +#ifdef HIOP_USE_AXOM + +bool hiopAlgFilterIPMQuasiNewton::save_state_to_file(const ::std::string& path) noexcept +{ + try { + sidre::DataStore ds; + sidre::Group* group = ds.getRoot()->createGroup("HiOp quasi-Newton alg state"); + this->save_state_to_sidre_group(*group); + + sidre::IOManager writer(this->get_nlp()->get_comm()); + int n_files; + MPI_Comm_size(this->get_nlp()->get_comm(), &n_files); + writer.write(ds.getRoot(), n_files, path, sidre::Group::getDefaultIOProtocol()); + return true; + } catch(const std::exception& exp) { + nlp->log->printf(hovError, "Error when saving checkpoint to file '%s'\n", path.c_str()); + nlp->log->printf(hovError, " Addtl info: %s\n", exp.what()); + return false; + } +} + +bool hiopAlgFilterIPMQuasiNewton::load_state_from_file(const ::std::string& path) noexcept +{ + try { + sidre::DataStore ds; + + sidre::IOManager reader(this->get_nlp()->get_comm()); + auto path2 = SidreHelper::check_path(path); + reader.read(ds.getRoot(), path2, false); + nlp->log->printf(hovScalars, "Loaded checkpoint file [%s].\n", path2.c_str()); + + const sidre::Group* group = ds.getRoot()->getGroup("HiOp quasi-Newton alg state"); + this->load_state_from_sidre_group(*group); + return true; + } catch(const std::exception& exp) { + nlp->log->printf(hovError, "Error in loading checkpoint from file '%s'\n", path.c_str()); + nlp->log->printf(hovError, " Addtl info: %s\n", exp.what()); + return false; + } +} + +void hiopAlgFilterIPMQuasiNewton::save_state_to_sidre_group(::axom::sidre::Group& group) +{ + using IndType = sidre::IndexType; + + //iterate state + //create views for each member that needs to be saved + SidreHelper::copy_iterate_to_views(group, "alg_iterate_", *it_curr); + + //state of quasi-Newton Hessian approximation + hiopHessianLowRank& hqn = dynamic_cast(*_Hess_Lagr); + const double hqn_params[] = {(double)hqn.l_max, + (double)hqn.l_curr, + hqn.sigma, + hqn.sigma0, + (double)hqn.matrixChanged}; + const size_type nhqn_params = sizeof(hqn_params) / sizeof(double); + SidreHelper::copy_array_to_view(group, "Hess_quasiNewton_params", hqn_params, nhqn_params); + + //quasi-Newton Hessian stores the previous iterate and corresponding derivatives + SidreHelper::copy_iterate_to_views(group, "Hess_quasiNewton_prev_iter", *hqn.it_prev_); + SidreHelper::copy_vec_to_view(group, "Hess_quasiNewton_prev_grad", *hqn.grad_f_prev_); + + auto* Jac_c = hqn.Jac_c_prev_; + SidreHelper::copy_array_to_view(group, + "Hess_quasiNewton_prev_Jacc", + Jac_c->local_data_const(), + Jac_c->get_local_size_n() * Jac_c->get_local_size_m()); + + + auto* Jac_d = hqn.Jac_d_prev_; + SidreHelper::copy_array_to_view(group, + "Hess_quasiNewton_prev_Jacd", + Jac_d->local_data_const(), + Jac_d->get_local_size_n() * Jac_d->get_local_size_m()); + + //quasi-Newton Hessian internal states + //memory matrices and internal representation + SidreHelper::copy_array_to_view(group, + "Hess_quasiNewton_St", + hqn.St_->local_data_const(), + hqn.St_->get_local_size_n() * hqn.St_->get_local_size_m()); + SidreHelper::copy_array_to_view(group, + "Hess_quasiNewton_Yt", + hqn.Yt_->local_data_const(), + hqn.Yt_->get_local_size_n() * hqn.Yt_->get_local_size_m()); + SidreHelper::copy_vec_to_view(group, "Hess_quasiNewton_D", *hqn.D_); + SidreHelper::copy_array_to_view(group, + "Hess_quasiNewton_L", + hqn.L_->local_data_const(), + hqn.L_->get_local_size_n() * hqn.L_->get_local_size_m()); + + + //algorithmic parameters for this state + //mu, iteration number, num MPI ranks + int nranks=1; +#ifdef HIOP_USE_MPI + MPI_Comm_size(get_nlp()->get_comm(), &nranks); +#endif + constexpr double version = HIOP_VERSION_MAJOR*100 + HIOP_VERSION_MINOR*10 + HIOP_VERSION_PATCH; + const double alg_params[] = {_mu, (double)iter_num_total_, (double)nranks, version}; + const size_type nparams = sizeof(alg_params) / sizeof(double); + SidreHelper::copy_array_to_view(group, "alg_params", alg_params, nparams); + + nlp->log->printf(hovScalars, + "Saved checkpoint to sidre::Group : ver %d.%d.%d on %d MPI ranks at " + "mu=%12.5e from iter=%d.\n", + HIOP_VERSION_MAJOR, + HIOP_VERSION_MINOR, + HIOP_VERSION_PATCH, + nranks, + _mu, + iter_num_total_); + +} + +void hiopAlgFilterIPMQuasiNewton::load_state_from_sidre_group(const sidre::Group& group) +{ + load_state_api_called_ = true; + + // + //algorithmic parameters + // + //!!!note: nparams needs to match the nparams from save_state_to_data_store + const int nparams = 4; + double alg_params[nparams]; + SidreHelper::copy_array_from_view(group, "alg_params", alg_params, nparams); + //!!! dev note: match order in save_state_to_data_store + _mu = alg_params[0]; + iter_num_total_ = alg_params[1]; + + int nranks=1; +#ifdef HIOP_USE_MPI + MPI_Comm_size(get_nlp()->get_comm(), &nranks); +#endif + if( (int)alg_params[2] != nranks ) { + ::std::stringstream ss; + ss << "Mismatch in the number of MPI ranks used to checkpoint. Checkpointing was " << + "done on " << (int)alg_params[2] << " ranks while HiOp currently runs on " << + nranks << " ranks.\n"; + throw std::runtime_error(ss.str()); + } + const int ver_major = ((int)alg_params[3] / 100); + const int ver_minor = ((int)alg_params[3] - ver_major*100)/10; + const int ver_patch = (int)alg_params[3] - ver_major*100 - ver_minor*10; + nlp->log->printf(hovScalars, + "Loaded checkpoint from sidre::Group. Found ver %d.%d.%d on %d MPI ranks at " + "mu=%12.5e from iter=%d.\n", + ver_major, + ver_minor, + ver_patch, + nranks, + _mu, + iter_num_total_); + + // + //iterate states + // + SidreHelper::copy_iterate_from_views(group, "alg_iterate_", *it_curr); + + // + //state of quasi-Newton Hessian approximation + // + hiopHessianLowRank& hqn = dynamic_cast(*_Hess_Lagr); + //!!!note: nparams needs to match the # of params from save_state_to_sidre_group + const int nhqn_params = 5; + double hqn_params[nhqn_params]; + SidreHelper::copy_array_from_view(group, "Hess_quasiNewton_params", hqn_params, nhqn_params); + + const size_type lim_mem_length = (size_type) hqn_params[1]; + //ensure the internals are allocated for this mem length + hqn.alloc_for_limited_mem(lim_mem_length); + + hqn.l_max = (size_type) hqn_params[0]; + hqn.l_curr = lim_mem_length; + hqn.sigma = hqn_params[2]; + hqn.sigma0 = hqn_params[3]; + hqn.matrixChanged = hqn_params[4]; + + assert(hqn.it_prev_); + //quasi-Newton Hessian stores the previous iterate and corresponding derivatives + SidreHelper::copy_iterate_from_views(group, "Hess_quasiNewton_prev_iter", *hqn.it_prev_); + SidreHelper::copy_vec_from_view(group, "Hess_quasiNewton_prev_grad", *hqn.grad_f_prev_); + + auto* Jac_c = hqn.Jac_c_prev_; + SidreHelper::copy_array_from_view(group, + "Hess_quasiNewton_prev_Jacc", + Jac_c->local_data(), + Jac_c->get_local_size_n() * Jac_c->get_local_size_m()); + + auto* Jac_d = hqn.Jac_d_prev_; + SidreHelper::copy_array_from_view(group, + "Hess_quasiNewton_prev_Jacd", + Jac_d->local_data(), + Jac_d->get_local_size_n() * Jac_d->get_local_size_m()); + + //quasi-Newton Hessian internal states + //memory matrices + SidreHelper::copy_array_from_view(group, + "Hess_quasiNewton_St", + hqn.St_->local_data(), + hqn.St_->get_local_size_n() * hqn.St_->get_local_size_m()); + SidreHelper::copy_array_from_view(group, + "Hess_quasiNewton_Yt", + hqn.Yt_->local_data(), + hqn.Yt_->get_local_size_n() * hqn.Yt_->get_local_size_m()); + SidreHelper::copy_vec_from_view(group, "Hess_quasiNewton_D", *hqn.D_); + SidreHelper::copy_array_from_view(group, + "Hess_quasiNewton_L", + hqn.L_->local_data(), + hqn.L_->get_local_size_n() * hqn.L_->get_local_size_m()); +} + +void hiopAlgFilterIPMQuasiNewton::checkpointing_stuff() +{ + if(nlp->options->GetString("checkpoint_save")=="no") { + return; + } + int chk_every_N = nlp->options->GetInteger("checkpoint_save_every_N_iter"); + //check iteration + if(iter_num_>0 && iter_num_ % chk_every_N==0) { + using ::std::string; + // replace "#" in checkpointing file with iteration number + string path = nlp->options->GetString("checkpoint_file"); + const auto s_it_num = ::std::to_string(iter_num_); + auto pos = path.find("#"); + while(string::npos != pos) { + path.replace(pos, 1, s_it_num); + pos = path.find("#", pos); + } + + nlp->log->printf(hovSummary, "Saving checkpoint at iter %d in '%s'.\n", iter_num_, path.c_str()); + //actual checkpointing via axom::sidre + save_state_to_file(path); + } +} +#endif // HIOP_USE_AXOM /****************************************************************************************************** * FULL NEWTON IPM @@ -1850,7 +2150,9 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() nlp->log->write("First residual-------------", *resid, hovIteration); - iter_num=0; nlp->runStats.nIter=iter_num; + iter_num_ = 0; + iter_num_total_ = 0; + nlp->runStats.nIter = iter_num_; bool disableLS = nlp->options->GetString("accept_every_trial_step")=="yes"; theta_max = theta_max_fact_*fmax(1.0,resid->get_theta()); @@ -1952,7 +2254,8 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() } //user callback - if(!nlp->user_callback_iterate(iter_num, _f_nlp, + if(!nlp->user_callback_iterate(iter_num_, + _f_nlp, logbar->f_logbar, *it_curr->get_x(), *it_curr->get_zl(), @@ -1975,7 +2278,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() /************************************************* * Termination check ************************************************/ - if(checkTermination(_err_nlp, iter_num, solver_status_)) { + if(checkTermination(_err_nlp, iter_num_, solver_status_)) { break; } if(NlpSolve_Pending!=solver_status_) break; //failure of the line search or user stopped. @@ -1989,7 +2292,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() if(!mu_updated) { break; } - nlp->log->printf(hovScalars, "Iter[%d] barrier params reduced: mu=%g tau=%g\n", iter_num, _mu, _tau); + nlp->log->printf(hovScalars, "Iter[%d] barrier params reduced: mu=%g tau=%g\n", iter_num_, _mu, _tau); //update only logbar problem and residual (the NLP didn't change) logbar->updateWithNlpInfo(*it_curr, _mu, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); @@ -2023,7 +2326,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() break; } } - nlp->log->printf(hovScalars, "Iter[%d] logbarObj=%23.17e (mu=%12.5e)\n", iter_num, logbar->f_logbar,_mu); + nlp->log->printf(hovScalars, "Iter[%d] logbarObj=%23.17e (mu=%12.5e)\n", iter_num_, logbar->f_logbar,_mu); /**************************************************** * Search direction calculation ***************************************************/ @@ -2055,7 +2358,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() bool switched; kkt = switch_to_fast_KKT(kkt, _mu, - iter_num, + iter_num_, linsol_safe_mode_on, linsol_safe_mode_max_iters, linsol_safe_mode_last_iter_switched_on, @@ -2079,7 +2382,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() bool switched; kkt = switch_to_safer_KKT(kkt, _mu, - iter_num, + iter_num_, linsol_safe_mode_on, linsol_safe_mode_max_iters, linsol_safe_mode_last_iter_switched_on, @@ -2119,7 +2422,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() nlp->log->printf(hovWarning, "Requesting additional accuracy and stability from the KKT linear system " "at iteration %d (safe mode ON) [1]\n", - iter_num); + iter_num_); continue; } } // end of if(!kkt->update(it_curr, _grad_f, _Jac_c, _Jac_d, _Hess_Lagr)) @@ -2128,7 +2431,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() if(fact_acceptor_ic) { bool linsol_safe_mode_on_before = linsol_safe_mode_on; //compute_search_direction call below updates linsol safe mode flag and linsol_safe_mode_lastiter - if(!compute_search_direction(kkt, linsol_safe_mode_on, linsol_forcequick, iter_num)) { + if(!compute_search_direction(kkt, linsol_safe_mode_on, linsol_forcequick, iter_num_)) { if(linsol_safe_mode_on_before || linsol_forcequick) { //it fails under safe mode, this is fatal @@ -2143,7 +2446,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() assert(fact_acceptor_dwd); bool linsol_safe_mode_on_before = linsol_safe_mode_on; //compute_search_direction call below updates linsol safe mode flag and linsol_safe_mode_lastiter - if(!compute_search_direction_inertia_free(kkt, linsol_safe_mode_on, linsol_forcequick, iter_num)) { + if(!compute_search_direction_inertia_free(kkt, linsol_safe_mode_on, linsol_forcequick, iter_num_)) { if(linsol_safe_mode_on_before || linsol_forcequick) { //it failed under safe mode delete kkt; @@ -2160,7 +2463,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() nlp->log->printf(hovSummary, "%s", nlp->runStats.kkt.get_summary_last_iter().c_str()); } - nlp->log->printf(hovIteration, "Iter[%d] full search direction -------------\n", iter_num); + nlp->log->printf(hovIteration, "Iter[%d] full search direction -------------\n", iter_num_); nlp->log->write("", *dir, hovIteration); /*************************************************************** * backtracking line search @@ -2300,7 +2603,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() // post line-search: filter is augmented whenever the switching condition or Armijo rule do not // hold for the trial point that was just accepted - if(nlp->options->GetString("force_resto")=="yes" && !within_FR_ && iter_num == 1) { + if(nlp->options->GetString("force_resto")=="yes" && !within_FR_ && iter_num_ == 1) { use_fr = apply_feasibility_restoration(kkt); if(use_fr) { // continue iterations if FR is accepted @@ -2382,7 +2685,8 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() nlp->log->printf(hovWarning, "Requesting additional accuracy and stability from the KKT linear system " - "at iteration %d (safe mode ON) [2]\n", iter_num); + "at iteration %d (safe mode ON) [2]\n", + iter_num_); // repeat linear solve (compute_search_direction) in safe mode (meaning additional accuracy // and stability is requested) @@ -2402,10 +2706,12 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() nlp->log->printf(hovScalars, "Iter[%d] -> accepted step primal=[%17.11e] dual=[%17.11e]\n", - iter_num,_alpha_primal, + iter_num_, + _alpha_primal, _alpha_dual); - iter_num++; - nlp->runStats.nIter=iter_num; + iter_num_++; + iter_num_total_++; + nlp->runStats.nIter = iter_num_; // fr problem has already updated dual, slacks and NLP functions if(!use_fr) { @@ -2439,7 +2745,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() // hiopIterate* pit=it_curr; it_curr=it_trial; it_trial=pit; - nlp->log->printf(hovIteration, "Iter[%d] -> full iterate:", iter_num); + nlp->log->printf(hovIteration, "Iter[%d] -> full iterate:", iter_num_); nlp->log->write("", *it_curr, hovIteration); nlp->runStats.tmSolverInternal.stop(); //----- @@ -2450,7 +2756,7 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() //update residual resid->update(*it_curr,_f_nlp, *_c, *_d,*_grad_f,*_Jac_c,*_Jac_d, *logbar); - nlp->log->printf(hovIteration, "Iter[%d] full residual:-------------\n", iter_num); + nlp->log->printf(hovIteration, "Iter[%d] full residual:-------------\n", iter_num_); nlp->log->write("", *resid, hovIteration); } @@ -2476,12 +2782,12 @@ hiopSolveStatus hiopAlgFilterIPMNewton::run() void hiopAlgFilterIPMNewton::outputIteration(int lsStatus, int lsNum, int use_soc, int use_fr) { - if(iter_num/10*10==iter_num) + if(iter_num_/10*10==iter_num_) nlp->log->printf(hovSummary, "iter objective inf_pr inf_du lg(mu) alpha_du alpha_pr linesrch\n"); if(lsStatus==-1) nlp->log->printf(hovSummary, "%4d %14.7e %7.3e %7.3e %6.2f %7.3e %7.3e -(-)\n", - iter_num, _f_nlp/nlp->get_obj_scale(), _err_nlp_feas, _err_nlp_optim, + iter_num_total_, _f_nlp/nlp->get_obj_scale(), _err_nlp_feas, _err_nlp_optim, log10(_mu), _alpha_dual, _alpha_primal); else { char stepType[2]; @@ -2500,7 +2806,7 @@ void hiopAlgFilterIPMNewton::outputIteration(int lsStatus, int lsNum, int use_so } nlp->log->printf(hovSummary, "%4d %14.7e %7.3e %7.3e %6.2f %7.3e %7.3e %d(%s)\n", - iter_num, _f_nlp/nlp->get_obj_scale(), _err_nlp_feas, _err_nlp_optim, + iter_num_total_, _f_nlp/nlp->get_obj_scale(), _err_nlp_feas, _err_nlp_optim, log10(_mu), _alpha_dual, _alpha_primal, lsNum, stepType); } } @@ -2794,7 +3100,7 @@ bool hiopAlgFilterIPMBase::apply_feasibility_restoration(hiopKKTLinSys* kkt) // FR problem inside a FR problem, see equation (33) // use wildcard function to update primal variables x it_trial->copyFrom(*it_curr); - if(!nlp->user_force_update(iter_num, + if(!nlp->user_force_update(iter_num_, _f_nlp, *it_trial->get_x(), *it_trial->get_zl(), diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index 11902144b..bf90fa97c 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -67,6 +67,14 @@ #include "hiopPDPerturbation.hpp" #include "hiopFactAcceptor.hpp" +#ifdef HIOP_USE_AXOM +namespace axom { +namespace sidre { +class Group; // forward declaration +} +} +#endif + #include "hiopTimer.hpp" namespace hiop @@ -117,6 +125,8 @@ class hiopAlgFilterIPMBase { { return filter.contains(theta, logbar_obj); } + + /// Setter for the primal steplength. inline void set_alpha_primal(const double alpha_primal) { _alpha_primal = alpha_primal; } protected: @@ -255,7 +265,11 @@ class hiopAlgFilterIPMBase { hiopResidual* resid, *resid_trial; - int iter_num; + /// Iteration number maintained internally by the algorithm and reset at each solve/run + int iter_num_; + /// Total iteration number over multiple solves/restarts using checkpoints. + int iter_num_total_; + double _err_nlp_optim, _err_nlp_feas, _err_nlp_complem;//not scaled by sd, sc, and sc double _err_nlp_optim0,_err_nlp_feas0,_err_nlp_complem0;//initial errors, not scaled by sd, sc, and sc double _err_log_optim, _err_log_feas, _err_log_complem;//not scaled by sd, sc, and sc @@ -339,10 +353,86 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase virtual ~hiopAlgFilterIPMQuasiNewton(); virtual hiopSolveStatus run(); + + // note that checkpointing is only available with a axom-enabled build +#ifdef HIOP_USE_AXOM + /** + * @brief Save state of HiOp algorithm to a sidre::Group as a checkpoint. + * + * @param group a reference to the group where state will be saved to + * + * @exception std::runtime indicates the group contains a view whose size does not match + * the size of the corresponding HiOp algorithm state variable of parameter. + * + * @details + * Each state variable of each parameter of HiOp algorithm will be saved in a named + * view within the group. A new view will be created within the group if it does not + * already exist. If it exists, the view must have same number of elements as the + * as the size of the corresponding state variable. This means that this method will + * throw an exception if an existing group is reused to save a problem that changed + * sizes since the group was created. + */ + virtual void save_state_to_sidre_group(::axom::sidre::Group& group); + + /** + * @brief Load state of HiOp algorithm from a sidre::Group checkpoint. + * + * @param group a pointer to group containing the a prevously saved HiOp algorithm state. + * + * @exception std::runtime indicates the group does not contain a view expected by this + * method or the view's number of elements mismatches the size of the corresponding HiOp + * state. The latter can occur if the file was saved with a different number of MPI ranks. + * + * @details + * Copies views from the sidre::Group passed as argument to HiOp algorithm's state variables + * and parameters. The group should be created by first calling save_state_to_sidre_group + * for a problem/NLP of the same sizes as the problem for which this method is called. + * The method expects views within the group with certain names. If one such view is not + * found or has a number of elements different than the size of the corresponding HiOp state, + * then a std::runtime_error exception is thrown. The latter can occur when the loading + * occurs for a instance of HiOp that is not ran on the same number of MPI ranks used to + * save the file. + */ + virtual void load_state_from_sidre_group(const ::axom::sidre::Group& group); + + /** + * @brief Save the state of the algorithm to the file for checkpointing. + * + * @param path the name of the file + * @return true if successful, false otherwise + * + * @details + * Internally, HiOp uses axom::sidre::DataStore and sidre's scalable IO. A detailed + * error description is sent to the log if an error or exception is caught. + */ + bool save_state_to_file(const ::std::string& path) noexcept; + + /** + * @brief Load the state of the algorithm from checkpoint file. + * + * @param path the name of the file to load from + * @return true if successful, false otherwise + * + * @details + * The file should contains a axom::sidre::DataStore that was previously saved using + * save_state_to_file(). A detailed error description is sent to the log if an error + * or exception is caught. + */ + bool load_state_from_file(const ::std::string& path) noexcept; +#endif // HIOP_USE_AXOM private: virtual void outputIteration(int lsStatus, int lsNum, int use_soc = 0, int use_fr = 0); + +#ifdef HIOP_USE_AXOM + ///@brief The options-based logic for saving checkpoint and the call to save_state(). + void checkpointing_stuff(); +#endif // HIOP_USE_AXOM + private: hiopNlpDenseConstraints* nlpdc; + ///@brief Indicates whether load checkpoint API was called previous to run method. + bool load_state_api_called_; + private: hiopAlgFilterIPMQuasiNewton() : hiopAlgFilterIPMBase(NULL) {}; hiopAlgFilterIPMQuasiNewton(const hiopAlgFilterIPMQuasiNewton& ) : hiopAlgFilterIPMBase(NULL){}; diff --git a/src/Optimization/hiopHessianLowRank.cpp b/src/Optimization/hiopHessianLowRank.cpp index a18af6d9b..cba17d07a 100644 --- a/src/Optimization/hiopHessianLowRank.cpp +++ b/src/Optimization/hiopHessianLowRank.cpp @@ -75,19 +75,22 @@ using namespace std; namespace hiop { -hiopHessianLowRank::hiopHessianLowRank(hiopNlpDenseConstraints* nlp_, int max_mem_len) - : l_max(max_mem_len), l_curr(-1), sigma(1.), sigma0(1.), nlp(nlp_), matrixChanged(false) + hiopHessianLowRank::hiopHessianLowRank(hiopNlpDenseConstraints* nlp_in, int max_mem_len) + : l_max(max_mem_len), l_curr(-1), sigma(1.), sigma0(1.), nlp(nlp_in), matrixChanged(false) { DhInv = nlp->alloc_primal_vec(); - St = nlp->alloc_multivector_primal(0,l_max); - Yt = St->alloc_clone(); //faster than nlp->alloc_multivector_primal(...); + St_ = nlp->alloc_multivector_primal(0, l_max); + Yt_ = St_->alloc_clone(); //faster than nlp->alloc_multivector_primal(...); //these are local - L = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 0, 0); - D = LinearAlgebraFactory::create_vector("DEFAULT", 0); - V = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 0, 0); + L_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 0, 0); + D_ = LinearAlgebraFactory::create_vector("DEFAULT", 0); + V_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 0, 0); //the previous iteration's objects are set to NULL - _it_prev=NULL; _grad_f_prev=NULL; _Jac_c_prev=NULL; _Jac_d_prev=NULL; + it_prev_ = new hiopIterate(nlp); + grad_f_prev_ = nlp->alloc_primal_vec(); + Jac_c_prev_ = nlp->alloc_Jac_c(); + Jac_d_prev_ = nlp->alloc_Jac_d(); //internal buffers for memory pool (none of them should be in n) #ifdef HIOP_USE_MPI @@ -138,7 +141,7 @@ hiopHessianLowRank::hiopHessianLowRank(hiopNlpDenseConstraints* nlp_, int max_me _Dx = DhInv->alloc_clone(); #ifdef HIOP_DEEPCHECKS - _Vmat = V->alloc_clone(); + _Vmat = V_->alloc_clone(); #endif yk = nlp->alloc_primal_vec(); @@ -151,11 +154,11 @@ hiopHessianLowRank::~hiopHessianLowRank() if(DhInv) delete DhInv; delete _Dx; - if(St) delete St; - if(Yt) delete Yt; - if(L) delete L; - if(D) delete D; - if(V) delete V; + delete St_; + delete Yt_; + delete L_; + delete D_; + delete V_; if(yk) delete yk; if(sk) delete sk; #ifdef HIOP_DEEPCHECKS @@ -163,10 +166,10 @@ hiopHessianLowRank::~hiopHessianLowRank() #endif - if(_it_prev) delete _it_prev; - if(_grad_f_prev)delete _grad_f_prev; - if(_Jac_c_prev) delete _Jac_c_prev; - if(_Jac_d_prev) delete _Jac_d_prev; + delete it_prev_; + delete grad_f_prev_; + delete Jac_c_prev_; + delete Jac_d_prev_; if(_buff_kxk) delete[] _buff_kxk; if(_buff_2lxk) delete[] _buff_2lxk; @@ -196,6 +199,24 @@ hiopHessianLowRank::~hiopHessianLowRank() } } +void hiopHessianLowRank::alloc_for_limited_mem(const size_type& mem_length) +{ + //note: St_ and Yt_ always have l_curr rows + if(l_curr == mem_length) { + assert(D_->get_size() == l_curr); + return; + } + delete D_; + delete L_; + delete Yt_; + delete St_; + St_ = nlp->alloc_multivector_primal(mem_length, l_max); + Yt_ = St_->alloc_clone(); + + //these are local + L_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", mem_length, mem_length); + D_ = LinearAlgebraFactory::create_vector("DEFAULT", mem_length); +} bool hiopHessianLowRank::updateLogBarrierDiagonal(const hiopVector& Dx) { @@ -222,8 +243,8 @@ void hiopHessianLowRank::print(FILE* f, hiopOutVerbosity v, const char* msg) con #endif nlp->log->printf(v, "sigma=%22.16f;\n", sigma); nlp->log->write("DhInv", *DhInv, v); - nlp->log->write("S_trans", *St, v); - nlp->log->write("Y_trans", *Yt, v); + nlp->log->write("S_trans", *St_, v); + nlp->log->write("Y_trans", *Yt_, v); fprintf(f, " [[Internal representation]]\n"); #ifdef HIOP_DEEPCHECKS @@ -231,8 +252,8 @@ void hiopHessianLowRank::print(FILE* f, hiopOutVerbosity v, const char* msg) con #else fprintf(f, "V matrix is available at this point (only its LAPACK factorization). Print it in updateInternalBFGSRepresentation() instead, before factorizeV()\n"); #endif - nlp->log->write("L", *L, v); - nlp->log->write("D", *D, v); + nlp->log->write("L", *L_, v); + nlp->log->write("D", *D_, v); } #endif @@ -256,7 +277,9 @@ bool hiopHessianLowRank::update(const hiopIterate& it_curr, const hiopVector& gr if(l_curr>=0) { size_type n=grad_f_curr_.get_size(); //compute s_new = x_curr-x_prev - hiopVector& s_new = new_n_vec1(n); s_new.copyFrom(*it_curr.x); s_new.axpy(-1.,*_it_prev->x); + hiopVector& s_new = new_n_vec1(n); + s_new.copyFrom(*it_curr.x); + s_new.axpy(-1.,*it_prev_->x); double s_infnorm=s_new.infnorm(); if(s_infnorm>=100*std::numeric_limits::epsilon()) { //norm of s not too small @@ -264,11 +287,13 @@ bool hiopHessianLowRank::update(const hiopIterate& it_curr, const hiopVector& gr // = graf_f_curr-grad_f_prev + (Jac_c_curr-Jac_c_prev)yc_curr+ (Jac_d_curr-Jac_c_prev)yd_curr - zl_curr*s_new + zu_curr*s_new hiopVector& y_new = new_n_vec2(n); y_new.copyFrom(grad_f_curr_); - y_new.axpy(-1., *_grad_f_prev); + y_new.axpy(-1., *grad_f_prev_); Jac_c_curr.transTimesVec (1.0, y_new, 1.0, *it_curr.yc); - _Jac_c_prev->transTimesVec(1.0, y_new,-1.0, *it_curr.yc); //!opt if nlp->Jac_c_isLinear no need for the multiplications - Jac_d_curr.transTimesVec (1.0, y_new, 1.0, *it_curr.yd); //!opt same here - _Jac_d_prev->transTimesVec(1.0, y_new,-1.0, *it_curr.yd); + //!opt if nlp->Jac_c_isLinear no need for the multiplications + Jac_c_prev_->transTimesVec(1.0, y_new,-1.0, *it_curr.yc); + //!opt same here + Jac_d_curr.transTimesVec (1.0, y_new, 1.0, *it_curr.yd); + Jac_d_prev_->transTimesVec(1.0, y_new,-1.0, *it_curr.yd); double sTy = s_new.dotProductWith(y_new), s_nrm2=s_new.twonorm(), y_nrm2=y_new.twonorm(); @@ -282,21 +307,21 @@ bool hiopHessianLowRank::update(const hiopIterate& it_curr, const hiopVector& gr if(l_max>0) { //compute the new row in L, update S and Y (either augment them or shift cols and add s_new and y_new) hiopVector& YTs = new_l_vec1(l_curr); - Yt->timesVec(0.0, YTs, 1.0, s_new); + Yt_->timesVec(0.0, YTs, 1.0, s_new); //update representation if(l_currappendRow(s_new); - Yt->appendRow(y_new); + St_->appendRow(s_new); + Yt_->appendRow(y_new); growL(l_curr, l_max, YTs); growD(l_curr, l_max, sTy); l_curr++; } else { //shift - St->shiftRows(-1); - Yt->shiftRows(-1); - St->replaceRow(l_max-1, s_new); - Yt->replaceRow(l_max-1, y_new); + St_->shiftRows(-1); + Yt_->shiftRows(-1); + St_->replaceRow(l_max-1, s_new); + Yt_->replaceRow(l_max-1, y_new); updateL(YTs,sTy); updateD(sTy); l_curr=l_max; @@ -304,8 +329,8 @@ bool hiopHessianLowRank::update(const hiopIterate& it_curr, const hiopVector& gr } //end of l_max>0 #ifdef HIOP_DEEPCHECKS nlp->log->printf(hovMatrices, "\nhiopHessianLowRank: these are L and D from the BFGS compact representation\n"); - nlp->log->write("L", *L, hovMatrices); - nlp->log->write("D", *D, hovMatrices); + nlp->log->write("L", *L_, hovMatrices); + nlp->log->write("D", *D_, hovMatrices); nlp->log->printf(hovMatrices, "\n"); #endif //update B0 (i.e., sigma) @@ -340,16 +365,18 @@ bool hiopHessianLowRank::update(const hiopIterate& it_curr, const hiopVector& gr } //save this stuff for next update - _it_prev->copyFrom(it_curr); _grad_f_prev->copyFrom(grad_f_curr_); - _Jac_c_prev->copyFrom(Jac_c_curr); _Jac_d_prev->copyFrom(Jac_d_curr); + it_prev_->copyFrom(it_curr); + grad_f_prev_->copyFrom(grad_f_curr_); + Jac_c_prev_->copyFrom(Jac_c_curr); + Jac_d_prev_->copyFrom(Jac_d_curr); nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianLowRank: storing the iteration info as 'previous'\n", s_infnorm); } else { //this is the first optimization iterate, just save the iterate and exit - if(NULL==_it_prev) _it_prev = it_curr.new_copy(); - if(NULL==_grad_f_prev) _grad_f_prev = grad_f_curr_.new_copy(); - if(NULL==_Jac_c_prev) _Jac_c_prev = Jac_c_curr.new_copy(); - if(NULL==_Jac_d_prev) _Jac_d_prev = Jac_d_curr.new_copy(); + it_prev_->copyFrom(it_curr); + grad_f_prev_->copyFrom(grad_f_curr_); + Jac_c_prev_->copyFrom(Jac_c_curr); + Jac_d_prev_->copyFrom(Jac_d_curr); nlp->log->printf(hovLinAlgScalarsVerb, "HessianLowRank on first update, just saving iteration\n"); @@ -372,36 +399,46 @@ bool hiopHessianLowRank::update(const hiopIterate& it_curr, const hiopVector& gr */ void hiopHessianLowRank::updateInternalBFGSRepresentation() { - size_type n=St->n(), l=St->m(); + size_type n=St_->n(); + size_type l=St_->m(); //grow L,D, andV if needed - if(L->m()!=l) { delete L; L=LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, l);} - if(D->get_size()!=l) { delete D; D=LinearAlgebraFactory::create_vector("DEFAULT", l); } - if(V->m()!=2*l) {delete V; V=LinearAlgebraFactory::create_matrix_dense("DEFAULT", 2*l, 2*l); } + if(L_->m()!=l) { + delete L_; + L_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, l); + } + if(D_->get_size()!=l) { + delete D_; + D_ = LinearAlgebraFactory::create_vector("DEFAULT", l); + } + if(V_->m()!=2*l) { + delete V_; + V_ = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 2*l, 2*l); + } //-- block (2,2) hiopMatrixDense& DpYtDhInvY = new_lxl_mat1(l); - symmMatTimesDiagTimesMatTrans_local(0.0, DpYtDhInvY, 1.0,*Yt,*DhInv); + symmMatTimesDiagTimesMatTrans_local(0.0, DpYtDhInvY, 1.0,*Yt_,*DhInv); #ifdef HIOP_USE_MPI const size_t buffsize=l*l*sizeof(double); memcpy(_buff1_lxlx3, DpYtDhInvY.local_data(), buffsize); #else - DpYtDhInvY.addDiagonal(1., *D); - V->copyBlockFromMatrix(l,l,DpYtDhInvY); + DpYtDhInvY.addDiagonal(1., *D_); + V_->copyBlockFromMatrix(l,l,DpYtDhInvY); #endif //-- block (1,2) hiopMatrixDense& StB0DhInvYmL = DpYtDhInvY; //just a rename hiopVector& B0DhInv = new_n_vec1(n); B0DhInv.copyFrom(*DhInv); B0DhInv.scale(sigma); - matTimesDiagTimesMatTrans_local(StB0DhInvYmL, *St, B0DhInv, *Yt); + matTimesDiagTimesMatTrans_local(StB0DhInvYmL, *St_, B0DhInv, *Yt_); #ifdef HIOP_USE_MPI memcpy(_buff1_lxlx3+l*l, StB0DhInvYmL.local_data(), buffsize); #else //substract L - StB0DhInvYmL.addMatrix(-1.0, *L); + StB0DhInvYmL.addMatrix(-1.0, *L_); // (1,2) block in V - V->copyBlockFromMatrix(0,l,StB0DhInvYmL); + V_->copyBlockFromMatrix(0,l,StB0DhInvYmL); #endif //-- block (2,2) @@ -409,11 +446,11 @@ void hiopHessianLowRank::updateInternalBFGSRepresentation() theDiag.addConstant(-1.0); //at this point theDiag=DhInv*B0-I theDiag.scale(sigma); hiopMatrixDense& StDS = DpYtDhInvY; //a rename - symmMatTimesDiagTimesMatTrans_local(0.0, StDS, 1.0, *St, theDiag); + symmMatTimesDiagTimesMatTrans_local(0.0, StDS, 1.0, *St_, theDiag); #ifdef HIOP_USE_MPI memcpy(_buff1_lxlx3+2*l*l, DpYtDhInvY.local_data(), buffsize); #else - V->copyBlockFromMatrix(0,0,StDS); + V_->copyBlockFromMatrix(0,0,StDS); #endif @@ -423,21 +460,21 @@ void hiopHessianLowRank::updateInternalBFGSRepresentation() // - block (2,2) DpYtDhInvY.copyFrom(_buff2_lxlx3); - DpYtDhInvY.addDiagonal(1., *D); - V->copyBlockFromMatrix(l,l,DpYtDhInvY); + DpYtDhInvY.addDiagonal(1., *D_); + V_->copyBlockFromMatrix(l,l,DpYtDhInvY); // - block (1,2) StB0DhInvYmL.copyFrom(_buff2_lxlx3+l*l); - StB0DhInvYmL.addMatrix(-1.0, *L); - V->copyBlockFromMatrix(0,l,StB0DhInvYmL); + StB0DhInvYmL.addMatrix(-1.0, *L_); + V_->copyBlockFromMatrix(0,l,StB0DhInvYmL); // - block (1,1) StDS.copyFrom(_buff2_lxlx3+2*l*l); - V->copyBlockFromMatrix(0,0,StDS); + V_->copyBlockFromMatrix(0,0,StDS); #endif #ifdef HIOP_DEEPCHECKS delete _Vmat; - _Vmat = V->new_copy(); + _Vmat = V_->new_copy(); _Vmat->overwriteLowerTriangleWithUpper(); #endif @@ -459,7 +496,7 @@ void hiopHessianLowRank::solve(const hiopVector& rhsx, hiopVector& x) { if(matrixChanged) updateInternalBFGSRepresentation(); - size_type n=St->n(), l=St->m(); + size_type n=St_->n(), l=St_->m(); #ifdef HIOP_DEEPCHECKS assert(rhsx.get_size()==n); assert(x.get_size()==n); @@ -475,12 +512,12 @@ void hiopHessianLowRank::solve(const hiopVector& rhsx, hiopVector& x) //2. stx= S^T*B0*DhInv*res and ytx=Y^T*DhInv*res hiopVector&stx=new_l_vec1(l), &ytx=new_l_vec2(l); stx.setToZero(); ytx.setToZero(); - Yt->timesVec(0.0,ytx,1.0,x); + Yt_->timesVec(0.0,ytx,1.0,x); hiopVector& B0DhInvx = new_n_vec1(n); B0DhInvx.copyFrom(x); //it contains DhInv*res B0DhInvx.scale(sigma); //B0*(DhInv*res) - St->timesVec(0.0,stx,1.0,B0DhInvx); + St_->timesVec(0.0, stx, 1.0, B0DhInvx); //3. solve with V hiopVector& spart=stx; hiopVector& ypart=ytx; @@ -489,9 +526,9 @@ void hiopHessianLowRank::solve(const hiopVector& rhsx, hiopVector& x) //4. multiply with DhInv*[B0*S Y], namely // result = DhInv*(B0*S*spart + Y*ypart) hiopVector& result = new_n_vec1(n); - St->transTimesVec(0.0, result, 1.0, spart); + St_->transTimesVec(0.0, result, 1.0, spart); result.scale(sigma); - Yt->transTimesVec(1.0, result, 1.0, ypart); + Yt_->transTimesVec(1.0, result, 1.0, ypart); result.componentMult(*DhInv); //5. x = first term - second term = x_computed_in_1 - result @@ -510,12 +547,14 @@ void hiopHessianLowRank::solve(const hiopVector& rhsx, hiopVector& x) * X is kxn */ void hiopHessianLowRank:: -symMatTimesInverseTimesMatTrans(double beta, hiopMatrixDense& W, - double alpha, const hiopMatrixDense& X) +symMatTimesInverseTimesMatTrans(double beta, + hiopMatrixDense& W, + double alpha, + const hiopMatrixDense& X) { if(matrixChanged) updateInternalBFGSRepresentation(); - size_type n=St->n(), l=St->m(); + size_type n=St_->n(), l=St_->m(); size_type k=W.m(); assert(X.m()==k); assert(X.n()==n); @@ -535,11 +574,12 @@ symMatTimesInverseTimesMatTrans(double beta, hiopMatrixDense& W, symmMatTimesDiagTimesMatTrans_local(beta,W,alpha,X,*DhInv); #endif //2. compute S1=X*DhInv*B0*S and Y1=X*DhInv*Y - hiopMatrixDense &S1=new_S1(X,*St), &Y1=new_Y1(X,*Yt); //both are kxl + auto& S1=new_S1(X,*St_); + auto& Y1=new_Y1(X,*Yt_); //both are kxl hiopVector& B0DhInv = new_n_vec1(n); B0DhInv.copyFrom(*DhInv); B0DhInv.scale(sigma); - matTimesDiagTimesMatTrans_local(S1, X, B0DhInv, *St); - matTimesDiagTimesMatTrans_local(Y1, X, *DhInv, *Yt); + matTimesDiagTimesMatTrans_local(S1, X, B0DhInv, *St_); + matTimesDiagTimesMatTrans_local(Y1, X, *DhInv, *Yt_); //3. reduce W, S1, and Y1 (dimensions: kxk, kxl, kxl) hiopMatrixDense& S2Y2 = new_kx2l_mat1(k,l); //Initialy S2Y2 = [Y1 S1] @@ -592,11 +632,15 @@ symMatTimesInverseTimesMatTrans(double beta, hiopMatrixDense& W, void hiopHessianLowRank::factorizeV() { - int N=V->n(), lda=N, info; - if(N==0) return; + int N = V_->n(); + int lda = N; + int info; + if(N==0) { + return; + } #ifdef HIOP_DEEPCHECKS - nlp->log->write("factorizeV: V is ", *V, hovMatrices); + nlp->log->write("factorizeV: V is ", *V_, hovMatrices); #endif char uplo='L'; //V is upper in C++ so it's lower in fortran @@ -606,7 +650,7 @@ void hiopHessianLowRank::factorizeV() int lwork=-1;//inquire sizes double Vwork_tmp; - DSYTRF(&uplo, &N, V->local_data(), &lda, _V_ipiv_vec, &Vwork_tmp, &lwork, &info); + DSYTRF(&uplo, &N, V_->local_data(), &lda, _V_ipiv_vec, &Vwork_tmp, &lwork, &info); assert(info==0); lwork=(int)Vwork_tmp; @@ -617,7 +661,7 @@ void hiopHessianLowRank::factorizeV() _V_work_vec = LinearAlgebraFactory::create_vector("DEFAULT", lwork); } else assert(_V_work_vec); - DSYTRF(&uplo, &N, V->local_data(), &lda, _V_ipiv_vec, _V_work_vec->local_data(), &lwork, &info); + DSYTRF(&uplo, &N, V_->local_data(), &lda, _V_ipiv_vec, _V_work_vec->local_data(), &lwork, &info); if(info<0) nlp->log->printf(hovError, "hiopHessianLowRank::factorizeV error: %d argument to dsytrf has an illegal value\n", -info); @@ -625,17 +669,19 @@ void hiopHessianLowRank::factorizeV() nlp->log->printf(hovError, "hiopHessianLowRank::factorizeV error: %d entry in the factorization's diagonal is exactly zero. Division by zero will occur if it a solve is attempted.\n", info); assert(info==0); #ifdef HIOP_DEEPCHECKS - nlp->log->write("factorizeV: factors of V: ", *V, hovMatrices); + nlp->log->write("factorizeV: factors of V: ", *V_, hovMatrices); #endif } void hiopHessianLowRank::solveWithV(hiopVector& rhs_s, hiopVector& rhs_y) { - int N=V->n(); - if(N==0) return; + int N = V_->n(); + if(N==0) { + return; + } - int l=rhs_s.get_size(); + int l = rhs_s.get_size(); #ifdef HIOP_DEEPCHECKS nlp->log->write("hiopHessianLowRank::solveWithV: RHS IN 's' part: ", rhs_s, hovMatrices); @@ -654,7 +700,7 @@ void hiopHessianLowRank::solveWithV(hiopVector& rhs_s, hiopVector& rhs_y) rhs.copyFromStarting(0, rhs_s); rhs.copyFromStarting(l, rhs_y); - DSYTRS(&uplo, &N, &one, V->local_data(), &lda, _V_ipiv_vec, rhs.local_data(), &N, &info); + DSYTRS(&uplo, &N, &one, V_->local_data(), &lda, _V_ipiv_vec, rhs.local_data(), &N, &info); if(info<0) nlp->log->printf(hovError, "hiopHessianLowRank::solveWithV error: %d argument to dsytrf has an illegal value\n", -info); assert(info==0); @@ -682,8 +728,10 @@ void hiopHessianLowRank::solveWithV(hiopVector& rhs_s, hiopVector& rhs_y) void hiopHessianLowRank::solveWithV(hiopMatrixDense& rhs) { - int N=V->n(); - if(0==N) return; + int N = V_->n(); + if(0==N) { + return; + } #ifdef HIOP_DEEPCHECKS nlp->log->write("solveWithV: RHS IN: ", rhs, hovMatrices); @@ -697,7 +745,7 @@ void hiopHessianLowRank::solveWithV(hiopMatrixDense& rhs) #ifdef HIOP_DEEPCHECKS assert(N==rhs.n()); #endif - DSYTRS(&uplo, &N, &nrhs, V->local_data(), &lda, _V_ipiv_vec, rhs.local_data(), &ldb, &info); + DSYTRS(&uplo, &N, &nrhs, V_->local_data(), &lda, _V_ipiv_vec, rhs.local_data(), &ldb, &info); if(info<0) nlp->log->printf(hovError, "hiopHessianLowRank::solveWithV error: %d argument to dsytrf has an illegal value\n", -info); assert(info==0); @@ -730,9 +778,9 @@ void hiopHessianLowRank::solveWithV(hiopMatrixDense& rhs) void hiopHessianLowRank::growL(const int& lmem_curr, const int& lmem_max, const hiopVector& YTs) { - int l=L->m(); + int l = L_->m(); #ifdef HIOP_DEEPCHECKS - assert(l==L->n()); + assert(l==L_->n()); assert(lmem_curr==l); assert(lmem_max>=l); #endif @@ -741,7 +789,7 @@ void hiopHessianLowRank::growL(const int& lmem_curr, const int& lmem_max, const hiopMatrixDense* newL = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l+1, l+1); assert(newL); //copy from L to newL - newL->copyBlockFromMatrix(0,0, *L); + newL->copyBlockFromMatrix(0,0, *L_); double* newL_mat=newL->local_data(); //doing the rest here const double* YTs_vec=YTs.local_data_const(); @@ -753,23 +801,23 @@ void hiopHessianLowRank::growL(const int& lmem_curr, const int& lmem_max, const for(int i=0; iget_size(); + int l = D_->get_size(); assert(l==lmem_curr); assert(lmem_max>=l); hiopVector* Dnew=LinearAlgebraFactory::create_vector("DEFAULT", l+1); double* Dnew_vec=Dnew->local_data(); - memcpy(Dnew_vec, D->local_data_const(), l*sizeof(double)); + memcpy(Dnew_vec, D_->local_data_const(), l*sizeof(double)); Dnew_vec[l]=sTy; - delete D; - D=Dnew; + delete D_; + D_ = Dnew; } /* L_{ij} = s_{i-1}^T y_{j-1}, if i>j, otherwise zero. Here i,j = 0,1,...,l_curr-1 @@ -778,14 +826,14 @@ void hiopHessianLowRank::growD(const int& lmem_curr, const int& lmem_max, const void hiopHessianLowRank::updateL(const hiopVector& YTs, const double& sTy) { int l=YTs.get_size(); - assert(l==L->m()); - assert(l==L->n()); + assert(l==L_->m()); + assert(l==L_->n()); #ifdef HIOP_DEEPCHECKS assert(l_curr==l); assert(l_curr==l_max); #endif const int lm1=l-1; - double* L_mat=L->local_data(); + double* L_mat=L_->local_data(); const double* yts_vec=YTs.local_data_const(); for(int i=1; iget_size(); - double* D_vec = D->local_data(); + int l=D_->get_size(); + double* D_vec = D_->local_data(); for(int i=0; in(); - assert(l_curr==St->m()); + size_type n=St_->n(); + assert(l_curr==St_->m()); assert(y.get_size()==n); - assert(St->get_local_size_n() == Yt->get_local_size_n()); + assert(St_->get_local_size_n() == Yt_->get_local_size_n()); //we have B+=B-B*s*B*s'/(s'*B*s)+yy'/(y'*s) //B0 is sigma*I. There is an additional diagonal log-barrier term _Dx @@ -936,8 +985,8 @@ void hiopHessianLowRank::timesVecCmn(double beta, hiopVector& y, double alpha, c bool print=false; if(print) { nlp->log->printf(hovMatrices, "---hiopHessianLowRank::timesVec \n"); - nlp->log->write("S=", *St, hovMatrices); - nlp->log->write("Y=", *Yt, hovMatrices); + nlp->log->write("S=", *St_, hovMatrices); + nlp->log->write("Y=", *Yt_, hovMatrices); nlp->log->write("DhInv=", *DhInv, hovMatrices); nlp->log->printf(hovMatrices, "sigma=%22.16e; addLogTerm=%d;\n", sigma, addLogTerm); if(addLogTerm) @@ -948,13 +997,14 @@ void hiopHessianLowRank::timesVecCmn(double beta, hiopVector& y, double alpha, c } //allocate and compute a_k and b_k + //! make sure the pointers within these std::vectors are deallocated a.resize(l_curr, nullptr); b.resize(l_curr, nullptr); - int n_local = Yt->get_local_size_n(); + int n_local = Yt_->get_local_size_n(); for(int k=0; kcopyFrom(Yt->local_data() + k*n_local); - sk->copyFrom(St->local_data() + k*n_local); + yk->copyFrom(Yt_->local_data() + k*n_local); + sk->copyFrom(St_->local_data() + k*n_local); double skTyk=yk->dotProductWith(*sk); if(skTyk < std::numeric_limits::epsilon()) { @@ -1102,754 +1152,4 @@ matTimesDiagTimesMatTrans_local(hiopMatrixDense& W, const hiopMatrixDense& S, co } } } - -/************************************************************************** - * this code is going to be removed - *************************************************************************/ -hiopHessianInvLowRank_obsolette::hiopHessianInvLowRank_obsolette(hiopNlpDenseConstraints* nlp_, int max_mem_len) - : hiopHessianLowRank(nlp_,max_mem_len) -{ - const hiopNlpDenseConstraints* nlp = dynamic_cast(nlp_); - // assert(nlp==NULL && "only NLPs with a small number of constraints are supported by HessianLowRank"); - - H0 = nlp->alloc_primal_vec(); - St = nlp->alloc_multivector_primal(0,l_max); - Yt = St->alloc_clone(); //faster than nlp->alloc_multivector_primal(...); - //these are local - R = LinearAlgebraFactory::create_matrix_dense("DEFAULT", 0, 0); - D = LinearAlgebraFactory::create_vector("DEFAULT", 0); - - //the previous iteration's objects are set to NULL - _it_prev=NULL; _grad_f_prev=NULL; _Jac_c_prev=NULL; _Jac_d_prev=NULL; - - - //internal buffers for memory pool (none of them should be in n) -#ifdef HIOP_USE_MPI - _buff_kxk = new double[nlp->m() * nlp->m()]; - _buff_lxk = new double[nlp->m() * l_max]; - _buff_lxl = new double[l_max*l_max]; -#else - //not needed in non-MPI mode - _buff_kxk = NULL; - _buff_lxk = NULL; - _buff_lxl = NULL; -#endif - - //auxiliary objects/buffers - _S1=_Y1=_S3=NULL; - _DpYtH0Y=NULL; - _l_vec1 = _l_vec2 = NULL; - _n_vec1 = H0->alloc_clone(); - _n_vec2 = H0->alloc_clone(); - //H0->setToConstant(sigma); - - sigma=sigma0; - sigma_update_strategy = SIGMA_STRATEGY1; - sigma_safe_min=1e-8; - sigma_safe_max=1e+8; - nlp->log->printf(hovScalars, "Hessian Low Rank: initial sigma is %g\n", sigma); -} -hiopHessianInvLowRank_obsolette::~hiopHessianInvLowRank_obsolette() -{ - if(H0) delete H0; - if(St) delete St; - if(Yt) delete Yt; - if(R) delete R; - if(D) delete D; - - if(_it_prev) delete _it_prev; - if(_grad_f_prev)delete _grad_f_prev; - if(_Jac_c_prev) delete _Jac_c_prev; - if(_Jac_d_prev) delete _Jac_d_prev; - - if(_buff_kxk) delete[] _buff_kxk; - if(_buff_lxk) delete[] _buff_lxk; - if(_buff_lxl) delete[] _buff_lxl; - if(_S1) delete _S1; - if(_Y1) delete _Y1; - if(_DpYtH0Y) delete _DpYtH0Y; - if(_S3) delete _S3; - if(_l_vec1) delete _l_vec1; - if(_l_vec2) delete _l_vec2; - if(_n_vec1) delete _n_vec1; - if(_n_vec2) delete _n_vec2; -} - -#include - -bool hiopHessianInvLowRank_obsolette:: -update(const hiopIterate& it_curr, const hiopVector& grad_f_curr_, - const hiopMatrix& Jac_c_curr_, const hiopMatrix& Jac_d_curr_) -{ - const hiopMatrixDense& Jac_c_curr = dynamic_cast(Jac_c_curr_); - const hiopMatrixDense& Jac_d_curr = dynamic_cast(Jac_d_curr_); - -#ifdef HIOP_DEEPCHECKS - assert(it_curr.zl->matchesPattern(nlp->get_ixl())); - assert(it_curr.zu->matchesPattern(nlp->get_ixu())); - assert(it_curr.sxl->matchesPattern(nlp->get_ixl())); - assert(it_curr.sxu->matchesPattern(nlp->get_ixu())); -#endif - - if(l_curr>0) { - size_type n=grad_f_curr_.get_size(); - //compute s_new = x_curr-x_prev - hiopVector& s_new = new_n_vec1(n); s_new.copyFrom(*it_curr.x); s_new.axpy(-1.,*_it_prev->x); - double s_infnorm=s_new.infnorm(); - if(s_infnorm>=100*std::numeric_limits::epsilon()) { //norm of s not too small - - //compute y_new = \grad J(x_curr,\lambda_curr) - \grad J(x_prev, \lambda_curr) (yes, J(x_prev, \lambda_curr)) - // = graf_f_curr-grad_f_prev + (Jac_c_curr-Jac_c_prev)yc_curr+ (Jac_d_curr-Jac_c_prev)yd_curr - zl_curr*s_new + zu_curr*s_new - hiopVector& y_new = new_n_vec2(n); - y_new.copyFrom(grad_f_curr_); - y_new.axpy(-1., *_grad_f_prev); - Jac_c_curr.transTimesVec (1.0, y_new, 1.0, *it_curr.yc); - _Jac_c_prev->transTimesVec(1.0, y_new,-1.0, *it_curr.yc); //!opt if nlp->Jac_c_isLinear no need for the multiplications - Jac_d_curr.transTimesVec (1.0, y_new, 1.0, *it_curr.yd); //!opt same here - _Jac_d_prev->transTimesVec(1.0, y_new,-1.0, *it_curr.yd); - //y_new.axzpy(-1.0, s_new, *it_curr.zl); - //y_new.axzpy( 1.0, s_new, *it_curr.zu); - - double sTy = s_new.dotProductWith(y_new), s_nrm2=s_new.twonorm(), y_nrm2=y_new.twonorm(); - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: s^T*y=%20.14e ||s||=%20.14e ||y||=%20.14e\n", sTy, s_nrm2, y_nrm2); - nlp->log->write("hiopHessianInvLowRank_obsolette s_new",s_new, hovIteration); - nlp->log->write("hiopHessianInvLowRank_obsolette y_new",y_new, hovIteration); - - if(sTy>s_nrm2*y_nrm2*std::numeric_limits::epsilon()) { //sTy far away from zero - //compute the new column in R, update S and Y (either augment them or shift cols and add s_new and y_new) - hiopVector& STy = new_l_vec1(l_curr-1); - St->timesVec(0.0, STy, 1.0, y_new); - //update representation - if(l_currappendRow(s_new); - Yt->appendRow(y_new); - growR(l_curr, l_max, STy, sTy); - growD(l_curr, l_max, sTy); - l_curr++; - } else { - //shift - St->shiftRows(-1); - Yt->shiftRows(-1); - St->replaceRow(l_max-1, s_new); - Yt->replaceRow(l_max-1, y_new); - updateR(STy,sTy); - updateD(sTy); - l_curr=l_max; - } - - //update B0 (i.e., sigma) - switch (sigma_update_strategy ) { - case SIGMA_STRATEGY1: - sigma=sTy/(s_nrm2*s_nrm2); - break; - case SIGMA_STRATEGY2: - sigma=y_nrm2*y_nrm2/sTy; - break; - case SIGMA_STRATEGY3: - sigma=sqrt(s_nrm2*s_nrm2 / y_nrm2 / y_nrm2); - break; - case SIGMA_STRATEGY4: - sigma=0.5*(sTy/(s_nrm2*s_nrm2)+y_nrm2*y_nrm2/sTy); - break; - case SIGMA_CONSTANT: - sigma=sigma0; - break; - default: - assert(false && "Option value for sigma_update_strategy was not recognized."); - break; - } // else of the switch - //safe guard it - sigma=fmax(fmin(sigma_safe_max, sigma), sigma_safe_min); - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: sigma was updated to %16.10e\n", sigma); - } else { //sTy is too small or negative -> skip - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: s^T*y=%12.6e not positive enough... skipping the Hessian update\n", sTy); - } - } else {// norm of s_new is too small -> skip - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: ||s_new||=%12.6e too small... skipping the Hessian update\n", s_infnorm); - } - - //save this stuff for next update - _it_prev->copyFrom(it_curr); _grad_f_prev->copyFrom(grad_f_curr_); - _Jac_c_prev->copyFrom(Jac_c_curr); _Jac_d_prev->copyFrom(Jac_d_curr); - nlp->log->printf(hovLinAlgScalarsVerb, "hiopHessianInvLowRank_obsolette: storing the iteration info as 'previous'\n", s_infnorm); - - } else { - //this is the first optimization iterate, just save the iterate and exit - if(NULL==_it_prev) _it_prev = it_curr.new_copy(); - if(NULL==_grad_f_prev) _grad_f_prev = grad_f_curr_.new_copy(); - if(NULL==_Jac_c_prev) _Jac_c_prev = Jac_c_curr.new_copy(); - if(NULL==_Jac_d_prev) _Jac_d_prev = Jac_d_curr.new_copy(); - - nlp->log->printf(hovLinAlgScalarsVerb, "HessianInvLowRank on first update, just saving iteration\n"); - - l_curr++; - } - return true; -} - -bool hiopHessianInvLowRank_obsolette::updateLogBarrierDiagonal(const hiopVector& Dx) -{ - H0->setToConstant(sigma); - H0->axpy(1.0,Dx); -#ifdef HIOP_DEEPCHECKS - assert(H0->allPositive()); -#endif - H0->invert(); - nlp->log->write("H0 InvHes", *H0, hovMatrices); - return true; -} - - -/* Y = beta*Y + alpha*this*X - * where 'this' is - * M^{-1} = H0 + [S HO*Y] [ R^{-T}*(D+Y^T*H0*Y)*R^{-1} -R^{-T} ] [ S^T ] - * [ -R^{-1} 0 ] [ Y^T*H0] - * - * M is is nxn, S,Y are nxl, R is upper triangular lxl, and X is nx1 - * Remember we store Yt=Y^T and St=S^T - */ -void hiopHessianInvLowRank_obsolette::apply(double beta, hiopVector& y, double alpha, const hiopVector& x) -{ - size_type n=St->n(), l=St->m(); -#ifdef HIOP_DEEPCHECKS - assert(y.get_size()==n); - assert(x.get_size()==n); - assert(H0->get_size()==n); -#endif - //0. y = beta*y + alpha*H0*y - y.scale(beta); - y.axzpy(alpha,*H0,x); - - //1. stx = S^T*x and ytx = Y^T*H0*x - hiopVector &stx=new_l_vec1(l), &ytx=new_l_vec2(l), &H0x=new_n_vec1(n); - St->timesVec(0.0,stx,1.0,x); - H0x.copyFrom(x); H0x.componentMult(*H0); - Yt->timesVec(0.0,ytx,1.0,H0x); - //2.ytx = R^{-T}* [(D+Y^T*H0*Y)*R^{-1} stx - ytx ] - // stx = -R^{-1}*stx - //2.1. stx = R^{-1}*stx - triangularSolve(*R, stx); - //2.2 ytx = -ytx + (D+Y^T*H0*Y)*R^{-1} stx - hiopMatrixDense& DpYtH0Y = *_DpYtH0Y; // ! this is computed in symmetricTimesMat - DpYtH0Y.timesVec(-1.0,ytx, 1.0, stx); - //2.3 ytx = R^{-T}* [(D+Y^T*H0*Y)*R^{-1} stx - ytx ] and stx = -R^{-1}*stx - triangularSolveTrans(*R,ytx); - - //3. add alpha(S*ytx + H0*Y*stx) to y - St->transTimesVec(1.0,y,alpha,ytx); - hiopVector& H0Ystx=new_n_vec1(n); - //H0Ystx=H0*Y*stx - Yt->transTimesVec(0.0, H0Ystx, -1.0, stx); //-1.0 since we haven't negated stx in 2.3 - H0Ystx.componentMult(*H0); - y.axpy(alpha,H0Ystx); -} -void hiopHessianInvLowRank_obsolette::apply(double beta, hiopMatrix& Y, double alpha, const hiopMatrix& X) -{ - assert(false && "not yet implemented"); -} - - -/* W = beta*W + alpha*X*this*X^T - * where 'this' is - * M^{-1} = H0 + [S HO*Y] [ R^{-T}*(D+Y^T*H0*Y)*R^{-1} -R^{-T} ] [ S^T ] - * [ -R^{-1} 0 ] [ Y^T*H0] - * - * W is kxk, H0 is nxn, S,Y are nxl, R is upper triangular lxl, and X is kxn - * Remember we store Yt=Y^T and St=S^T - */ -void hiopHessianInvLowRank_obsolette:: -symmetricTimesMat(double beta, hiopMatrixDense& W, - double alpha, const hiopMatrixDense& X) -{ - size_type n=St->n(), l=St->m(), k=W.m(); - assert(n==H0->get_size()); - assert(k==W.n()); - assert(l==Yt->m()); - assert(n==Yt->n()); assert(n==St->n()); - assert(k==X.m()); assert(n==X.n()); - - //1.--compute W = beta*W + alpha*X*HO*X^T by calling symmMatTransTimesDiagTimesMat_local -#ifdef HIOP_USE_MPI - int myrank, ierr; - ierr=MPI_Comm_rank(nlp->get_comm(),&myrank); assert(MPI_SUCCESS==ierr); - if(0==myrank) - symmMatTimesDiagTimesMatTrans_local(beta,W,alpha,X,*H0); - else - symmMatTimesDiagTimesMatTrans_local(0.0,W,alpha,X,*H0); - //W will be MPI_All_reduced later -#else - symmMatTimesDiagTimesMatTrans_local(beta,W,alpha,X,*H0); -#endif - - //2.--compute S1=S^T*X^T=St*X^T and Y1=Y^T*H0*X^T=Yt*H0*X^T - hiopMatrixDense& S1 = new_S1(*St,X); - St->timesMatTrans_local(0.0,S1,1.0,X); - hiopMatrixDense& Y1 = new_Y1(*Yt,X); - matTimesDiagTimesMatTrans_local(Y1,*Yt,*H0,X); - - //3.--compute Y^T*H0*Y from D+Y^T*H0*Y - hiopMatrixDense& DpYtH0Y = new_DpYtH0Y(*Yt); - symmMatTimesDiagTimesMatTrans_local(0.0,DpYtH0Y, 1.0,*Yt,*H0); -#ifdef HIOP_USE_MPI - //!opt - use one buffer and one reduce call - ierr=MPI_Allreduce(S1.local_data(), _buff_lxk,l*k, MPI_DOUBLE,MPI_SUM,nlp->get_comm()); assert(ierr==MPI_SUCCESS); - S1.copyFrom(_buff_lxk); - ierr=MPI_Allreduce(Y1.local_data(), _buff_lxk,l*k, MPI_DOUBLE,MPI_SUM,nlp->get_comm()); assert(ierr==MPI_SUCCESS); - Y1.copyFrom(_buff_lxk); - ierr=MPI_Allreduce(W.local_data(), _buff_kxk,k*k, MPI_DOUBLE,MPI_SUM,nlp->get_comm()); assert(ierr==MPI_SUCCESS); - W.copyFrom(_buff_kxk); - - ierr=MPI_Allreduce(DpYtH0Y.local_data(), _buff_lxl,l*l, MPI_DOUBLE,MPI_SUM,nlp->get_comm()); assert(ierr==MPI_SUCCESS); - DpYtH0Y.copyFrom(_buff_lxl); -#endif - //add D to finish calculating D+Y^T*H0*Y - DpYtH0Y.addDiagonal(1., *D); - - //now we have W = beta*W+alpha*X*HO*X^T and the remaining term in the form - // [S1^T Y1^T] [ R^{-T}*DpYtH0Y*R^{-1} -R^{-T} ] [ S1 ] - // [ -R^{-1} 0 ] [ Y1 ] - // So that W is updated with - // W += (S1^T*R^{-T})*DpYtH0Y*(R^{-1}*S1) - (S1^T*R^{-T})*Y1 - Y1^T*(R^{-1}*S1) - // or W += S2^T *DpYtH0Y* S2 - S2^T *Y1 - Y1^T*S2 - - //4.-- compute S1 = R\S1 - triangularSolve(*R,S1); - - //5.-- W += S2^T *DpYtH0Y* S2 - //S3=DpYtH0Y*S2 - hiopMatrix& S3=new_S3(DpYtH0Y, S1); - DpYtH0Y.timesMat(0.0,S3,1.0,S1); - // W += S2^T * S3 - S1.transTimesMat(1.0,W,1.0,S3); - - //6.-- W += - S2^T *Y1 - (S2^T *Y1)^T - // W = W - S2^T*Y1 - S1.transTimesMat(1.0,W,-1.0,Y1); - // W = W - Y1*S2^T - Y1.transTimesMat(1.0,W,-1.0,S1); - - // -- done -- -#ifdef HIOP_DEEPCHECKS - W.print(); - assert(W.assertSymmetry(1e-14)); -#endif - -} - -/* symmetric multiplication W = beta*W + alpha*X*Diag*X^T - * W is kxk local, X is kxn distributed and Diag is n, distributed - * The ops are perform locally. The reduce is done separately/externally to decrease comm - */ -void hiopHessianInvLowRank_obsolette:: -symmMatTimesDiagTimesMatTrans_local(double beta, hiopMatrixDense& W, - double alpha, const hiopMatrixDense& X, - const hiopVector& d) -{ - size_type k=W.m(); - size_type n_local=X.get_local_size_n(); - - assert(X.m()==k); - -#ifdef HIOP_DEEPCHECKS - assert(W.n()==k); - assert(d.get_size()==X.n()); - assert(d.get_local_size()==n_local); -#endif - //#define chunk 512; //!opt - const double *xi, *xj; - double acc; - double* Wdata=W.local_data(); - const double* Xdata=X.local_data_const(); - const double* dd=d.local_data_const(); - for(int i=0; im(); -#ifdef HIOP_DEEPCHECKS - assert(l==R->n()); - assert(lmem_curr==l); - assert(lmem_max>l); -#endif - //newR = [ R S^T*y ] - // [ 0 s^T*y ] - hiopMatrixDense* newR = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l+1, l+1); - assert(newR); - //copy from R to newR - newR->copyBlockFromMatrix(0,0, *R); - - double* newR_mat=newR->local_data(); //doing the rest here - const double* STy_vec=STy.local_data_const(); - //for(int j=0; jget_size(); - assert(l==lmem_curr); - assert(lmem_max>l); - - hiopVector* Dnew=LinearAlgebraFactory::create_vector("DEFAULT", l+1); - double* Dnew_vec=Dnew->local_data(); - memcpy(Dnew_vec, D->local_data_const(), l*sizeof(double)); - Dnew_vec[l]=sTy; - - delete D; - D=Dnew; -} - -void hiopHessianInvLowRank_obsolette::updateR(const hiopVector& STy, const double& sTy) -{ - int l=STy.get_size(); - assert(l==R->m()); -#ifdef HIOP_DEEPCHECKS - assert(l==R->n()); -#endif - const int lm1=l-1; - double* R_mat=R->local_data(); - const double* sTy_vec=STy.local_data_const(); - for(int i=0; iget_size(); - double* D_vec = D->local_data(); - for(int i=0; in()==k); -#endif - if(NULL!=_S1 && _S1->n()!=l) { delete _S1; _S1=NULL; } - - if(NULL==_S1) _S1=LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, k); - - return *_S1; -} - -hiopMatrixDense& hiopHessianInvLowRank_obsolette::new_Y1(const hiopMatrixDense& Yt, const hiopMatrixDense& X) -{ - //Y1 is Yt*H0*X^T = Y^T*H0*X^T, where Y^T is lxn, H0 is diag nxn, X is kxn - size_type k=X.m(), l=Yt.m(); -#ifdef HIOP_DEEPCHECKS - const size_type n=Yt.n(); - assert(X.n()==n); - if(_Y1!=NULL) assert(_Y1->n()==k); -#endif - - if(NULL!=_Y1 && _Y1->n()!=l) { delete _Y1; _Y1=NULL; } - - if(NULL==_Y1) _Y1=LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, k); - return *_Y1; -} -hiopMatrixDense& hiopHessianInvLowRank_obsolette::new_DpYtH0Y(const hiopMatrixDense& Yt) -{ - size_type l = Yt.m(); -#ifdef HIOP_DEEPCHECKS - if(_DpYtH0Y!=NULL) assert(_DpYtH0Y->m()==_DpYtH0Y->n()); -#endif - if(_DpYtH0Y!=NULL && _DpYtH0Y->m()!=l) {delete _DpYtH0Y; _DpYtH0Y=NULL;} - if(_DpYtH0Y==NULL) _DpYtH0Y=LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, l); - return *_DpYtH0Y; -} - -/* S3 = DpYtH0H * S2, where S2=R\S1. DpYtH0H is symmetric (!opt) lxl and S2 is lxk */ -hiopMatrixDense& hiopHessianInvLowRank_obsolette::new_S3(const hiopMatrixDense& Left, const hiopMatrixDense& Right) -{ - int l=Left.m(), k=Right.n(); -#ifdef HIOP_DEEPCHECKS - assert(Right.m()==l); - assert(Left.n()==l); - if(_S3!=NULL) assert(_S3->m()<=l); //< when the representation grows, = when it doesn't -#endif - if(_S3!=NULL && _S3->m()!=l) { - delete _S3; - _S3=NULL; - } - - if(_S3==NULL) { - _S3 = LinearAlgebraFactory::create_matrix_dense("DEFAULT", l, k); - } - return *_S3; -} -hiopVector& hiopHessianInvLowRank_obsolette::new_l_vec1(int l) -{ - if(_l_vec1!=NULL && _l_vec1->get_size()==l) { - return *_l_vec1; - } - if(_l_vec1!=NULL) { - delete _l_vec1; - } - _l_vec1 = LinearAlgebraFactory::create_vector("DEFAULT", l); - return *_l_vec1; -} -hiopVector& hiopHessianInvLowRank_obsolette::new_l_vec2(int l) -{ - if(_l_vec2!=NULL && _l_vec2->get_size()==l) { - return *_l_vec2; - } - if(_l_vec2!=NULL) { - delete _l_vec2; - } - _l_vec2 = LinearAlgebraFactory::create_vector("DEFAULT", l); - return *_l_vec2; -} -hiopVector& hiopHessianInvLowRank_obsolette::new_l_vec3(int l) -{ - if(_l_vec3!=NULL && _l_vec3->get_size()==l) return *_l_vec3; - - if(_l_vec3 != NULL) { - delete _l_vec3; - } - _l_vec3 = LinearAlgebraFactory::create_vector("DEFAULT", l); - return *_l_vec3; -} - - //#ifdef HIOP_DEEPCHECKS -// #include -// using namespace std; -// void hiopHessianInvLowRank_obsolette::timesVecCmn(double beta, hiopVector& y, double alpha, const hiopVector& x, bool addLogTerm) -// { -// long long n=St->n(); -// assert(l_curr-1==St->m()); -// assert(y.get_size()==n); -// //we have B+=B-B*s*B*s'/(s'*B*s)+yy'/(y'*s) -// //B0 is sigma*I (and is NOT this->H0, since this->H0=(B0+Dx)^{-1}) - -// bool print=true; -// if(print) { -// nlp->log->printf(hovMatrices, "---hiopHessianInvLowRank_obsolette::timesVec \n"); -// nlp->log->write("S':", *St, hovMatrices); -// nlp->log->write("Y':", *Yt, hovMatrices); -// nlp->log->write("H0:", *H0, hovMatrices); -// nlp->log->printf(hovMatrices, "sigma=%22.16e addLogTerm=%d\n", sigma, addLogTerm); -// nlp->log->printf(hovMatrices, "y=beta*y + alpha*this*x : beta=%g alpha=%g\n", beta, alpha); -// nlp->log->write("x_in:", x, hovMatrices); -// nlp->log->write("y_in:", y, hovMatrices); -// } - -// hiopVectorPar *yk=dynamic_cast(nlp->alloc_primal_vec()); -// hiopVectorPar *sk=dynamic_cast(nlp->alloc_primal_vec()); -// //allocate and compute a_k and b_k -// vector a(l_curr),b(l_curr); -// for(int k=0; kcopyFrom(Yt->local_data()[k]); -// sk->copyFrom(St->local_data()[k]); -// double skTyk=yk->dotProductWith(*sk); -// assert(skTyk>0); -// b[k]=dynamic_cast(nlp->alloc_primal_vec()); -// b[k]->copyFrom(*yk); -// b[k]->scale(1/sqrt(skTyk)); - -// a[k]=dynamic_cast(nlp->alloc_primal_vec()); - -// //compute ak by an inner loop -// a[k]->copyFrom(*sk); -// if(addLogTerm) -// a[k]->componentDiv(*H0); -// else -// a[k]->scale(sigma); -// for(int i=0; idotProductWith(*sk); -// a[k]->axpy(biTsk, *b[i]); -// double aiTsk = a[i]->dotProductWith(*sk); -// a[k]->axpy(aiTsk, *a[i]); -// } -// double skTak = a[k]->dotProductWith(*sk); -// a[k]->scale(1/sqrt(skTak)); -// } - -// //new we have B= Dx+B_0 + sum{ bk bk' - ak ak' : k=0,1,...,l_curr-1} (H0=(Dx+B0)^{-1}) -// //compute the product with x -// //y = beta*y+alpha*H0_inv*x + alpha* sum { bk'x bk - ak'x ak : k=0,1,...,l_curr-1} -// y.scale(beta); -// if(addLogTerm) -// y.axdzpy(alpha,x,*H0); -// else -// y.axpy(alpha*sigma, x); -// for(int k=0; kdotProductWith(x); -// double akTx = a[k]->dotProductWith(x); - -// y.axpy( alpha*bkTx, *b[k]); -// y.axpy(-alpha*akTx, *a[k]); -// } - -// if(print) { -// nlp->log->write("y_out:", y, hovMatrices); -// } - -// for(vector::iterator it=a.begin(); it!=a.end(); ++it) -// delete *it; -// for(vector::iterator it=b.begin(); it!=b.end(); ++it) -// delete *it; - -// delete yk; -// delete sk; -// } - -// void hiopHessianInvLowRank_obsolette::timesVec(double beta, hiopVector& y, double alpha, const hiopVector&x) -// { -// this->timesVecCmn(beta, y, alpha, x, true); -// } - -// void hiopHessianInvLowRank_obsolette::timesVec_noLogBarrierTerm(double beta, hiopVector& y, double alpha, const hiopVector&x) -// { -// this->timesVecCmn(beta, y, alpha, x, false); -// } - -//#endif - }; diff --git a/src/Optimization/hiopHessianLowRank.hpp b/src/Optimization/hiopHessianLowRank.hpp index b487ea086..05457c1e4 100644 --- a/src/Optimization/hiopHessianLowRank.hpp +++ b/src/Optimization/hiopHessianLowRank.hpp @@ -91,7 +91,7 @@ namespace hiop class hiopHessianLowRank : public hiopMatrix { public: - hiopHessianLowRank(hiopNlpDenseConstraints* nlp_, int max_memory_length); + hiopHessianLowRank(hiopNlpDenseConstraints* nlp_in, int max_memory_length); virtual ~hiopHessianLowRank(); /* return false if the update destroys hereditary positive definitness and the BFGS update is not taken*/ @@ -124,6 +124,7 @@ class hiopHessianLowRank : public hiopMatrix virtual void timesVecCmn(double beta, hiopVector& y, double alpha, const hiopVector&x, bool addLogBarTerm = false) const; protected: + friend class hiopAlgFilterIPMQuasiNewton; int l_max; //max memory size int l_curr; //number of pairs currently stored double sigma; //initial scaling factor of identity @@ -143,12 +144,14 @@ class hiopHessianLowRank : public hiopMatrix bool matrixChanged; //these are matrices from the compact representation; they are updated at each iteration. // more exactly Bk=B0-[B0*St' Yt']*[St*B0*St' L]*[St*B0] - // [ L' -D] [Yt ] - hiopMatrixDense *St,*Yt; //we store the transpose to easily access columns in S and T - hiopMatrixDense *L; //lower triangular from the compact representation - hiopVector* D; //diag + // [ L' -D] [Yt ] + //transpose of S and T are store to easily access columns + hiopMatrixDense* St_; + hiopMatrixDense* Yt_; + hiopMatrixDense* L_; //lower triangular from the compact representation + hiopVector* D_; //diag //these are matrices from the representation of the inverse - hiopMatrixDense* V; + hiopMatrixDense* V_; #ifdef HIOP_DEEPCHECKS //copy of the matrix - needed to check the residual hiopMatrixDense* _Vmat; @@ -158,9 +161,10 @@ class hiopHessianLowRank : public hiopMatrix void updateL(const hiopVector& STy, const double& sTy); void updateD(const double& sTy); //also stored are the iterate, gradient obj, and Jacobians at the previous optimization iteration - hiopIterate *_it_prev; - hiopVector *_grad_f_prev; - hiopMatrixDense *_Jac_c_prev, *_Jac_d_prev; + hiopIterate *it_prev_; + hiopVector *grad_f_prev_; + hiopMatrixDense *Jac_c_prev_; + hiopMatrixDense *Jac_d_prev_; //internal helpers void updateInternalBFGSRepresentation(); @@ -206,6 +210,10 @@ class hiopHessianLowRank : public hiopMatrix } private: //utilities + + /// @brief Ensures the internal containers are ready to work with "limited memory" mem_length + void alloc_for_limited_mem(const size_type& mem_length); + /* symmetric multiplication W = beta*W + alpha*X*Diag*X^T */ static void symmMatTimesDiagTimesMatTrans_local(double beta, hiopMatrixDense& W_, double alpha, const hiopMatrixDense& X_, @@ -223,7 +231,7 @@ class hiopHessianLowRank : public hiopMatrix hiopHessianLowRank() {}; hiopHessianLowRank(const hiopHessianLowRank&) {}; hiopHessianLowRank& operator=(const hiopHessianLowRank&) {return *this;}; -public: + /* methods that need to be implemented as the class inherits from hiopMatrix*/ public: virtual hiopMatrix* alloc_clone() const @@ -363,110 +371,6 @@ class hiopHessianLowRank : public hiopMatrix /* check symmetry */ virtual bool assertSymmetry(double tol=1e-16) const { return true; } #endif - -}; - -/* Low-rank representation for the inverse of a low-rank matrix (Hessian). - * In HIOP we need to solve with the BFGS quasi-Newton Hessian given by - * the matrix M=B0+[B0*S Y] [S^T*B0*S L] [ S^T*B0 ] - * [ L^T -D] [ Y ] - * Reference: Byrd, Nocedal, Schnabel, "Representations of quasi-Newton matrices and - * and there use in limited memory methods", Math. Programming 63 (1994), p. 129-156. - * - * To save on computations, we maintain a direct representaton of its inverse - * M^{-1} = H0 + [S HO*Y] [ R^{-T}*(D+Y^T*H0*Y)*R^{-1} -R^{-T} ] [ S^T ] - * [ -R^{-1} 0 ] [ Y^T*H0] - * Here, H0=inv(B0). Check the above reference to see what each matrix represents. - */ -class hiopHessianInvLowRank_obsolette : public hiopHessianLowRank -{ -public: - hiopHessianInvLowRank_obsolette(hiopNlpDenseConstraints* nlp, int max_memory_length); - virtual bool update(const hiopIterate& x_curr, const hiopVector& grad_f_curr, - const hiopMatrix& Jac_c_curr, const hiopMatrix& Jac_d_curr); - - virtual bool updateLogBarrierDiagonal(const hiopVector& Dx); - - /* ! - * these methods use quantities computed in symmetricTimesMat. They should be called - * after symmetricTimesMat - */ - virtual void apply(double beta, hiopVector& y, double alpha, const hiopVector& x); - virtual void apply(double beta, hiopMatrix& Y, double alpha, const hiopMatrix& X); - - /* W = beta*W + alpha*X*this*X^T - * ! make sure this is called before 'apply' - */ - virtual void symmetricTimesMat(double beta, hiopMatrixDense& W, - double alpha, const hiopMatrixDense& X); - - virtual ~hiopHessianInvLowRank_obsolette(); - -private: //internal methods - - /* symmetric multiplication W = beta*W + alpha*X*Diag*X^T */ - static void symmMatTimesDiagTimesMatTrans_local(double beta, hiopMatrixDense& W_, - double alpha, const hiopMatrixDense& X_, - const hiopVector& d); - /* W=S*Diag*X^T */ - static void matTimesDiagTimesMatTrans_local(hiopMatrixDense& W, const hiopMatrixDense& S, const hiopVector& d, const hiopMatrixDense& X); - - /* rhs = R \ rhs, where R is upper triangular lxl and rhs is lx */ - static void triangularSolve(const hiopMatrixDense& R, hiopMatrixDense& rhs); - static void triangularSolve(const hiopMatrixDense& R, hiopVector& rhs); - static void triangularSolveTrans(const hiopMatrixDense& R, hiopVector& rhs); - - //grows R when the number of BFGS updates is less than the max memory - void growR(const int& l_curr, const int& l_max, const hiopVector& STy, const double& sTy); - void growD(const int& l_curr, const int& l_max, const double& sTy); - void updateR(const hiopVector& STy, const double& sTy); - void updateD(const double& sTy); -private: - hiopVector* H0; - hiopMatrixDense *St,*Yt; //we store the transpose to easily access columns in S and T - hiopMatrixDense *R; - hiopVector* D; - - int sigma_update_strategy; - double sigma_safe_min, sigma_safe_max; - - //also stored are the iterate, gradient obj, and Jacobians at the previous iterations - hiopIterate *_it_prev; - hiopVector *_grad_f_prev; - hiopMatrixDense *_Jac_c_prev, *_Jac_d_prev; - - //internals buffers - double* _buff_kxk; // size = num_constraints^2 - double* _buff_lxk; // size = q-Newton mem size x num_constraints - double* _buff_lxl; - //auxiliary objects - hiopMatrixDense *_S1, *_Y1, *_DpYtH0Y; //aux matrices to hold St*X, Yt*H0*X, and D+Y^T*H0*Y in symmetricTimesMat - hiopMatrixDense& new_S1(const hiopMatrixDense& St, const hiopMatrixDense& X); - hiopMatrixDense& new_Y1(const hiopMatrixDense& Yt, const hiopMatrixDense& X); - hiopMatrixDense& new_DpYtH0Y(const hiopMatrixDense& Yt); - //similar for S3=DpYtH0Y*S2 - hiopMatrixDense *_S3; - hiopMatrixDense& new_S3(const hiopMatrixDense& Left, const hiopMatrixDense& Right); - hiopVector *_l_vec1, *_l_vec2, *_l_vec3, *_n_vec1, *_n_vec2; - hiopVector& new_l_vec1(int l); - hiopVector& new_l_vec2(int l); - hiopVector& new_l_vec3(int l); - inline hiopVector& new_n_vec1(size_type n) - { -#ifdef HIOP_DEEPCHECKS - assert(_n_vec1!=NULL); - assert(_n_vec1->get_size()==n); -#endif - return *_n_vec1; - } - inline hiopVector& new_n_vec2(size_type n) - { -#ifdef HIOP_DEEPCHECKS - assert(_n_vec2!=NULL); - assert(_n_vec2->get_size()==n); -#endif - return *_n_vec2; - } }; } //~namespace diff --git a/src/Optimization/hiopKKTLinSys.cpp b/src/Optimization/hiopKKTLinSys.cpp index 12e265003..d91a3d64d 100644 --- a/src/Optimization/hiopKKTLinSys.cpp +++ b/src/Optimization/hiopKKTLinSys.cpp @@ -947,7 +947,6 @@ bool hiopKKTLinSys::compute_directions_w_IR(const hiopResidual* resid, hiopItera bool bret = bicgIR_->solve(dir, resid); nlp_->runStats.kkt.nIterRefinInner += bicgIR_->get_sol_num_iter(); - nlp_->runStats.kkt.tmSolveInner.stop(); if(!bret) { nlp_->log->printf(hovWarning, "%s", bicgIR_->get_convergence_info().c_str()); diff --git a/src/Utils/SidreHelper.hpp b/src/Utils/SidreHelper.hpp new file mode 100644 index 000000000..69bca66bf --- /dev/null +++ b/src/Utils/SidreHelper.hpp @@ -0,0 +1,285 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC. +// Produced at the Lawrence Livermore National Laboratory (LLNL). +// LLNL-CODE-742473. All rights reserved. +// +// This file is part of HiOp. For details, see https://github.com/LLNL/hiop. HiOp +// is released under the BSD 3-clause license (https://opensource.org/licenses/BSD-3-Clause). +// Please also read "Additional BSD Notice" below. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// i. Redistributions of source code must retain the above copyright notice, this list +// of conditions and the disclaimer below. +// ii. Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the disclaimer (as noted below) in the documentation and/or +// other materials provided with the distribution. +// iii. Neither the name of the LLNS/LLNL nor the names of its contributors may be used to +// endorse or promote products derived from this software without specific prior written +// permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY +// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES +// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +// SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS +// OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED +// AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, +// EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +// +// Additional BSD Notice +// 1. This notice is required to be provided under our contract with the U.S. Department +// of Energy (DOE). This work was produced at Lawrence Livermore National Laboratory under +// Contract No. DE-AC52-07NA27344 with the DOE. +// 2. Neither the United States Government nor Lawrence Livermore National Security, LLC +// nor any of their employees, makes any warranty, express or implied, or assumes any +// liability or responsibility for the accuracy, completeness, or usefulness of any +// information, apparatus, product, or process disclosed, or represents that its use would +// not infringe privately-owned rights. +// 3. Also, reference herein to any specific commercial products, process, or services by +// trade name, trademark, manufacturer or otherwise does not necessarily constitute or +// imply its endorsement, recommendation, or favoring by the United States Government or +// Lawrence Livermore National Security, LLC. The views and opinions of authors expressed +// herein do not necessarily state or reflect those of the United States Government or +// Lawrence Livermore National Security, LLC, and shall not be used for advertising or +// product endorsement purposes. + +/** + * @file SidreHelper.hpp + * + * @author Cosmin G. Petra , LLNL + * @author Tom Epperly , LLNL + * + */ +#ifndef HIOP_SIDRE_HELP +#define HIOP_SIDRE_HELP + +#include "hiopVector.hpp" + +#include +#include +#include +#include + +#include +#include + +namespace hiop +{ +/** + * @brief Holder of functionality needed by HiOp for checkpointing based on axom::sidre + */ +class SidreHelper +{ +public: + /** + * @brief Copy raw array to sidre::View within specified sidre::Group. + * + * @params group contains the view where the copy should be made to. + * @params view_name is the name of the view where to copy + * @params arr_src is the source double array + * @params size is the number of elements of the array + * + * @exception std::runtime indicates the group contains a view with a number of elements + * different than expected size. + * + * @details A view with the specified name will be created if does not already exist. If + * exists, the view should have the same number of elements as the argument `size`. + */ + + static void copy_array_to_view(::axom::sidre::Group& group, + const ::std::string& view_name, + const double* arr_src, + const hiop::size_type& size) + { + auto view = get_or_create_view(group, view_name, size); + if(view->getNumElements() != size) { + ::std::stringstream ss; + ss << "Size mismatch between HiOp state and existing sidre::view '" << view_name << + "' when copying to view. HiOp state has " << size << " doubles, while the view " << + "has " << view->getNumElements() << " double elements.\n"; + throw ::std::runtime_error(ss.str()); + } + + const auto stride(view->getStride()); + double *const arr_dest(view->getArray()); + if(1==stride) { + ::std::copy(arr_src, arr_src+size, arr_dest); + } else { + for(::axom::sidre::IndexType i=0; igetNumElements() != size) { + ::std::stringstream ss; + ss << "Size mismatch between HiOp state and sidre::View '" << view_name << + "' when copying from the view. HiOp state is " << size << " doubles, "<< + "while the view has " << view_const->getNumElements() << " double elements.\n"; + throw ::std::runtime_error(ss.str()); + } + + // const_cast becase View does not have a const getArray() + auto view = const_cast<::axom::sidre::View*>(view_const); + assert(view); + const double *const arr_src = view->getArray(); + const auto stride(view->getStride()); + + if(1==stride) { + ::std::copy(arr_src, arr_src+size, arr_dest); + } else { + for(hiop::index_type i=0; i range = {"yes", "no"}; + constexpr char msgcs[] = "Save state of NLP solver to file indicated by 'checkpoint_file'."; + register_str_option("checkpoint_save", range[1], range, msgcs); + + constexpr char msgcsN[] = "Iteration frequency of saving checkpoints to disk."; + register_int_option("checkpoint_save_every_N_iter", 10, 1, 1e+6, msgcsN); + constexpr char msgcf[] = "Path to checkpoint file to load from or save to."; + register_str_option("checkpoint_file", "hiop_state_chk", msgcf); + + constexpr char msgclos[] = "On (re)start the NLP solver will load checkpoint file " + "specified by 'checkpoint_file' option."; + register_str_option("checkpoint_load_on_start", range[1], range, msgclos); + } +} void hiopOptionsNLP::ensure_consistence() { //check that the values of different options are consistent @@ -1551,6 +1568,26 @@ void hiopOptionsNLP::ensure_consistence() } set_val("moving_lim_rel", 0.); } + +#ifndef HIOP_USE_AXOM + const vector chkpnt_opts = {"checkpoint_save", + "checkpoint_save_every_N_iter", + "checkpoint_file", + "checkpoint_load_on_start"}; + for(string opt : chkpnt_opts) { + if(is_user_defined(opt.c_str())) { + log_printf(hovWarning, + "Checkpointing not available since HiOp was not built with AXOM. All checkpointing options " + "are ignored.\n"); + //reset them to as not being user defined to avoid triggering the message. + for(auto opt2 : chkpnt_opts) { + mOptions_[opt2]->specifiedInFile = false; + mOptions_[opt2]->specifiedAtRuntime = false; + } + break; + } + } +#endif } /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////