added api for lumped weight; some methods renaming

Author: cnpetra

src/Drivers/Dense/NlpDenseConsEx1.hpp

@@ -112,10 +112,10 @@ class DiscretizedFunction : public hiop::hiopVectorPar
 class DenseConsEx1 : public hiop::hiopInterfaceDenseConstraints
 {
 public:
-  DenseConsEx1(int n_mesh_elem, double mesh_ratio=1.0, bool use_weighted_vars=false)
+  DenseConsEx1(int n_mesh_elem, double mesh_ratio=1.0, int type_weighted_vars=0)
       : n_vars(n_mesh_elem),
         comm(MPI_COMM_WORLD),
-        use_weighted_space_(use_weighted_vars),
+        type_weighted_space_((hiopInterfaceBase::WeightedSpaceType)type_weighted_vars),
         solver_(nullptr)
   {
     // create the members
@@ -259,7 +259,7 @@ class DenseConsEx1 : public hiop::hiopInterfaceDenseConstraints
                         double alpha_pr,
                         int ls_trials);
 
-  bool applyM(const size_type& n, const double* x_in, double* y_out)
+  bool apply_M(const size_type& n, const double* x_in, double* y_out)
   {
     x->copyFrom(x_in);
     _mesh->applyM(*x);
@@ -270,7 +270,7 @@ class DenseConsEx1 : public hiop::hiopInterfaceDenseConstraints
   /**
    * Computes action y of the inner product weight matrix H on a vector x, namely y=H*x
    */
-  virtual bool applyH(const size_type& n, const double* x_in, double* y_out)
+  virtual bool apply_H(const size_type& n, const double* x_in, double* y_out)
   {
     x->copyFrom(x_in);
     _mesh->applyM(*x);
@@ -281,7 +281,7 @@ class DenseConsEx1 : public hiop::hiopInterfaceDenseConstraints
   /**
    * Computes action y of the inverse of H on a vector x, namely y=H^{-1}*x
    */
-  virtual bool applyHinv(const size_type& n, const double* x_in, double* y_out)
+  virtual bool apply_Hinv(const size_type& n, const double* x_in, double* y_out)
   {
     x->copyFrom(x_in);
     _mesh->applyMinv(*x);
@@ -290,17 +290,20 @@ class DenseConsEx1 : public hiop::hiopInterfaceDenseConstraints
   }
 
   /**
-   * Enables the use of weighted inner products via @applyM, @applyH, and @applyHinv
+   * Enables the use of weighted inner products via @apply_M, @apply_H, and @apply_Hinv
    */
-  virtual bool useWeightedVectorSpace()
+  virtual hiopInterfaceBase::WeightedSpaceType get_weighted_space_type()
   {
-    return use_weighted_space_;
+    // return 
+    // - hiopInterfaceBase::Euclidean for no weighted space 
+    // - hiopInterfaceBase::HilbertLumped for space weighted by lumped mass matrix
+    return type_weighted_space_;
   }
 
 private:
   int n_vars;
   MPI_Comm comm;
-  bool use_weighted_space_;
+  hiopInterfaceBase::WeightedSpaceType type_weighted_space_;
   Ex1Meshing1D* _mesh;
 
   int n_local;

src/Interface/hiopInterface.hpp

@@ -61,7 +61,7 @@
 
 namespace hiop
 {
-/** Solver status codes. */
+// Solver status codes
 enum hiopSolveStatus
 {
   //(partial) success
@@ -142,6 +142,24 @@ class hiopInterfaceBase
     hiopQuadratic,
     hiopNonlinear
   };
+  /**
+   * Type of weighting to be used internally for computing norm and inner products, for determining
+   * representations of dual bound variables, and initial quasi-Newton Hessian approximations.
+   */
+  enum WeightedSpaceType
+  {
+    // no weighted space.
+    Euclidean = 0,
+    
+    // lumps mass matrix and uses lumped matrix everywhere
+    HilbertLumped,
+    
+    /**
+     * Experimental: lumps mass matrix and uses for dual bound representers and initial Hessian
+     * approximations, however, uses H and H-inverse norms in termination criteria
+     */
+    Hilbert
+  };
 
 public:
   hiopInterfaceBase() {};
@@ -288,17 +306,17 @@ class hiopInterfaceBase
    * mesh independent performance of the algorithm. Currently M is lumped into a diagonal
    * and used internally for the abovementioned dual variables.
    *
-   * @note If the variables are not coming from discretizations (specified by a false 
-   * value returned by @useWeightedVectorSpace), the methods should return false. HiOp 
-   * will use internally M=I. Otherwise, should return true or false depending on whether 
-   * the user code succesfully applied M. 
+   * @note If the variables are not coming from discretizations (specified by 
+   * a hiopInterfaceBase::Euclidean return value from @get_weighted_space_type), the method is
+   * not called and HiOp will use internally M=I. Otherwise, should true or false depending 
+   * on whether the user code succesfully applied M. 
    *
    * @param[in] n the (global) number of variables
    * @param[in] x the array to which mass action is applied
    * @param[out] y the result of apply
    * @return see note above.
    */
-  virtual bool applyM(const size_type& n, const double* x, double* y)
+  virtual bool apply_M(const size_type& n, const double* x, double* y)
   {
     return true;
   }
@@ -316,10 +334,9 @@ class hiopInterfaceBase
    * interior-point method for PDE-constrained optimization using finite element 
    * discretizations, Optimiz. Meth. and Software, Vol. 38, 2023.
    *
-   * @note If the variables are not coming from discretizations (specified by a false 
-   * value returned by @useWeightedVectorSpace), the methods should return false. HiOp 
-   * will use internally H=I. Otherwise, should return true or false depending on whether 
-   * the user code succesfully applied H.
+   * @note The method is called only when hiopInterfaceBase::Hilbert is returned by 
+   * @get_weighted_space_type. In this case, it should return true or false depending on 
+   * whether the user code succesfully applied H.
    * 
    * @note Currently HiOp only uses H for computing the inner products and norms (including
    * for vector representations of duals discretizations). The lumped mass matrix M is used
@@ -331,39 +348,46 @@ class hiopInterfaceBase
    * @param[out] y the result of apply
    * @return see note above
    */
-  virtual bool applyH(const size_type& n, const double* x, double* y)
+  virtual bool apply_H(const size_type& n, const double* x, double* y)
   {
     return true;
   }
 
   /**
    * Computes action y of the inverse of H on a vector x, namely y=H^{-1}*x
    *
-   * See @applyH for a discussion of H and additional notes. The inverse of H plays an
+   * See @apply_H for a discussion of H and additional notes. The inverse of H plays an
    * important role in computing the "dual" norms and, in turn, to ensure mesh
-   * independent behavior of the IPM solver(s). Also see notes from @applyH.
+   * independent behavior of the IPM solver(s). Also see notes from @apply_H.
    *
    * @param[in] n the (global) number of variables
    * @param[in] x the array to which the inverse is applied
    * @param[out] y the result of apply
-   * @return see return of @applyH
+   * @return see return of @apply_H
    */
-  virtual bool applyHinv(const size_type& n, const double* x, double* y)
+  virtual bool apply_Hinv(const size_type& n, const double* x, double* y)
   {
     return true;
   }
 
   /**
-   * Enables the use of weighted inner products via @applyM, @applyH, and @applyHinv
+   * Enables the use of weighted inner products via @apply_M, @apply_H, and @apply_Hinv
+   * 
+   * See @WeightedSpaceType for the return values of this function and their meaning. If the
+   * return value is
+   *  - hiopInterfaceBase::Euclidean: @apply_M, @apply_H, and @apply_Hinv need not be overriden by user.
+   * If provided, they are not called.
+   *  - hiopInterfaceBase::HilbertLumped: only @apply_M is called. @apply_H and @apply_Hinv are not.
+   *  - hiopInterfaceBase::Hilbert: @apply_M, @apply_H, and @apply_Hinv are called.
    *
-   * See also notes for @applyM, @applyH, and @applyHinv.
+   * See also notes for @apply_M, @apply_H, and @apply_Hinv.
    *
-   * @return true if enabled, false to default to Euclidean space with <u,v>=u^T*v
+   * @return WeightedSpaceType (see @WeightedSpaceType and notes above
    */
-  virtual bool useWeightedVectorSpace()
+  virtual WeightedSpaceType get_weighted_space_type()
   {
     //the default impl. instructs HiOp to use Euclidean/l^2 variables (H=I) internally
-    return false;
+    return Euclidean;
   }
 
   /**

src/Optimization/VectorSpace.cpp

@@ -63,7 +63,7 @@ VectorSpace::VectorSpace(hiopNlpFormulation* nlp)
 {
   assert(nlp);
   M_lump_ = nullptr;
-  if(nlp->useWeightedVectorSpace()) {
+  if(nlp->get_weighted_space_type()) {
     vec_n_ = nlp_->alloc_primal_vec();
     vec_n2_ = nlp_->alloc_primal_vec();
   } else {
@@ -81,7 +81,7 @@ VectorSpace::~VectorSpace()
 
 //bool VectorSpace::apply_M(const hiopVector& x, hiopVector& y) const
 //{
-//  if(nlp_->useWeightedVectorSpace()) {
+//  if(nlp_->get_weighted_space_type()) {
 //    return nlp_->eval_M(x, y);
 //  } else {
 //    y.copyFrom(x);
@@ -100,7 +100,7 @@ bool VectorSpace::apply_M_lumped(const hiopVector& x, hiopVector& y) const
 // Computes H primal norm
 double VectorSpace::norm_H_primal(const hiopVector& x) const
 {
-  if(nlp_->useWeightedVectorSpace()) {      
+  if(nlp_->get_weighted_space_type()) {      
     nlp_->eval_H(x, *vec_n_);
     auto dp = x.dotProductWith(*vec_n_);
     return ::std::sqrt(dp);
@@ -111,7 +111,7 @@ double VectorSpace::norm_H_primal(const hiopVector& x) const
 // Computes H dual norm
 double VectorSpace::norm_H_dual(const hiopVector& x) const
 {
-  if(nlp_->useWeightedVectorSpace()) {      
+  if(nlp_->get_weighted_space_type()) {      
     nlp_->eval_H_inv(x, *vec_n_);
     auto dp = x.dotProductWith(*vec_n_);
     return ::std::sqrt(dp);
@@ -122,7 +122,7 @@ double VectorSpace::norm_H_dual(const hiopVector& x) const
 
 double VectorSpace::norm_stationarity(const hiopVector& x) const
 {
-  if(nlp_->useWeightedVectorSpace()) {
+  if(nlp_->get_weighted_space_type()) {
     vec_n_->copyFrom(x);
     vec_n_->componentDiv(*M_lumped());
     return vec_n_->infnorm();
@@ -134,7 +134,7 @@ double VectorSpace::norm_stationarity(const hiopVector& x) const
 // Compute norm one weighted by M, i.e., 1^T*M*|x|
 double VectorSpace::norm_M_one(const hiopVector& x) const
 {
-  if(nlp_->useWeightedVectorSpace()) {
+  if(nlp_->get_weighted_space_type()) {
     //use vec_n2_ since vec_n_ may be changed in M_lumped_
     vec_n2_->copyFrom(x);
     vec_n2_->component_abs();
@@ -147,7 +147,7 @@ double VectorSpace::norm_M_one(const hiopVector& x) const
 
 double VectorSpace::norm_complementarity(const hiopVector& x) const
 {
-  if(nlp_->useWeightedVectorSpace()) {
+  if(nlp_->get_weighted_space_type()) {
     // since both x (slacks) and the bound duals are in the same (primal) space, inf norm "is
     // mesh independent".
     return x.infnorm();
@@ -159,7 +159,7 @@ double VectorSpace::norm_complementarity(const hiopVector& x) const
 // Computes the "volume" of the space, 1^T M*1 
 double VectorSpace::volume() const
 {
-  if(nlp_->useWeightedVectorSpace()) {
+  if(nlp_->get_weighted_space_type()) {
     double vol_total = nlp_->m_ineq_low() + nlp_->m_ineq_upp();
     if(nlp_->n_low() > 0 || nlp_->n_upp() > 0) {
       //compute ||1||_M
@@ -184,7 +184,7 @@ const hiopVector* VectorSpace::M_lumped() const
 {
   if(M_lump_ == nullptr) {
     M_lump_ = nlp_->alloc_primal_vec();
-    if(nlp_->useWeightedVectorSpace()) {    
+    if(nlp_->get_weighted_space_type()) {    
       vec_n_->setToConstant(1.);
       nlp_->eval_M(*vec_n_, *M_lump_);
     } else {
@@ -197,7 +197,7 @@ const hiopVector* VectorSpace::M_lumped() const
 void VectorSpace::
 add_linear_damping_term(const hiopVector& ixl, const hiopVector& ixu, const double& ct, hiopVector& x) const
 {
-  if(nlp_->useWeightedVectorSpace()) {
+  if(nlp_->get_weighted_space_type()) {
     vec_n_->copyFrom(ixl);
     vec_n_->axpy(-1.0, ixu);
     vec_n_->componentMult(*M_lumped());
@@ -209,7 +209,7 @@ add_linear_damping_term(const hiopVector& ixl, const hiopVector& ixu, const doub
 
 double VectorSpace::log_barrier_eval_local(const hiopVector& x, const hiopVector& ix) const
 {
-  if(nlp_->useWeightedVectorSpace()) {
+  if(nlp_->get_weighted_space_type()) {
     return x.logBarrierWeighted_local(ix, *M_lumped());
   } else {
     return x.logBarrier_local(ix);
@@ -219,7 +219,7 @@ double VectorSpace::log_barrier_eval_local(const hiopVector& x, const hiopVector
 // Adds (to `gradx`) the gradient of the weighted log, namely gradx = gradx - mu * M_lumped * ix/s
 void VectorSpace::log_barrier_grad_add(const double& mu, const hiopVector& s, const hiopVector& ix, hiopVector& gradx) const
 {
-  if(nlp_->useWeightedVectorSpace()) {
+  if(nlp_->get_weighted_space_type()) {
     vec_n_->copyFrom(ix);
     vec_n_->componentDiv(s);
     vec_n_->componentMult(*M_lumped());
@@ -235,7 +235,7 @@ double VectorSpace::linear_damping_term_local(const hiopVector& s,
                                               const double& mu,
                                               const double& kappa_d) const
 {
-  if(nlp_->useWeightedVectorSpace()) {
+  if(nlp_->get_weighted_space_type()) {
     vec_n_->copyFrom(s);
     vec_n_->componentMult(*M_lumped());
     return vec_n_->linearDampingTerm_local(ixl, ixr, mu, kappa_d);

src/Optimization/hiopNlpFormulation.cpp

@@ -779,9 +779,9 @@ bool hiopNlpFormulation::eval_M(const hiopVector& x, hiopVector& y)
          "weighted vector spaces not supported when fixed variables are removed");
 
   runStats.tm_eval_M_apply.start();
-  bool bret = interface_base.applyM(nlp_transformations_.n_pre(),
-                                    x_full->local_data_const(),
-                                    y_full->local_data());
+  bool bret = interface_base.apply_M(nlp_transformations_.n_pre(),
+                                     x_full->local_data_const(),
+                                     y_full->local_data());
   runStats.tm_eval_M_apply.stop();
   runStats.n_eval_M_apply++;
 
@@ -802,9 +802,9 @@ bool hiopNlpFormulation::eval_H(const hiopVector& x, hiopVector& y)
          "weighted vector spaces not supported when fixed variables are removed");
 
   runStats.tm_eval_H_apply.start();
-  bool bret = interface_base.applyH(nlp_transformations_.n_pre(),
-                                    x_full->local_data_const(),
-                                    y_full->local_data());
+  bool bret = interface_base.apply_H(nlp_transformations_.n_pre(),
+                                     x_full->local_data_const(),
+                                     y_full->local_data());
   runStats.tm_eval_H_apply.stop();
   runStats.n_eval_H_apply++;
 
@@ -823,9 +823,9 @@ bool hiopNlpFormulation::eval_H_inv(const hiopVector& x, hiopVector& y)
          "weighted vector spaces not supported when fixed variables are removed");
 
   runStats.tm_eval_Hinv_apply.start();
-  bool bret = interface_base.applyHinv(nlp_transformations_.n_pre(),
-                                       x_full->local_data_const(),
-                                       y_full->local_data());
+  bool bret = interface_base.apply_Hinv(nlp_transformations_.n_pre(),
+                                        x_full->local_data_const(),
+                                        y_full->local_data());
   runStats.tm_eval_Hinv_apply.stop();
   runStats.n_eval_Hinv_apply++;
 

src/Optimization/hiopNlpFormulation.hpp

@@ -118,9 +118,9 @@ class hiopNlpFormulation
   virtual bool eval_Jac_d(hiopVector& x, bool new_x, hiopMatrix& Jac_d) = 0;
   virtual bool eval_Jac_c_d(hiopVector& x, bool new_x, hiopMatrix& Jac_c, hiopMatrix& Jac_d);
 
-  inline bool useWeightedVectorSpace()
+  inline bool get_weighted_space_type()
   {
-    return interface_base.useWeightedVectorSpace();
+    return interface_base.get_weighted_space_type();
   }
   /* Wrapper over user-defined M apply that also applies the NLP transformations. */
   bool eval_M(const hiopVector& x, hiopVector& y);