Skip to content

Commit

Permalink
process constraints with RAJA
Browse files Browse the repository at this point in the history
  • Loading branch information
@nychiang
nychiang committed May 4, 2023
1 parent a8d4183 commit edc4e7b
Show file tree
Hide file tree
Showing 11 changed files with 431 additions and 91 deletions.
35 changes: 35 additions & 0 deletions src/LinAlg/VectorCudaKernels.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
}

}

}
42 changes: 40 additions & 2 deletions src/LinAlg/hiopVector.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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:
Expand Down
18 changes: 18 additions & 0 deletions src/LinAlg/hiopVectorCuda.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

15 changes: 15 additions & 0 deletions src/LinAlg/hiopVectorCuda.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_; }

Expand Down
18 changes: 18 additions & 0 deletions src/LinAlg/hiopVectorHip.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

15 changes: 15 additions & 0 deletions src/LinAlg/hiopVectorHip.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_; }

Expand Down
88 changes: 88 additions & 0 deletions src/LinAlg/hiopVectorPar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.;
}
}
}


};
15 changes: 15 additions & 0 deletions src/LinAlg/hiopVectorPar.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
*/
Expand Down
15 changes: 15 additions & 0 deletions src/LinAlg/hiopVectorRaja.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_;
Expand Down
Loading

0 comments on commit edc4e7b

Please sign in to comment.