Update IPM NLP solvers to be aware of finite element 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, bool use_weighted_vars=false)
       : n_vars(n_mesh_elem),
         comm(MPI_COMM_WORLD),
+        use_weighted_space_(use_weighted_vars),
         solver_(nullptr)
   {
     // create the members
@@ -257,10 +258,49 @@ class DenseConsEx1 : public hiop::hiopInterfaceDenseConstraints
                         double alpha_du,
                         double alpha_pr,
                         int ls_trials);
+
+  bool applyM(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 applyH(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 applyHinv(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 @applyM, @applyH, and @applyHinv
+   */
+  virtual bool useWeightedInnerProducts()
+  {
+    return use_weighted_space_;
+  }
 
 private:
   int n_vars;
   MPI_Comm comm;
+  bool use_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

@@ -279,6 +279,90 @@ 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.
+   *
+   * @note If the variables are not coming from discretizations (specified by a false 
+   * value returned by @useWeightedInnerProducts), 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. 
+   *
+   * @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)
+  {
+    //the default impl. instructs HiOp to use Euclidean/l^2 variables (M=I) internally
+    return false;
+  }
+
+  /**
+   * 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 If the variables are not coming from discretizations (specified by a false 
+   * value returned by @useWeightedInnerProducts), 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 M.
+   *
+   * @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 applyH(const size_type& n, const double* x, double* y)
+  {
+    //the default impl. instructs HiOp to use Euclidean/l^2 variables (H=I) internally
+    return false;
+  }
+
+  /**
+   * 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
+   * important role in computing the "dual" norms and, in turn, to ensure mesh
+   * independent behavior of the IPM solver(s). The inverse apply is also needed by the
+   * specialized solves HiOp performs when the quasi-Newton solver is used. 
+   *
+   * @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
+   */
+  virtual bool applyHinv(const size_type& n, const double* x, double* y)
+  {
+    //the default impl. instructs HiOp to use Euclidean/l^2 variables (H=I) internally
+    return false;
+  }
+
+  /**
+   * Enables the use of weighted inner products via @applyM, @applyH, and @applyHinv
+   *
+   * See also notes for @applyM, @applyH, and @applyHinv.
+   *
+   * @return true if enabled, false to default to Euclidean space with <u,v>=u^T*v
+   */
+  virtual bool useWeightedInnerProducts()
+  {
+    return false;
+  }
+
   /**
    * 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)
 {

src/LinAlg/VectorHipKernels.hpp

@@ -187,10 +187,17 @@ 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 */
-double linear_damping_term_kernel(int n, const double* vd, const double* ld, const double* rd, double mu, double kappa_d);
+double linear_damping_term_kernel(int n,
+                                  const double* vd,
+                                  const double* ld,
+                                  const double* rd,
+                                  double mu,
+                                  double kappa_d);
 /** @brief compute min(d1) */
 double min_local_kernel(int n, double* d1);
 /** @brief Checks if selected elements of `d1` are positive */

src/LinAlg/hiopVector.hpp

@@ -643,6 +643,21 @@ class hiopVector
    */
   virtual double logBarrier_local(const hiopVector& select) const = 0;
 
+  /**
+   * @brief Weighted sum of selected log(this[i])
+   *
+   * @param[in] weights - weights vector
+   * @param[in] select - pattern selection
+   *
+   * @pre `this`, `weights`, and `select` have same partitioning.
+   * @pre Selected elements of `this` are > 0.
+   * @post `this`, `weights`, and `select` are not modified
+   *
+   * @warning This is local method only!
+   */
+  virtual double logBarrierWeighted_local(const hiopVector& select, const hiopVector& weights) const = 0;
+
+
   /**
    * @brief adds the gradient of the log barrier, namely this[i]=this[i]+alpha*1/select(x[i])
    *