Update IPM NLP solvers to be aware of weighted spaces

src/Drivers/Dense/CMakeLists.txt

@@ -32,6 +32,7 @@ install(
 ##########################################################
 add_test(NAME NlpDenseCons1_5H  COMMAND ${RUNCMD} "$<TARGET_FILE:NlpDenseConsEx1.exe>"  "500" "1.0" "-selfcheck")
 add_test(NAME NlpDenseCons1_5K  COMMAND ${RUNCMD} "$<TARGET_FILE:NlpDenseConsEx1.exe>" "5000" "1.0" "-selfcheck")
+add_test(NAME NlpDenseCons1_25K_InfDim COMMAND ${RUNCMD} "$<TARGET_FILE:NlpDenseConsEx1.exe>" "25000" "1.0" "-selfcheck" "-use_L2")
 add_test(NAME NlpDenseCons1_50K COMMAND ${RUNCMD} "$<TARGET_FILE:NlpDenseConsEx1.exe>" "50000" "1.0" "-selfcheck")
 if(HIOP_USE_MPI)
   add_test(NAME NlpDenseCons1_50K_mpi COMMAND ${MPICMD} -n 2 "$<TARGET_FILE:NlpDenseConsEx1.exe>" "50000" "1.0" "-selfcheck")

src/Drivers/Dense/NlpDenseConsEx1.cpp

@@ -76,7 +76,15 @@ bool Ex1Meshing1D::get_vecdistrib_info(size_type global_n, index_type* cols)
   for(int i = 0; i <= comm_size; i++) cols[i] = col_partition[i];
   return true;
 }
-void Ex1Meshing1D::applyM(DiscretizedFunction& f) { f.componentMult(*this->_mass); }
+void Ex1Meshing1D::applyM(DiscretizedFunction& f)
+{
+  f.componentMult(*this->_mass);
+}
+
+void Ex1Meshing1D::applyMinv(DiscretizedFunction& f)
+{
+  f.componentDiv(*this->_mass);
+}
 
 // converts the local indexes to global indexes
 index_type Ex1Meshing1D::getGlobalIndex(index_type i_local) const

src/Drivers/Dense/NlpDenseConsEx1.hpp

@@ -68,8 +68,8 @@ class Ex1Meshing1D
   index_type* get_col_partition() const { return col_partition; }
   MPI_Comm get_comm() const { return comm; }
 
-  virtual void applyM(DiscretizedFunction& f);
-
+  void applyM(DiscretizedFunction& f);
+  void applyMinv(DiscretizedFunction& f);
 protected:
   hiop::hiopVector* _mass;  // the length or the mass of the elements
   double _a, _b;            // end points
@@ -112,9 +112,10 @@ class DiscretizedFunction : public hiop::hiopVectorPar
 class DenseConsEx1 : public hiop::hiopInterfaceDenseConstraints
 {
 public:
-  DenseConsEx1(int n_mesh_elem = 100, double mesh_ratio = 1.0)
+  DenseConsEx1(int n_mesh_elem, double mesh_ratio=1.0, int type_weighted_vars=0)
       : n_vars(n_mesh_elem),
         comm(MPI_COMM_WORLD),
+        type_weighted_space_((hiopInterfaceBase::WeightedSpaceType)type_weighted_vars),
         solver_(nullptr)
   {
     // create the members
@@ -257,10 +258,52 @@ class DenseConsEx1 : public hiop::hiopInterfaceDenseConstraints
                         double alpha_du,
                         double alpha_pr,
                         int ls_trials);
+
+  bool apply_M(const size_type& n, const double* x_in, double* y_out)
+  {
+    x->copyFrom(x_in);
+    _mesh->applyM(*x);
+    x->copyTo(y_out);
+    return true;
+  }
+
+  /**
+   * Computes action y of the inner product weight matrix H on a vector x, namely y=H*x
+   */
+  virtual bool apply_H(const size_type& n, const double* x_in, double* y_out)
+  {
+    x->copyFrom(x_in);
+    _mesh->applyM(*x);
+    x->copyTo(y_out);
+    return true;
+  }
+
+  /**
+   * Computes action y of the inverse of H on a vector x, namely y=H^{-1}*x
+   */
+  virtual bool apply_Hinv(const size_type& n, const double* x_in, double* y_out)
+  {
+    x->copyFrom(x_in);
+    _mesh->applyMinv(*x);
+    x->copyTo(y_out);
+    return true;
+  }
+
+  /**
+   * Enables the use of weighted inner products via @apply_M, @apply_H, and @apply_Hinv
+   */
+  virtual hiopInterfaceBase::WeightedSpaceType get_weighted_space_type()
+  {
+    // 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;
+  hiopInterfaceBase::WeightedSpaceType type_weighted_space_;
   Ex1Meshing1D* _mesh;
 
   int n_local;

src/Drivers/Dense/NlpDenseConsEx1Driver.cpp

@@ -21,17 +21,26 @@ static bool self_check(size_type n, double obj_value);
 #ifdef HIOP_USE_AXOM
 static bool do_load_checkpoint_test(const size_type& mesh_size, const double& ratio, const double& obj_val_expected);
 #endif
-static bool parse_arguments(int argc, char** argv, size_type& n, double& distortion_ratio, bool& self_check)
+static bool parse_arguments(int argc, char** argv, size_type& n, double& distortion_ratio, bool& wei_vars, bool& self_check)
 {
   n = 20000;
   distortion_ratio = 1.;
+  wei_vars = false;
   self_check = false;  // default options
 
   switch(argc) {
     case 1:
       // no arguments
       return true;
       break;
+    case 5: // 4 arguments - -use_L2 expected
+    {
+      if(std::string(argv[4]) == "-use_L2") {
+        wei_vars = true;
+      } else {
+        return false;
+      }
+    }
     case 4:  // 3 arguments - -selfcheck expected
     {
       if(std::string(argv[3]) == "-selfcheck") {
@@ -91,15 +100,17 @@ int main(int argc, char** argv)
   }
 #endif
   bool selfCheck;
+  //flags the use of Ex1 with the mass weighted variables (L2 inner product) or without (Euclidean inner product) 
+  bool wei_vars;
   size_type mesh_size;
   double ratio;
   double objective = 0.;
-  if(!parse_arguments(argc, argv, mesh_size, ratio, selfCheck)) {
+  if(!parse_arguments(argc, argv, mesh_size, ratio, wei_vars, selfCheck)) {
     usage(argv[0]);
     return 1;
   }
 
-  DenseConsEx1 problem(mesh_size, ratio);
+  DenseConsEx1 problem(mesh_size, ratio, wei_vars);
   hiop::hiopNlpDenseConstraints nlp(problem);
 
   hiop::hiopAlgFilterIPM solver(&nlp);
@@ -110,7 +121,9 @@ int main(int argc, char** argv)
 
   // this is used for testing when the driver is called with -selfcheck
   if(selfCheck) {
-    if(!self_check(mesh_size, objective)) return -1;
+    if(!self_check(mesh_size, objective)) {
+      return -1;
+    }
   } else {
     if(rank == 0) {
       printf("Optimal objective: %22.14e. Solver status: %d. Number of iterations: %d\n",
@@ -142,9 +155,9 @@ int main(int argc, char** argv)
 
 static bool self_check(size_type n, double objval)
 {
-#define num_n_saved 3  // keep this is sync with n_saved and objval_saved
-  const size_type n_saved[] = {500, 5000, 50000};
-  const double objval_saved[] = {8.6156700e-2, 8.6156106e-02, 8.6161001e-02};
+#define num_n_saved 4  // keep this is sync with n_saved and objval_saved
+  const size_type n_saved[] = {500, 5000, 50000, 25000};
+  const double objval_saved[] = {8.6156700e-2, 8.6156106e-02, 8.6161001e-02, 8.6155777e-02};
 
 #define relerr 1e-6
   bool found = false;

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() {};
@@ -282,6 +300,99 @@ class hiopInterfaceBase
     return false;  // defaults to serial
   }
 
+  /**
+   * Computes action y of the mass matrix on a vector x, namely y=M*x
+   *
+   * The mass matrix is generally the weight matrix associated with L^2 finite element 
+   * discretizations. This matrix is used internally to build and maintain Riesz-like 
+   * representations of the dual variables (associated with bound constraints) to ensure
+   * 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 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 apply_M(const size_type& n, const double* x, double* y)
+  {
+    return true;
+  }
+
+  /**
+   * Computes action y of the inner product weight matrix H on a vector x, namely y=H*x
+   *
+   * The weight matrix H is generally associated with L^2 or H^1 finite element 
+   * discretizations and corresponding weighted inner products: <u_h,v_h> = u^T H v. For L^2, 
+   * H is the mass matrix, while for H^1 is the mass plus stiffness. This matrix is applied 
+   * internally to obtain termination criteria and derive Hessian representations that are 
+   * mesh independent and consistent with the function space nature of the problem.
+   *
+   * Additional info: C. G. Petra et. al., On the implementation of a quasi-Newton 
+   * interior-point method for PDE-constrained optimization using finite element 
+   * discretizations, Optimiz. Meth. and Software, Vol. 38, 2023.
+   *
+   * @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
+   * in the quasi-Newton Hessian approximations. This can be revisited, but will require
+   * the user to provide a method for solve with H plus a diagonal positive definite matrix.
+   *
+   * @param[in] n the (global) number of variables
+   * @param[in] x the array to which the H matrix is applied
+   * @param[out] y the result of apply
+   * @return see note above
+   */
+  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 @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 @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 @apply_H
+   */
+  virtual bool apply_Hinv(const size_type& n, const double* x, double* y)
+  {
+    return true;
+  }
+
+  /**
+   * 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 @apply_M, @apply_H, and @apply_Hinv.
+   *
+   * @return WeightedSpaceType (see @WeightedSpaceType and notes above
+   */
+  virtual WeightedSpaceType get_weighted_space_type()
+  {
+    //the default impl. instructs HiOp to use Euclidean/l^2 variables (H=I) internally
+    return Euclidean;
+  }
+
   /**
    * Method provides a primal or starting point. This point is subject to internal adjustments.
    *

src/LinAlg/VectorCudaKernels.cu

@@ -1042,6 +1042,35 @@ double log_barr_obj_kernel(int n, double* d1, const double* id)
   return sum;
 }
 
+/** @brief compute sum(w[i]*log(d1[i])) forall i where id[i]=1*/
+double log_barr_wei_obj_kernel(int n, double* d1, const double* id, const double* w)
+{
+  thrust::device_ptr<double> dev_v = thrust::device_pointer_cast(d1);
+  thrust::device_ptr<const double> id_v = thrust::device_pointer_cast(id);
+  thrust::device_ptr<const double> wei_v = thrust::device_pointer_cast(w);
+
+  // wrap raw pointer with a device_ptr 
+  thrust_log_select<double> log_select_op;
+  thrust::plus<double> plus_op;
+  thrust::multiplies<double> mult_op;
+
+  // TODO: how to avoid this temp vec?
+  thrust::device_ptr<double> v_temp = thrust::device_malloc(n*sizeof(double));
+
+  // compute v=x*id
+  thrust::transform(thrust::device, dev_v, dev_v+n, id_v, v_temp, mult_op);
+  // compute v[i] = ( log(v[i]) if v[i]>0, otherwise 0 )
+  thrust::transform(thrust::device, v_temp, v_temp+n, v_temp, log_select_op); 
+  // compute v[i] = w[i]*v[i] 
+  thrust::transform(thrust::device, v_temp, v_temp+n, wei_v, v_temp, mult_op); 
+  // sum up
+  const double sum = thrust::reduce(thrust::device, v_temp, v_temp+n, 0.0, plus_op);
+
+  thrust::device_free(v_temp);
+
+  return sum;
+}
+
 /** @brief compute sum(d1[i]) */
 double thrust_sum_kernel(int n, double* d1)
 {

src/LinAlg/VectorCudaKernels.hpp

@@ -187,6 +187,8 @@ void thrust_component_sqrt_kernel(int n, double* d1);
 void thrust_negate_kernel(int n, double* d1);
 /** @brief compute sum(log(d1[i])) forall i where id[i]=1*/
 double log_barr_obj_kernel(int n, double* d1, const double* id);
+/** @brief compute sum(w[i]*log(d1[i])) forall i where id[i]=1*/
+double log_barr_wei_obj_kernel(int n, double* d1, const double* id, const double* w);
 /** @brief compute sum(d1[i]) */
 double thrust_sum_kernel(int n, double* d1);
 /** @brief Linear damping term */

src/LinAlg/VectorHipKernels.cpp

@@ -937,6 +937,36 @@ double log_barr_obj_kernel(int n, double* d1, const double* id)
   return sum;
 }
 
+/** @brief compute sum(w[i]*log(d1[i])) forall i where id[i]=1*/
+double log_barr_wei_obj_kernel(int n, double* d1, const double* id, const double* w)
+{
+  thrust::device_ptr<double> dev_v = thrust::device_pointer_cast(d1);
+  thrust::device_ptr<const double> id_v = thrust::device_pointer_cast(id);
+  thrust::device_ptr<const double> wei_v = thrust::device_pointer_cast(w);
+
+  // wrap raw pointer with a device_ptr 
+  thrust_log_select<double> log_select_op;
+  thrust::plus<double> plus_op;
+  thrust::multiplies<double> mult_op;
+
+  // TODO: how to avoid this temp vec?
+  thrust::device_ptr<double> v_temp = thrust::device_malloc(n*sizeof(double));
+
+  // compute v=x*id
+  thrust::transform(thrust::device, dev_v, dev_v+n, id_v, v_temp, mult_op);
+  // compute v[i] = ( log(v[i]) if v[i]>0, otherwise 0 )
+  thrust::transform(thrust::device, v_temp, v_temp+n, v_temp, log_select_op); 
+  // compute v[i] = w[i]*v[i] 
+  thrust::transform(thrust::device, v_temp, v_temp+n, wei_v, v_temp, mult_op); 
+  // sum up
+  const double sum = thrust::reduce(thrust::device, v_temp, v_temp+n, 0.0, plus_op);
+
+  thrust::device_free(v_temp);
+
+  return sum;
+}
+
+
 /** @brief compute sum(d1[i]) */
 double thrust_sum_kernel(int n, double* d1)
 {