hiopInterface.hpp

// Copyright (c) 2017, Lawrence Livermore National Security, LLC.
// Produced at the Lawrence Livermore National Laboratory (LLNL).
// LLNL-CODE-742473. All rights reserved.
//
// This file is part of HiOp. For details, see https://github.com/LLNL/hiop. HiOp
// is released under the BSD 3-clause license (https://opensource.org/licenses/BSD-3-Clause).
// Please also read "Additional BSD Notice" below.
//
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
// i. Redistributions of source code must retain the above copyright notice, this list
// of conditions and the disclaimer below.
// ii. Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the disclaimer (as noted below) in the documentation and/or
// other materials provided with the distribution.
// iii. Neither the name of the LLNS/LLNL nor the names of its contributors may be used to
// endorse or promote products derived from this software without specific prior written
// permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
// SHALL LAWRENCE LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
// OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
// AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
// EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Additional BSD Notice
// 1. This notice is required to be provided under our contract with the U.S. Department
// of Energy (DOE). This work was produced at Lawrence Livermore National Laboratory under
// Contract No. DE-AC52-07NA27344 with the DOE.
// 2. Neither the United States Government nor Lawrence Livermore National Security, LLC
// nor any of their employees, makes any warranty, express or implied, or assumes any
// liability or responsibility for the accuracy, completeness, or usefulness of any
// information, apparatus, product, or process disclosed, or represents that its use would
// not infringe privately-owned rights.
// 3. Also, reference herein to any specific commercial products, process, or services by
// trade name, trademark, manufacturer or otherwise does not necessarily constitute or
// imply its endorsement, recommendation, or favoring by the United States Government or
// Lawrence Livermore National Security, LLC. The views and opinions of authors expressed
// herein do not necessarily state or reflect those of the United States Government or
// Lawrence Livermore National Security, LLC, and shall not be used for advertising or
// product endorsement purposes.

/**
 * @file hiopInterface.hpp
 *
 * @author Cosmin G. Petra <petra1@llnl.gov>, LLNL
 * @author Nai-Yuan Chiang <chiang7@lnnl.gov>, LLNL
 *
 */

#ifndef HIOP_INTERFACE_BASE
#define HIOP_INTERFACE_BASE

#include "hiop_defs.hpp"
#include "hiopMPI.hpp"

namespace hiop
{
// Solver status codes
enum hiopSolveStatus
{
  //(partial) success
  Solve_Success = 0,
  Solve_Success_RelTol = 1,
  Solve_Acceptable_Level = 2,
  Infeasible_Problem = 5,
  Iterates_Diverging = 6,
  Feasible_Not_Optimal = 7,
  // solver stopped based on user-defined criteria that are not related to optimality
  Max_Iter_Exceeded = 10,
  Max_CpuTime_Exceeded = 11,
  User_Stopped = 12,

  // NLP algorithm/solver reports issues in solving the problem and stops without being certain
  // that is solved the problem to optimality or that the problem is infeasible.
  // Feasible_Point_Found,
  NlpAlgorithm_failure = -1,
  Diverging_Iterates = -2,
  Search_Dir_Too_Small = -3,
  Steplength_Too_Small = -4,
  Err_Step_Computation = -5,
  // errors related to user-provided data (e.g., inconsistent problem specification, 'nans' in the
  // function/sensitivity evaluations, invalid options)
  Invalid_Problem_Definition = -11,
  Invalid_Parallelization = -12,
  Invalid_UserOption = -13,
  Invalid_Number = -14,
  Error_In_User_Function = -15,
  Error_In_FR = -16,

  // ungraceful errors and returns
  Exception_Unrecoverable = -100,
  Memory_Alloc_Problem = -101,
  SolverInternal_Error = -199,

  // unknown NLP solver errors or return codes
  UnknownNLPSolveStatus = -1000,
  SolveInitializationError = -1001,

  // intermediary statuses for the solver
  NlpSolve_IncompleteInit = -10001,
  NlpSolve_SolveNotCalled = -10002,
  NlpSolve_Pending = -10003
};

/** Base class for the solver's interface that has no assumptions how the
 *  matrices are stored. The vectors are dense and distributed row-wise.
 *  The data distribution is decided by the calling code (that implements
 *  this interface) and specified to the optimization via 'get_vecdistrib_info'
 *
 *  Three possible implementations are for sparse NLPs (hiopInterfaceSparse),
 *  mixed dense-sparse NLPs (hiopInterfaceMDS), and NLPs with small
 *  number of global constraints (hiopInterfaceDenseConstraints).
 *
 * @note Please take notice of the following notes regarding the implementation of
 * hiop::hiopInterfaceMDS on the device. All pointers marked as "managed by Umpire"
 * are allocated by HiOp using the Umpire's API. They all are addressed in the
 * same memory space; however, the memory space can be host (typically CPU),
 * device (typically GPU), or unified memory (um) spaces as per Umpire
 * specification. The selection of the memory space is done via the option
 * "mem_space" of HiOp. It is the responsibility of the implementers of the
 * HiOp's interfaces to work with the "managed by Umpire" pointers in the same
 * memory space as the one specified by the "mem_space" option.
 *
 * @note The above note does not currently apply to the NLP interfaces
 * hiop::hiopInterfaceDenseConstraints and hiop::hiopInterfaceSparse) and the pointers
 * marked as "managed by Umpire" are in the host/CPU memory space (subject to change
 * in future versions of HiOp).
 */
class hiopInterfaceBase
{
public:
  // Types indicating linearity or nonlinearity.
  enum NonlinearityType
  {
    hiopLinear = 0,
    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() {};
  virtual ~hiopInterfaceBase() {};

  /** Specifies the problem dimensions.
   *
   * @param n global number of variables
   * @param m number of constraints
   */
  virtual bool get_prob_sizes(size_type& n, size_type& m) = 0;

  /** Specifies the type of optimization problem
   * @param[out] type indicating whether the optimization problem is
   *                  linearily, quadratically, or general nonlinearily.
   * TODO: need to `deepcheck` is this return value matches the returned type array from
   * `get_vars_info` and `get_cons_info`
   */
  virtual bool get_prob_info(NonlinearityType& type)
  {
    type = hiopInterfaceBase::hiopNonlinear;
    return true;
  }

  /** Specifies bounds on the variables.
   *
   * @param[in] n global number of constraints
   * @param[out] xlow array of lower bounds. A value of -1e20 or less means no lower
   *                  bound is present (managed by Umpire)
   * @param[out] xupp array of upper bounds. A value of 1e20 or more means no upper
   *                  bound is present  (managed by Umpire)
   * @param[out] type array of indicating whether the variables enters the objective
   *                  linearily, quadratically, or general nonlinearily. Momentarily
   *                  all bounds should be marked as nonlinear (allocated on host).
   */
  virtual bool get_vars_info(const size_type& n, double* xlow, double* xupp, NonlinearityType* type) = 0;

  /** Specififes the bounds on the constraints.
   *
   * @param[in] m number of constraints
   * @param[out] clow array of lower bounds for constraints. A value of -1e20 or less means no lower
   *                  bound is present (managed by Umpire)
   * @param[out] cupp array of upper bounds for constraints. A value of 1e20 or more means no upper
   *                  bound is present (managed by Umpire)
   * @param[out] type array of indicating whether the constraint is linear, quadratic, or general
   *                  nonlinear. Momentarily all bounds should be marked as nonlinear (allocated on host).
   */
  virtual bool get_cons_info(const size_type& m, double* clow, double* cupp, NonlinearityType* type) = 0;

  /** Method the evaluation of the objective function.
   *
   * @param[in] n global size of the problem
   * @param[in] x array with the local entries of the primal variable (managed by Umpire)
   * @param[in] new_x whether x has been changed from the previous calls to other evaluation methods
   *                  (gradient, constraints, Jacobian, and Hessian).
   * @param[out] obj_value the value of the objective function at @p x
   *
   * @note When MPI is enabled, each rank returns the objective value in @p obj_value. @p x points to
   * the local entries and the function is responsible for knowing the local buffer size.
   */
  virtual bool eval_f(const size_type& n, const double* x, bool new_x, double& obj_value) = 0;

  /** Method for the evaluation of the gradient of objective.
   *
   * @param[in] n global size of the problem
   * @param[in] x array with the local entries of the primal variable  (managed by Umpire)
   * @param[in] new_x whether x has been changed from the previous calls to other evaluation methods
   * (                function, constraints, Jacobian, and Hessian)
   * @param[out] gradf the entries of the gradient of the objective function at @p x, local to the
   * MPI rank (managed by Umpire)
   *
   *  @note When MPI is enabled, each rank should access only the local buffers @p x and @p gradf.
   */
  virtual bool eval_grad_f(const size_type& n, const double* x, bool new_x, double* gradf) = 0;

  /** Evaluates a subset of the constraints @p cons(@p x). The subset is of size
   *  @p num_cons and is described by indexes in the @p idx_cons array. The method will be called at each
   *  iteration separately for the equality constraints subset and for the inequality constraints subset.
   *  This is done for performance considerations, to avoid auxiliary/temporary storage and copying.
   *
   * @param[in] n the global number of variables
   * @param[in] m the number of constraints
   * @param[in] num_cons the number constraints/size of subset to be evaluated
   * @param[in] idx_cons: indexes in {1,2,...,m} of the constraints to be evaluated  (managed by Umpire)
   * @param[in] x the point where the constraints need to be evaluated  (managed by Umpire)
   * @param[in] new_x whether x has been changed from the previous call to f, grad_f, or Jac
   * @param[out] cons array of size num_cons containing the value of the  constraints indicated by
   *                    @p idx_cons  (managed by Umpire)
   *
   *  @note When MPI is enabled, every rank populates @p cons since the constraints are not distributed.
   */
  virtual bool eval_cons(const size_type& n,
                         const size_type& m,
                         const size_type& num_cons,
                         const index_type* idx_cons,
                         const double* x,
                         bool new_x,
                         double* cons) = 0;

  /** Evaluates the constraints body @p cons(@p x), both equalities and inequalities, in one call.
   *
   * @param[in] n the global number of variables
   * @param[in] m the number of constraints
   * @param[in] x the point where the constraints need to be evaluated  (managed by Umpire)
   * @param[in] new_x whether x has been changed from the previous call to f, grad_f, or Jac
   * @param[out] cons array of size num_cons containing the value of the  constraints indicated by
   *                     @p idx_cons  (managed by Umpire)
   *
   * HiOp will first call the other hiopInterfaceBase::eval_cons() twice. If the implementer/user wants the
   * functionality  of this "one-call" overload, he should return false from the other
   * hiopInterfaceBase::eval_cons() (during both calls).
   *
   *  @note When MPI is enabled, every rank populates @p cons since the constraints are not distributed.
   */
  virtual bool eval_cons(const size_type& n, const size_type& m, const double* x, bool new_x, double* cons)
  {
    return false;
  }

  /** Passes the communicator, defaults to MPI_COMM_WORLD (dummy for non-MPI builds)  */
  virtual bool get_MPI_comm(MPI_Comm& comm_out)
  {
    comm_out = MPI_COMM_WORLD;
    return true;
  }

  /**
   * Method for column partitioning specification for distributed memory vectors. Process P owns
   * cols[P], cols[P]+1, ..., cols[P+1]-1, P={0,1,...,NumRanks}.
   *
   * Example: for a vector x of @p global_n=6 elements on 3 ranks, the column partitioning is
   * @p cols=[0,2,4,6].
   *
   * The caller manages memory associated with @p cols, which is an array of size NumRanks+1
   */
  virtual bool get_vecdistrib_info(size_type global_n, index_type* cols)
  {
    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.
   *
   * @note Avoid using this method since it will be removed in a future release and replaced with
   * the same-name method below.
   *
   * The method returns true (and populates @p x0) or returns false, in which case HiOp will
   * internally set @p x0 to all zero (still subject to internal adjustements).
   *
   * By default, HiOp first calls the overloaded primal-dual starting point specification
   * (overloaded) method get_starting_point() (see below). If the above returns false, HiOp will then call
   * this method.
   *
   * @param[in] n the global number of variables
   * @param[out] x0 the user-defined initial values for the primal variablers (managed by Umpire)
   *
   */
  virtual bool get_starting_point(const size_type& n, double* x0) { return false; }

  /**
   * Method provides a primal or a primal-dual starting point. This point is subject
   * to internal adjustments in HiOp.
   *
   * If the user (implementer of this method) has good estimates only of the primal variables,
   * the method should populate @p x0 with these values and return true. The @p duals_avail
   * should be set to false; internally, HiOp will not access @p z_bndL0, @p z_bndU0, and
   * @p lambda0 in this case.
   *
   * If the user (implementer of this method) has good estimates of the duals of bound constraints
   * and of inequality and equality constraints, @p duals_avail boolean argument should
   * be set to true and the respective duals should be provided (in @p z_bndL0 and @p z_bndU0 and
   * @p lambda0, respectively). In this case, the user should also set @p x0 to his/her estimate
   * of primal variables and return true.
   *
   * If user does not have high-quality (primal or primal-dual) starting points, the method should
   * return false (see note below).
   *
   * @note When this method returns false, HiOp will call the overload
   * get_starting_point() for only primal variables (see the above function). This behaviour is for backward compatibility
   * and will be removed in a future release.
   *
   * @param[in] n the global number of variables
   * @param[in] m the number of constraints
   * @param[out] x0 the user-defined initial values for the primal variablers (managed by Umpire)
   * @param[out] duals_avail a boolean argument which indicates whether the initial values of duals are given by the user
   * @param[out] z_bndL0 the user-defined initial values for the duals of the variable lower bounds (managed by Umpire)
   * @param[out] z_bndU0 the user-defined initial values for the duals of the variable upper bounds (managed by Umpire)
   * @param[out] lambda0 the user-defined initial values for the duals of the constraints (managed by Umpire)
   * @param[out] slacks_avail a boolean argument which indicates whether the initial values for the inequality slacks
   *                            (added by HiOp internally) are given by the user
   * @param[out] ineq_slack the user-defined initial values for the slacks added to transfer inequalities to equalities
   *                          (managed by Umpire)
   *
   */
  virtual bool get_starting_point(const size_type& n,
                                  const size_type& m,
                                  double* x0,
                                  bool& duals_avail,
                                  double* z_bndL0,
                                  double* z_bndU0,
                                  double* lambda0,
                                  bool& slacks_avail,
                                  double* ineq_slack)
  {
    duals_avail = false;
    slacks_avail = false;
    return false;
  }

  /**
   * Method provides a primal-dual starting point for warm start. This point is subject
   * to internal adjustments in HiOp.
   *
   * User provides starting point for all the iterate variable used in HiOp.
   * This method is for advanced users, as it will skip all the other safeguard in HiOp, e.g., project x into bounds.
   *
   * @param[in] n the global number of variables
   * @param[in] m the number of constraints
   * @param[out] x0 the user-defined initial values for the primal variablers (managed by Umpire)
   * @param[out] z_bndL0 the user-defined initial values for the duals of the variable lower bounds (managed by Umpire)
   * @param[out] z_bndU0 the user-defined initial values for the duals of the variable upper bounds (managed by Umpire)
   * @param[out] lambda0 the user-defined initial values for the duals of the constraints (managed by Umpire)
   * @param[out] ineq_slack the user-defined initial values for the slacks added to transfer inequalities to equalities
   *                          (managed by Umpire)
   * @param[out] vl0 the user-defined initial values for the duals of the constraint lower bounds (managed by Umpire)
   * @param[out] vu0 the user-defined initial values for the duals of the constraint upper bounds (managed by Umpire)
   *
   */
  virtual bool get_warmstart_point(const size_type& n,
                                   const size_type& m,
                                   double* x0,
                                   double* z_bndL0,
                                   double* z_bndU0,
                                   double* lambda0,
                                   double* ineq_slack,
                                   double* vl0,
                                   double* vu0)
  {
    return false;
  }

  /**
   * Callback method called by HiOp when the optimal solution is reached. User should use it
   * to retrieve primal-dual optimal solution.
   *
   * @param[in] status status of the solution process
   * @param[in] n global number of variables
   * @param[in] x array of (local) entries of the primal variables at solution (managed by Umpire, see note below)
   * @param[in] z_L array of (local) entries of the dual variables for lower bounds at solution (managed by Umpire,
   * see note below)
   * @param[in] z_U array of (local) entries of the dual variables for upper bounds at solution (managed by Umpire,
   * see note below)
   * @param[in] g array of the values of the constraints body at solution (managed by Umpire, see note below)
   * @param[in] lambda array of (local) entries of the dual variables for constraints at solution (managed by Umpire,
   * see note below)
   * @param[in] obj_value objective value at solution
   *
   * @note HiOp's option `callback_mem_space` can be used to change the memory location of array parameters managaged by
   * Umpire. More specifically, when `callback_mem_space` is set to `host` (and `mem_space` is `device`), HiOp transfers the
   * arrays from device to host first, and then passes/returns pointers on host for the arrays managed by Umpire. These
   * pointers can be then used in host memory space (without the need to rely on or use Umpire).
   *
   */
  virtual void solution_callback(hiopSolveStatus status,
                                 size_type n,
                                 const double* x,
                                 const double* z_L,
                                 const double* z_U,
                                 size_type m,
                                 const double* g,
                                 const double* lambda,
                                 double obj_value)
  {}

  /**
   * Callback for the (end of) iteration. This method is not called during the line-searche
   * procedure. @see solution_callback() for an explanation of the parameters.
   *
   * @note If the user (implementer) of this methods returns false, HiOp will stop the
   * the optimization with hiop::hiopSolveStatus ::User_Stopped return code.
   *
   * @param[in] iter the current iteration number
   * @param[in] obj_value objective value
   * @param[in] logbar_obj_value log barrier objective value
   * @param[in] n global number of variables
   * @param[in] x array of (local) entries of the primal variables (managed by Umpire, see note below)
   * @param[in] z_L array of (local) entries of the dual variables for lower bounds (managed by Umpire, see note below)
   * @param[in] z_U array of (local) entries of the dual variables for upper bounds (managed by Umpire, see note below)
   * @param[in] m_ineq the number of inequality constraints
   * @param[in] s array of the slacks added to transfer inequalities to equalities (managed by Umpire, see note below)
   * @param[in] m the number of constraints
   * @param[in] g array of the values of the constraints body (managed by Umpire, see note below)
   * @param[in] lambda array of (local) entries of the dual variables for constraints (managed by Umpire, see note below)
   * @param[in] inf_pr inf norm of the primal infeasibilities
   * @param[in] inf_du inf norm of the dual infeasibilities
   * @param[in] onenorm_pr one norm of the primal infeasibilities
   * @param[in] mu the log barrier parameter
   * @param[in] alpha_du dual step size
   * @param[in] alpha_pr primal step size
   * @param[in] ls_trials the number of line search iterations
   *
   * @note HiOp's option `callback_mem_space` can be used to change the memory location of array parameters managaged by
   * Umpire. More specifically, when `callback_mem_space` is set to `host` (and `mem_space` is `device`), HiOp transfers the
   * arrays from device to host first, and then passes/returns pointers on host for the arrays managed by Umpire. These
   * pointers can be then used in host memory space (without the need to rely on or use Umpire).
   *
   */
  virtual bool iterate_callback(int iter,
                                double obj_value,
                                double logbar_obj_value,
                                int n,
                                const double* x,
                                const double* z_L,
                                const double* z_U,
                                int m_ineq,
                                const double* s,
                                int m,
                                const double* g,
                                const double* lambda,
                                double inf_pr,
                                double inf_du,
                                double onenorm_pr,
                                double mu,
                                double alpha_du,
                                double alpha_pr,
                                int ls_trials)
  {
    return true;
  }

  /**
   * This method is used to provide user all the internal hiop iterates. @see solution_callback()
   * for an explanation of the parameters.
   *
   * @param[in] x array of (local) entries of the primal variables (managed by Umpire, see note below)
   * @param[in] z_L array of (local) entries of the dual variables for lower bounds (managed by Umpire, see note below)
   * @param[in] z_U array of (local) entries of the dual variables for upper bounds (managed by Umpire, see note below)
   * @param[in] yc array of (local) entries of the dual variables for equality constraints (managed by Umpire, see note
   * below)
   * @param[in] yd array of (local) entries of the dual variables for inequality constraints (managed by Umpire, see note
   * below)
   * @param[in] s array of the slacks added to transfer inequalities to equalities (managed by Umpire, see note below)
   * @param[in] v_L array of (local) entries of the dual variables for constraint lower bounds (managed by Umpire, see note
   * below)
   * @param[in] v_U array of (local) entries of the dual variables for constraint upper bounds (managed by Umpire, see note
   * below)
   *
   * @note HiOp's option `callback_mem_space` can be used to change the memory location of array parameters managaged by
   * Umpire. More specifically, when `callback_mem_space` is set to `host` (and `mem_space` is `device`), HiOp transfers the
   * arrays from device to host first, and then passes/returns pointers on host for the arrays managed by Umpire. These
   * pointers can be then used in host memory space (without the need to rely on or use Umpire).
   *
   */
  virtual bool iterate_full_callback(const double* x,
                                     const double* z_L,
                                     const double* z_U,
                                     const double* yc,
                                     const double* yd,
                                     const double* s,
                                     const double* v_L,
                                     const double* v_U)
  {
    return true;
  }

  /**
   * A wildcard function used to change the primal variables.
   *
   * @note If the user (implementer) of this methods returns false, HiOp will stop the
   * the optimization with hiop::hiopSolveStatus::User_Stopped return code.
   */
  virtual bool force_update_x(const int n, double* x) { return true; }

private:
  hiopInterfaceBase(const hiopInterfaceBase&) {};
  void operator=(const hiopInterfaceBase&) {};
};

/** Specialized interface for NLPs with 'global' but few constraints.
 */
class hiopInterfaceDenseConstraints : public hiopInterfaceBase
{
public:
  hiopInterfaceDenseConstraints() {};
  virtual ~hiopInterfaceDenseConstraints() {};
  /**
   * Evaluates the Jacobian of the subset of constraints indicated by idx_cons and of size num_cons.
   * Example: Assuming idx_cons[k]=i, which means that the gradient of the (i+1)th constraint is
   * to be evaluated, one needs to do Jac[k][0]=d/dx_0 con_i(x), Jac[k][1]=d/dx_1 con_i(x), ...
   * When MPI enabled, each rank computes only the local columns of the Jacobian, that is the partials
   * with respect to local variables.
   *
   * The parameter 'Jac' is passed as as a contiguous array storing the dense Jacobian matrix by rows.
   *
   * Parameters: see eval_cons
   */
  virtual bool eval_Jac_cons(const size_type& n,
                             const size_type& m,
                             const size_type& num_cons,
                             const index_type* idx_cons,
                             const double* x,
                             bool new_x,
                             double* Jac) = 0;

  /**
   * Evaluates the Jacobian of equality and inequality constraints in one call.
   *
   * The main difference from the above 'eval_Jac_cons' is that the implementer/user of this
   * method does not have to split the constraints into equalities and inequalities; instead,
   * HiOp does this internally.
   *
   * The parameter 'Jac' is passed as as a contiguous array storing the dense Jacobian matrix by rows.
   *
   * TODO: build an example (new one-call Nlp formulation derived from ex2) to illustrate this
   * feature and to test HiOp's internal implementation of eq.-ineq. spliting.
   */
  virtual bool eval_Jac_cons(const size_type& n, const size_type& m, const double* x, bool new_x, double* Jac)
  {
    return false;
  }
};

/**
 * Specialized interface for NLPs having mixed DENSE and sparse (MDS) blocks in the
 * Jacobian and Hessian.
 *
 * More specifically, this interface is for specifying optimization problem in x
 * split as (xs,xd), the rule of thumb being that xs have sparse derivatives and
 * xd have dense derivatives
 *
 * min f(x) s.t. g(x) <= or = 0, lb<=x<=ub
 * such that
 *  - Jacobian w.r.t. xs and LagrHessian w.r.t. (xs,xs) are sparse
 *  - Jacobian w.r.t. xd and LagrHessian w.r.t. (xd,xd) are dense
 *  - LagrHessian w.r.t (xs,xd) is zero (later this assumption will be relaxed)
 *
 * @note HiOp expects the sparse variables first and then the dense variables. In many cases,
 * the implementer has to (inconviniently) keep a map between his internal variables
 * indexes and the indexes HiOp.
 *
 * @note This interface is 'local' in the sense that data is not assumed to be
 * distributed across MPI ranks ('get_vecdistrib_info' should return 'false')
 *
 */
class hiopInterfaceMDS : public hiopInterfaceBase
{
public:
  hiopInterfaceMDS() {};
  virtual ~hiopInterfaceMDS() {};

  /**
   *  Returns the sizes and number of nonzeros of the sparse and dense blocks within MDS
   *
   * @param[out] nx_sparse number of sparse variables
   * @param[out] nx_ense number of dense variables
   * @param[out] nnz_sparse_Jace number of nonzeros in the Jacobian of the equalities w.r.t.
   *             sparse variables
   * @param[out] nnz_sparse_Jaci number of nonzeros in the Jacobian of the inequalities w.r.t.
   *             sparse variables
   * @param[out] nnz_sparse_Hess_Lagr_SS number of nonzeros in the (sparse) Hessian w.r.t.
   *             sparse variables
   * @param[out] nnz_sparse_Hess_Lagr_SD reserved, should be set to 0
   */
  virtual bool get_sparse_dense_blocks_info(int& nx_sparse,
                                            int& nx_dense,
                                            int& nnz_sparse_Jaceq,
                                            int& nnz_sparse_Jacineq,
                                            int& nnz_sparse_Hess_Lagr_SS,
                                            int& nnz_sparse_Hess_Lagr_SD) = 0;

  /**
   * Evaluates the Jacobian of constraints split in the sparse (triplet format) and
   * dense matrices (rows storage)
   *
   * This method is called twice per Jacobian evaluation, once for equalities and once for
   * inequalities (see 'eval_cons' for more information). It is advantageous to provide
   * this method when the underlying NLP's constraints come naturally split in equalities
   * and inequalities. When it is not convenient to do so, use 'eval_Jac_cons' below.
   *
   * @param[in] n number of variables
   * @param[in] m Number of constraints
   * @param[in] num_cons number of constraints to evaluate (size of idx_cons array)
   * @param[in] idx_cons indexes of the constraints to evaluate (managed by Umpire)
   * @param[in] x the point at which to evaluate (managed by Umpire)
   * @param[in] new_x indicates whether any of the other eval functions have been evaluated
   *                  previously (false) or not (true) at x
   * @param[in] nsparse number of sparse variables
   * @param[in] ndense number of dense variables
   * @param[in] nnzJacS number of nonzeros in the sparse Jacobian
   * @param[out] iJacS array of row indexes in the sparse Jacobian (managed by Umpire)
   * @param[out] jJacS array of column indexes in the sparse Jacobian (managed by Umpire)
   * @param[out] MJacS array of nonzero values in the sparse Jacobian (managed by Umpire)
   * @param[out] JacD array with the values of the dense Jacobian (managed by Umpire)
   *
   * The implementer of this method should be aware of the following observations.
   * 1) 'JacD' parameter will be always non-null
   * 2) When 'iJacS' and 'jJacS' are non-null, the implementer should provide the (i,j)
   * indexes.
   * 3) When 'MJacS' is non-null, the implementer should provide the values corresponding to
   * entries specified by 'iJacS' and 'jJacS'
   * 4) 'iJacS' and 'jJacS' are both either non-null or null during a call.
   * 5) Both 'iJacS'/'jJacS' and 'MJacS' can be non-null during the same call or only one of
   * them non-null; but they will not be both null.
   */
  virtual bool eval_Jac_cons(const size_type& n,
                             const size_type& m,
                             const size_type& num_cons,
                             const index_type* idx_cons,
                             const double* x,
                             bool new_x,
                             const size_type& nsparse,
                             const size_type& ndense,
                             const size_type& nnzJacS,
                             index_type* iJacS,
                             index_type* jJacS,
                             double* MJacS,
                             double* JacD) = 0;

  /**
   * Evaluates the Jacobian of equality and inequality constraints in one call. This Jacobian is
   * mixed dense-sparse (MDS), which means is structurally split in the sparse (triplet format) and
   * dense matrices (rows storage)
   *
   * The main difference from the above 'eval_Jac_cons' is that the implementer/user of this
   * method does not have to split the constraints into equalities and inequalities; instead,
   * HiOp does this internally.  HiOp will call this method whenever the implementer/user returns
   * false from the 'eval_Jac_cons' above (which is called for equalities and inequalities separately).
   *
   * @param[in] n number of variables
   * @param[in] m Number of constraints
   * @param[in] x the point at which to evaluate (managed by Umpire)
   * @param[in] new_x indicates whether any of the other eval functions have been evaluated previously
   *                  (false) or not (true) at x
   * @param[in] nsparse number of sparse variables
   * @param[in] ndense number of dense variables
   * @param[in] nnzJacS number of nonzeros in the sparse Jacobian
   * @param[out] iJacS array of row indexes in the sparse Jacobian (managed by Umpire)
   * @param[out] jJacS array of column indexes in the sparse Jacobian (managed by Umpire)
   * @param[out] MJacS array of nonzero values in the sparse Jacobian (managed by Umpire)
   * @param[out] JacD array with the values of the dense Jacobian (managed by Umpire)
   *
   * Notes for implementer of this method:
   * 1) 'JacD' parameter will be always non-null.
   * 2) When 'iJacS' and 'jJacS' are non-null, the implementer should provide the (i,j) indexes.
   * 3) When 'MJacS' is non-null, the implementer should provide the values corresponding to
   * entries specified by 'iJacS' and 'jJacS' (managed by Umpire).
   * 4) 'iJacS' and 'jJacS' are both either non-null or null during a call.
   * 5) Both 'iJacS'/'jJacS' and 'MJacS' can be non-null during the same call or only one of them
   * non-null; but they will not be both null.
   */
  virtual bool eval_Jac_cons(const size_type& n,
                             const size_type& m,
                             const double* x,
                             bool new_x,
                             const size_type& nsparse,
                             const size_type& ndense,
                             const size_type& nnzJacS,
                             index_type* iJacS,
                             index_type* jJacS,
                             double* MJacS,
                             double* JacD)
  {
    return false;
  }

  /**
   * Evaluates the Hessian of the Lagrangian function in 3 structural blocks: HSS is the Hessian
   * w.r.t. (xs,xs),  HDD is the Hessian w.r.t. (xd,xd), and HSD is the Hessian w.r.t (xs,xd).
   * Please consult the user manual for a details on the form the Lagrangian function takes.
   *
   * @note HSD is for now assumed to be zero. The implementer should return nnzHSD=0
   * during the first call to 'eval_Hess_Lagr'. On subsequent calls, HiOp will pass the
   * triplet arrays for HSD set to NULL and the implementer (obviously) should not use them.
   *
   * @param[in] n number of variables
   * @param[in] m Number of constraints
   * @param[in] x the point at which to evaluate (managed by Umpire)
   * @param[in] new_x indicates whether any of the other eval functions have been evaluated
   *                  previously (false) or not (true) at x
   * @param[in] obj_factor scalar that multiplies the objective term in the Lagrangian function
   * @param[in] lambda array with values of the multipliers used by the Lagrangian function
   * @param[in] new_lambda indicates whether lambda  values changed since last call
   * @param[in] nsparse number of sparse variables
   * @param[in] ndense number of dense variables
   * @param[in] nnzHSS number of nonzeros in the (sparse) Hessian w.r.t. sparse variables
   * @param[out] iHSS array of row indexes in the Hessian w.r.t. sparse variables (managed by
   *                  Umpire)
   * @param[out] jHSS array of column indexes in the Hessian w.r.t. sparse variables
   *                  (managed by Umpire)
   * @param[out] MHSS array of nonzero values in the Hessian w.r.t. sparse variables
   *                  (managed by Umpire)
   * @param[out] HDDD array with the values of the Hessian w.r.t. to dense variables
   *                  (managed by Umpire)
   * @param[out] iHSD is reserved and should not be accessed
   * @param[out] jHSD is reserved and should not be accessed
   * @param[out] MHSD is reserved and should not be accessed
   * @param[out] HHSD is reserved and should not be accessed
   *
   * Notes
   * 1)-5) from 'eval_Jac_cons' apply to xxxHSS and HDD arrays
   * 6) The order is multipliers is: lambda=[lambda_eq, lambda_ineq]
   */
  virtual bool eval_Hess_Lagr(const size_type& n,
                              const size_type& m,
                              const double* x,
                              bool new_x,
                              const double& obj_factor,
                              const double* lambda,
                              bool new_lambda,
                              const size_type& nsparse,
                              const size_type& ndense,
                              const size_type& nnzHSS,
                              index_type* iHSS,
                              index_type* jHSS,
                              double* MHSS,
                              double* HDD,
                              size_type& nnzHSD,
                              index_type* iHSD,
                              index_type* jHSD,
                              double* MHSD) = 0;
};

/** Specialized interface for NLPs with sparse Jacobian and Hessian matrices.
 *
 * More specifically, this interface is for specifying optimization problem:
 *
 * min f(x) s.t. g(x) <=, =, or >= 0, lb<=x<=ub
 *
 * such that Jacobian w.r.t. x and Hessian of the Lagrangian w.r.t. x are sparse
 *
 * @note this interface is 'local' in the sense that data is not assumed to be
 * distributed across MPI ranks ('get_vecdistrib_info' should return 'false').
 * Acceleration can be however obtained using OpenMP and CUDA via Raja
 * abstraction layer that HiOp uses and via linear solver.
 *
 */
class hiopInterfaceSparse : public hiopInterfaceBase
{
public:
  hiopInterfaceSparse() {};
  virtual ~hiopInterfaceSparse() {};

  /** Get the number of variables and constraints, nonzeros
   * and get the number of nonzeros in Jacobian and Heesian
   */
  virtual bool get_sparse_blocks_info(size_type& nx,
                                      size_type& nnz_sparse_Jaceq,
                                      size_type& nnz_sparse_Jacineq,
                                      size_type& nnz_sparse_Hess_Lagr) = 0;

  /** Evaluates the sparse Jacobian of constraints.
   *
   * This method is called twice per Jacobian evaluation, once for equalities and once for
   * inequalities (see 'eval_cons' for more information). It is advantageous to provide
   * this method when the underlying NLP's constraints come naturally split in equalities
   * and inequalities. When it is not convenient to do so, see the overloaded method.
   *
   * Parameters:
   *  - first six: see eval_cons (in parent class)
   *  - nnzJacS, iJacS, jJacS, MJacS: number of nonzeros, (i,j) indexes, and values of
   * the sparse Jacobian.
   *
   */
  virtual bool eval_Jac_cons(const size_type& n,
                             const size_type& m,
                             const size_type& num_cons,
                             const index_type* idx_cons,
                             const double* x,
                             bool new_x,
                             const size_type& nnzJacS,
                             index_type* iJacS,
                             index_type* jJacS,
                             double* MJacS) = 0;

  /** Evaluates the sparse Jacobian of equality and inequality constraints in one call.
   *
   * The main difference from the overloaded counterpart is that the implementer/user of this
   * method does not have to split the constraints into equalities and inequalities; instead,
   * HiOp does this internally.
   *
   * Parameters:
   *  - first four: number of variables, number of constraints, (primal) variables at which the
   * Jacobian should be evaluated, and boolean flag indicating whether the variables 'x' have
   * changed since a previous call to ny of the function and derivative evaluations.
   *  - nnzJacS, iJacS, jJacS, MJacS: number of nonzeros, (i,j) indexes, and values of
   * the sparse Jacobian block; indexes are within the sparse Jacobian block
   *
   * Notes for implementer of this method:
   * 1) When 'iJacS' and 'jJacS' are non-null, the implementer should provide the (i,j)
   * indexes.
   * 2) When 'MJacS' is non-null, the implementer should provide the values corresponding to
   * entries specified by 'iJacS' and 'jJacS'
   * 3) 'iJacS' and 'jJacS' are both either non-null or null during a call.
   * 4) Both 'iJacS'/'jJacS' and 'MJacS' can be non-null during the same call or only one of them
   * non-null; but they will not be both null.
   *
   * HiOp will call this method whenever the implementer/user returns false from the 'eval_Jac_cons'
   * (which is called for equalities and inequalities separately) above.
   */
  virtual bool eval_Jac_cons(const size_type& n,
                             const size_type& m,
                             const double* x,
                             bool new_x,
                             const size_type& nnzJacS,
                             index_type* iJacS,
                             index_type* jJacS,
                             double* MJacS)
  {
    return false;
  }

  /** Evaluates the sparse Hessian of the Lagrangian function.
   *
   * @note 1)-4) from 'eval_Jac_cons' applies to xxxHSS
   * @note 5) The order of multipliers is: lambda=[lambda_eq, lambda_ineq]
   */
  virtual bool eval_Hess_Lagr(const size_type& n,
                              const size_type& m,
                              const double* x,
                              bool new_x,
                              const double& obj_factor,
                              const double* lambda,
                              bool new_lambda,
                              const size_type& nnzHSS,
                              index_type* iHSS,
                              index_type* jHSS,
                              double* MHSS) = 0;

  /** Specifying the get_MPI_comm code defined in the base class
   */
  virtual bool get_MPI_comm(MPI_Comm& comm_out)
  {
    comm_out = MPI_COMM_SELF;
    return true;
  }
};

}  // namespace hiop
#endif