From 7a939b7dc25fe9463339afd4824228d2373552e4 Mon Sep 17 00:00:00 2001 From: "Cosmin G. Petra" Date: Mon, 26 Aug 2024 15:38:42 -0700 Subject: [PATCH 01/25] added notes and some drafty interface --- src/Interface/hiopInterface.hpp | 23 ++++++++- src/Optimization/hiopAlgFilterIPM.hpp | 72 +++++++++++++++++++++++++++ 2 files changed, 93 insertions(+), 2 deletions(-) diff --git a/src/Interface/hiopInterface.hpp b/src/Interface/hiopInterface.hpp index f0b6cece..bea73076 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) @@ -497,6 +497,25 @@ class hiopInterfaceBase return true; } + /** + * This method is called after each iteration and should be implemented by the user to instruct HiOp + * to save a complete state of the algorithm to disk. An axom::sidre::DataStore will be used for IO, which + * can be set by the user upon creation of the HiOp algorithm class. If not set, HiOp will used one. + * + * This is feature is useful for IO checkpointing, for example; it allows the internal algorithm to be + * restarted at the same state as before saving. @see hiop::hiopAlgorithm::tbd for the method to be used + * to load a previously saved state. + * + * The method provided by the user should return true if HiOp should save the state of the algorithm at + * at current iteration. The argument is passed by HiOp to indicate the current iterate number. + * + * @param[in] iter the current iteration number + */ + virtual bool save_state_iterate_callback(int iter) + { + return false; + } + /** * A wildcard function used to change the primal variables. * diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index 11902144..f3b636bf 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -117,6 +117,75 @@ class hiopAlgFilterIPMBase { { return filter.contains(theta, logbar_obj); } + + //notes on checkpointing + // 1. need to allow user to pass axom::sidre::DataStore. HiOp will put all the info into a group + // Question: should a DataStore be passed any time restarting is invoked (this is a bit cumbersome), + // or just once, by calling the HiOp algorithm class + // + // + // + // + // (for when checkpointing is used without the user setting a data store, so HiOp will create it and do the IO) + // 2. need to allow user to pass a string with the file where the DataStore will be + // writting to/reading from. A default name will be used for empty filename. + + /** + * Setter for the axom::sidre DataStore used to manage the data associated with the state of + * NLP the algorithm. If the setter is not called by the user, the DataStore will be created + * internally by the iteration checkpointing object. + */ + inline void set_state_data_manager(axom::sidre::DataStore& mng) + { + iter_chkpnt_.set_data_manger(mng); + } + + /** + * The method saves the state of the algorithm in the axom::sidre::DataStore object that was + * previously provided by @set_checkpoint_data_manager. If this method has not been previously + * called, the HiOp will create such instance and will save it on disk under a default name. + */ + inline void save_state() + { + iter_chkpnt_.save(this); + } + + /** + * The method saves the state of the algorithm in the file specified by the string argument. If + * the string is empty, a file with a default name will be created. + * + * Internally, HiOp uses axom::sidre::DataStore object that is created internally and shares the + * IO code with @save_state. This method disregards previous calls to @set_checkpoint_data_manager. + */ + inline void save_state(const std::string& filename) + { + iter_chkpnt_.save(filename, this); + } + + + /** + * This method loads the state of the algorithm from a axom::sidre::DataStore that was previously + * provided by @set_checkpoint_data_manager. This DataStore instance needs to be properly + * initialized and have a group called "HiOpState". + */ + inline void load_state() + { + iter_chkpnt_.load(this); + } + + /** + * This method loads the state of the algorithm from the file whose name is passed as a string + * argument. HiOp expected that the file contains a axom::sidre::DataStore that was previously saved + * using one of the @save_state methods above. + * + */ + inline void load_state(const std::string& filename) + { + iter_chkpnt_.load(filename, this); + } + + + /// Setter for the primal steplength. inline void set_alpha_primal(const double alpha_primal) { _alpha_primal = alpha_primal; } protected: @@ -245,6 +314,9 @@ class hiopAlgFilterIPMBase { hiopNlpFormulation* nlp; hiopFilter filter; + /// Helper for saving/loading algorithm state to disk. + Checkpointing iter_chkpnt_; + hiopLogBarProblem* logbar; /* Iterate, search directions (managed by this (algorithm) class) */ From 181f7f9f735c60fecdd40b9e4706d2559c939cd7 Mon Sep 17 00:00:00 2001 From: "Cosmin G. Petra" Date: Wed, 28 Aug 2024 11:26:24 -0700 Subject: [PATCH 02/25] added draft of the api for checkpointing --- CMakeLists.txt | 12 ++- src/Interface/hiopInterface.hpp | 7 +- src/Optimization/hiopAlgFilterIPM.cpp | 45 ++++++++++++ src/Optimization/hiopAlgFilterIPM.hpp | 102 ++++++++------------------ 4 files changed, 92 insertions(+), 74 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 8b2db654..06787ee6 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,15 @@ if(HIOP_USE_RAJA) message(STATUS "Found umpire pkg-config: ${umpire_CONFIG}") endif() +if(HIOP_USE_AXOM) + 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}") +endif() + + if(HIOP_WITH_KRON_REDUCTION) set(HIOP_UMFPACK_DIR CACHE PATH "Path to UMFPACK directory") include(FindUMFPACK) diff --git a/src/Interface/hiopInterface.hpp b/src/Interface/hiopInterface.hpp index bea73076..0d7b4ef2 100644 --- a/src/Interface/hiopInterface.hpp +++ b/src/Interface/hiopInterface.hpp @@ -500,14 +500,15 @@ class hiopInterfaceBase /** * This method is called after each iteration and should be implemented by the user to instruct HiOp * to save a complete state of the algorithm to disk. An axom::sidre::DataStore will be used for IO, which - * can be set by the user upon creation of the HiOp algorithm class. If not set, HiOp will used one. + * can be set by the user upon creation of the HiOp algorithm class. If not set, HiOp will create one + * internally. TODO: mention API methods from the Alg class. * - * This is feature is useful for IO checkpointing, for example; it allows the internal algorithm to be + * This feature is useful for IO checkpointing, for example; it allows the internal algorithm to be * restarted at the same state as before saving. @see hiop::hiopAlgorithm::tbd for the method to be used * to load a previously saved state. * * The method provided by the user should return true if HiOp should save the state of the algorithm at - * at current iteration. The argument is passed by HiOp to indicate the current iterate number. + * at current iteration. The argument is passed HiOp to indicate the current iterate number. * * @param[in] iter the current iteration number */ diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index 2ac5fcc8..3af3107c 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -1485,6 +1485,51 @@ void hiopAlgFilterIPMQuasiNewton::outputIteration(int lsStatus, int lsNum, int u } } +//declaration required by C++14, not anymore by C++17 or after +constexpr char hiopAlgFilterIPMQuasiNewton::default_state_filename[]; + +void hiopAlgFilterIPMQuasiNewton::save_state_to_file(const ::std::string& path_in) +{ + auto path = path_in=="" ? default_state_filename : path_in; + + axom_sidre_DataStore* ds = new axom_sidre_DataStore(0); + + this->save_state_to_data_store(ds); + + //::axom::sidre::IOManager writer(this->get_nlp()->get_comm()); + //int num_files; + //MPI_Comm_size(this->get_nlp()->get_comm(), &num_files); + //writer.write(ds_->getRoot(), num_files, path.str(), ::axom::sidre::Group::getDefaultIOProtocol()); + + delete ds; +} + +void hiopAlgFilterIPMQuasiNewton::load_state_from_file(const ::std::string& path_in) +{ + auto path = path_in=="" ? default_state_filename : path_in; + //todo +} + +void hiopAlgFilterIPMQuasiNewton::save_state_to_data_store(void* ds) +{ + //Group* nlp_group = ds->getRoot()->createGroup("hiop solver"); + + //create views for each member that needs to be saved + + const double* x = it_curr->get_x()->local_data_host(); + //destination = nlp_group->createViewAndAllocate("x", ::axom::sidre::DOUBLE_ID, size); +} + +void hiopAlgFilterIPMQuasiNewton::load_state_from_data_store(const void* ds) +{ + //Group* nlp_group = ds->getRoot()->createGroup("hiop solver"); + + //create views for each member that needs to be saved + + const double* x = it_curr->get_x()->local_data_host(); + //destination = nlp_group->createViewAndAllocate("x", ::axom::sidre::DOUBLE_ID, size); +} + /****************************************************************************************************** * FULL NEWTON IPM diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index f3b636bf..60a0f669 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -67,6 +67,8 @@ #include "hiopPDPerturbation.hpp" #include "hiopFactAcceptor.hpp" +#include "Checkpointing.hpp" + #include "hiopTimer.hpp" namespace hiop @@ -118,73 +120,6 @@ class hiopAlgFilterIPMBase { return filter.contains(theta, logbar_obj); } - //notes on checkpointing - // 1. need to allow user to pass axom::sidre::DataStore. HiOp will put all the info into a group - // Question: should a DataStore be passed any time restarting is invoked (this is a bit cumbersome), - // or just once, by calling the HiOp algorithm class - // - // - // - // - // (for when checkpointing is used without the user setting a data store, so HiOp will create it and do the IO) - // 2. need to allow user to pass a string with the file where the DataStore will be - // writting to/reading from. A default name will be used for empty filename. - - /** - * Setter for the axom::sidre DataStore used to manage the data associated with the state of - * NLP the algorithm. If the setter is not called by the user, the DataStore will be created - * internally by the iteration checkpointing object. - */ - inline void set_state_data_manager(axom::sidre::DataStore& mng) - { - iter_chkpnt_.set_data_manger(mng); - } - - /** - * The method saves the state of the algorithm in the axom::sidre::DataStore object that was - * previously provided by @set_checkpoint_data_manager. If this method has not been previously - * called, the HiOp will create such instance and will save it on disk under a default name. - */ - inline void save_state() - { - iter_chkpnt_.save(this); - } - - /** - * The method saves the state of the algorithm in the file specified by the string argument. If - * the string is empty, a file with a default name will be created. - * - * Internally, HiOp uses axom::sidre::DataStore object that is created internally and shares the - * IO code with @save_state. This method disregards previous calls to @set_checkpoint_data_manager. - */ - inline void save_state(const std::string& filename) - { - iter_chkpnt_.save(filename, this); - } - - - /** - * This method loads the state of the algorithm from a axom::sidre::DataStore that was previously - * provided by @set_checkpoint_data_manager. This DataStore instance needs to be properly - * initialized and have a group called "HiOpState". - */ - inline void load_state() - { - iter_chkpnt_.load(this); - } - - /** - * This method loads the state of the algorithm from the file whose name is passed as a string - * argument. HiOp expected that the file contains a axom::sidre::DataStore that was previously saved - * using one of the @save_state methods above. - * - */ - inline void load_state(const std::string& filename) - { - iter_chkpnt_.load(filename, this); - } - - /// Setter for the primal steplength. inline void set_alpha_primal(const double alpha_primal) { _alpha_primal = alpha_primal; } @@ -314,9 +249,6 @@ class hiopAlgFilterIPMBase { hiopNlpFormulation* nlp; hiopFilter filter; - /// Helper for saving/loading algorithm state to disk. - Checkpointing iter_chkpnt_; - hiopLogBarProblem* logbar; /* Iterate, search directions (managed by this (algorithm) class) */ @@ -411,6 +343,36 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase virtual ~hiopAlgFilterIPMQuasiNewton(); virtual hiopSolveStatus run(); + + //work in progress + virtual void save_state_to_data_store(void* sidre_data_store); + virtual void load_state_from_data_store(const void* sidre_data_store); + + static constexpr char default_state_filename[] = "hiop_qn_state.sidre"; + + /** + * @brief save the state of the algorithm to the file + * @param path the name of the file + * + * @details + * Internally, HiOp uses axom::sidre::DataStore, which is saved to the file. If argument is the + * empty string, HiOp will attempt saving the state to the path specified by default_state_filename + * static member. + */ + void save_state_to_file(const ::std::string& path=""); + + /** + * @brief load the state of the algorithm from file + * @param path the name of the file to load from + * + * @details + * The file should contains a axom::sidre::DataStore that was previously saved using save_state_to_file(). + * If argument is the empty string, HiOp will attempt loading state from the path specified by + * default_state_filename static member. + * + */ + void load_state_from_file(const ::std::string& path=""); + private: virtual void outputIteration(int lsStatus, int lsNum, int use_soc = 0, int use_fr = 0); private: From 0e428f532991f81ea4477e4684aced50b861cfe3 Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Wed, 28 Aug 2024 11:30:47 -0700 Subject: [PATCH 03/25] fixed compilation issues --- src/Optimization/hiopAlgFilterIPM.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index 60a0f669..4be26885 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -67,13 +67,13 @@ #include "hiopPDPerturbation.hpp" #include "hiopFactAcceptor.hpp" -#include "Checkpointing.hpp" - #include "hiopTimer.hpp" namespace hiop { - +//temporary dummy definition for axom::sidre::DataStore until axom is going to be included +using axom_sidre_DataStore=double; + class hiopAlgFilterIPMBase { public: hiopAlgFilterIPMBase(hiopNlpFormulation* nlp_, const bool within_FR = false); From 354b82be87d067595160a605b768553391a55d96 Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Wed, 28 Aug 2024 14:08:28 -0700 Subject: [PATCH 04/25] integrated AXOM code compiles; does nothing --- src/Interface/hiop_defs.hpp.in | 1 + src/Optimization/hiopAlgFilterIPM.cpp | 44 +++++++++++++++++++++------ src/Optimization/hiopAlgFilterIPM.hpp | 16 ++++++---- src/Utils/hiopOptions.cpp | 2 ++ 4 files changed, 47 insertions(+), 16 deletions(-) diff --git a/src/Interface/hiop_defs.hpp.in b/src/Interface/hiop_defs.hpp.in index 0af21826..68570c8b 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 3af3107c..a179341a 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -66,6 +66,14 @@ #include "hiopCppStdUtils.hpp" +#ifdef HIOP_USE_AXOM +#include +#include +#include +#include +using namespace axom; +#endif + #include #include #include @@ -74,6 +82,9 @@ namespace hiop { +//#ifdef HIOP_USE_AXOM + +//#ifdef HIOP_USE_AXOM hiopAlgFilterIPMBase::hiopAlgFilterIPMBase(hiopNlpFormulation* nlp_in, const bool within_FR) : soc_dir(nullptr), @@ -1485,6 +1496,8 @@ void hiopAlgFilterIPMQuasiNewton::outputIteration(int lsStatus, int lsNum, int u } } +#ifdef HIOP_USE_AXOM + //declaration required by C++14, not anymore by C++17 or after constexpr char hiopAlgFilterIPMQuasiNewton::default_state_filename[]; @@ -1492,14 +1505,14 @@ void hiopAlgFilterIPMQuasiNewton::save_state_to_file(const ::std::string& path_i { auto path = path_in=="" ? default_state_filename : path_in; - axom_sidre_DataStore* ds = new axom_sidre_DataStore(0); + sidre::DataStore* ds = new sidre::DataStore(); this->save_state_to_data_store(ds); - //::axom::sidre::IOManager writer(this->get_nlp()->get_comm()); - //int num_files; - //MPI_Comm_size(this->get_nlp()->get_comm(), &num_files); - //writer.write(ds_->getRoot(), num_files, path.str(), ::axom::sidre::Group::getDefaultIOProtocol()); + 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.c_str(), sidre::Group::getDefaultIOProtocol()); delete ds; } @@ -1510,17 +1523,28 @@ void hiopAlgFilterIPMQuasiNewton::load_state_from_file(const ::std::string& path //todo } -void hiopAlgFilterIPMQuasiNewton::save_state_to_data_store(void* ds) +void hiopAlgFilterIPMQuasiNewton::save_state_to_data_store(::axom::sidre::DataStore* ds) { - //Group* nlp_group = ds->getRoot()->createGroup("hiop solver"); + using IndType = sidre::IndexType; + sidre::Group* nlp_group = ds->getRoot()->createGroup("hiop solver"); //create views for each member that needs to be saved const double* x = it_curr->get_x()->local_data_host(); - //destination = nlp_group->createViewAndAllocate("x", ::axom::sidre::DOUBLE_ID, size); + const IndType size = it_curr->get_x()->get_local_size(); + sidre::View* dest = nlp_group->createViewAndAllocate("x", axom::sidre::DOUBLE_ID, size); + double *const dest_ptr(dest->getArray()); + const auto stride(dest->getStride()); + if(1==stride) { + ::std::copy(x, x+size, dest_ptr); + } else { + for(IndType i=0; igetRoot()->createGroup("hiop solver"); @@ -1529,7 +1553,7 @@ void hiopAlgFilterIPMQuasiNewton::load_state_from_data_store(const void* ds) const double* x = it_curr->get_x()->local_data_host(); //destination = nlp_group->createViewAndAllocate("x", ::axom::sidre::DOUBLE_ID, size); } - +#endif // HIOP_USE_AXOM /****************************************************************************************************** * FULL NEWTON IPM diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index 4be26885..a44208b6 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -67,12 +67,14 @@ #include "hiopPDPerturbation.hpp" #include "hiopFactAcceptor.hpp" +#ifdef HIOP_USE_AXOM +#include +#endif + #include "hiopTimer.hpp" namespace hiop { -//temporary dummy definition for axom::sidre::DataStore until axom is going to be included -using axom_sidre_DataStore=double; class hiopAlgFilterIPMBase { public: @@ -344,9 +346,11 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase virtual hiopSolveStatus run(); - //work in progress - virtual void save_state_to_data_store(void* sidre_data_store); - virtual void load_state_from_data_store(const void* sidre_data_store); + // note that checkpointing is only available with a axom-enabled build +#ifdef HIOP_USE_AXOM + + virtual void save_state_to_data_store(::axom::sidre::DataStore* data_store); + virtual void load_state_from_data_store(const ::axom::sidre::DataStore* data_store); static constexpr char default_state_filename[] = "hiop_qn_state.sidre"; @@ -372,7 +376,7 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase * */ void load_state_from_file(const ::std::string& path=""); - +#endif // HIOP_USE_AXOM private: virtual void outputIteration(int lsStatus, int lsNum, int use_soc = 0, int use_fr = 0); private: diff --git a/src/Utils/hiopOptions.cpp b/src/Utils/hiopOptions.cpp index c8529645..f2d1cbb7 100644 --- a/src/Utils/hiopOptions.cpp +++ b/src/Utils/hiopOptions.cpp @@ -1289,6 +1289,8 @@ void hiopOptionsNLP::register_options() ""); } + // checkpointing and restarting + } void hiopOptionsNLP::ensure_consistence() From a14a445007461f5186eefebe92b8c5fb56765f66 Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Wed, 28 Aug 2024 15:48:41 -0700 Subject: [PATCH 05/25] added user options for checkpointing --- src/Optimization/hiopAlgFilterIPM.cpp | 20 ++++++++++++++ src/Optimization/hiopAlgFilterIPM.hpp | 5 ++++ src/Utils/hiopOptions.cpp | 40 +++++++++++++++++++++++++-- 3 files changed, 63 insertions(+), 2 deletions(-) diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index a179341a..1f6b52a7 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -1106,6 +1106,11 @@ 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 ************************************************/ @@ -1553,6 +1558,21 @@ void hiopAlgFilterIPMQuasiNewton::load_state_from_data_store(const sidre::DataSt const double* x = it_curr->get_x()->local_data_host(); //destination = nlp_group->createViewAndAllocate("x", ::axom::sidre::DOUBLE_ID, size); } + +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 + ::std::string path = nlp->options->GetString("checkpoint_file"); + + + + // replace # + +} #endif // HIOP_USE_AXOM /****************************************************************************************************** diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index a44208b6..d155ffed 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -379,6 +379,11 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase #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 private: hiopNlpDenseConstraints* nlpdc; private: diff --git a/src/Utils/hiopOptions.cpp b/src/Utils/hiopOptions.cpp index f2d1cbb7..0f5d7eef 100644 --- a/src/Utils/hiopOptions.cpp +++ b/src/Utils/hiopOptions.cpp @@ -1290,9 +1290,25 @@ void hiopOptionsNLP::register_options() } // checkpointing and restarting - -} + // - currently only for IPM Quasi-Newton solver + // - only available with HIOP_USE_AXOM + { + vector range = {"yes", "no"}; + constexpr char msgcs[] = "Save state of NLP solver to file specified 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. If present character '#' will be replaced " + "with iteration number."; + 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'."; + register_str_option("checkpoint_load_on_start", range[1], range, msgclos); + } +} void hiopOptionsNLP::ensure_consistence() { //check that the values of different options are consistent @@ -1553,6 +1569,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 } ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// From d4979557717ca5736112bf491d9c2ccb9efa7bca Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Tue, 3 Sep 2024 14:13:56 -0700 Subject: [PATCH 06/25] more work on load checkpoint EOD --- src/Optimization/hiopAlgFilterIPM.cpp | 133 ++++++++++++++++++++------ src/Optimization/hiopAlgFilterIPM.hpp | 31 +++--- src/Utils/hiopOptions.cpp | 2 +- 3 files changed, 122 insertions(+), 44 deletions(-) diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index 1f6b52a7..7c514169 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -987,8 +987,29 @@ 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; + // + // starting point: + // - user provided (with slack adjustments and lsq eq. duals initialization + // or + // - loaded checkpoint + // + if(nlp->options->GetString("checkpoint_load_on_start") != "yes") { + //this also evaluates the nlp + startingProcedure(*it_curr, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); + _mu=mu0; + } else { + // + //checkpoint load + // + //load from file: will populate it_curr, _Hess_lagr, and algorithmic parameters + load_state_from_file(nlp->options->GetString("checkpoint_file")); + //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); @@ -1503,13 +1524,8 @@ void hiopAlgFilterIPMQuasiNewton::outputIteration(int lsStatus, int lsNum, int u #ifdef HIOP_USE_AXOM -//declaration required by C++14, not anymore by C++17 or after -constexpr char hiopAlgFilterIPMQuasiNewton::default_state_filename[]; - -void hiopAlgFilterIPMQuasiNewton::save_state_to_file(const ::std::string& path_in) +void hiopAlgFilterIPMQuasiNewton::save_state_to_file(const ::std::string& path) { - auto path = path_in=="" ? default_state_filename : path_in; - sidre::DataStore* ds = new sidre::DataStore(); this->save_state_to_data_store(ds); @@ -1522,41 +1538,90 @@ void hiopAlgFilterIPMQuasiNewton::save_state_to_file(const ::std::string& path_i delete ds; } -void hiopAlgFilterIPMQuasiNewton::load_state_from_file(const ::std::string& path_in) +void hiopAlgFilterIPMQuasiNewton::load_state_from_file(const ::std::string& path) { - auto path = path_in=="" ? default_state_filename : path_in; - //todo + sidre::DataStore* ds = new sidre::DataStore(); + sidre::IOManager reader(this->get_nlp()->get_comm()); + reader.read(ds->getRoot(), path, false); + + this->load_state_from_data_store(ds); + + delete ds; } -void hiopAlgFilterIPMQuasiNewton::save_state_to_data_store(::axom::sidre::DataStore* ds) +void hiopAlgFilterIPMQuasiNewton:: +copy_vec_to_new_view(const ::std::string& name, const hiopVector* vec, sidre::Group* nlp_group) { using IndType = sidre::IndexType; - sidre::Group* nlp_group = ds->getRoot()->createGroup("hiop solver"); + const IndType size = vec->get_local_size(); + sidre::View* dest = nlp_group->createViewAndAllocate(name, sidre::DOUBLE_ID, size); - //create views for each member that needs to be saved - - const double* x = it_curr->get_x()->local_data_host(); - const IndType size = it_curr->get_x()->get_local_size(); - sidre::View* dest = nlp_group->createViewAndAllocate("x", axom::sidre::DOUBLE_ID, size); - double *const dest_ptr(dest->getArray()); const auto stride(dest->getStride()); + double *const dest_ptr(dest->getArray()); + const double* arr = vec->local_data_host_const(); if(1==stride) { - ::std::copy(x, x+size, dest_ptr); + ::std::copy(arr, arr+size, dest_ptr); } else { for(IndType i=0; igetView(name); + if(view) { + + } + else { + nlp->log->printf(hovWarning, "Could not find view '%s' in checkpointing file.\n", name.c_str()); } } -void hiopAlgFilterIPMQuasiNewton::load_state_from_data_store(const sidre::DataStore* ds) +void hiopAlgFilterIPMQuasiNewton::save_state_to_data_store(::axom::sidre::DataStore* ds) { - //Group* nlp_group = ds->getRoot()->createGroup("hiop solver"); + using IndType = sidre::IndexType; + sidre::Group* nlp_group = ds->getRoot()->createGroup("hiop solver"); //create views for each member that needs to be saved + copy_vec_to_new_view("x", it_curr->get_x(), nlp_group); + copy_vec_to_new_view("d", it_curr->get_d(), nlp_group); + copy_vec_to_new_view("sxl", it_curr->get_sxl(), nlp_group); + copy_vec_to_new_view("sxu", it_curr->get_sxu(), nlp_group); + copy_vec_to_new_view("sdl", it_curr->get_sdl(), nlp_group); + copy_vec_to_new_view("sdu", it_curr->get_sdu(), nlp_group); + copy_vec_to_new_view("yc", it_curr->get_yc(), nlp_group); + copy_vec_to_new_view("zl", it_curr->get_zl(), nlp_group); + copy_vec_to_new_view("zu", it_curr->get_zu(), nlp_group); + copy_vec_to_new_view("vl", it_curr->get_vl(), nlp_group); + copy_vec_to_new_view("vu", it_curr->get_vu(), nlp_group); + + //quasi-Newton Hessian approximation + + //algorithmic parameters for this state + //mu, iteration number + const double alg_params[] = {_mu}; + +} + +void hiopAlgFilterIPMQuasiNewton::load_state_from_data_store(const sidre::DataStore* ds) +{ + const sidre::Group* nlp_group = ds->getRoot()->getGroup("hiop solver"); + + copy_vec_from_view("x", it_curr->get_x(), nlp_group); + copy_vec_from_view("d", it_curr->get_d(), nlp_group); + copy_vec_from_view("sxl", it_curr->get_sxl(), nlp_group); + copy_vec_from_view("sxu", it_curr->get_sxu(), nlp_group); + copy_vec_from_view("sdl", it_curr->get_sdl(), nlp_group); + copy_vec_from_view("sdu", it_curr->get_sdu(), nlp_group); + copy_vec_from_view("yc", it_curr->get_yc(), nlp_group); + copy_vec_from_view("zl", it_curr->get_zl(), nlp_group); + copy_vec_from_view("zu", it_curr->get_zu(), nlp_group); + copy_vec_from_view("vl", it_curr->get_vl(), nlp_group); + copy_vec_from_view("vu", it_curr->get_vu(), nlp_group); + - const double* x = it_curr->get_x()->local_data_host(); - //destination = nlp_group->createViewAndAllocate("x", ::axom::sidre::DOUBLE_ID, size); } void hiopAlgFilterIPMQuasiNewton::checkpointing_stuff() @@ -1566,12 +1631,20 @@ void hiopAlgFilterIPMQuasiNewton::checkpointing_stuff() } int chk_every_N = nlp->options->GetInteger("checkpoint_save_every_N_iter"); //check iteration - ::std::string path = nlp->options->GetString("checkpoint_file"); - + 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"); + auto pos = path.find("#"); + if(string::npos != pos) { + auto s_it_num = ::std::to_string(iter_num); + path.replace(pos, 1, s_it_num); + } - - // replace # - + 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 diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index d155ffed..7784ea85 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -352,30 +352,24 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase virtual void save_state_to_data_store(::axom::sidre::DataStore* data_store); virtual void load_state_from_data_store(const ::axom::sidre::DataStore* data_store); - static constexpr char default_state_filename[] = "hiop_qn_state.sidre"; - /** - * @brief save the state of the algorithm to the file + * @brief Save the state of the algorithm to the file * @param path the name of the file * * @details - * Internally, HiOp uses axom::sidre::DataStore, which is saved to the file. If argument is the - * empty string, HiOp will attempt saving the state to the path specified by default_state_filename - * static member. + * Internally, HiOp uses axom::sidre::DataStore and sidre's scalable IO. */ - void save_state_to_file(const ::std::string& path=""); + void save_state_to_file(const ::std::string& path); /** * @brief load the state of the algorithm from file * @param path the name of the file to load from * * @details - * The file should contains a axom::sidre::DataStore that was previously saved using save_state_to_file(). - * If argument is the empty string, HiOp will attempt loading state from the path specified by - * default_state_filename static member. - * + * The file should contains a axom::sidre::DataStore that was previously saved using + * save_state_to_file(). */ - void load_state_from_file(const ::std::string& path=""); + void load_state_from_file(const ::std::string& path); #endif // HIOP_USE_AXOM private: virtual void outputIteration(int lsStatus, int lsNum, int use_soc = 0, int use_fr = 0); @@ -383,7 +377,18 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase #ifdef HIOP_USE_AXOM ///@brief The options-based logic for saving checkpoint and the call to save_state(). void checkpointing_stuff(); -#endif + + /** + * @brief Copy HiOp vector to a (new) axom::sidre::View. + * + * @details A new view is created/allocated within the sidre view. Pointer is managed by the DataStore. + */ + void copy_vec_to_new_view(const ::std::string& name, const hiopVector* vec, ::axom::sidre::Group* nlp_group); + + /// Copy content of the named sidre view into HiOp Vector. + void copy_vec_from_view(const ::std::string& name, hiopVector* vec, const axom::sidre::Group* nlp_group); +#endif // HIOP_USE_AXOM + private: hiopNlpDenseConstraints* nlpdc; private: diff --git a/src/Utils/hiopOptions.cpp b/src/Utils/hiopOptions.cpp index 0f5d7eef..dde6d3a5 100644 --- a/src/Utils/hiopOptions.cpp +++ b/src/Utils/hiopOptions.cpp @@ -1301,7 +1301,7 @@ void hiopOptionsNLP::register_options() register_int_option("checkpoint_save_every_N_iter", 10, 1, 1e+6, msgcsN); constexpr char msgcf[] = "Path to checkpoint file. If present character '#' will be replaced " - "with iteration number."; + "with the iteration number."; register_str_option("checkpoint_file", "hiop_state_#.chk", msgcf); constexpr char msgclos[] = "On (re)start the NLP solver will load checkpoint file " From d4900a973ab7dd9542118ac4e9b8ee9fa6f8cae4 Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Wed, 4 Sep 2024 09:17:33 -0700 Subject: [PATCH 07/25] semi-operation checkpointing checkpoint interface complete additional states need to be saved --- src/Optimization/hiopAlgFilterIPM.cpp | 103 +++++++++++++++++++------- src/Optimization/hiopAlgFilterIPM.hpp | 34 +++++++-- 2 files changed, 105 insertions(+), 32 deletions(-) diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index 7c514169..d978066a 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -79,6 +79,8 @@ using namespace axom; #include #include #include +#include +#include namespace hiop { @@ -1002,12 +1004,20 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() //checkpoint load // //load from file: will populate it_curr, _Hess_lagr, and algorithmic parameters - load_state_from_file(nlp->options->GetString("checkpoint_file")); - //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; + auto chkpnt_ok = load_state_from_file(nlp->options->GetString("checkpoint_file")); + if(chkpnt_ok) { + //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; + } + } else { + nlp->log->printf(hovWarning, "Using default starting procedure (no checkpoint load!).\n"); + //this also evaluates the nlp + startingProcedure(*it_curr, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); + _mu=mu0; } + solver_status_ = NlpSolve_SolveNotCalled; } @@ -1524,29 +1534,45 @@ void hiopAlgFilterIPMQuasiNewton::outputIteration(int lsStatus, int lsNum, int u #ifdef HIOP_USE_AXOM -void hiopAlgFilterIPMQuasiNewton::save_state_to_file(const ::std::string& path) +bool hiopAlgFilterIPMQuasiNewton::save_state_to_file(const ::std::string& path) noexcept { - sidre::DataStore* ds = new sidre::DataStore(); - - this->save_state_to_data_store(ds); - - 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.c_str(), sidre::Group::getDefaultIOProtocol()); - + auto success = true; + sidre::DataStore* ds = nullptr; + try { + ds = new sidre::DataStore(); + this->save_state_to_data_store(ds); + + 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.c_str(), sidre::Group::getDefaultIOProtocol()); + } catch(const std::exception& exp) { + nlp->log->printf(hovError, "Error in saving checkpoint to file '%s'\n", path.c_str()); + nlp->log->printf(hovError, " Addtl info: %s\n", exp.what()); + success = false; + } delete ds; + return success; } -void hiopAlgFilterIPMQuasiNewton::load_state_from_file(const ::std::string& path) +bool hiopAlgFilterIPMQuasiNewton::load_state_from_file(const ::std::string& path) noexcept { - sidre::DataStore* ds = new sidre::DataStore(); - sidre::IOManager reader(this->get_nlp()->get_comm()); - reader.read(ds->getRoot(), path, false); - - this->load_state_from_data_store(ds); - + sidre::DataStore* ds = nullptr; + auto success = true; + try { + ds = new sidre::DataStore(); + sidre::IOManager reader(this->get_nlp()->get_comm()); + reader.read(ds->getRoot(), path, false); + + this->load_state_from_data_store(ds); + + } 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()); + success = false; + } delete ds; + return success; } void hiopAlgFilterIPMQuasiNewton:: @@ -1570,12 +1596,35 @@ copy_vec_to_new_view(const ::std::string& name, const hiopVector* vec, sidre::Gr void hiopAlgFilterIPMQuasiNewton:: copy_vec_from_view(const ::std::string& name, hiopVector* vec, const sidre::Group* nlp_group) { - const sidre::View* view = nlp_group->getView(name); + const sidre::View* view_const = nlp_group->getView(name); + if(!view_const) { + ::std::stringstream ss; + ss << "Could not find view '" << name << "in sidre::Group.\n"; + throw ::std::runtime_error(ss.str()); + } + // const_cast becase View does not have a const getArray() + sidre::View* view = const_cast(view_const); if(view) { - - } - else { - nlp->log->printf(hovWarning, "Could not find view '%s' in checkpointing file.\n", name.c_str()); + const hiop::size_type size = vec->get_local_size(); + if(view->getNumElements() != size) { + ::std::stringstream ss; + ss << "Size mismatch for state/view '" << name << "' between hiop and sidre::View: " << + " HiOp state is " << size << " doubles, while the view has " << view->getNumElements() << + " double elements.\n"; + throw ::std::runtime_error(ss.str()); + } + double* arr_dest = vec->local_data_host(); + 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 Date: Wed, 4 Sep 2024 09:19:48 -0700 Subject: [PATCH 08/25] removed save checkpoint callback from the interface --- src/Interface/hiopInterface.hpp | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/src/Interface/hiopInterface.hpp b/src/Interface/hiopInterface.hpp index 0d7b4ef2..47c54d3c 100644 --- a/src/Interface/hiopInterface.hpp +++ b/src/Interface/hiopInterface.hpp @@ -496,26 +496,6 @@ class hiopInterfaceBase { return true; } - - /** - * This method is called after each iteration and should be implemented by the user to instruct HiOp - * to save a complete state of the algorithm to disk. An axom::sidre::DataStore will be used for IO, which - * can be set by the user upon creation of the HiOp algorithm class. If not set, HiOp will create one - * internally. TODO: mention API methods from the Alg class. - * - * This feature is useful for IO checkpointing, for example; it allows the internal algorithm to be - * restarted at the same state as before saving. @see hiop::hiopAlgorithm::tbd for the method to be used - * to load a previously saved state. - * - * The method provided by the user should return true if HiOp should save the state of the algorithm at - * at current iteration. The argument is passed HiOp to indicate the current iterate number. - * - * @param[in] iter the current iteration number - */ - virtual bool save_state_iterate_callback(int iter) - { - return false; - } /** * A wildcard function used to change the primal variables. From d21377479c288af1b664a09ee43e1b39adfb7e9c Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Wed, 4 Sep 2024 09:26:34 -0700 Subject: [PATCH 09/25] fixed typos in comments --- src/Optimization/hiopAlgFilterIPM.cpp | 3 --- src/Optimization/hiopAlgFilterIPM.hpp | 10 +++++----- 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index d978066a..f9c45acd 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -84,9 +84,6 @@ using namespace axom; namespace hiop { -//#ifdef HIOP_USE_AXOM - -//#ifdef HIOP_USE_AXOM hiopAlgFilterIPMBase::hiopAlgFilterIPMBase(hiopNlpFormulation* nlp_in, const bool within_FR) : soc_dir(nullptr), diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index e0a4eb87..4437f231 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -353,9 +353,9 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase * @param data_store a pointer to DataStore * * @details - * A new sidre::group "hiop solver" is create within data_store. Then a sidre::View is - * created for each of the algorithm's states and the states are copied in the - * corresponding views. + * A new sidre::group "hiop solver" is created within data_store. Then a sidre::View + * is created within the group for each of the algorithm's states and the states are + * copied in the corresponding views. */ virtual void save_state_to_data_store(::axom::sidre::DataStore* data_store); @@ -405,11 +405,11 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase /** * @brief Copy HiOp vector to a (new) axom::sidre::View. * - * @details A new view is created/allocated within the sidre view. Pointer is managed by the DataStore. + * @details A new view is created/allocated within the sidre group. */ void copy_vec_to_new_view(const ::std::string& name, const hiopVector* vec, ::axom::sidre::Group* nlp_group); - /// Copy content of the named sidre view into HiOp Vector. + /// Copy content of the named sidre view from nlp_group into HiOp Vector. void copy_vec_from_view(const ::std::string& name, hiopVector* vec, const axom::sidre::Group* nlp_group); #endif // HIOP_USE_AXOM From 4a9fbf1fc6eb1b07d2a81f118909d6e3e589bf7b Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Fri, 6 Sep 2024 17:03:19 -0700 Subject: [PATCH 10/25] moved sidre-related code from Algorithm class to a "utils" helper --- src/Optimization/hiopAlgFilterIPM.cpp | 134 ++++++----------- src/Utils/SidreHelper.hpp | 206 ++++++++++++++++++++++++++ 2 files changed, 252 insertions(+), 88 deletions(-) create mode 100644 src/Utils/SidreHelper.hpp diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index f9c45acd..d0c72b9d 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -67,10 +67,7 @@ #include "hiopCppStdUtils.hpp" #ifdef HIOP_USE_AXOM -#include -#include -#include -#include +#include "SidreHelper.hpp" using namespace axom; #endif @@ -79,8 +76,6 @@ using namespace axom; #include #include #include -#include -#include namespace hiop { @@ -988,7 +983,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() // // starting point: - // - user provided (with slack adjustments and lsq eq. duals initialization + // - user provided (with slack adjustments and lsq eq. duals initialization) // or // - loaded checkpoint // @@ -1014,7 +1009,6 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() startingProcedure(*it_curr, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); _mu=mu0; } - solver_status_ = NlpSolve_SolveNotCalled; } @@ -1572,102 +1566,66 @@ bool hiopAlgFilterIPMQuasiNewton::load_state_from_file(const ::std::string& path return success; } -void hiopAlgFilterIPMQuasiNewton:: -copy_vec_to_new_view(const ::std::string& name, const hiopVector* vec, sidre::Group* nlp_group) -{ - using IndType = sidre::IndexType; - const IndType size = vec->get_local_size(); - sidre::View* dest = nlp_group->createViewAndAllocate(name, sidre::DOUBLE_ID, size); - - const auto stride(dest->getStride()); - double *const dest_ptr(dest->getArray()); - const double* arr = vec->local_data_host_const(); - if(1==stride) { - ::std::copy(arr, arr+size, dest_ptr); - } else { - for(IndType i=0; igetView(name); - if(!view_const) { - ::std::stringstream ss; - ss << "Could not find view '" << name << "in sidre::Group.\n"; - throw ::std::runtime_error(ss.str()); - } - // const_cast becase View does not have a const getArray() - sidre::View* view = const_cast(view_const); - if(view) { - const hiop::size_type size = vec->get_local_size(); - if(view->getNumElements() != size) { - ::std::stringstream ss; - ss << "Size mismatch for state/view '" << name << "' between hiop and sidre::View: " << - " HiOp state is " << size << " doubles, while the view has " << view->getNumElements() << - " double elements.\n"; - throw ::std::runtime_error(ss.str()); - } - double* arr_dest = vec->local_data_host(); - 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; igetRoot()->createGroup("hiop solver"); - - //create views for each member that needs to be saved - copy_vec_to_new_view("x", it_curr->get_x(), nlp_group); - copy_vec_to_new_view("d", it_curr->get_d(), nlp_group); - copy_vec_to_new_view("sxl", it_curr->get_sxl(), nlp_group); - copy_vec_to_new_view("sxu", it_curr->get_sxu(), nlp_group); - copy_vec_to_new_view("sdl", it_curr->get_sdl(), nlp_group); - copy_vec_to_new_view("sdu", it_curr->get_sdu(), nlp_group); - copy_vec_to_new_view("yc", it_curr->get_yc(), nlp_group); - copy_vec_to_new_view("zl", it_curr->get_zl(), nlp_group); - copy_vec_to_new_view("zu", it_curr->get_zu(), nlp_group); - copy_vec_to_new_view("vl", it_curr->get_vl(), nlp_group); - copy_vec_to_new_view("vu", it_curr->get_vu(), nlp_group); + //metadata + + //iterate states + //create views for each member that needs to be saved + SidreHelper::copy_vec_to_view(*nlp_group, "x", *it_curr->get_x()); + SidreHelper::copy_vec_to_view(*nlp_group, "d", *it_curr->get_d()); + SidreHelper::copy_vec_to_view(*nlp_group, "sxl", *it_curr->get_sxl()); + SidreHelper::copy_vec_to_view(*nlp_group, "sxu", *it_curr->get_sxu()); + SidreHelper::copy_vec_to_view(*nlp_group, "sdl", *it_curr->get_sdl()); + SidreHelper::copy_vec_to_view(*nlp_group, "sdu", *it_curr->get_sdu()); + SidreHelper::copy_vec_to_view(*nlp_group, "yc", *it_curr->get_yc()); + SidreHelper::copy_vec_to_view(*nlp_group, "zl", *it_curr->get_zl()); + SidreHelper::copy_vec_to_view(*nlp_group, "zu", *it_curr->get_zu()); + SidreHelper::copy_vec_to_view(*nlp_group, "vl", *it_curr->get_vl()); + SidreHelper::copy_vec_to_view(*nlp_group, "vu", *it_curr->get_vu()); + //quasi-Newton Hessian approximation //algorithmic parameters for this state //mu, iteration number - const double alg_params[] = {_mu}; + const double alg_params[] = {_mu, (double)iter_num}; + const size_type nparams = sizeof(alg_params) / sizeof(double); + SidreHelper::copy_array_to_view(*nlp_group, "alg_params", alg_params, nparams); } void hiopAlgFilterIPMQuasiNewton::load_state_from_data_store(const sidre::DataStore* ds) { const sidre::Group* nlp_group = ds->getRoot()->getGroup("hiop solver"); - copy_vec_from_view("x", it_curr->get_x(), nlp_group); - copy_vec_from_view("d", it_curr->get_d(), nlp_group); - copy_vec_from_view("sxl", it_curr->get_sxl(), nlp_group); - copy_vec_from_view("sxu", it_curr->get_sxu(), nlp_group); - copy_vec_from_view("sdl", it_curr->get_sdl(), nlp_group); - copy_vec_from_view("sdu", it_curr->get_sdu(), nlp_group); - copy_vec_from_view("yc", it_curr->get_yc(), nlp_group); - copy_vec_from_view("zl", it_curr->get_zl(), nlp_group); - copy_vec_from_view("zu", it_curr->get_zu(), nlp_group); - copy_vec_from_view("vl", it_curr->get_vl(), nlp_group); - copy_vec_from_view("vu", it_curr->get_vu(), nlp_group); - - + //metadata + + //iterate states + SidreHelper::copy_vec_from_view(*nlp_group, "x", *it_curr->get_x()); + SidreHelper::copy_vec_from_view(*nlp_group, "d", *it_curr->get_d()); + SidreHelper::copy_vec_from_view(*nlp_group, "sxl", *it_curr->get_sxl()); + SidreHelper::copy_vec_from_view(*nlp_group, "sxu", *it_curr->get_sxu()); + SidreHelper::copy_vec_from_view(*nlp_group, "sdl", *it_curr->get_sdl()); + SidreHelper::copy_vec_from_view(*nlp_group, "sdu", *it_curr->get_sdu()); + SidreHelper::copy_vec_from_view(*nlp_group, "yc", *it_curr->get_yc()); + SidreHelper::copy_vec_from_view(*nlp_group, "zl", *it_curr->get_zl()); + SidreHelper::copy_vec_from_view(*nlp_group, "zu", *it_curr->get_zu()); + SidreHelper::copy_vec_from_view(*nlp_group, "vl", *it_curr->get_vl()); + SidreHelper::copy_vec_from_view(*nlp_group, "vu", *it_curr->get_vu()); + + //quasi-Newton Hessian approximation + + //algorithmic parameters + //!!! dev note: nparams needs to match the nparams from save_state_to_data_store + const int nparams = 2; + double alg_params[nparams]; + SidreHelper::copy_array_from_view(*nlp_group, "alg_params", alg_params, nparams); + //!!! dev note: match order in save_state_to_data_store + _mu = alg_params[0]; + iter_num = alg_params[1]; } void hiopAlgFilterIPMQuasiNewton::checkpointing_stuff() diff --git a/src/Utils/SidreHelper.hpp b/src/Utils/SidreHelper.hpp new file mode 100644 index 00000000..1c1f37c7 --- /dev/null +++ b/src/Utils/SidreHelper.hpp @@ -0,0 +1,206 @@ +// 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 Date: Sun, 8 Sep 2024 16:21:29 -0700 Subject: [PATCH 11/25] switched to refs; some testing of options-based checkpointing --- src/Optimization/hiopAlgFilterIPM.cpp | 133 +++++++++++++++----------- src/Optimization/hiopAlgFilterIPM.hpp | 65 +++++++------ src/Utils/SidreHelper.hpp | 13 ++- 3 files changed, 125 insertions(+), 86 deletions(-) diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index d0c72b9d..edb43da0 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -991,6 +991,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() //this also evaluates the nlp startingProcedure(*it_curr, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); _mu=mu0; + iter_num = 0; } else { // //checkpoint load @@ -1007,7 +1008,8 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() nlp->log->printf(hovWarning, "Using default starting procedure (no checkpoint load!).\n"); //this also evaluates the nlp startingProcedure(*it_curr, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); - _mu=mu0; + _mu=mu0; + iter_num = 0; } solver_status_ = NlpSolve_SolveNotCalled; } @@ -1020,7 +1022,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()); @@ -1527,105 +1529,120 @@ void hiopAlgFilterIPMQuasiNewton::outputIteration(int lsStatus, int lsNum, int u bool hiopAlgFilterIPMQuasiNewton::save_state_to_file(const ::std::string& path) noexcept { - auto success = true; - sidre::DataStore* ds = nullptr; try { - ds = new sidre::DataStore(); - this->save_state_to_data_store(ds); + 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.c_str(), sidre::Group::getDefaultIOProtocol()); + writer.write(ds.getRoot(), n_files, path, sidre::Group::getDefaultIOProtocol()); + return true; } catch(const std::exception& exp) { - nlp->log->printf(hovError, "Error in saving checkpoint to file '%s'\n", path.c_str()); + 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()); - success = false; + return false; } - delete ds; - return success; } bool hiopAlgFilterIPMQuasiNewton::load_state_from_file(const ::std::string& path) noexcept { - sidre::DataStore* ds = nullptr; - auto success = true; try { - ds = new sidre::DataStore(); + sidre::DataStore ds; + sidre::IOManager reader(this->get_nlp()->get_comm()); - reader.read(ds->getRoot(), path, false); - - this->load_state_from_data_store(ds); - + 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()); - success = false; + return false; } - delete ds; - return success; } -void hiopAlgFilterIPMQuasiNewton::save_state_to_data_store(::axom::sidre::DataStore* ds) +void hiopAlgFilterIPMQuasiNewton::save_state_to_sidre_group(::axom::sidre::Group& group) { using IndType = sidre::IndexType; - sidre::Group* nlp_group = ds->getRoot()->createGroup("hiop solver"); //metadata //iterate states //create views for each member that needs to be saved - SidreHelper::copy_vec_to_view(*nlp_group, "x", *it_curr->get_x()); - SidreHelper::copy_vec_to_view(*nlp_group, "d", *it_curr->get_d()); - SidreHelper::copy_vec_to_view(*nlp_group, "sxl", *it_curr->get_sxl()); - SidreHelper::copy_vec_to_view(*nlp_group, "sxu", *it_curr->get_sxu()); - SidreHelper::copy_vec_to_view(*nlp_group, "sdl", *it_curr->get_sdl()); - SidreHelper::copy_vec_to_view(*nlp_group, "sdu", *it_curr->get_sdu()); - SidreHelper::copy_vec_to_view(*nlp_group, "yc", *it_curr->get_yc()); - SidreHelper::copy_vec_to_view(*nlp_group, "zl", *it_curr->get_zl()); - SidreHelper::copy_vec_to_view(*nlp_group, "zu", *it_curr->get_zu()); - SidreHelper::copy_vec_to_view(*nlp_group, "vl", *it_curr->get_vl()); - SidreHelper::copy_vec_to_view(*nlp_group, "vu", *it_curr->get_vu()); + SidreHelper::copy_vec_to_view(group, "x", *it_curr->get_x()); + SidreHelper::copy_vec_to_view(group, "d", *it_curr->get_d()); + SidreHelper::copy_vec_to_view(group, "sxl", *it_curr->get_sxl()); + SidreHelper::copy_vec_to_view(group, "sxu", *it_curr->get_sxu()); + SidreHelper::copy_vec_to_view(group, "sdl", *it_curr->get_sdl()); + SidreHelper::copy_vec_to_view(group, "sdu", *it_curr->get_sdu()); + SidreHelper::copy_vec_to_view(group, "yc", *it_curr->get_yc()); + SidreHelper::copy_vec_to_view(group, "yd", *it_curr->get_yd()); + SidreHelper::copy_vec_to_view(group, "zl", *it_curr->get_zl()); + SidreHelper::copy_vec_to_view(group, "zu", *it_curr->get_zu()); + SidreHelper::copy_vec_to_view(group, "vl", *it_curr->get_vl()); + SidreHelper::copy_vec_to_view(group, "vu", *it_curr->get_vu()); - //quasi-Newton Hessian approximation + //state of quasi-Newton Hessian approximation //algorithmic parameters for this state - //mu, iteration number - const double alg_params[] = {_mu, (double)iter_num}; + //mu, iteration number, num MPI ranks + int nranks=1; +#ifdef HIOP_USE_MPI + MPI_Comm_size(get_nlp()->get_comm(), &nranks); +#endif + + const double alg_params[] = {_mu, (double)iter_num, (double)nranks}; const size_type nparams = sizeof(alg_params) / sizeof(double); - SidreHelper::copy_array_to_view(*nlp_group, "alg_params", alg_params, nparams); + SidreHelper::copy_array_to_view(group, "alg_params", alg_params, nparams); } -void hiopAlgFilterIPMQuasiNewton::load_state_from_data_store(const sidre::DataStore* ds) +void hiopAlgFilterIPMQuasiNewton::load_state_from_sidre_group(const sidre::Group& group) { - const sidre::Group* nlp_group = ds->getRoot()->getGroup("hiop solver"); - //metadata - //iterate states - SidreHelper::copy_vec_from_view(*nlp_group, "x", *it_curr->get_x()); - SidreHelper::copy_vec_from_view(*nlp_group, "d", *it_curr->get_d()); - SidreHelper::copy_vec_from_view(*nlp_group, "sxl", *it_curr->get_sxl()); - SidreHelper::copy_vec_from_view(*nlp_group, "sxu", *it_curr->get_sxu()); - SidreHelper::copy_vec_from_view(*nlp_group, "sdl", *it_curr->get_sdl()); - SidreHelper::copy_vec_from_view(*nlp_group, "sdu", *it_curr->get_sdu()); - SidreHelper::copy_vec_from_view(*nlp_group, "yc", *it_curr->get_yc()); - SidreHelper::copy_vec_from_view(*nlp_group, "zl", *it_curr->get_zl()); - SidreHelper::copy_vec_from_view(*nlp_group, "zu", *it_curr->get_zu()); - SidreHelper::copy_vec_from_view(*nlp_group, "vl", *it_curr->get_vl()); - SidreHelper::copy_vec_from_view(*nlp_group, "vu", *it_curr->get_vu()); - - //quasi-Newton Hessian approximation - //algorithmic parameters //!!! dev note: nparams needs to match the nparams from save_state_to_data_store - const int nparams = 2; + const int nparams = 3; double alg_params[nparams]; - SidreHelper::copy_array_from_view(*nlp_group, "alg_params", 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 = 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()); + } + + //iterate states + SidreHelper::copy_vec_from_view(group, "x", *it_curr->get_x()); + SidreHelper::copy_vec_from_view(group, "d", *it_curr->get_d()); + SidreHelper::copy_vec_from_view(group, "sxl", *it_curr->get_sxl()); + SidreHelper::copy_vec_from_view(group, "sxu", *it_curr->get_sxu()); + SidreHelper::copy_vec_from_view(group, "sdl", *it_curr->get_sdl()); + SidreHelper::copy_vec_from_view(group, "sdu", *it_curr->get_sdu()); + SidreHelper::copy_vec_from_view(group, "yc", *it_curr->get_yc()); + SidreHelper::copy_vec_from_view(group, "yd", *it_curr->get_yd()); + SidreHelper::copy_vec_from_view(group, "zl", *it_curr->get_zl()); + SidreHelper::copy_vec_from_view(group, "zu", *it_curr->get_zu()); + SidreHelper::copy_vec_from_view(group, "vl", *it_curr->get_vl()); + SidreHelper::copy_vec_from_view(group, "vu", *it_curr->get_vu()); + + //state of quasi-Newton Hessian approximation + } void hiopAlgFilterIPMQuasiNewton::checkpointing_stuff() diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index 4437f231..fe150bc5 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -68,7 +68,11 @@ #include "hiopFactAcceptor.hpp" #ifdef HIOP_USE_AXOM -#include +namespace axom { +namespace sidre { +class Group; // forward declaration +} +} #endif #include "hiopTimer.hpp" @@ -349,31 +353,47 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase // note that checkpointing is only available with a axom-enabled build #ifdef HIOP_USE_AXOM /** - * @brief Save state of HiOp algorithm to a sidre data store. - * @param data_store a pointer to DataStore + * @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 - * A new sidre::group "hiop solver" is created within data_store. Then a sidre::View - * is created within the group for each of the algorithm's states and the states are - * copied in the corresponding views. + * 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_data_store(::axom::sidre::DataStore* data_store); + virtual void save_state_to_sidre_group(::axom::sidre::Group& group); /** - * @brief Load state of HiOp algorithm from a sidre data store. - * @param data_store a pointer to DataStore + * @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 data store sidre::group (named "hiop solver") to HiOp algorithm's - * state variables. The group should be created by first calling save_state_to_data_store - * for a problem/NLP of the same sizes as the problem for which load_state_from_data_store - * is called. This ensures the views have the names and sizes expected by this method. - * Otherwise a std::runtime_error exception is thrown. + * 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_data_store(const ::axom::sidre::DataStore* data_store); + virtual void load_state_from_sidre_group(const ::axom::sidre::Group& group); /** - * @brief Save the state of the algorithm to the file + * @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 * @@ -384,7 +404,8 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase bool save_state_to_file(const ::std::string& path) noexcept; /** - * @brief load the state of the algorithm from file + * @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 * @@ -401,16 +422,6 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase #ifdef HIOP_USE_AXOM ///@brief The options-based logic for saving checkpoint and the call to save_state(). void checkpointing_stuff(); - - /** - * @brief Copy HiOp vector to a (new) axom::sidre::View. - * - * @details A new view is created/allocated within the sidre group. - */ - void copy_vec_to_new_view(const ::std::string& name, const hiopVector* vec, ::axom::sidre::Group* nlp_group); - - /// Copy content of the named sidre view from nlp_group into HiOp Vector. - void copy_vec_from_view(const ::std::string& name, hiopVector* vec, const axom::sidre::Group* nlp_group); #endif // HIOP_USE_AXOM private: diff --git a/src/Utils/SidreHelper.hpp b/src/Utils/SidreHelper.hpp index 1c1f37c7..428b0257 100644 --- a/src/Utils/SidreHelper.hpp +++ b/src/Utils/SidreHelper.hpp @@ -101,6 +101,7 @@ class SidreHelper "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) { @@ -153,7 +154,7 @@ class SidreHelper ::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"; + "while the view has " << view_const->getNumElements() << " double elements.\n"; throw ::std::runtime_error(ss.str()); } @@ -182,6 +183,16 @@ class SidreHelper copy_array_from_view(group, view_name, arr, size); } + + /// Add '.root' extension if path is not a valid file + static ::std::string check_path(::std::string path) + { + ::std::ifstream f(path, ::std::ifstream::in); + //this hack is to trigger a failure (f.good() returns false) if 'path' exists but it + //is a directory. + f.seekg(0, ::std::ios::end); + return f.good() ? path : (path + ".root"); + } private: /** * @brief Get or create new view within a sidre::Group From 085eb8849fb848e838f2f04552e95b97c728685f Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Wed, 11 Sep 2024 14:09:14 -0700 Subject: [PATCH 12/25] added sidre copy to/from dense matrices --- src/Optimization/hiopAlgFilterIPM.cpp | 127 +++++++++++++++++++------- src/Optimization/hiopAlgFilterIPM.hpp | 4 +- src/Utils/SidreHelper.hpp | 70 +++++++++++++- 3 files changed, 167 insertions(+), 34 deletions(-) diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index edb43da0..6ed3e69f 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -1192,7 +1192,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.", @@ -1572,22 +1572,53 @@ void hiopAlgFilterIPMQuasiNewton::save_state_to_sidre_group(::axom::sidre::Group //metadata - //iterate states + //iterate state //create views for each member that needs to be saved - SidreHelper::copy_vec_to_view(group, "x", *it_curr->get_x()); - SidreHelper::copy_vec_to_view(group, "d", *it_curr->get_d()); - SidreHelper::copy_vec_to_view(group, "sxl", *it_curr->get_sxl()); - SidreHelper::copy_vec_to_view(group, "sxu", *it_curr->get_sxu()); - SidreHelper::copy_vec_to_view(group, "sdl", *it_curr->get_sdl()); - SidreHelper::copy_vec_to_view(group, "sdu", *it_curr->get_sdu()); - SidreHelper::copy_vec_to_view(group, "yc", *it_curr->get_yc()); - SidreHelper::copy_vec_to_view(group, "yd", *it_curr->get_yd()); - SidreHelper::copy_vec_to_view(group, "zl", *it_curr->get_zl()); - SidreHelper::copy_vec_to_view(group, "zu", *it_curr->get_zu()); - SidreHelper::copy_vec_to_view(group, "vl", *it_curr->get_vl()); - SidreHelper::copy_vec_to_view(group, "vu", *it_curr->get_vu()); + 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()); - //state of quasi-Newton Hessian approximation //algorithmic parameters for this state //mu, iteration number, num MPI ranks @@ -1595,10 +1626,8 @@ void hiopAlgFilterIPMQuasiNewton::save_state_to_sidre_group(::axom::sidre::Group #ifdef HIOP_USE_MPI MPI_Comm_size(get_nlp()->get_comm(), &nranks); #endif - const double alg_params[] = {_mu, (double)iter_num, (double)nranks}; const size_type nparams = sizeof(alg_params) / sizeof(double); - SidreHelper::copy_array_to_view(group, "alg_params", alg_params, nparams); } @@ -1607,7 +1636,7 @@ void hiopAlgFilterIPMQuasiNewton::load_state_from_sidre_group(const sidre::Group //metadata //algorithmic parameters - //!!! dev note: nparams needs to match the nparams from save_state_to_data_store + //!!!note: nparams needs to match the nparams from save_state_to_data_store const int nparams = 3; double alg_params[nparams]; SidreHelper::copy_array_from_view(group, "alg_params", alg_params, nparams); @@ -1628,21 +1657,57 @@ void hiopAlgFilterIPMQuasiNewton::load_state_from_sidre_group(const sidre::Group } //iterate states - SidreHelper::copy_vec_from_view(group, "x", *it_curr->get_x()); - SidreHelper::copy_vec_from_view(group, "d", *it_curr->get_d()); - SidreHelper::copy_vec_from_view(group, "sxl", *it_curr->get_sxl()); - SidreHelper::copy_vec_from_view(group, "sxu", *it_curr->get_sxu()); - SidreHelper::copy_vec_from_view(group, "sdl", *it_curr->get_sdl()); - SidreHelper::copy_vec_from_view(group, "sdu", *it_curr->get_sdu()); - SidreHelper::copy_vec_from_view(group, "yc", *it_curr->get_yc()); - SidreHelper::copy_vec_from_view(group, "yd", *it_curr->get_yd()); - SidreHelper::copy_vec_from_view(group, "zl", *it_curr->get_zl()); - SidreHelper::copy_vec_from_view(group, "zu", *it_curr->get_zu()); - SidreHelper::copy_vec_from_view(group, "vl", *it_curr->get_vl()); - SidreHelper::copy_vec_from_view(group, "vu", *it_curr->get_vu()); + 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() diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index fe150bc5..72d11d5a 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -394,7 +394,7 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase /** * @brief Save the state of the algorithm to the file for checkpointing. * - * @param path the name of the file + * @param path the name of the file * @return true if successful, false otherwise * * @details @@ -406,7 +406,7 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase /** * @brief Load the state of the algorithm from checkpoint file. * - * @param path the name of the file to load from + * @param path the name of the file to load from * @return true if successful, false otherwise * * @details diff --git a/src/Utils/SidreHelper.hpp b/src/Utils/SidreHelper.hpp index 428b0257..69bca66b 100644 --- a/src/Utils/SidreHelper.hpp +++ b/src/Utils/SidreHelper.hpp @@ -122,11 +122,45 @@ class SidreHelper const double* arr = vec.local_data_host_const(); copy_array_to_view(group, view_name, arr, size); } + + /** + * @brief Copy vectors of HiOp iterate to multiple sidre::View(s) in specified sidre::Group + * + * @params group contains the views where the vectors members of iterate will be copied to. + * @params view_name_prefix is the string prefixing the views' names for each vector + * @params it is the HiOp iterate object + * + * @exception std::runtime indicates the group contains a view with a number of elements + * different than size apparent from the iterate's vector. + * + * @details For each vector member of the iterate, a view will be created if does not + * already exist. If exists, the view should have the same number of elements as the + * corresponding (HiOp) vector member of iterate passed as argument. The name of view + * is formed by appending the iterate's name ("x", "s", "y", etc.) to view_name_prefix. + */ + static void copy_iterate_to_views(::axom::sidre::Group& group, + const ::std::string& view_name_prefix, + const hiopIterate& it) + { + copy_vec_to_view(group, view_name_prefix+"x", *it.get_x()); + copy_vec_to_view(group, view_name_prefix+"d", *it.get_d()); + copy_vec_to_view(group, view_name_prefix+"sxl", *it.get_sxl()); + copy_vec_to_view(group, view_name_prefix+"sxu", *it.get_sxu()); + copy_vec_to_view(group, view_name_prefix+"sdl", *it.get_sdl()); + copy_vec_to_view(group, view_name_prefix+"sdu", *it.get_sdu()); + copy_vec_to_view(group, view_name_prefix+"yc", *it.get_yc()); + copy_vec_to_view(group, view_name_prefix+"yd", *it.get_yd()); + copy_vec_to_view(group, view_name_prefix+"zl", *it.get_zl()); + copy_vec_to_view(group, view_name_prefix+"zu", *it.get_zu()); + copy_vec_to_view(group, view_name_prefix+"vl", *it.get_vl()); + copy_vec_to_view(group, view_name_prefix+"vu", *it.get_vu()); + } + /** * @brief Copy raw array from sidre::View within specified sidre::Group. * - * @params group contains the view where the copy should be made to. + * @params group contains the view that should be copied from * @params view_name is the name of the view where to copy * @params arr_dest is the source double array * @params size is the number of elements of the array @@ -184,6 +218,40 @@ class SidreHelper } + /** + * @brief Copy iterate from multiple sidre::View(s) within specified sidre::Group. + * + * @params group contains the views where the copy should be made to. + * @params view_name_prefix is a string prefixing the views' name + * @params it is the HiOp iterate object where the views will be copied to. + * + * @exception std::runtime indicates the group contains a view with a number of elements + * different than expected size or that a view with the expected name does not exist. + * + * @details All views must exist and have a number of elements identical to the size of + * corresponding vector from `it`. The views are located by their names, which are + * expected to be view_name_prefix concatenated with the iterate's name, see method + * copy_vec_to_views. + */ + + static void copy_iterate_from_views(const ::axom::sidre::Group& group, + const ::std::string& view_name_prefix, + hiopIterate& it) + { + copy_vec_from_view(group, view_name_prefix+"x", *it.get_x()); + copy_vec_from_view(group, view_name_prefix+"d", *it.get_d()); + copy_vec_from_view(group, view_name_prefix+"sxl", *it.get_sxl()); + copy_vec_from_view(group, view_name_prefix+"sxu", *it.get_sxu()); + copy_vec_from_view(group, view_name_prefix+"sdl", *it.get_sdl()); + copy_vec_from_view(group, view_name_prefix+"sdu", *it.get_sdu()); + copy_vec_from_view(group, view_name_prefix+"yc", *it.get_yc()); + copy_vec_from_view(group, view_name_prefix+"yd", *it.get_yd()); + copy_vec_from_view(group, view_name_prefix+"zl", *it.get_zl()); + copy_vec_from_view(group, view_name_prefix+"zu", *it.get_zu()); + copy_vec_from_view(group, view_name_prefix+"vl", *it.get_vl()); + copy_vec_from_view(group, view_name_prefix+"vu", *it.get_vu()); + } + /// Add '.root' extension if path is not a valid file static ::std::string check_path(::std::string path) { From b21c3c59cada2abae9510380f165b34f078eb6ae Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Wed, 11 Sep 2024 14:09:49 -0700 Subject: [PATCH 13/25] instrumentation for saving quasi-Newton internals to sidre --- src/Optimization/hiopHessianLowRank.cpp | 1004 ++++------------------- src/Optimization/hiopHessianLowRank.hpp | 132 +-- 2 files changed, 170 insertions(+), 966 deletions(-) diff --git a/src/Optimization/hiopHessianLowRank.cpp b/src/Optimization/hiopHessianLowRank.cpp index a18af6d9..1612f2fd 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,16 +399,26 @@ 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); @@ -394,7 +431,7 @@ void hiopHessianLowRank::updateInternalBFGSRepresentation() 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 @@ -409,7 +446,7 @@ 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 @@ -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 b487ea08..05457c1e 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 From bfe1b407af32c15d7e42ad8d22ba05bf8911e7d6 Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Thu, 12 Sep 2024 06:39:33 -0700 Subject: [PATCH 14/25] updated iteration counter to keep track of total number over restarts --- src/Optimization/hiopAlgFilterIPM.cpp | 110 +++++++++++++++----------- src/Optimization/hiopAlgFilterIPM.hpp | 6 +- 2 files changed, 67 insertions(+), 49 deletions(-) diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index 6ed3e69f..6e32fa6b 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -87,7 +87,9 @@ hiopAlgFilterIPMBase::hiopAlgFilterIPMBase(hiopNlpFormulation* nlp_in, const boo d_soc(nullptr), within_FR_{within_FR}, pd_perturb_{nullptr}, - fact_acceptor_{nullptr} + fact_acceptor_{nullptr}, + iter_num_(0), + iter_num_total_(0) { nlp = nlp_in; //force completion of the nlp's initialization @@ -981,6 +983,9 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() nlp->runStats.tmOptimizTotal.start(); + iter_num_ = 0; + iter_num_total_ = 0; + // // starting point: // - user provided (with slack adjustments and lsq eq. duals initialization) @@ -991,7 +996,6 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() //this also evaluates the nlp startingProcedure(*it_curr, _f_nlp, *_c, *_d, *_grad_f, *_Jac_c, *_Jac_d); _mu=mu0; - iter_num = 0; } else { // //checkpoint load @@ -1006,10 +1010,10 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() } } else { nlp->log->printf(hovWarning, "Using default starting procedure (no checkpoint load!).\n"); - //this also evaluates the nlp + 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 = 0; } solver_status_ = NlpSolve_SolveNotCalled; } @@ -1022,7 +1026,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() nlp->log->write("First residual-------------", *resid, hovIteration); - 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()); @@ -1109,7 +1113,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(), @@ -1138,7 +1142,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::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. @@ -1152,7 +1156,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); @@ -1183,7 +1187,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 ***************************************************/ @@ -1204,7 +1208,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; } @@ -1212,7 +1216,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; @@ -1224,7 +1228,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 @@ -1365,7 +1369,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 @@ -1425,10 +1429,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) { @@ -1460,7 +1468,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(); //----- @@ -1471,7 +1479,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); } @@ -1497,12 +1505,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]; @@ -1520,7 +1528,7 @@ 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); } } @@ -1626,7 +1634,7 @@ void hiopAlgFilterIPMQuasiNewton::save_state_to_sidre_group(::axom::sidre::Group #ifdef HIOP_USE_MPI MPI_Comm_size(get_nlp()->get_comm(), &nranks); #endif - const double alg_params[] = {_mu, (double)iter_num, (double)nranks}; + const double alg_params[] = {_mu, (double)iter_num_total_, (double)nranks}; const size_type nparams = sizeof(alg_params) / sizeof(double); SidreHelper::copy_array_to_view(group, "alg_params", alg_params, nparams); } @@ -1642,7 +1650,7 @@ void hiopAlgFilterIPMQuasiNewton::load_state_from_sidre_group(const sidre::Group 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 = alg_params[1]; + iter_num_total_ = alg_params[1]; int nranks=1; #ifdef HIOP_USE_MPI @@ -1717,17 +1725,17 @@ void hiopAlgFilterIPMQuasiNewton::checkpointing_stuff() } int chk_every_N = nlp->options->GetInteger("checkpoint_save_every_N_iter"); //check iteration - if(iter_num>0 && iter_num % chk_every_N==0) { + 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"); auto pos = path.find("#"); if(string::npos != pos) { - auto s_it_num = ::std::to_string(iter_num); + auto s_it_num = ::std::to_string(iter_num_); path.replace(pos, 1, s_it_num); } - nlp->log->printf(hovSummary, "Saving checkpoint at iter %d in '%s'.\n", iter_num, path.c_str()); + 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); } @@ -2098,7 +2106,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()); @@ -2200,7 +2210,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(), @@ -2223,7 +2234,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. @@ -2237,7 +2248,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); @@ -2271,7 +2282,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 ***************************************************/ @@ -2303,7 +2314,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, @@ -2327,7 +2338,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, @@ -2367,7 +2378,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)) @@ -2376,7 +2387,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 @@ -2391,7 +2402,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; @@ -2408,7 +2419,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 @@ -2548,7 +2559,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 @@ -2630,7 +2641,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) @@ -2650,10 +2662,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) { @@ -2687,7 +2701,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(); //----- @@ -2698,7 +2712,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); } @@ -2724,12 +2738,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]; @@ -2748,7 +2762,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); } } @@ -3042,7 +3056,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 72d11d5a..a2f69201 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -265,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 From 5bc2af68dcff63eecc465719ba9e7cab78b24bab Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Fri, 13 Sep 2024 08:01:20 -0700 Subject: [PATCH 15/25] updated doc; replace all # --- doc/src/sections/solver_options.tex | 10 ++++++++++ src/Optimization/hiopAlgFilterIPM.cpp | 5 +++-- src/Utils/hiopOptions.cpp | 9 ++++----- 3 files changed, 17 insertions(+), 7 deletions(-) diff --git a/doc/src/sections/solver_options.tex b/doc/src/sections/solver_options.tex index 64057a22..c34e1e29 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} +\Hi can save/load its internal state to/from disk. This can be helphul when running a job on a cluster that enforces limits on the job's running time. This functionality is currently available only for the quasi-Newton algorithm. The checkpointing is done using Axom's scalable Sidre data manager and IO (see \url{https://axom.readthedocs.io/en/develop/axom/sidre/docs/sphinx/index.html}) and requires an Axom-enabled build (use ``-DHIOP_USE_AXOM=ON'' with cmake). + +\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/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index 6e32fa6b..9f964e21 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -1729,10 +1729,11 @@ void hiopAlgFilterIPMQuasiNewton::checkpointing_stuff() 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("#"); - if(string::npos != pos) { - auto s_it_num = ::std::to_string(iter_num_); + 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()); diff --git a/src/Utils/hiopOptions.cpp b/src/Utils/hiopOptions.cpp index dde6d3a5..7ce51f91 100644 --- a/src/Utils/hiopOptions.cpp +++ b/src/Utils/hiopOptions.cpp @@ -1294,18 +1294,17 @@ void hiopOptionsNLP::register_options() // - only available with HIOP_USE_AXOM { vector range = {"yes", "no"}; - constexpr char msgcs[] = "Save state of NLP solver to file specified by 'checkpoint_file'."; + 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. If present character '#' will be replaced " - "with the iteration number."; - register_str_option("checkpoint_file", "hiop_state_#.chk", msgcf); + 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'."; + "specified by 'checkpoint_file' option."; register_str_option("checkpoint_load_on_start", range[1], range, msgclos); } } From a664e35d635299155f63b6577ea7cf7d42b0f380 Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Fri, 13 Sep 2024 10:24:47 -0700 Subject: [PATCH 16/25] added example on how to use checkpoint API --- src/Drivers/Dense/NlpDenseConsEx1.cpp | 61 ++++++++++- src/Drivers/Dense/NlpDenseConsEx1.hpp | 46 ++++++++- src/Drivers/Dense/NlpDenseConsEx1Driver.cpp | 108 +++++++++++++++++--- src/Optimization/hiopAlgFilterIPM.cpp | 47 +++++---- src/Optimization/hiopAlgFilterIPM.hpp | 5 + 5 files changed, 229 insertions(+), 38 deletions(-) diff --git a/src/Drivers/Dense/NlpDenseConsEx1.cpp b/src/Drivers/Dense/NlpDenseConsEx1.cpp index df97ee21..7c2a2a23 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 state ex1"); + + //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 d9b0235f..b7ef19f6 100644 --- a/src/Drivers/Dense/NlpDenseConsEx1.hpp +++ b/src/Drivers/Dense/NlpDenseConsEx1.hpp @@ -16,6 +16,17 @@ #define MPI_Comm int #endif +#ifdef HIOP_USE_AXOM +namespace axom { +namespace sidre { +// forward declarations +class DataStore; +class Group; +} +} +#endif + + #include /* Example 1: a simple infinite-dimensional QP in the optimiz. function variable x:[0,1]->R @@ -78,7 +89,7 @@ class Ex1Meshing1D MPI_Comm comm; int my_rank, comm_size; index_type* col_partition; - + friend class DiscretizedFunction; private: @@ -112,7 +123,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 +231,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 +266,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 4fbebc4b..ce3c5455 100644 --- a/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp +++ b/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp @@ -7,9 +7,21 @@ #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); +static bool do_load_checkpoint_test(const size_type& mesh_size, + const double& ratio, + const double& obj_val_expected); static bool parse_arguments(int argc, char **argv, size_type& n, double& distortion_ratio, bool& self_check) { @@ -67,24 +79,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 +112,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 +161,58 @@ static bool self_check(size_type n, double objval) return true; } + +/** + * 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) +{ +#ifdef HIOP_USE_AXOM + //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 state ex11"); + 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; + } + +#endif + + return true; +} diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index 9f964e21..e71d1e4d 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -922,7 +922,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(); @@ -984,36 +985,38 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() nlp->runStats.tmOptimizTotal.start(); iter_num_ = 0; - iter_num_total_ = 0; // // starting point: // - user provided (with slack adjustments and lsq eq. duals initialization) - // or - // - loaded checkpoint + // - 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") { + 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 { - // - //checkpoint load - // - //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) { - //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; + if(!load_state_api_called_) { + // + //checkpoint load from file + // + //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, "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; + } + //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; } @@ -1641,6 +1644,8 @@ void hiopAlgFilterIPMQuasiNewton::save_state_to_sidre_group(::axom::sidre::Group void hiopAlgFilterIPMQuasiNewton::load_state_from_sidre_group(const sidre::Group& group) { + load_state_api_called_ = true; + //metadata //algorithmic parameters diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index a2f69201..1a9db4d2 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -430,6 +430,11 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase private: hiopNlpDenseConstraints* nlpdc; +#ifdef HIOP_USE_AXOM + ///@brief Indicates whether load checkpoint API was called previous to run method. + bool load_state_api_called_; +#endif // HIOP_USE_AXOM + private: hiopAlgFilterIPMQuasiNewton() : hiopAlgFilterIPMBase(NULL) {}; hiopAlgFilterIPMQuasiNewton(const hiopAlgFilterIPMQuasiNewton& ) : hiopAlgFilterIPMBase(NULL){}; From 261ccf9d89e8b6571eb5262e1d085d4150ebce14 Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Fri, 13 Sep 2024 10:30:55 -0700 Subject: [PATCH 17/25] clean up --- src/Drivers/Dense/NlpDenseConsEx1.cpp | 4 ++-- src/Drivers/Dense/NlpDenseConsEx1.hpp | 11 ----------- src/Drivers/Dense/NlpDenseConsEx1Driver.cpp | 6 +++--- src/Optimization/hiopAlgFilterIPM.hpp | 2 -- 4 files changed, 5 insertions(+), 18 deletions(-) diff --git a/src/Drivers/Dense/NlpDenseConsEx1.cpp b/src/Drivers/Dense/NlpDenseConsEx1.cpp index 7c2a2a23..f66d0b7c 100644 --- a/src/Drivers/Dense/NlpDenseConsEx1.cpp +++ b/src/Drivers/Dense/NlpDenseConsEx1.cpp @@ -210,7 +210,7 @@ bool DenseConsEx1::iterate_callback(int iter, { #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_) { + if(iter > 0 && (iter % 5 == 0) && nullptr!=solver_) { // //Example of how to save HiOp state to axom::sidre::Group // @@ -222,7 +222,7 @@ bool DenseConsEx1::iterate_callback(int iter, //the actual saving of state to group try { solver_->save_state_to_sidre_group(*group); - } catch(::std::runtime_error& e) { + } catch(std::runtime_error& e) { //user chooses action when an error occured in saving the state... //we choose to stop HiOp return false; diff --git a/src/Drivers/Dense/NlpDenseConsEx1.hpp b/src/Drivers/Dense/NlpDenseConsEx1.hpp index b7ef19f6..967255f0 100644 --- a/src/Drivers/Dense/NlpDenseConsEx1.hpp +++ b/src/Drivers/Dense/NlpDenseConsEx1.hpp @@ -16,17 +16,6 @@ #define MPI_Comm int #endif -#ifdef HIOP_USE_AXOM -namespace axom { -namespace sidre { -// forward declarations -class DataStore; -class Group; -} -} -#endif - - #include /* Example 1: a simple infinite-dimensional QP in the optimiz. function variable x:[0,1]->R diff --git a/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp b/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp index ce3c5455..3bfe567e 100644 --- a/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp +++ b/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp @@ -184,8 +184,8 @@ static bool do_load_checkpoint_test(const size_type& mesh_size, //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) + //"hiop_state_ex1.root" was created from the interface class' iterate_callback method + //(saved every 5 iterations) sidre::DataStore ds; try { @@ -199,7 +199,7 @@ static bool do_load_checkpoint_test(const size_type& mesh_size, //the actual API call try { - const sidre::Group* group = ds.getRoot()->getGroup("hiop state ex11"); + const sidre::Group* group = ds.getRoot()->getGroup("hiop state ex1"); solver.load_state_from_sidre_group(*group); } catch(std::runtime_error& e) { printf("Failed to load from sidre::group. Error: [%s]", e.what()); diff --git a/src/Optimization/hiopAlgFilterIPM.hpp b/src/Optimization/hiopAlgFilterIPM.hpp index 1a9db4d2..4195fc9d 100644 --- a/src/Optimization/hiopAlgFilterIPM.hpp +++ b/src/Optimization/hiopAlgFilterIPM.hpp @@ -430,10 +430,8 @@ class hiopAlgFilterIPMQuasiNewton : public hiopAlgFilterIPMBase private: hiopNlpDenseConstraints* nlpdc; -#ifdef HIOP_USE_AXOM ///@brief Indicates whether load checkpoint API was called previous to run method. bool load_state_api_called_; -#endif // HIOP_USE_AXOM private: hiopAlgFilterIPMQuasiNewton() : hiopAlgFilterIPMBase(NULL) {}; From 0506d3826cef646ce640257b584ace03cd1a3623 Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Fri, 13 Sep 2024 17:10:15 -0700 Subject: [PATCH 18/25] added metadata --- src/Drivers/Dense/NlpDenseConsEx1.cpp | 2 +- src/Drivers/Dense/NlpDenseConsEx1Driver.cpp | 2 +- src/Optimization/hiopAlgFilterIPM.cpp | 42 +++++++++++++++++---- 3 files changed, 36 insertions(+), 10 deletions(-) diff --git a/src/Drivers/Dense/NlpDenseConsEx1.cpp b/src/Drivers/Dense/NlpDenseConsEx1.cpp index f66d0b7c..cd0a3167 100644 --- a/src/Drivers/Dense/NlpDenseConsEx1.cpp +++ b/src/Drivers/Dense/NlpDenseConsEx1.cpp @@ -217,7 +217,7 @@ bool DenseConsEx1::iterate_callback(int iter, //We first manufacture a Group. User code supposedly already has one. sidre::DataStore ds; - sidre::Group* group = ds.getRoot()->createGroup("hiop state ex1"); + sidre::Group* group = ds.getRoot()->createGroup("HiOp quasi-Newton alg state"); //the actual saving of state to group try { diff --git a/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp b/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp index 3bfe567e..c2dcc108 100644 --- a/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp +++ b/src/Drivers/Dense/NlpDenseConsEx1Driver.cpp @@ -199,7 +199,7 @@ static bool do_load_checkpoint_test(const size_type& mesh_size, //the actual API call try { - const sidre::Group* group = ds.getRoot()->getGroup("hiop state ex1"); + 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()); diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index e71d1e4d..666aefec 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -1581,8 +1581,6 @@ void hiopAlgFilterIPMQuasiNewton::save_state_to_sidre_group(::axom::sidre::Group { using IndType = sidre::IndexType; - //metadata - //iterate state //create views for each member that needs to be saved SidreHelper::copy_iterate_to_views(group, "alg_iterate_", *it_curr); @@ -1637,20 +1635,32 @@ void hiopAlgFilterIPMQuasiNewton::save_state_to_sidre_group(::axom::sidre::Group #ifdef HIOP_USE_MPI MPI_Comm_size(get_nlp()->get_comm(), &nranks); #endif - const double alg_params[] = {_mu, (double)iter_num_total_, (double)nranks}; + 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; - - //metadata + // //algorithmic parameters + // //!!!note: nparams needs to match the nparams from save_state_to_data_store - const int nparams = 3; + 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 @@ -1666,13 +1676,29 @@ void hiopAlgFilterIPMQuasiNewton::load_state_from_sidre_group(const sidre::Group 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()); + 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; From b78547549bab638d4f69ba85aa0620a1f3f2594c Mon Sep 17 00:00:00 2001 From: Cosmin G Petra Date: Fri, 13 Sep 2024 17:18:06 -0700 Subject: [PATCH 19/25] testing and clean up --- src/Optimization/hiopAlgFilterIPM.cpp | 6 ++++++ src/Optimization/hiopHessianLowRank.cpp | 10 +++++----- src/Optimization/hiopKKTLinSys.cpp | 2 +- 3 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/Optimization/hiopAlgFilterIPM.cpp b/src/Optimization/hiopAlgFilterIPM.cpp index 666aefec..95745871 100644 --- a/src/Optimization/hiopAlgFilterIPM.cpp +++ b/src/Optimization/hiopAlgFilterIPM.cpp @@ -1002,6 +1002,7 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() // //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) { @@ -1012,6 +1013,11 @@ hiopSolveStatus hiopAlgFilterIPMQuasiNewton::run() _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)) { diff --git a/src/Optimization/hiopHessianLowRank.cpp b/src/Optimization/hiopHessianLowRank.cpp index 1612f2fd..cba17d07 100644 --- a/src/Optimization/hiopHessianLowRank.cpp +++ b/src/Optimization/hiopHessianLowRank.cpp @@ -423,8 +423,8 @@ void hiopHessianLowRank::updateInternalBFGSRepresentation() 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) @@ -436,9 +436,9 @@ void hiopHessianLowRank::updateInternalBFGSRepresentation() 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) @@ -450,7 +450,7 @@ void hiopHessianLowRank::updateInternalBFGSRepresentation() #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 diff --git a/src/Optimization/hiopKKTLinSys.cpp b/src/Optimization/hiopKKTLinSys.cpp index 12e26500..2be80927 100644 --- a/src/Optimization/hiopKKTLinSys.cpp +++ b/src/Optimization/hiopKKTLinSys.cpp @@ -947,7 +947,7 @@ 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(); + //nlp_->runStats.kkt.tmSolveInner.stop(); if(!bret) { nlp_->log->printf(hovWarning, "%s", bicgIR_->get_convergence_info().c_str()); From 4e673d5817429cd730096fa231e0195af03e9db1 Mon Sep 17 00:00:00 2001 From: "Cosmin G. Petra" Date: Sun, 22 Sep 2024 16:44:09 -0700 Subject: [PATCH 20/25] update user manual with checkpointing --- doc/src/sections/solver_options.tex | 10 +++++----- doc/src/techrep_main.tex | 27 +++++++++++++++++++++++++-- 2 files changed, 30 insertions(+), 7 deletions(-) diff --git a/doc/src/sections/solver_options.tex b/doc/src/sections/solver_options.tex index c34e1e29..a2a73be4 100755 --- a/doc/src/sections/solver_options.tex +++ b/doc/src/sections/solver_options.tex @@ -423,16 +423,16 @@ \subsubsection{Problem preprocessing} \medskip -\subsubsection{Checkpointing of the solver state and restarting} -\Hi can save/load its internal state to/from disk. This can be helphul when running a job on a cluster that enforces limits on the job's running time. This functionality is currently available only for the quasi-Newton algorithm. The checkpointing is done using Axom's scalable Sidre data manager and IO (see \url{https://axom.readthedocs.io/en/develop/axom/sidre/docs/sphinx/index.html}) and requires an Axom-enabled build (use ``-DHIOP_USE_AXOM=ON'' with cmake). +\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\_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\_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$. +\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 2edb3f11..f16f36fe 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 From 77c88a2276a03989bb3bf37d904f3d4cc40b35f7 Mon Sep 17 00:00:00 2001 From: "Cosmin G. Petra" Date: Sun, 22 Sep 2024 16:56:30 -0700 Subject: [PATCH 21/25] updated pdf user manual --- doc/hiop_usermanual.pdf | Bin 455210 -> 461514 bytes 1 file changed, 0 insertions(+), 0 deletions(-) diff --git a/doc/hiop_usermanual.pdf b/doc/hiop_usermanual.pdf index f69a6e790444af9cc9c0bf17caead72d76b8d07e..ded0954614efc35ba8f828f5f792a77eb8d8e96a 100644 GIT binary patch delta 187520 zcmZs?b8Ig_*Qi_Dwr$(CZQJ&*w%c9Xwrz9Qwr%aU_x zsp>m!GNJjOX}&W{9j~qtHhxbijjJ*9{3)>Bd`%2b)uE9_6QOJo6vX94E0ggbZ2b_5#)G!Pt zgQ#vQ=PCilsWn-&@2VNpoa7V~Wt?RUwn5O~d;kEyPCXfJ?Mess(X z4U2)@4IkeY;r;IP7!-$|RHBQ2vHlNCgq<~=APS4(#z4v=R8^aXJ%kXi%OfV6ZH#T3 zp$Y2Nq;o$-R8wB17Hb4PxMb2(!e$Is?ffbbkKzZuKW^&!8inS@N2#9jn7=w0z1dz5 zD8Qq|v`p4rMUGvjJ#y>n?GS)#7Fc_|W8M`hE}R9NEZOKym=m$rbb*VNgv^ZAg8;LRM@wH zBCFfCDw^;0JX>9akjp9#1R5MIy0jR#7k2ev*^o7aA700&ML`uEbWj*`2QybU7jt9# z|B3!5uUMEko0CB#z*8&{sS&vTcaAZ0uyUn5KvDr}we%IYq|p55>pKz(7)Dr+6rt6N zl(b6z7KPV&v_B87gJzD#?=Dq7Bqkq!bx=Ku`Vp$^ zdxHTAh@dCms7CPp@F6COAMbQ7bd>v%WNWbT%tnA! z;^bI7excL9x#A`iXjhYBBQr!0v7nu1Lx#jE-b6w|OyDRK!SILT@fs>3_e6t%uHjZP zPIF4~(iEso)Y$h#LxJLXOjvX{Mh97&) zBQ9}?dx%?7iuMV!h4+f-Y*$yEo$LzyY<1br^b2N+^;A{U?Pls3$1;1YcPqfP8gn!2 z+UD*2aCmoRZcA`3ejcGb;%r>AyBw3yDsS_jt<6hX*!Q@2VxObE3@490Z+B7p z;w!~0*K@V?*_N%kY!%~K)^1V}@Kn9dXUQlFL5obUxJhcEypp-F50=zJ7b?5*^ixiRFmS;zA9G>Fv}SaqeAdbZ#H zZ)Hvo$m*fqkVfX4hCRFEhB~ktOC?Dn>hc z!zsB(Ls9+((;|v>_vq3C5S{ zrPG>4AMHX!Q=he|_;@M)lCy!h^5$yVa8hyP%ly|jOu zX&b`~1v=>Xfkk!xa^~P_8F0MDo3Oj|JkpZSb-ZuH((6~jF3>tp$S0B}ZvK!H0x7tn*e|8Kt0MJEXt>v2&bsC?w38vHsD#Cg zPkEcmfo-*Kh|1^qJ|IZuW%2?4%R`*Lx4*M{l`evyYKOTY*ZzzQOPIg#j8bpdLSssr z|2L-ZOw~tl<_y5Lp_k5V0E&vq!p`>JWe+oVI&dun4d6yc-ZO^-ZRbpLCYhTW%b!?_ z!Xu}vT*1$#Ji7#TW}*!Tt1Q_0@aI#2ILa8ygrW}B^122HSftD6uJ0=kiZVUL>G4zB zDuHo6?*^YlK$wk)Z78KV2HH}X?SL4TN<5``^u9^V*RXNeL>6o z)^Jr_wZxbHG~N4DI%DZ%0wLn;Nd$F7RM>d+uKD4?%w`~>G~jheV0wx@(l!5{IszO$ zF(OH3v0Dbu#$GktXwaljU~~% z>kU093Vv3ZCwxxt26?q!ww#%a{TVsaDIM!&5Yu7&W5JeuuRn0;OT3vtED#NXviLF~ zxk{3+glTOd6d&~_y{Z`^hdu^OT1*_%1R#>$Oq^^-)GmS7w$v1fI2`eP`ikPdh>ddd zl{o}*Jr&?UeUgFmPYe|fZSbW7)dAv7cD)&OH}7*&t`ae1dHz8eKv*5Y%dMN+km{_Bi6jay-30DYq>Qj8n;sJ)ab zj-dA`r!b0{(sCVV;XjMA#W`JQjiM)s6(_-ghv@UJ_Hg%m~bxutJ))!dQG!ak?U`yMCC@&0#8hB=<1a!g=Gb1bLME@#d z!C$*82{iz@KAUrZ@;~k@j%JgDfQRl2Ow?I<>QtMTRl7F3qeBdtP+y^uSZF$&Q_;Yd z?^=jxI>u4emBh5PI6t70BZ{>>KGH0$>qUvA17j!L$@(2K4-RaLMH~myD`~h_*(g1= z;s0e(?M+*-82O1l2XLB0@Tpm4T~)_%FfNlq9o-#fg&-{=N0PuUZawk1hxSqUW|Pzf zd#8;W6%-?#zhm@a)o~~iyo{-Z_!ycIA=nT;#=Rs?(Z7D}OS2eamu(~!5wQ8&={mL` z{*hf_W|q?W!N_c#FMHm^IY$~_Y@n5HE1uQMt5!>AKcw;j0^ohzkh|ig5J^n>%;sce zd6RBIL2Ot|oy4_YCf6(MWSg*<7e*W%T!v7Odg_O@ncubhf1jZ+OnWd;FZnPJd95w| zbl2&#Fr-_xhw_Grqf{k6A*5Kyh+w5zJ65;pt;UeA`zIIZ1KHwwt=ew2sOIvgb~11f z6(N?37dX`b&5a+cEt$C4mJOXCRZHEX~IV68z(P%|ip0Cg2W<5%dfE>d(hqP)gROw*N6QHm3mj&Hsvs!} z8=5X-^RTkJag!xyusZ=AFkh2Y+5ax9B=l>DilbFU&$%IaH#j<$qiBAjs2I zU9mMq9RD`Rh-Qc2;1={q5hqK-v-H370>RZ~XM`jWk_iJDN3_Kh_*w)IgK5d3)*>Gy z(_r~dQ5A1<3xW2X+ntQDsXL7~HI}WvYaMQVH?Z55*Jb~xaeT|t2AIbYt?Op_uS$W% zIdsi-S%nBamFpi6CfbTus$$LY_GI7rY@h?U2sK%+dX zauAybjz5boj$tu&g>t^*UBQ85P6RK4W@GjCK~ufCsaT+mj0A_HB_pRZ=O)Hc zfx(PI0Vs3ZSy2Y1q8B3rTc@<86)?k{igyN;s1cY<^l>ti_QRj!5~XWq&^V~ix5+ce zc<=*dsn4@bR8~7yE}VHmK&h(Nfxx6%A2 zq!tXKh){AlHX}5GgveCyFvYG(Nom#g9v1Nl7G;9E9nO0g%lw| zoFQ|rTgwY5dB3v)Q}K`2F`)=)fn0{LfZe1s=M2~(nj-3bQiweQA(8@Hge#mQ#EtOL z0MjuVB{!za^_vVBu~lZ^CZjz3$#Tt*h!$nOEj;@62y0I%yG*H1r56GDj_uzT*5{88 zljAq_5nCT^oPB*6?uwezwchswZi?!VOa}IEYUWe=9-k-a+HYSL?k?+&v$jfp|CZhA zKNEDHEW#DB+I?TFcdnG~KCy}D+P__V0SL~poBpf4*Ylc~*(pCgQ)d1zy1yRGfR0Zq z1Ku8_br`xCTsYRf0!9f3eZ;hBbgACDQ}M zdKM|h{GX*!Ssz1s?f(2NUk0O#iDsX_wOk}aEo9vhejlfVYrycG^&Dh;3?l5<@^k1bCj;6`jPz7y? z87iGB-ToX(5OU10lnRQDK_pbEC(H#A42=Sa%)~YlUKu75J+VS7 zv=)?nOkdng?WNuc{A)~S{=$v7X18P4iq?S(L9sYjFaTs6e7f_ z*1#{2PH_XI@i#6)*FvA8AIUHu3m9z @=j97h;QlGDM)Ck#)* z@3D|#isPO9ynar6%g7;HjtI~!rP+aw7Ikm&uI`x^~FpTXO7drZQs zksm0@M7-A9jj^=JU=|WeC*Q~_p%!PD_|j!`9-L8!IaQulNwqOi4laiLU!WKdLtghC zApK@lL7Tv)ZsS0eA7#)~ofl4a(A4vQ@G)vofqHEsC)6r5aN<0`S_r+tn}u8`(*ypH zU(@>Q>oDH;-di>)W&+05ijQwSl?5i+@Fh3|)wy@FkKg!d4LVKJ1;S)q0YR_fO#0n~ zH#TCdCsfpa#s+V-Bu}FJMLgYC@I>6*SkrBLHg`d4^N(HY$8)P*Im@!@T}w?Kjt@s} z{_e}z+BX|-JbHO%2$&RMXAb(XvQG)I*eJ85$=a#=y**n_@S-+_`` z#x{StqHC-n!KVOLAvGi;kdET~ArMYKfR*2FHodp&}ZzrTR1KvYJbr$ z{0!+_kS~kkG`L5#g+v|@C;y0WS%0R}5LpB{swW}_A*TaCqX2MLMC5Zg8rCR$wm&Tm zr@T#nI0Jmth86*gr)M#k-w`)9A5DgipHC0}oI?F0U*rncXZlv7ZKc`F!DxO2#G|_3 zkNa!zL$i9Pla<_Z-&F0C!yh z)`E0DS~Y+bsP&ao{c^<}QK!3OD6-Fw)ARjjeLznpVbG&bvrl_we(vwh#>QuhV2FPP z?aos^z)OUF6IMsQl9IP~B}Je3vNP`Smoq=|9Pa5L_MGl?{;C+gm1s{e*E8nkE5_8f ziRUU$5b)6B+pkr7H)}_4%zVv+SW-C(6@pUUFMn${du$_`Ti;VjXdyoyQ5E0g1s=am8+-{g zG4dq+G~%TD-1k8Xp!}rTi^hS_QUc2 z+zaq?xqok`lsv<<7@ob(8Mj;oMF1j+~&U0F0N7@~mo+B;B z1Q!N2v`=6M)jMkphzv0jW&ReAeI>?fbpeV3M^Kwego)0kNs7zxcUjd~ZbexR31T$( zv}Bnyr4rlA`Ou`^X<*F-wc5YeaByL@$5{%@2Y75~US`vj743@Y=_ER4O}780#e^RZ zpS%3kvnA4}Pn95GVofj`0E8LwIAXPy^^601C>H;)lI@hIzs=LE_YF6e3_}ap6W0%l0^UODBr3V#sX?|k4)s!p+WXCbRWOC zLi&7&V&OZTq20vn2kM&esxCG+JMi=iU5kq{-L6jrc7PVVi35c10=e6PH)>-c1H#03 z{+0aGcGaCOjTN-iPoE?a2$~6aJ5{nHrS&s@0gbdbYpo9o(zZK+bA`Hz5JAN2O656C zFc0YczBK{mvoGroIT<;1$nFjmcnu*EdG7y*LC3N1ND8OSVUVNRKO~P$_c(g>U}tI) zl?X+k6(nXfA_4de;*g8{#Vf?9Y_>rZsTVm>|9X}3;(aeY{PI^Jm|VAQ=8)}Xp!{dw zP7S(jYFOglrr@$>xn?2$YWdefeUiP~{%^+*S+<*)PF#99NaDVbTSV$-xk8GZCo##M z)Y8i*iu2g8iPdQUHW4Nj#zoicoYSAWr48g|ijTY!qeg(5%12==rtHG8@q||j3R)+K z_xQ@K&Bz6tG@7scpSvwMDD^*F{ZTz-@;XJtEhs@oWcAQMM2W!PbUh%U^H&#O^pI}B zS^DrK@U;5n@<&c@vO(5g`uj9Sb2Zk@?n|{ZURNc@66@S@Q@C(mmYa{Omo2r;sm#O{ zrM)RQoT`8h)g;t1??9>ll70I|J@Qi8RZPnc%OWTI z8kzg&yMv4{L?77KL0>8eQ&rk4hRneduBoz-KQy50RiM=;)W;XOQ2P*V1zw&x22`Ms znBr*OvU$c4>&`zE6=DXaESBV?vX2|E$~KxiHRup%8HpCFq{?;M&z& zH7_#!3S^dW&;r$OfdXXM$rT@0>gZ}9c0(v(UC6ioqP=J@C*ZAk)oxj5jyW)JT1$>= ze+uB5F0vC=Iq00l#GYmL?H82g3TdP)JpX~|I{m6Fj#%9nuRFM)!)mHA(~_U9JLsgE z6aeG^Z@({D$SQjVWgpZ!^dkXA5fy1nVv;ejwRm@>EU(aS+1qxZqx&f@ai15k_{S6y z492=XV@zt8HZh#E`hCjrt30W^A&rYya$lBQW%(NxWVNLPx< z$z;>kJ>wqo&CdwpDE?*i+zGZykiZ$GSxo-L6k)m|{G+KSnRy_tQRBz{-Q;%U4`Jkc zpB#r6W7x-YfF2f(VL0>B3w#$2&c#x2d*ZRWwE2boC|VE zI@|6Ds%mFTk#{!h;kDATArpkzCVnh+z{ zLLbX(&td5$DS4N-L+U_v(Bfdp;^bmGcAq@t4MQgaicj^-vbjFvFQhWDJ1YPa16XaD z(zOo|ykaRyyq&F`*K?zkqxM2aO>#4SfIo+=)VD_|-{zN;*Z#-#V=;zc6nfyRULWkh zFZ%`I&Z=>a7=o%->+EZ+66P_zu&K7InY}OQ%h9Dxl+MJ$5QSk;m}_B0kYR9EP2o#HFWixw25SvXtDjOK)XzuT zv^wY~R0InpowLXee>M18BW3Kn@9rZ=PM=SZ2kMAEmOJ1^HYMHdX#$8_In&N#4JUJ? z35PHd0uGUxNH2KdIEqwSS6!dvtXT*`lc;pb7;QumFAl{XR|tJgJ{2N;hJ$8eC3_yS z-DAeRNENjRnR=>mr^?n>;lyvtH*kOOJvlT2Dkz{cT-jK1laB@#4c5rK@Lv2=(IWQ1 zzS~y3O*9%_y6&a|rw>TbWW**49)7bQs$PvoM6SxyQyays>(r*qO%%OwZ>z-bklqZ9 zK$p(UK~XNmxmYM^CMixlgQVVTyQ%Z;(IIEbzAH*~KBh&A3DGiW0^W42azCs+))T<|X`_GSo>aw=LuF6qD*&WiJ3xlAm^}uAB_Tr4 z9vJnx*&aWnkPmehJ)`_mo!l0^IT3u-P*^JuNRV@qIHy}jvF$j`ZuzHmSL7a%)cvYc z_4<7dKTK|gnO1z?+ z17yJq?n^Gju2gePBel&}9KGa?7Eq+FsW|5#6kRN;(}26lYSZz0bCrVR=s#yPS}VYf z-5rgC9S{0?h^2!PX>X@KR7OEUV)^VkH3Jecq+dWkY*hr#S56O-5fG5P;KNkRAH!za zKt06op24`E_L8-5}aW8(TyN7LGpbXu&IRc~+symPX>K14x=4KaF!OF7Sh zQqj>BnNE-zG&_8o^aU7mhPinVw}!v3u-a?5OVJ%&)L|RzzgjZoJ$aNfGdK*F4E?WA z4Uhs(m4Htma)C1;c5i(usA0m_DoS9%;+s>s_D(8cD38~~TYN#pg|Xn`Jhr%mCXxl} zQNUD(>te$>IN5C9!~1nlDZ@O z!KDey^eA_AC1Bj0R}24xE5pRPgS%4kO?--?9#`90)_}b3E?Q7N*uv!(hwK!S$^j={ za0RMM`@>=l*`F(Sn5h~E2c_V4rV8H5+tnvo zPBU0ipbE~^VgSw%C}XQ>(?l|jx@K#VsS1yR(w!u$x^maiGRZWQ4RDG_(X$BRi*)Vm z2rU{Q0015Km}^nE7g1FF0)v~|-(9C(I-0z93Q znBWfHL#o!TZ9&njnj#K#=<*PIiS`fM&tr2Q=|JKJePKNUnWWoRpL^Z|b*ymt2AySf z^A0vVBr)aB))5-WvOR0GmF~@Ms~i-;IeN{G`+N#?mU2$A4`L6J50<;$rGTFt{EnLo zV$Jii_@l6d=mfXBuJwSiN)hx01GvQCt{HrSjq*}sm3*q=?!gCbgH?k)!P;3-gW0w^ zJhN_>>e4~x9)!&>Y7t#+yMKz!J?vIJ54}#sU{p8goSsN{Zg6K#lx2<*(pB8=RB(3EQzcBPOCul^DDH?#UpZ8GC;|Ifk)-q%_W0Mv|S2-$U zeDfG5nz`SZGDrMJ8+R4WFwvFA3V%`U!pUR(8SgQkTMj9OPcE)l-fkZ0m(bY>3VqGe zX_pQ-md&1j4ZGsmljus20Am-wYCr?&b#|xc3>u=RVt%*g67S3o+p~#|r?y(@HWI(W zuW9Y{bk7*(^D&hW<_)Zo!mo^Z?L{3FTFh|M9%pRKPG3K_(!$b+fD%=cQG~1^m+h_i zN?F3$HhArMK2vaB&U=|Dd~8zp8~*wGzG#NNO8bTnB17<(CPV-|AO!n zj4jmxr}zR-Dk=$C#dTU@2Iy?CQ0S>y6N1l{SBbLy`Ad_O0(xP|ggur{i>?SSWHODo{w*GCMh`#pN7oyJKmto!*Bk|SeVm4ZQPYU$SFk_+g!n38twLO1~9 z4Qkj`^LhntujrP)#8=_y`N>FB7TO{Q;Q*y60B;=f7kMUkz~Kva{*{>Rzf}MmOG*$H zItUl{|A+U}$y*$cxAhCV717W~(t)q1i6a{$l|3}Px*IgSkntm##u`NwWPe6}K5rm| zP>zKX8pl(1Vx)rQ-TnjjijfR;m~Ia%sF(HmhJ4vXEXh&T3&Ev4PRfS)Zw-Ja>>vVZo5?CSo#|4t6GEq|5R zaA6! z%+1fkhUxaY=V30i#_swmG%7`TnXk9_d19A&qL#l7GP4+$xUO zT!{8J;K&<~S?wbyToVJyqH)?I?N3HWXGKCzl_=>UroptWNh%C)@)6`3faFkeFAM}b zKoi-3fW*_a5${tZQ-_`E+)sg8rZ94xZj5|P+$wnKBHr+fD(^@1M%Cy@VIcUL0Gw)z z>w3vfT3ER?)&lUDEE(ZJcE++)yroT`6SHa{Hw}v_b5TH-bY8MND8p0m1E=35)mS?F zAg&(G1j3_?};W8EnopmRkNf`B)ZV zEJHA%{*lC!95`JQ{Lw!bR1Cu9!4}bB6T3?q+6g{M`kxHGaTk& z&Xid3bK<-)3qrSS}lRq&RysA(Ro|ayEJO1E+jnhwDE6F5*!8_ohJCPWw-;hL;TzNn@l=ASOyY|OeK<; zUa9qb*&A0=A;?JY45Y z(gK*BS{>y%V1UWSXC?_sG;|_ z>)KH>4`Kb!75^yl6=h_EQ35)M%|1ELX|CS|;(xRc^R{E9x8t^6@G=?qkK62(B}J-ZfYZZcZ~KvmvpKFQAawc_E=}=gYo-glftg*sWsq5N;eB z!UEQA~Y*C4b--s086$-TK2Qg81GYn3N zt*TV19Jnh`Xawhh=AN8B(_0Xs`wMivEQ-05!c5c#Xs80!8C-nd_?r^;CdA_@Ht>>GLk^ZeaWox;G~@at)Ebg%dM;>(_8P* z6MYFGC-L&z-~eaIwd%YRsM`YzM5w~8Py&@uf%kZuSw_0VC+O5cw#L^2FpYg_2G5g) zA;Nc~LSK)0i-0kl8XiPpmDaU)zT))T;`3&ajLq052DI(OT8>y}8Zx zJ~GwtJ99J{Q@wQH?5AQcRzZiT$+J8eSU%i39sSI_C#u>*W*op8Sscr#>hG=B(DY&( zw5{;|dYsd!tb4j>T;hQM^1YE|cc}}-m;BMCqjt}nDqsIZ16K(e(^8J53g-!<=eaN= zeODdm@ee+A%Pk(qC<+XURPg)IW34jRN2_0i+m*qsd9-84x8ZWQvPV zt_e5^b*Y3P^GYxaD+`fLe{_|J&-Tb{EL!ag$d$kPzGb~?lF;;E) zar2{LBiEaj{=KmT)VS!DVQcAZq?tyMC|*9!|7TZYXH)DqcCKs0R%>yct-*rT0?0}R zxL~$nx5!#U+_hi#x^d#6OC3N2itbkf8a&pBt=kBH61N85$>Cw zHeO`Yg1;Jz`?g#1epClz<~eY$Yv1#LcVquh3dXA> zvB^3#Fj=&2F=82FGW(0?9vvl(PVX>Wa&%AxC!j|*ZE?0r+Pv-9GrHp0IC^}|ZtnOl zU|Ap>q!IW5X4D&_?qCk)WW5>$9b^(Z?$?^EJ>vE-xhYlV6LPfvYE=6#3A}u;Fg1J> z+xSNzPl5;7ulip+5sku##So;(HnTo~K?U{`;kCtokU_Nr}R z*4-*OLA(`MzJ8+{UNqYR$}y|eJE^lx{;~e8;@k?l0p>RJ*ngjPH?+G=U_V?!G%&ro z`|r&u78Ujs39rPSNY%&)LUhwx>!0m?if7$oW5) zW4LD0arh;iyz*khv1n-u#raSszkdPFK4r%p8L8>{(;w&nJ}-6yeN3%|Bb!DV$0^4D z5Y*|xN>Lqy#pvus{h45&H>DhaNPs6LCvzrL`G+}sP_ZMXSNNbzOJB3x@^4)Ew34bm zhy0XQGaZF9lh=gpXDL1aNrrLN{XxE>pZ;ovCt)Z3i3ywjDfPOv^GuBid0C?+VpK7(5ScaVTK-_J@a^i0$}G+%{14zSj<$KIj3vO zQ!Ff0Ir3Ikl1|W6ujmYH{c=VoFd|hvxM<74SUP3Do;g!O@P!C_X{r!qJedXI=*u{& z*APPvU*Zb#PRQ#2v{?rHdwZkY%gz56Kw?DgM+O3*)zZ`F$lA^^Ds9X8NHUq~R@RSY zQ2beW{5y=#S$q4E$9%@;@jG4b#}>vZIHag{Dw$ugBafGyU{U6E;4ES6$g;TeAN`* z&TtRF+H?*xWboT~pz}b306fh-?#`E*PPO$%fL_rzrE(sj{}EL27CyBJI4wDw4roCGQahi{>34=S%P}mzFjD&avZRA zQ$vfTSRI1zac&TCoA0%*dTB%DCptdFQd`vZ*_o0$O{Jtb$^|9GPvPCsH#kpactyl5 z-7XYDecooSM6KlOYeRSbhXs-z!`GJFgVt?J0RQD;EDjFAzv`hRt5|_5gem^p%@P#g zCg~_XMU{2|pXF$UH%a#8HM8&SyE7MOwjd#C}91 zlu+9xXXu`pcgdpczLP96RjDW0`kz-G@G`av`RU5NXERB7OjFu`<4`39G z=(WO!7v=Tj0fol^Bk_Wt))asABCGIZwaB6s>8#Eob{Q4oXtl`gI+2;;XhlBRhf>X3 zejsNpvb7eMcb%hj;;Eq}IzMz=T%%>Dt{p8xE=_4CBXL-STZQ?Wfx3g{qIpF?6uK06 zPi1^7L_D$?#y2Q1)W7nM*O-Qrk&X z{y7^dn7f;evKZU`l1f~$DvE+%U=2ebeBdjqf1av86(flmJnt@L>6;VB^-)P`8NUfr zB&iuBk*)B_a$_oxTk0kz9cu$1QwMJbV)(Z9{wu;9BGv}hE-l!6RqP*Fn0DC0#gcW% zSOHBS9aRlx3{d61S^Uv7P*c^5UzeN-vXXj}=j*L0@-*4Tfss$wf=${GNj!9&M6U){ z^++H}L)s((?gG-gvNP|hbxkQB!g?DC*5+On0Z(hiSFH6SH7`Fa4el@Ctt-;;=_^hI z-ljvl5yk*|9a?KYT+>9J3}l+oF$eBGFh#0nSL&1}re3bbfEjYDSgI#f)Gb$spGs2e zO}k#7i;L^kU4a#(P$?`L+gUZiGLZO^^8>HohmEf-2i5Sc<5ai8*=h*v}PX?odFYDNSY*uKp=QKE0y$gV;-;njg?(>Bn&(r)juvKAxp!J0y_hN27d zrEtpiUjB+kY-lrf((pMO;Co@S|4y535~5UeevE&IrtNGgJ-qRI`!9It1Zj3EQR$u&H41%l)Aeew6ob z0)|0DpXq-?v!cbD;YYS}D><;aB|9Wad~$5N3qJ2prPhM<$B=)*mX7adB8vd@S)?Q# z>6T1Zs2>p4zX)W*tYYWMy#X~7dHz}c^dDdq^V6*RJ%T#6^98NLM{f&Q>-LeGoehjn zKK$l)VRge)2I-t7jH4pOF{Iomr*gNYoKI*r@|jU)Yee#iB7kV-r;drX4b?A>q3u*| zN)o0OXn!0W$6y6Sbtv_NK29P8E-y_oRZAumjAUG>Ld2M~#mZeZHtkPZx|%6m{b zF@wE)eSlwq#4gRlQN<+A@o%yD==+9E^OVZIFXz*azj`|~17~k%*2ba&{a~~|{dTr7 z&W|gbjkpHQ;Y)G^X?`L%L=JPcR^_G)CDtM$P=QfSWhDrtsaZFgOlfi4AdJLmFaJ`U z+}V8@7Q(y;DzqRhH2|E9yaK|dvnn&qxeiljF#PPfp~gNz=d`!%cLjRI^D`QP9Btsl zwh)SuV-bdTwVV5^jeQ!l&8M;jItpj&k=oCag1R7~z!rZq(7gJQ(DY6<(B2RxG3cKA zGQ+$dg?CbP22>PXDhPucxj7%;(pfLfJZb>IdrS&j<4<}epkmQga zQ`oKluStm{)HLk0c!3mhCu_Q*hRKio2q(_K(_fni=#KXL^CKeb4bk%t!9ae}Od6WK zbo{#gKxm~h*n|DwD8rD4=bU*8)_rus!y%5@TEvWC4tZB^3o%Op`FokmS}LWxMDDbW zQLc{o_AopBDvf-e8*=K+mb^Q@*2~TpSKyu!A<-;51sra*a#K3wUt>%G-CVpsnGo%9uw?x)qiDf2JwF-vTKT zn#z;KEq?610X7x=+()VgRRPl#Wj&`U`2Ub~kQP6y;fn^h0wH6Zkn++uR*ozAa#My# zh_BE*3X`*f#Q$3TOx+EB3hmvILMO?8I9i!1MfooJ_@RxPH2%T_0a1C_5bV;O8A8hn z%sTR1$4){V{^=Y%gmb1r2vTE`HknQMuy>>`5wt@t zE*L@i_S@Wlo)er!poa(?Be}AE%J^1O||cq*kgZvPzWmtHuOSOu5vqh<*EuTAxvMrV5|KCadKqVQDTr}*B%P{5~GLmN$8u>M=SYdl; z!`XGQ6{|2_1T0%F&nO%5epo7TSBz6?VyOJh$li#PhpkA&`DI1_3&IKjmM|?MX7*`- zeBd%f0N1krj7sDkWtQ-7>lui9+AQEutN;bqyswvsLJSokYRB%Jx9@)4<35PJTt7R} z)-M%Bg|PKp6T|PCT*?#$KRGzg#a4f%SDCU$A|jq;+F6%yHs8nmA=$raPDg0IaS2!> zs7{)}^zBlJ!7zPMW3d1a8!Tk?qc$u(8WPzNd%_MS@ch#QzLBJ<@V~em&qBb+$5ydJ z4vc4E{_ag}vYQ4UTq2;PoCTC+x%&BR9ze?;qv(fU<-q>4y~9x>qok94cEvc59cuY; zwH@?b-V3tbu7JCDn!7pZM3aJ~CgFm1r`EO=$ek@7$Mq~$l$8!B)~7D~eSets&lMiN z4k7aJ9P(e>3?T8;J2t8ezvp<^X^fT=dJZr7Jz2Gc$8J+PM0glyjk1(mGl8;%P<;Xr z8C5b%=w<62Oc>+gC(z8boQ`aKGh@7@1lqo|L}Bh)n0GStP=IHQ^A7VT@_S+lV}bb^ zY7~*Y^;lkKLoEU{<;)MSLY>cO5#89q?)rHUn!{pb*d{UgSOqu&mxys$i3-UuTk(Qc!EY@toLSB3j9 zW=w~D68)$bcT$C)zClHF?EVi`&i_rC!_M))Tj!`yc(_^qpWZoE9uC(3lk%I=+IGO_ zM(}$!$RBxZz7>zAL^Wm6=$W=fc*ix4EpbJ%DCFR){*n;8`^wJ%I7-<6=KwY*K9p=Uu=!*v9 zAAsedB8MP22nQMIlLn1h6$cE#BE!@9aL1r^6J<~vcr!0t{aDlQX$%>5XCPohe>B1L zSCE#(*-8_*Vw}){$-B1wM_4na{EOvRNg4|a zD(tWyLeS`vRLpU__sGx(DI#Pi%P=FV98QGg*S$n^I;+^+*l&MT7c3t_)8N4PgZ#XV zfW-p3e*P%B4#br(m>hqssFgrt|^l*3xREoG0_c!Y)3XQVt3CvJyS9^^bC z-H?6$bw?GR)z^jwiI-YJD0by7Z|LHC45ln#o7^G|%ECobWV5Z76b{%ZieXQTl)N4;Z_67;o}jhhwI{BhpEp9sgHg1M!0kB* zp{<*t?r`>+E&)`S(qxt5KZ6jkjI@JXx4bU){7fC^_^8!JO_LFx$$i7Qity66q;z=H zNgJ1T%7s-L;;sK7+1J&ZJCj+c(x^!{j2~f2AKSWdHXa+Z|8YjQL=j(JSRI>s>EtvD z|7Ou=SX&)@IoyACk94Pwye;5d{A96?L2*CY>E)0r!T`Ga3Hz=&YZ-}qk%)quz4nub zmXQrDiZU_f?}@94IKG3Mo62z6tBf+HH{M9;8EZ%U&->p**zNS;K|mq2D}tK2h|L%7 z6Wxpy>w8{BLaij+tR%W))H4RgOyPc(>OHCHalD@WwmC!MwpO$kV(d63 zuB8H^GJyYlmEb=w$!k0cJMz6h`c-{6{j1!iDzv^a7(--|ixhBvA_1ka47e{r*f45U zqb~8}{^a4jiNe0jZQUPQay9M)V_xE!DU`Y83fy1LfcZesnpLU(sd7BS===r@|Y~qPKMz>w22BXLDiMulm=B$O`hDTa^t|1 zt_$HpdXJ~n%E!YuOM1Aw%tSy4|(ylZ$1be-yR-&aOzt)Hx$`?+{> zw%gQDVP9r<>H2j0y$Xjn1EjLCn|+|2+7CPd)vgxMf3jYLw8xgQBT<*ai1;W0jPM&n zCc-Xds!UoWtjPOyr8ld;%~~w1D~^9H8I~kWR%_VN@=0s^chGUr7gF#Shr5bwt}Q0O+T@?=Fe zQwE!H6ine%8=g`}QYWi$r;asG{@lT^C?%K`;?-6&c7-%yhjY>lU0?MWw`e4wgE&CP zkVT7}C5K?do4*tDexTB@*@H&*_QT+;H%|~>CB#YKPGyacz|)#GSGeUoP)~T1UUa?S zO22Q*Vfy!^ub#n;a;R5;#nzZukJ-u(RN;JQSZ6xu32R9o9=zU|A{c%nOa z)`p-+l4XW&IIW$Laxhwnnf@=xFL3;x^Oye*e*7On9S0Zd|9~H}aVF;L(gK${S`HgL zXnr%bmKOsBi1O^6w3H6ulNa{>tMs5pq__~ytpxN7jGjeTx;=iC5lP}>);ew%IYaV9 z8wJu4QiKMWm6qKetviuD0=t)q#gR|7eyx8DzWH019h4v*_TcLsLsvug6JQIm&NtDk zmQor`N%Lw0V(*8B62CW`o`7ou0nfL>CinrBWY?W{4D>Y=?Q^II?489&d* zz0=qfW;ENy6=pNK-z-qSM`zT-wfyq7qx@QLY*HP{KgcQ!7OvOJbmtsodZKej|5%by z<0?t;FW}qrxOwfB2l5w!U0B_EQ6~5P4Ee>wASml2jgz+Zw5!_lItRr3>TgQC6u*%m zONIQXIK-wq{^cikD3C0q`Wi9*IynEkpJMs)50IT#AnfPnl7eocdDMa0&EI{Wl=*(G zOJ9=6vDIFXEw}c6sfm4l?tAy{%+`D8*)n4N=;ZN|!0=zW!sWIV0kiWCGK5A)f8fGR zMBS5s&K^g(6c`uhU9Rz`{w`;0Pk#f%i{U$+H$7Da3A1erILIcWjR;j5(OByo!9Vg| z?WY(bj6y#ZYLRl(y~OD;qiqZyH=b~YmoVBP5{D%d!D+d~vy_mi?L8R9+dAOH(M(%^$3U*GuFBJ+g!^dO=y+1U3?0 zRsCBWSuNq;C=EcynNPS2GdPnC6es_yf|*1CY9`te_MU^Zu3a z-l*D7POWOTc0M;nV>cCh+g{ZFTBeW>KftK4fo$9folm)tUkXv;zlfaM@@QIi7x7j$ z92r9Ro9k-%jedf0JIOoGCZ~N=;(4T5A@9n8WX`u($IJi4JMLC!;f&Kh6bGpxoNEX zq^4T1iBzlTGCZeB?92W20bsL&OB5_CKJF`OZi0}9s2qA=L$Fsn_?R2!d(^PM(qE-q zcpUCN>&oAq^^&K$r{Q{%b_!&$nGNz`8yi~D7A@K&rm{NSo(PeP(ig1{7!iQ5`kjK) z{H$5Y@u`K190D|5M*Q-=!8Bo?)Jd(3&fiu?Z}Sn3SfRCr)AZGxAX(Z~5DdztU-Uck zi=$vsjLj;6FDvocAbs%t!2 z1mW?@>iCW8YsKa=m(A0q`g!O2#b2>Xb;fe0$^o~18PBWoBYY$UIqq5JG{5>l%s#hZ z3vJ>3oOx{N)nMpxQBoB?5Bl-QiY>GEf-}00#`f_uDxZd0c;$~S;1>8?#V{}pEV@AU zU|E`=MI*I;E&#V?<4prOyx{A~v><2rnfFvyw#H{`Dxh z{exf3SNPVIB1aa~iR-nr)5?I#EYYNCWHLX%m9C)A-R~DkE%QMSVFM!;XL%3wT&ze) zn9bbCKVzN+mF*5JM@7JvSOy0EbTJ{c=l8bakkjDveTi`ZtQB00IGbs|P^B-wNgiw! zIhU{}2?j2-<%U5HEHdW}npU9&I5wDhP`m#wnt&5NS)J|5F((b9!Iw|ElD(P^)gS|s z&hNi*_}o&%h4E(_CU#oM^aX=Pe|uk`SSA`1~6!lquBsbZp@9JI^=y^IP1{L zAh(dNfRGhF2F=n6*Jwwe5z>m*3{hOnwd z72ch8R~8n!&>p&McM6pXl#VlONX*wZ_P zmdjKh67_b%(SHIw6!99;fb$3iD3&L&(+|n~w>t0leQYIk$PnamKsz`caNBT4v3lfS z)$^YS`K&2{LRmXSM0GX5x_rdoz6;r{DQN*Wm}YO9VH5=V;|IIwfGmCr#qN(gG5o~w zC%eIRcwzC^U>Y2`xpU*~p*YO!e_|gV8)c>Up4!ogEaOjowpp_ju2CiZZ}%xTnS*%4 z%1IFS8~9a`w(AFpA{dZ2tL1wVdY5mjv~BqkYl zO}+MROlFTfwT7v!)jjg68`(W{sW{N3+&@3xO?N%QDVG#fFO1zrP$VYhB0EU&&eA9g z?-xTb?Pb!!8r=}BU3oZTu|-DabncV7)bn=t#Ws)@C$DoELN`uFdJDek`ez{lbB?P2>_Z3%2dcort`JAOhCIddfgsuW1N|jjvqrQI zUAx{D$r`Ni^US720zGtYf&sWmnJ0$O5=4Zl2y%^Cz`{u4>Id16D}b8T#Gf=Q>ztnW zZFyf7RNQrKPhFDMEeV2w-kKKxZ*;N&jN55?L^pp;9IkG ztGwXBr)mvKREQ|_TeyTM2G@2V*qKNL(!t9XH{}kc1zZQ?cI&J!j=SCyv=pFQg5FNC zNH_$aHKCVr4A%NV>q|uln#9~6|AvjYeMLZ`Sf=bIJlJl#shai|^xXV`|B4%dY?#Ep zFIjFoi^TnAnUr~0KLqr;8Y*b0^uBhNA#rYbfkZur2M~TiX_omC1DF%R3}$+!$Lm=L z?P<4jfB1LV#UCm=R?OC)0u%JT9?Dq@mDpdd8|!9cQz7#wiUGvsegva+>j=(E#zvSM z4+IY%2#+?l94MFh$A^#V7EgX0FfWGXs8E-@cN=q*`Eg znp3rfZEbtzNcT(yJE+AH_C~rf@MWRzKrsg+xUXkyOcy>~Z(ue+Ug!{^e8o~n3 zU9?}Ef_o0^V`4{EylMO@jg)a~3|GX;J!EQp<-8tjBivWdC!=$KYNvNkp6&Rbzsz9N zxq!-KJEvjIwo#B;;?XJ(En@O5Ay>soO(ibVn~je8Ah2wAP)~16%Xaf5hvrxj)Tqqo z%v@Sioxk~kSj+ga)w+D;Y;*n9x@;bTA}q(V)elmtm}TCQnYx+61>mMt8s90B%&tK0*60V zo$;j-uHWDWZ3mff>KT@NJxqvyu-~m9$vo9y-0!73a`aI_LNamVHC@zh|ENIrM?@9>6`yGhG^~G*&M|EP4nJ6i0R!UyiN(n_h z4hd}o?|1afQmF8Z|K$2*3Gmfun3)QRxBHf(;Ui{fzQ~7qN_dvFRDBejn0zUHFofY(}tz--TS z8rxX~DIR5koc`q~px8s^mY7yw^#ePimPQ{5VN)`{*hc%~8B~)_cL^&kGJvZSriJ!< z=967%UH~}y2a_2i==6RvsAD*TocV$}dbaegH-0Arn|y&wGeJ#nyZD2$(Sp%EFRtl^ zkO0YGSA-(`UrcR&OVi?iWaN{fsu~27JdGQDn1=jaW=>=@oAHGlVwLbL2#nf*m!S4Y zlBMJL!XLFAkoHNjBm(iI`oao1CHc|!+Zl}S_&u?K`JX)AK3@EmoD8+0-j19*fGo1g zF9^X?9lB7gppg_&9H8hW7NaNMLu4EFa1rD9P3&XYlQdBfo>Gz)XmUEv8?+2w2rV`J zo=qwXR=+IE*UEK#^o7WYyrqBBKmK^Hbf}7ZDoV=eAaoXy)T ziNJo3`*f`h51D1@h6O+VRynj6D1`kH+WY@fxyVGBwm_`SXdA`dQP%Hi-ncp?F;XVD!aZJ^LHb=# zTMo+a?*P%O{N6yUZ;c2*VcY%H8wGA7bwA4fkS7QE0mbTm{EwY*Tz?i)m#4j(h5cvs z1TSeZ;B&6>@0ss8qNK=D*4j~a4d&F3S@jKg*`UpiP_{#X`g4a6d{%S=0+XS_lct<) zn||Sdnd7f7W-;U;>~ZE6%Uw5yqg!L>)J`JUrYoE@{4j}m=8uW54G%7?x^80>9KP!* z=Tz%SkHP&tAarEC&CI^pzJSRpa5sGK68FXdV5c3Y4 zUG)cSaIItfhGm$;;ja-&!U5$;Ie%=@bTsyCRDJ7ctBx45)7dp|Oj}$qQ^ahq4@^8) z`&lhFoq_$E^K>;gB#*Z!cphD)h!WsyFI zNZXSnLRq1@pQQ^rF!uR7-c+&vX|=KHZ)B{)XLfqOZp=1IkM^5Su@q@mHSAK6dSb{6 zCS`2^LjL1Zb@7UGgjWal&X`WmXk)$8;Gq1hz2rV-&$&f|RSlC}u;%6dWmwvOZaIr9NRxJ8E zt{CgdW7}LE5UZ?O>$w!zcjN#j)l*eaOze!6y2-uEwQ65YwB`#LBhtKoW`9Fg8(5a5 ze!Dkj>6P8uhIgK)n=DSO2Th30%d~qW zY(G8cI=Z-aM8M1X^0td<$B`qLEp*CRM<>)>3mYc@&T6Ko-o*!c=Ga%jP3nb8SmNm>3(&LS}?#??54(^MMsP_QGXXaaN6i3UA}XG%X(H=M)X=yz5#!Neh6 z{N#*{pqsQ$U26JbIX*B|uUc+bFwJ?eu*CMnPT~Ggv{!|RzX#0q)fZrGGgJ8_&YZT~ zLU=*Qw0mId_(7*Yvwhn29ghB7{0h)w7?F*con&u=ds>kV=>A72rzCtSKY%l zxF{aE+`|0RLBe<^rjQp1S@qXQ6gq4BmTZv?I$wajcT7QnqvRw=Na}g&>+W^6fp+PoCEx)Bh;hIL`ZC{`U zKNt3wl5Y1~O~ulhQy{{;Ok)zc#H4m`3>I38vNSL(B7?m;$x#MI_WZ3?VuV26RIi7l zHJ~RKFXb^JAnYCxjANebZIyVN%tRUpdninnCW6+)Q3Ch$3hl3*4seo2@M7HC__m@A z3_OIsz5vUql3@4>1}(GXXc$|0wJbW`<}F!99yJ$bM{`95-sXeFJ?<Vr|KW?q)2?}^8rTk%c0t__Nc)cUMX@{ZKB`?K{}4^Tjd0hIRwM8kNZ)ojn*x=G z11O@Y@7cUjnD0mS&FCOsf8om7bOes z72UqeH@UW=b`@>b{y<+nfbS1LR?%IU!VBhEsmfE)g!3%4h=PW}T1$Kj{GcXYzK zXY%#vxVw&7eX+ezh0paCU7QDES(vK&k+zo{&3!fQ!yU;{9W<1Z8J0f$Ov~WsJYBmv&A7G{0Qv*VWpj?{e1F z%N_XKSyH!rD}BLwD@JZJuhWOdTgro;a`tcLW9BX_tcSySJ{Qm^?5=JsJK80LzJtG3 z?8d>BB8mv)$j@0mZsy0qr$}e++@ZsVCsz1D_!~HZSHuc;e*ScQ1Mi?#Ex_5>1(7({ zprK7VbG<+>8*}UAM701BGxDsa$rKEAi+LNPP3U&R|1re>8G(OM3(ECog*&)qID<&E z?Dr8tzpCB^jbh^z{V{K{4?KH#=>uB@dg5rp4JqQo$O7Ye0g>&;8v$#4qN`%}??=6q z_8XmrXGUM{hwNwBN`Pptq3Fr|&8<6lswWhN6EDB2`xw7%dD!s}B&Gnl@adYa|2QRb z#doP5G9e*{{$+R_XuQ0YGG?#83Z9vcyqSpV@HUgBHSCQy+p*EyJd`=_ZukCpBc$Zv z%Y6UOwg?21jx3gvN>u1IE_Vv{Lj=o)t_2HPWIOZqoFmnXZ9wFAgaTE<BxGfxPxCDL3e}i9COvfEmIbaaXM;iC)0;KC-CF`3Y`$B{Et?+2U6Fhs)bpu!trbcwh1xm7T1h+D7AMj8vA=~CJZ+X9ip=ZXfbN2{6ImPO`x z>ZLvYDEfAxlfbg=-Fdbif1Z}W9x1%vt+Nk>M!Il^CLKvr_UNG9P-O-BWmGwGrz!V)(01Ezn1_U|KkN>$X#N`7d8&$9F3C zVph1`B1~=Y7s<=0$XH)+&>KM#!|I*4tSndFGTsiZEKPRo^J3grj|VQnN)8&F_-DtAZ3_*+Oj^j2TMPJkebM2A#D zAipxyv9oyL@wDX=-`g-geKF?PdScp`{FWT#W0HUL%OIfD;^)rqJ@5#2<%%TH*Nf6T zq|OqU!Akzt5R10j%Wky=mr$JqjXg}>B)<#Z3}EYLXWmQHUS=23NFMQH%?iK-7SjMZ zXE)Jo%XQ!@bxR)sTqG-c4R_NXZd`I+2COW4TIU9ziC_YE>9G+$l_HIQs^(GEZj12M z&LUb*WB>a=(j1jtSS~NmpNrxb;_N8ozD~>4M<|*mE$(HJwev&Hx5ho*NUvrO*4Ecc zfi{AesW56iB%#?zRq*s3D<;srYwS_+ic@E?EmQN(jYVwh&b7s?4=dFJ{Pz>7q{hO~ z$?S3J*&!59LGfg((QfLZ^5XtqH51=GNF@-IC~Rk*B@jqG@#@n~Qu>bF4FBd9v;Dr~RZ|o0sS9S4CIVKQe;{oR(>bxt1{eiw?6NzDZ zuKp2-iO-C=<>d>4=!&8qDu}m2gty}JX|?JCvi2Zb+BM%NMUot+PWS>Z9aY_n@~!? z2nS~w{E%&XUK5^yhSp;_5NhU(qWKp(aC=();^6Zv5f5!N2&*1lpe|yrvYhM<{PN(MtADT;H5OL& zE=Xv*3<38x!|2OIz0m0AV#*q~Qv=)k6ngHJ#Qgt*Kb)*A|5xfeF%}0BoP&iewSEwS z7TDC$O5K;l^vg9cx(j?~!D>#W3kU;|fx#s)rMii8glew7!nzHmm~p*7ZmrP8A*-C)aI@VGoQro*3X3hhJGFq6{?tr_{ z%C6eAO6v(qwSj9&S0ufbNBb@&6Zh+il}^HGNIEJImN@k|OwgCJO77ZP~`0G3zw{#U)Fnm$I@ z??_^|ud7LI;RBS3ZftdemRV=DVL-7bbqyoT?HLtAW?Iy=juA_*);@*1GS-sM`s>h5 zt*cT!%y^2EZ^g<@zyXg+5RIe2R}18AE%TMP7w-3KA?@;EgT&44xgO*_J&%I8(H_Uw z{6YG&z?YM(^6ij?gZ{^PNH7>{BDXbRZq+^+G^x zbJF#Q@ToPY%kMTd>5dk^VMFt5BJ)Q>^K$#(4yY^p2{OBfBkswR=}mxN%#iFk!_1*@ zRh7#?J(>9#chn8@5X7sx7QcXWyVZ zc(96)bY|q*?~gNiS>Eq+3t&4>{Z+q!_47=}Pb&NAp1I-@IrZ<#%colV;ROwBk|eC& z71z6%&K8s&bT3;jFN+P#O?Dc6?MRxEFoGbu`cDHkaWHQFgk7~)RHaZvI25&(EKBzh zEqzi~?&IGQ9i~k=-BsE|0tDPXi*UvA)nUN*8#b&UZT}mf7aJO z_l`|9Q5bmfiJ{#}v?1ywiQ4`wVY#8P&c7hvox8I6@y^%5g9d3V-$kIHCF=6n_xofT z;+|OfzF%HZG7H4qJxe}5;9eoY{yp$~Kb<=i1|PG*+Z&J0G1=@;SfW!wX+xxUVfV-0 z8;47O_!qKBu_lTKE&x(j)-JKX=rd(-%@yC|zr{-qXN6Ui?pmbxn_!ylCEvkU4*$>` z{9aTldZ#@8XVz$e-eXIT1^%FBa}sDPisiJ#HE{{CtW-}b{efir>tXB$*3V1L9{(#z zY5`OtQ8-`?2@LXuYBN}ULQF69;(D1|*L9_)L%T8ozuNsB@gD$&c?1@P&ocxHO2Op5 z1HmtZ-IV)hcnDW}$g^4RpIIMR@hp0+*F9vh9wSN~RCVCwu|bRO?Q9YvHr59i2TCsy z>L?zWaqPU*KFMq<#4B$R4PKya9j<3PU)nLT!AqrbbZfHwSw9 zU|awP>?&%z9wu;pwtX;8DtM9>%*Z97qp_T&pDmLT=DyT$P*D}3?AC0fJ#GBEy!Egn zyzt}6x|J!D*dFUxLUs5}+N#ujfTd{H0SS*!yTC5mh|r0^{3fBW`@|DM1c~{{_4g8V zULY?PA1h2AvS`NdiQ6&~EC#$Z`0{s53v~tuGrh|*$N``!($lH!@vsI8cKzn~mt=k9 zl%Jj(8Ei0h+kzH4xSF5Arzi5=LGqy?-aXZzAr_iRQe-^2DFRhsE(UaeX1LKpoz3_P zo*j5gX#3J%cD#V*EVy<#?o!&q-IX}qrKJR8DH@T}YDcyr0w+P(ejW&pJ4?Ojew+?Q z?jjDRNg=?lO|mb<(3|8-BA&IwrbhT+`J7PkRsdhS4A8|Cd2OKx543IYigQuCBH0_{i=g6XXPWWL@=PS!le{FqjU_-S%O9f@#9wK*`s1 zD@wB=cu{i%<$G@I$kOwNKEX19)swp2RuU1d&}zz6q?9;kz%bF~Gbn~0Q?oaPGVXG8tq z!2s5V;CPXD&{Kh8`4ZU(`3*Qjt!ENY8Yw6{r8k*0dS-wQx+8l6g zXd5}~=fL^VqtB=Q43>(M_F@Y$rD6L`lkv7p(2`Kr=K8+wiNNCq3&mv#kr_9ZmA{(b zu&Awd7m07YC}AOCD6Fnr&#U3^D*F+?-T<$>v#KWt+RvDI$l{-(U4^OSldOvY3=E+U zkE)`Uzm7=h04SK~iQu16alWz%O_M%#Ga_gVR^m136sak2>zO=@h;Teklwe9S2_({& z5HgWzBX-R6(AhZPqFdJ;a4ix9e+_JjV9TVyNDE=rGWqZm+gz1R62H<@AxfD#SrGT~pvk1uK>heQlsFP2ua*=H1Hw%o zY!t^|>T@ohfGt*;iXtRN@?lO+6F?G3o!au?BR5kqbIV$qcN~%PU346fgSdp-B0wB9 z^t+}061`#8CFur%VN&Gw4o6qcGccr`sDT9q9(O?a52PrAcvfnXpWKW>tUgJT*a2#4 z&O)MG?1~%b_R|MqkdxTK=A&U!WX=TXlk(tOD`WHn8@fI*;6OBR#Dt$l-h`)0qdm_LU#KJ8TwL5`>G>a&3K>Cz7I`v%5wOnUsk}wWE*&pzP12l9 z8zYTbCk>u5MH!ejxlkkFnMw*F*4wUsx zLBJ%H-C0ij13PWNq59)<)GvRU#r3P4_+cH%@jXI}4pPJ0KS5i!9Z=%lOPmJBAe4jA zw-Nvv1h2>R!^n2m5j=I8Na@e;zwpLBtS72f)kg~9>z8RM_lF&pc+ykMNW0Of)zqmLFC}Qi*t++KS**=zE=5>9Yn5bM5OlTA=hj+t z=?SRjNC;MyCXF7xG)pFcE5zrp{r>AxOiv7Tb2vHj4`M3ZJirF6gx_rcT!I0sKKPn$ zj%V}L@SaWo^|p{+2&?yJ*YCEN-*>nX<#}{SDS22_qHmJT>a^X5=^#f)*1H#Iznbsi=kLztYx;)pnsB@jCWh zjH;T}FCr=LvVqu2aI|Wv<9tYi)da9=HkjP8sz3g=6?TDZAz)}Eu38>&wNhol0qxxB z?j9UD0clYyd53x2-~~uAN1C*f1(FDK0i=92f;}=(X6z1q>So#PJ3$7Cyp}!6=`hwD zL8t;$6SpwLQH9a2Kmysw!V-KLR!8;|l}JpeF(>{nK0qR-E0zHHfJTTkY=jtoRb1#s zwgfUBgUJ<%gEZ1ZXog5bP`&7lEuF6>_iBzq2+Ct#32C0I46>gwTgeUf3mTiSNd8J> z{J@K>x0uGvRuMu*f|ufX=M`b}zuB!F)yBnpDA(vtj zE+WYbtZ}br2iJo|)|>a5Vg+e$!LJ9wfQZa}BZ|xksVk$Ey?X3l-2T~| zDiEt6hTjs6@kqiT$JI^wi~2SHLorIReqaNQ90BkJ3D%9#+k`0kA;m1A{~D)APcN>`G8dNxrcqYEpB-l%74TE{h8SU^+_^}SfOaI&pP!T; z(O&6kL`JP{XX!&4{bEi4FX%QE9~$UkK4&{6Cfo$D-661rnFnk1l4IXtEE&>d6~yvL z$pGa(XoON*oIzn#Ea7?KgC?BCcP(okSk&aW#M|M-BTxzKbn0vko~;6DC@9dEs_SAM zfTZIddj69qi+gMZTirDeSvxG=C3iT92?wDy`)rweI79WEj*-pmp8u0N;q?3YW|BCt zlmPVYFq}uV5+d)qR`H~yk1$bxvX$D^1wcvd-}70X@7EuRtWH{Kc2!lK;H`c{b~WJn zvMW8}yX`*`73;jXG;K9*vWZ)y45}!&$}YbOml-K;vq)F+O2<*Rxl78;q61jOxL+ex z%{L~Xqn}30;FvFDxAj&ozet<;Z8Z=rFr2n4Ae3LntMa3{eOmP@t#!W6m8&XSfF2}U zEB=b=$TkafZ+vml4R1Z?%K|H;%(mplI3AN~DsE?|^^vhJAziC3B&jq-DNk@h8K^zp zlWxNDBH0Rxim-2gIx9&$-t7MdOreG*)z@A-<_^y=-#jSmQ?*C4Y@KTQtDq7=uIvV* z4(P>wV-6~QpsEmj*FQy~8^=HI0#bNRCldVk0i?tb1t>q9^pmgFg7o)MB^R#@0>&ZI z2_y<|A6j`biOf4%274V?5ctCKJYzA{lDL4~E~d!PgigX}t-p4-=-8(M6uS9;ks)VX z*D;ZspzZ~mFX70@H};$NNL_^8Yey2xCYZ3?TEAN&zaLiA8QY~KzVl1O0j z9iowpa-&G~1p1w47AhYPd1253U5*LEI?`DFfE;WQh{--tU1kc_=;@v7S+vu(xP^8( zKVoRy+fQ_!AFg28RDeG$kVXN9lD7W7hsv)jrjm7f4YTp+#QMm?aFq#1QB1(NPdf>U zip!vVn)%%rL>nQqvx#B1q`Pr*Yy0CLgjFbHC1h-~zeLS7*J3gL(%Z#UVI#)u^d+qG{M8&o?3n&UmdZ&QL@bsW(%)UZacF+hk zSIE^-8X9JT^Lk_0=W;6?Y8!?oOtKeeYl&O~U^@6*PfaKHpZsyCH+Llytu|g5hjMs*> z3u{Aoxgw5VjoUT((*V(_D${pLSup*}p?x}hjTX~q5ysT%|6=Q$!UGAmEgajnJ+ZBc zF|jAs#5Ov%GqID2ZQD*JwrxAPbMMQ2IOnDNxxcEe+O=z~ehs!7Uji}aD0e!C*@t1O}0WRg&P#tC%g`r{}>IY#E{BLlkIHJ zw5w%)t!}Dn+Yx}SKapwu)`IIM`H%Mw+jXL(7QoH|bs0*sj_J)O@;%_Hy$Kdn~FPwB7PMFPuIB(>eRcp`jDF}ez6r^WKk*@ zre_N!I7B~dnURj?y^E}f0a7aTLEGwTexd6!f!8M1JagL6FbV;bCxQwF2uQDyUx6r6 zPg$RF=bFY6_JV#6X!-MPkvpFS?J;VcjJuzMJfPdWbsmOZVkoQ7tV}=Q5OoR`Lg>ut zUEgFQKC|+V-_Nu_p+A7;o!AP4?BljD?$AI`>H39{LeZ zZ@G$$5LmP#7_S4se3pD?i%_mLze``Y$13c7 zQ3a(tXi!<^wttbLUGELm)z_eJP&w_jRNzo;xQQiJ#I_UEp)uH?oZ1@1Lpbyt^R2H=Kq;$$^F z=ud~_SxvnAjsEFZw}&iQ@1hi0*-3FUPv6Oyl_~gg08%68?@0AEFBe!~;fq2)fCpIL zvxcyInPUwFQE{2OKibq@FoIaVI0k;!!r%N!>o>0#8w8Em1z*TAUVyG<9QY~K8#Ds- z{M__vh`q*}O%>n-F z=USEZH)3)fm8{5LJv7w^-zzl6bF8I9Ln?0(5rR%|cIH^#u&7q5i~vK8pM;Nw%BzP- zGvWWPF~SL>fPLRIyucOyO0VBe6+>||W%c6teLk07-N&45+Rv6~TKyW$TGt!;Urz|~ zEqzboK8lHvR-_j>bY5Dav-;(s^O;@Y3`ngFCG8}{(0=o%&wcc3EM`|Uh5Wy;Ot^|Y4>`MN`xFaS+T*yBbVjyk5c-|=qy zA;Y9LL}J^T3b?a&xRBz?RZtG^(O~lRyw?BSD`0>pJVODlfe>{-Uo!EpqQJw&q|Yaa zn(@L_qk9&fkQ5JzsHJOh933g+t9BAD30!enZXZii{a%*`VHgE9B7R8ErGSN58ii|U zyTAw*>3_P_USb`+a46{-O#QPnD{gMQYzb92FI*Ky27KWv{huvnU3sMpF}&C*pC>Z9Hbk`bJV|&1#$kd zR}pJ+Q#%>AKh8!lAyd`cbBwUf6)dQjwpFyzR3MA^SVdCe1v~pOeKLOl!frInfoU+a zK+-Kb4jeQyeaSJjw;iKG@_rdAOcO^mt3KCi>nGr`j~^ioRq31^MpsS?C^iSyO>_tb zu;`JwD5<-X_+e4{+=2-|jy-%w*^*4EzAF$(U zVTGQ}Si+?wYD$%3J!cJRYX1vXtA=~>daQU*fAPF*8obWTB~eG?eN z7vR?5G&19YF~y7{6d%gz_8DJtU~)-%!Hg8Xx-^2mtaflBlhVIFgk+)WIMQKapr!|* zJH$+^hO6k@^dznBezQ%`->u3fNEL2z_-IT?&Q6jG@+F_YL?I|U3WF$!)55rntp|5) zcFtRcfZ%>-xl1;HMcVs_BmC~zERwMT$n<$!*rkaN@=RrX}W@O9 z!su?j3QS}c+u$uYFm55g%(HB<3PM8|5~PPpxx3%Dqt#2?XdtIuue$T>z#WAlvx3@P zQy6k|pQlBwQlWOk+1{3X$*eqGpg=loGqMG*KSd-yK_ApwN#YzwS~`2{Iv1E`<)T=b zakfV1pMT8wHJy1JpV@a1uq7S30v*!mZf;8Zw0q1+(|s8j&s!!lDJDEi&D+u0A=rG_ z>aqi$BC_M7Un-r?+WR)WLX6wF{4wmn2OhQS?>UaVAphk{LXkouzK-d$?L zwa^5}X}*8;Dlf!qs|#L0xnE{2zgEA;vVu6hh~$Qf)iH$ z9DWSN#&&!0KZ*US;{CZdnWC_TouO<0!~I$w=Ew+yNcqj}W`&Am6OfU+$_fr#wC@RRrqF5AP^}!5NM;lZQ^2QUV%@PZdBgpGz_B;)^b=%*?Fn7uU6&=Wu}c4#>re85acy84`_)=;R+quWATA`d@8Yx z*={$2a=!98M|97fTF<$dly54;slLdmUwQwCsnJX(8P(vVPY^Kh5xAE z+q2u~c7f3nP$a$`-&XUhY;;}lIsE{rx11mKoHcpa@*|yZRu3n)DLy+lk(c`IG_Gfq zsTKPMDARk}EF)8xdl6gmf!+_ZkVun`p?rxQbLd3vJ+qmKz@Bv$S0sR*R@eX)pZbOA zm+ZPn>6kr<<(Q(gA4=OfI(YmP|4xeIk9Ot9YCJL7Yzp#8rMf-iWRtIJ6vCWX-{AaD zbddyx@}neX%OI5O>o@#4v0cbB11EXs)g(^q6QkSVmq&XZfVjxrVL{ny_V>?Sxw?Xw zoZ=`23gc1^(4;G5(kgb|)BxU5o&RK^48EAj-=Sc9$n_)P{fpAgWRafGb+5f%t!fVp z;k)^+evGNYX5vHqo2|S5o_@5_wPP&IjTB0c#|Q#l*SE1lg2V}m`pxd(klyV)=BS>p ze+@BO@Nb(RDTYQ$FPh$1pG^L+%Ca1cvfug%{(vUh0c$IYN4vXoi>rFiWiMMSyk+<3 z5Pd9S5DRBB#l2rR3$BY6P1fb(<-J7ATu+&Zmzo2O*`F7|SchP~YJ7h_bmUfv**HMN zgC#H^>p3Yz?NH}$6a>08q&cWr3e&g`7}Xwgi0+RVjSg*xOzH7kxltXCH56;d){$PK zqZkAn0PF~bwnn#A4l`KY2%U+7EruBHG8c>kGfnY)3G0uE4WW}&R=qTgZtAsr4J8{E z>CrTl6QjOQ5*k)GGT3Vl)k-tbLbNJsP}UAD=#Oc?laVheo{};4@PK7v>H1BsFjMOA z_*MYr4cvpH`#Z*4$HrG3#i1{klzl|Jm4J=ivTdiViY27wi91bU3;G&*?I#HEyfJg|e|#y