Process constraints on device

src/LinAlg/VectorCudaKernels.cu

@@ -622,6 +622,28 @@ __global__ void compute_cusum_cu(int n, int* vec, const double* id)
   }
 }
 
+/** @brief set veq[i]=1 and vineq[i]=0 if v1[i-1]=v2[i-1], and veq[i]=0 and vineq[i]=01if v1[i-1]!=v2[i-1]*/
+__global__ void mark_cons_cu(int n, int* veq, int* vineq, const double* v1, const double* v2)
+{
+  const int num_threads = blockDim.x * gridDim.x;
+  const int tid = blockIdx.x * blockDim.x + threadIdx.x;    
+  for (int i = tid; i < n; i += num_threads) {
+    if(i==0) {
+      veq[i] = 0;
+      vineq[i] = 0;
+    } else {
+      // from i=1..n
+      if(v1[i-1]==v2[i-1]){
+        veq[i] = 1;
+        vineq[i] = 0;
+      } else {
+        veq[i] = 0;   
+        vineq[i] = 1;     
+      }
+    }
+  }
+}
+
 /// @brief Copy the entries in 'dd' where corresponding 'ix' is nonzero, to vd starting at start_index_in_dest.
 __global__ void copyToStartingAt_w_pattern_cu(int n_src, 
                                               int n_dest,
@@ -1462,6 +1484,19 @@ void compute_cusum_kernel(int sz, int* buf, const double* id)
   thrust::inclusive_scan(dev_v, dev_v + sz, dev_v); // in-place scan
 }
 
+/** @brief compute cusum from the given pattern*/
+void compute_cusum_kernel(int sz, int* ebuf, int* ibuf, const double* v1, const double* v2)
+{
+  int num_blocks = (sz+block_size-1)/block_size;
+  mark_cons_cu<<<num_blocks,block_size>>>(sz, ebuf, ibuf, v1, v2);
+
+  thrust::device_ptr<int> dev_e = thrust::device_pointer_cast(ebuf);
+  thrust::device_ptr<int> dev_i = thrust::device_pointer_cast(ibuf);
+
+  thrust::inclusive_scan(dev_e, dev_e + sz, dev_e); // in-place scan
+  thrust::inclusive_scan(dev_i, dev_i + sz, dev_i); // in-place scan
+}
+
 }
 
 }

src/LinAlg/hiopVector.hpp

@@ -1019,8 +1019,8 @@ class hiopVector
    * @param[in] xu - lower bound of primal variable `x`
    * @param[in] ixl - index of the variables with lower bounds
    * @param[in] ixu - index of the variables with upper bounds
-   * @param[out] n_bnds_low - number of variables with lower bounds only
-   * @param[out] n_bnds_upp - number of variables with upper bounds only
+   * @param[out] n_bnds_low - number of variables with lower bounds
+   * @param[out] n_bnds_upp - number of variables with upper bounds
    * @param[out] n_bnds_lu - number of variables with both lower and upper bounds
    * @param[out] n_fixed_vars - number of fixed variables
    * @param[in] fixed_var_tol - tolerance used to define fixed variables
@@ -1050,6 +1050,44 @@ class hiopVector
                                 const double& fixed_var_tol,
                                 const double& fixed_var_perturb) = 0;
 
+  /**
+   * @brief process constraints. Firstly the constraints are split to equalities and 
+   *        inequalities. Then it preprocesses inequality bounds and returns counts of
+   *        the constraints with lower, upper, and lower and upper bounds.
+   *
+   * @param[in] this - crhs, the right hand side of equality constraints
+   * @param[in] gl_vec - lower bounds for all the constraints
+   * @param[in] gu_vec - upper bounds for all the constraints
+   * @param[out] dl_vec - lower bounds for inequality constraints
+   * @param[out] du_vec - upper bounds for inequality constraints
+   * @param[out] idl_vec - index of the inequality constraints with lower bounds
+   * @param[out] idu_vec - index of the inequality constraints with upper bounds
+   * @param[out] n_ineq_low - number of inequality constraints with lower bounds
+   * @param[out] n_ineq_upp - number of inequality constraints with lower bounds
+   * @param[out] n_ineq_lu - number of inequality constraints with both lower and upper bounds
+   * @param[out] cons_eq_mapping - a map between equality constaints and full list of constraints
+   * @param[out] cons_ineq_mapping - a map between inequality constaints and full list of constraints
+   * @param[out] eqcon_type - types of all the equality constraints
+   * @param[out] incon_type - types of all the inequality constraints
+   * @param[in] cons_type - types of all the constraints
+   * 
+   * @pre this is a local method
+   */
+  virtual void process_constraints_local(const hiopVector& gl_vec,
+                                         const hiopVector& gu_vec,
+                                         hiopVector& dl_vec,
+                                         hiopVector& du_vec,
+                                         hiopVector& idl_vec,
+                                         hiopVector& idu_vec,
+                                         size_type& n_ineq_low,
+                                         size_type& n_ineq_upp,
+                                         size_type& n_ineq_lu,
+                                         hiopVectorInt& cons_eq_mapping,
+                                         hiopVectorInt& cons_ineq_mapping,
+                                         hiopInterfaceBase::NonlinearityType* eqcon_type,
+                                         hiopInterfaceBase::NonlinearityType* incon_type,
+                                         hiopInterfaceBase::NonlinearityType* cons_type) = 0;
+
 protected:
   size_type n_; //we assume sequential data
 protected:

src/LinAlg/hiopVectorCuda.cpp

@@ -1151,5 +1151,23 @@ void hiopVectorCuda::relax_bounds_vec(hiopVector& xu,
                                   fixed_var_perturb);
 }
 
+void hiopVectorCuda::process_constraints_local(const hiopVector& gl_vec,
+                                               const hiopVector& gu_vec,
+                                               hiopVector& dl_vec,
+                                               hiopVector& du_vec,
+                                               hiopVector& idl_vec,
+                                               hiopVector& idu_vec,
+                                               size_type& n_ineq_low,
+                                               size_type& n_ineq_upp,
+                                               size_type& n_ineq_lu,
+                                               hiopVectorInt& cons_eq_mapping,
+                                               hiopVectorInt& cons_ineq_mapping,
+                                               hiopInterfaceBase::NonlinearityType* eqcon_type,
+                                               hiopInterfaceBase::NonlinearityType* incon_type,
+                                               hiopInterfaceBase::NonlinearityType* cons_type)
+{
+  assert(0&&"not yet");
+}
+
 } // namespace hiop
 

src/LinAlg/hiopVectorCuda.hpp

@@ -362,6 +362,21 @@ class hiopVectorCuda : public hiopVector
                                 const double& fixed_var_tol,
                                 const double& fixed_var_perturb);
 
+  virtual void process_constraints_local(const hiopVector& gl_vec,
+                                         const hiopVector& gu_vec,
+                                         hiopVector& dl_vec,
+                                         hiopVector& du_vec,
+                                         hiopVector& idl_vec,
+                                         hiopVector& idu_vec,
+                                         size_type& n_ineq_low,
+                                         size_type& n_ineq_upp,
+                                         size_type& n_ineq_lu,
+                                         hiopVectorInt& cons_eq_mapping,
+                                         hiopVectorInt& cons_ineq_mapping,
+                                         hiopInterfaceBase::NonlinearityType* eqcon_type,
+                                         hiopInterfaceBase::NonlinearityType* incon_type,
+                                         hiopInterfaceBase::NonlinearityType* cons_type);
+
   /* functions for this class */
   inline MPI_Comm get_mpi_comm() const { return comm_; }
 

src/LinAlg/hiopVectorHip.cpp

@@ -1155,5 +1155,23 @@ void hiopVectorHip::relax_bounds_vec(hiopVector& xu,
                                          fixed_var_perturb);
 }
 
+void hiopVectorHip::process_constraints_local(const hiopVector& gl_vec,
+                                              const hiopVector& gu_vec,
+                                              hiopVector& dl_vec,
+                                              hiopVector& du_vec,
+                                              hiopVector& idl_vec,
+                                              hiopVector& idu_vec,
+                                              size_type& n_ineq_low,
+                                              size_type& n_ineq_upp,
+                                              size_type& n_ineq_lu,
+                                              hiopVectorInt& cons_eq_mapping,
+                                              hiopVectorInt& cons_ineq_mapping,
+                                              hiopInterfaceBase::NonlinearityType* eqcon_type,
+                                              hiopInterfaceBase::NonlinearityType* incon_type,
+                                              hiopInterfaceBase::NonlinearityType* cons_type)
+{
+  assert(0&&"not yet");
+}
+
 } // namespace hiop
 

src/LinAlg/hiopVectorHip.hpp

@@ -364,6 +364,21 @@ class hiopVectorHip : public hiopVector
                                 const double& fixed_var_tol,
                                 const double& fixed_var_perturb);
 
+  virtual void process_constraints_local(const hiopVector& gl_vec,
+                                         const hiopVector& gu_vec,
+                                         hiopVector& dl_vec,
+                                         hiopVector& du_vec,
+                                         hiopVector& idl_vec,
+                                         hiopVector& idu_vec,
+                                         size_type& n_ineq_low,
+                                         size_type& n_ineq_upp,
+                                         size_type& n_ineq_lu,
+                                         hiopVectorInt& cons_eq_mapping,
+                                         hiopVectorInt& cons_ineq_mapping,
+                                         hiopInterfaceBase::NonlinearityType* eqcon_type,
+                                         hiopInterfaceBase::NonlinearityType* incon_type,
+                                         hiopInterfaceBase::NonlinearityType* cons_type);
+
   /* functions for this class */
   inline MPI_Comm get_mpi_comm() const { return comm_; }
 

src/LinAlg/hiopVectorPar.cpp

@@ -1406,4 +1406,92 @@ void hiopVectorPar::relax_bounds_vec(hiopVector& xu,
   }
 }
 
+void hiopVectorPar::process_constraints_local(const hiopVector& gl_vec,
+                                              const hiopVector& gu_vec,
+                                              hiopVector& dl_vec,
+                                              hiopVector& du_vec,
+                                              hiopVector& idl_vec,
+                                              hiopVector& idu_vec,
+                                              size_type& n_ineq_low,
+                                              size_type& n_ineq_upp,
+                                              size_type& n_ineq_lu,
+                                              hiopVectorInt& cons_eq_mapping,
+                                              hiopVectorInt& cons_ineq_mapping,
+                                              hiopInterfaceBase::NonlinearityType* eqcon_type,
+                                              hiopInterfaceBase::NonlinearityType* incon_type,
+                                              hiopInterfaceBase::NonlinearityType* cons_type)
+{
+#ifdef HIOP_DEEPCHECKS
+  assert(gl_vec.get_local_size() == gu_vec.get_local_size());
+  assert(dl_vec.get_local_size() == du_vec.get_local_size());
+  assert(dl_vec.get_local_size() == idl_vec.get_local_size());
+  assert(dl_vec.get_local_size() == idu_vec.get_local_size());
+  assert(dl_vec.get_local_size() + this->get_local_size() == gu_vec.get_local_size());
+#endif
+
+  double *dl = dl_vec.local_data();
+  double *du = du_vec.local_data();
+  const double *gl = gl_vec.local_data_const();
+  const double *gu = gu_vec.local_data_const();
+  double* idl = idl_vec.local_data();
+  double* idu = idu_vec.local_data();
+
+  double *c_rhs = this->local_data();
+  index_type *eqcon_map = cons_eq_mapping.local_data();
+  index_type *incon_map = cons_ineq_mapping.local_data();
+
+  size_type n_eq = this->get_local_size();
+  size_type n_in = dl_vec.get_local_size();
+  size_type n_cons = n_eq + n_in;
+
+  /* splitting (preprocessing) step done on the CPU */
+  size_type it_eq = 0;
+  size_type it_in = 0;
+  for(index_type i = 0; i < n_cons; i++) {
+    if(gl[i] == gu[i]) {
+      eqcon_type[it_eq] = cons_type[i]; 
+      c_rhs[it_eq] = gl[i]; 
+      eqcon_map[it_eq] = i;
+      it_eq++;
+    } else {
+#ifdef HIOP_DEEPCHECKS
+      assert(gl[i] <= gu[i] && "please fix the inconsistent inequality constraints, otherwise the problem is infeasible");
+#endif
+      incon_map[it_in] = cons_type[i];
+      dl[it_in] = gl[i]; 
+      du[it_in] = gu[i]; 
+      incon_map[it_in] = i;
+      it_in++;
+    }
+  }
+  assert(it_eq == n_eq);
+  assert(it_in == n_in);
+
+  /* iterate over the inequalities and build the idl(ow) and idu(pp) vectors */
+  n_ineq_low = 0;
+  n_ineq_upp = 0;
+  n_ineq_lu = 0;
+
+  for(index_type i=0; i<n_in; i++) {
+    if(dl[i]>-1e20) { 
+      idl[i] = 1.;
+      n_ineq_low++; 
+      if(du[i]< 1e20) {
+        n_ineq_lu++;
+      }
+    }
+    else {
+      idl[i]=0.;
+    }
+
+    if(du[i]< 1e20) { 
+      idu[i] = 1.;
+      n_ineq_upp++; 
+    } else {
+      idu[i]=0.;
+    }
+  }
+}
+
+
 };

src/LinAlg/hiopVectorPar.hpp

@@ -324,6 +324,21 @@ class hiopVectorPar : public hiopVector
                                 const double& fixed_var_tol,
                                 const double& fixed_var_perturb);
 
+  virtual void process_constraints_local(const hiopVector& gl_vec,
+                                         const hiopVector& gu_vec,
+                                         hiopVector& dl_vec,
+                                         hiopVector& du_vec,
+                                         hiopVector& idl_vec,
+                                         hiopVector& idu_vec,
+                                         size_type& n_ineq_low,
+                                         size_type& n_ineq_upp,
+                                         size_type& n_ineq_lu,
+                                         hiopVectorInt& cons_eq_mapping,
+                                         hiopVectorInt& cons_ineq_mapping,
+                                         hiopInterfaceBase::NonlinearityType* eqcon_type,
+                                         hiopInterfaceBase::NonlinearityType* incon_type,
+                                         hiopInterfaceBase::NonlinearityType* cons_type);
+
   /**
    * @brief accessor to the execution policy
    */

src/LinAlg/hiopVectorRaja.hpp

@@ -327,6 +327,21 @@ class hiopVectorRaja : public hiopVector
                                 const double& fixed_var_tol,
                                 const double& fixed_var_perturb);
 
+  virtual void process_constraints_local(const hiopVector& gl_vec,
+                                         const hiopVector& gu_vec,
+                                         hiopVector& dl_vec,
+                                         hiopVector& du_vec,
+                                         hiopVector& idl_vec,
+                                         hiopVector& idu_vec,
+                                         size_type& n_ineq_low,
+                                         size_type& n_ineq_upp,
+                                         size_type& n_ineq_lu,
+                                         hiopVectorInt& cons_eq_mapping,
+                                         hiopVectorInt& cons_ineq_mapping,
+                                         hiopInterfaceBase::NonlinearityType* eqcon_type,
+                                         hiopInterfaceBase::NonlinearityType* incon_type,
+                                         hiopInterfaceBase::NonlinearityType* cons_type);
+
   const ExecSpace<MEMBACKEND, EXECPOLICYRAJA>& exec_space() const
   {
     return exec_space_;