1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #ifndef MFEM_HIOP
13 #define MFEM_HIOP
14 
15 #include "linalg.hpp"
16 #include "../config/config.hpp"
17 #include "../general/globals.hpp"
18 
19 #ifdef MFEM_USE_MPI
20 #include "operator.hpp"
21 #endif
22 
23 #ifdef MFEM_USE_HIOP
24 
25 #include "hiopInterface.hpp"
26 #include "hiopNlpFormulation.hpp"
27 
28 namespace mfem
29 {
30 
31 /// Internal class - adapts the OptimizationProblem class to HiOp's interface.
32 class HiopOptimizationProblem : public hiop::hiopInterfaceDenseConstraints
33 {
34 private:
35 
36 #ifdef MFEM_USE_MPI
37    MPI_Comm comm;
38 #endif
39 
40    // Problem info.
41    const OptimizationProblem &problem;
42 
43    // Local and global number of variables and constraints.
44    const long long ntdofs_loc, m_total;
45    long long ntdofs_glob;
46 
47    // Initial guess.
48    const Vector *x_start;
49 
50    Vector constr_vals;
51    DenseMatrix constr_grads;
52    bool constr_info_is_current;
53    void UpdateConstrValsGrads(const Vector x);
54 
55 public:
HiopOptimizationProblem(const OptimizationProblem & prob)56    HiopOptimizationProblem(const OptimizationProblem &prob)
57       : problem(prob),
58         ntdofs_loc(prob.input_size), m_total(prob.GetNumConstraints()),
59         ntdofs_glob(ntdofs_loc),
60         x_start(NULL),
61         constr_vals(m_total), constr_grads(m_total, ntdofs_loc),
62         constr_info_is_current(false)
63    {
64 #ifdef MFEM_USE_MPI
65       // Used when HiOp with MPI support is called by a serial driver.
66       comm = MPI_COMM_WORLD;
67 #endif
68    }
69 
70 #ifdef MFEM_USE_MPI
HiopOptimizationProblem(const MPI_Comm & comm_,const OptimizationProblem & prob)71    HiopOptimizationProblem(const MPI_Comm& comm_,
72                            const OptimizationProblem &prob)
73       : comm(comm_),
74         problem(prob),
75         ntdofs_loc(prob.input_size), m_total(prob.GetNumConstraints()),
76         ntdofs_glob(0),
77         x_start(NULL),
78         constr_vals(m_total), constr_grads(m_total, ntdofs_loc),
79         constr_info_is_current(false)
80    {
81       MPI_Allreduce(&ntdofs_loc, &ntdofs_glob, 1, MPI_LONG_LONG_INT,
82                     MPI_SUM, comm);
83    }
84 #endif
85 
setStartingPoint(const Vector & x0)86    void setStartingPoint(const Vector &x0) { x_start = &x0; }
87 
88    /** Extraction of problem dimensions:
89     *  n is the number of variables, m is the number of constraints. */
90    virtual bool get_prob_sizes(long long int& n, long long int& m);
91 
92    /** Provide an primal starting point. This point is subject to adjustments
93     *  internally in HiOp. */
94    virtual bool get_starting_point(const long long &n, double *x0);
95 
96    virtual bool get_vars_info(const long long& n, double *xlow, double* xupp,
97                               NonlinearityType* type);
98 
99    /** bounds on the constraints
100     *  (clow<=-1e20 means no lower bound, cupp>=1e20 means no upper bound) */
101    virtual bool get_cons_info(const long long &m, double *clow, double *cupp,
102                               NonlinearityType* type);
103 
104    /** Objective function evaluation.
105     *  Each rank returns the global objective value. */
106    virtual bool eval_f(const long long& n, const double *x, bool new_x,
107                        double& obj_value);
108 
109    /** Gradient of the objective function (local chunk). */
110    virtual bool eval_grad_f(const long long &n, const double *x, bool new_x,
111                             double *gradf);
112 
113    /** Evaluates a subset of the constraints cons(x). The subset is of size
114     *  num_cons and is described by indexes in the idx_cons array,
115     *  i.e. cons[c] = C(x)[idx_cons[c]] where c = 0 .. num_cons-1.
116     *  The methods may be called multiple times, each time for a subset of the
117     *  constraints, for example, for the subset containing the equalities and
118     *  for the subset containing the inequalities. However, each constraint will
119     *  be inquired EXACTLY once. This is done for performance considerations,
120     *  to avoid temporary holders and memory copying.
121     *
122     *  Parameters:
123     *   - n, m: the global number of variables and constraints
124     *   - num_cons, idx_cons (array of size num_cons): the number and indexes of
125     *     constraints to be evaluated
126     *   - x: the point where the constraints are to be evaluated
127     *   - new_x: whether x has been changed from the previous call to f, grad_f,
128     *     or Jac
129     *   - cons: array of size num_cons containing the value of the  constraints
130     *     indicated by idx_cons
131     *
132     *  When MPI enabled, every rank populates cons, since the constraints are
133     *  not distributed.
134     */
135    virtual bool eval_cons(const long long &n, const long long &m,
136                           const long long &num_cons, const long long *idx_cons,
137                           const double *x, bool new_x, double *cons);
138 
139    /** Evaluates the Jacobian of the subset of constraints indicated by
140     *  idx_cons. The idx_cons is assumed to be of size num_cons.
141     *  Example: if cons[c] = C(x)[idx_cons[c]] where c = 0 .. num_cons-1, then
142     *  one needs to do Jac[c][j] = d cons[c] / dx_j, j = 1 .. n_loc.
143     *  Jac is computed and stored in a contiguous vector (offset by rows).
144     *
145     *  Parameters: see eval_cons().
146     *
147     *  When MPI enabled, each rank computes only the local columns of the
148     *  Jacobian, that is the partials with respect to local variables.
149     */
150    virtual bool eval_Jac_cons(const long long &n, const long long &m,
151                               const long long &num_cons,
152                               const long long *idx_cons,
153                               const double *x, bool new_x, double *Jac);
154 
155    /** Specifies column partitioning for distributed memory vectors.
156     *  Process p owns vector entries with indices cols[p] to cols[p+1]-1,
157     *  where p = 0 .. nranks-1. The cols array is of size nranks + 1.
158     *  Example: for a vector x of 6 entries (globally) on 3 ranks, the uniform
159     *  column partitioning is cols=[0,2,4,6].
160     */
161    virtual bool get_vecdistrib_info(long long global_n, long long *cols);
162 
163 #ifdef MFEM_USE_MPI
get_MPI_comm(MPI_Comm & comm_out)164    virtual bool get_MPI_comm(MPI_Comm &comm_out)
165    {
166       comm_out = comm;
167       return true;
168    }
169 #endif
170 };
171 
172 /// Adapts the HiOp functionality to the MFEM OptimizationSolver interface.
173 class HiopNlpOptimizer : public OptimizationSolver
174 {
175 protected:
176    HiopOptimizationProblem *hiop_problem;
177 
178 #ifdef MFEM_USE_MPI
179    MPI_Comm comm;
180 #endif
181 
182 public:
183    HiopNlpOptimizer();
184 #ifdef MFEM_USE_MPI
185    HiopNlpOptimizer(MPI_Comm comm_);
186 #endif
187    virtual ~HiopNlpOptimizer();
188 
189    virtual void SetOptimizationProblem(const OptimizationProblem &prob);
190 
191    /// Solves the optimization problem with xt as initial guess.
192    virtual void Mult(const Vector &xt, Vector &x) const;
193 };
194 
195 } // mfem namespace
196 
197 #endif //MFEM_USE_HIOP
198 #endif //MFEM_HIOP guard
199