MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
hiop.cpp
Go to the documentation of this file.
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 #include "../config/config.hpp"
13 #include "hiop.hpp"
14 
15 #ifdef MFEM_USE_HIOP
16 #include <iostream>
17 
18 #include "hiopAlgFilterIPM.hpp"
19 
20 using namespace hiop;
21 
22 namespace mfem
23 {
24 
25 bool HiopOptimizationProblem::get_prob_sizes(long long &n, long long &m)
26 {
27  n = ntdofs_glob;
28  m = problem.GetNumConstraints();
29 
30  return true;
31 }
32 
33 bool HiopOptimizationProblem::get_starting_point(const long long &n, double *x0)
34 {
35  MFEM_ASSERT(x_start != NULL && ntdofs_loc == x_start->Size(),
36  "Starting point is not set properly.");
37 
38  memcpy(x0, x_start->GetData(), ntdofs_loc * sizeof(double));
39 
40  return true;
41 }
42 
43 bool HiopOptimizationProblem::get_vars_info(const long long &n,
44  double *xlow, double *xupp,
45  NonlinearityType *type)
46 {
47  MFEM_ASSERT(n == ntdofs_glob, "Global input mismatch.");
48  MFEM_ASSERT(problem.GetBoundsVec_Lo() && problem.GetBoundsVec_Hi(),
49  "Solution bounds are not set!");
50 
51  const int s = ntdofs_loc * sizeof(double);
52  std::memcpy(xlow, problem.GetBoundsVec_Lo()->GetData(), s);
53  std::memcpy(xupp, problem.GetBoundsVec_Hi()->GetData(), s);
54 
55  return true;
56 }
57 
58 bool HiopOptimizationProblem::get_cons_info(const long long &m,
59  double *clow, double *cupp,
60  NonlinearityType *type)
61 {
62  MFEM_ASSERT(m == m_total, "Global constraint size mismatch.");
63 
64  int csize = 0;
65  if (problem.GetC())
66  {
67  csize = problem.GetEqualityVec()->Size();
68  const int s = csize * sizeof(double);
69  std::memcpy(clow, problem.GetEqualityVec()->GetData(), s);
70  std::memcpy(cupp, problem.GetEqualityVec()->GetData(), s);
71  }
72  if (problem.GetD())
73  {
74  const int s = problem.GetInequalityVec_Lo()->Size() * sizeof(double);
75  std::memcpy(clow + csize, problem.GetInequalityVec_Lo()->GetData(), s);
76  std::memcpy(cupp + csize, problem.GetInequalityVec_Hi()->GetData(), s);
77  }
78 
79  return true;
80 }
81 
82 bool HiopOptimizationProblem::eval_f(const long long &n, const double *x,
83  bool new_x, double &obj_value)
84 {
85  MFEM_ASSERT(n == ntdofs_glob, "Global input mismatch.");
86 
87  if (new_x) { constr_info_is_current = false; }
88 
89  Vector x_vec(ntdofs_loc);
90  x_vec = x;
91  obj_value = problem.CalcObjective(x_vec);
92 
93  return true;
94 }
95 
96 bool HiopOptimizationProblem::eval_grad_f(const long long &n, const double *x,
97  bool new_x, double *gradf)
98 {
99  MFEM_ASSERT(n == ntdofs_glob, "Global input mismatch.");
100 
101  if (new_x) { constr_info_is_current = false; }
102 
103  Vector x_vec(ntdofs_loc), gradf_vec(ntdofs_loc);
104  x_vec = x;
105  problem.CalcObjectiveGrad(x_vec, gradf_vec);
106  std::memcpy(gradf, gradf_vec.GetData(), ntdofs_loc * sizeof(double));
107 
108  return true;
109 }
110 
111 bool HiopOptimizationProblem::eval_cons(const long long &n, const long long &m,
112  const long long &num_cons,
113  const long long *idx_cons,
114  const double *x, bool new_x,
115  double *cons)
116 {
117  MFEM_ASSERT(n == ntdofs_glob, "Global input mismatch.");
118  MFEM_ASSERT(m == m_total, "Constraint size mismatch.");
119  MFEM_ASSERT(num_cons <= m, "num_cons should be at most m = " << m);
120 
121  if (num_cons == 0) { return true; }
122 
123  if (new_x) { constr_info_is_current = false; }
124  Vector x_vec(ntdofs_loc);
125  x_vec = x;
126  UpdateConstrValsGrads(x_vec);
127 
128  for (int c = 0; c < num_cons; c++)
129  {
130  MFEM_ASSERT(idx_cons[c] < m_total, "Constraint index is out of bounds.");
131  cons[c] = constr_vals(idx_cons[c]);
132  }
133 
134  return true;
135 }
136 
137 bool HiopOptimizationProblem::eval_Jac_cons(const long long &n,
138  const long long &m,
139  const long long &num_cons,
140  const long long *idx_cons,
141  const double *x, bool new_x,
142  double *Jac)
143 {
144  MFEM_ASSERT(n == ntdofs_glob, "Global input mismatch.");
145  MFEM_ASSERT(m == m_total, "Constraint size mismatch.");
146  MFEM_ASSERT(num_cons <= m, "num_cons should be at most m = " << m);
147 
148  if (num_cons == 0) { return true; }
149 
150  if (new_x) { constr_info_is_current = false; }
151  Vector x_vec(ntdofs_loc);
152  x_vec = x;
153  UpdateConstrValsGrads(x_vec);
154 
155  for (int c = 0; c < num_cons; c++)
156  {
157  MFEM_ASSERT(idx_cons[c] < m_total, "Constraint index is out of bounds.");
158  for (int j = 0; j < ntdofs_loc; j++)
159  {
160  // The matrix is stored by rows.
161  Jac[c * ntdofs_loc + j] = constr_grads(idx_cons[c], j);
162  }
163  }
164 
165  return true;
166 }
167 
168 bool HiopOptimizationProblem::get_vecdistrib_info(long long global_n,
169  long long *cols)
170 {
171 #ifdef MFEM_USE_MPI
172  int nranks;
173  MPI_Comm_size(comm, &nranks);
174 
175  long long *sizes = new long long[nranks];
176  MPI_Allgather(&ntdofs_loc, 1, MPI_LONG_LONG_INT, sizes, 1,
177  MPI_LONG_LONG_INT, comm);
178  cols[0] = 0;
179  for (int r = 1; r <= nranks; r++)
180  {
181  cols[r] = sizes[r-1] + cols[r-1];
182  }
183 
184  delete [] sizes;
185  return true;
186 #else
187  // Returning false means that Hiop runs in non-distributed mode.
188  return false;
189 #endif
190 }
191 
192 void HiopOptimizationProblem::UpdateConstrValsGrads(const Vector x)
193 {
194  if (constr_info_is_current) { return; }
195 
196  // If needed (e.g. for CG spaces), communication should be handled by the
197  // operators' Mult() and GetGradient() methods.
198 
199  int cheight = 0;
200  if (problem.GetC())
201  {
202  cheight = problem.GetC()->Height();
203 
204  // Values of C.
205  Vector vals_C(constr_vals.GetData(), cheight);
206  problem.GetC()->Mult(x, vals_C);
207 
208  // Gradients C.
209  const Operator &oper_C = problem.GetC()->GetGradient(x);
210  const DenseMatrix *grad_C = dynamic_cast<const DenseMatrix *>(&oper_C);
211  MFEM_VERIFY(grad_C, "Hiop expects DenseMatrices as operator gradients.");
212  MFEM_ASSERT(grad_C->Height() == cheight && grad_C->Width() == ntdofs_loc,
213  "Incorrect dimensions of the C constraint gradient.");
214  for (int i = 0; i < cheight; i++)
215  {
216  for (int j = 0; j < ntdofs_loc; j++)
217  {
218  constr_grads(i, j) = (*grad_C)(i, j);
219  }
220  }
221  }
222 
223  if (problem.GetD())
224  {
225  const int dheight = problem.GetD()->Height();
226 
227  // Values of D.
228  Vector vals_D(constr_vals.GetData() + cheight, dheight);
229  problem.GetD()->Mult(x, vals_D);
230 
231  // Gradients of D.
232  const Operator &oper_D = problem.GetD()->GetGradient(x);
233  const DenseMatrix *grad_D = dynamic_cast<const DenseMatrix *>(&oper_D);
234  MFEM_VERIFY(grad_D, "Hiop expects DenseMatrices as operator gradients.");
235  MFEM_ASSERT(grad_D->Height() == dheight && grad_D->Width() == ntdofs_loc,
236  "Incorrect dimensions of the D constraint gradient.");
237  for (int i = 0; i < dheight; i++)
238  {
239  for (int j = 0; j < ntdofs_loc; j++)
240  {
241  constr_grads(i + cheight, j) = (*grad_D)(i, j);
242  }
243  }
244  }
245 
246  constr_info_is_current = true;
247 }
248 
249 HiopNlpOptimizer::HiopNlpOptimizer() : OptimizationSolver(), hiop_problem(NULL)
250 {
251 #ifdef MFEM_USE_MPI
252  // Set in case a serial driver uses a parallel MFEM build.
253  comm = MPI_COMM_WORLD;
254  int initialized, nret = MPI_Initialized(&initialized);
255  MFEM_ASSERT(MPI_SUCCESS == nret, "Failure in calling MPI_Initialized!");
256  if (!initialized)
257  {
258  nret = MPI_Init(NULL, NULL);
259  MFEM_ASSERT(MPI_SUCCESS == nret, "Failure in calling MPI_Init!");
260  }
261 #endif
262 }
263 
264 #ifdef MFEM_USE_MPI
266  : OptimizationSolver(comm_), hiop_problem(NULL), comm(comm_) { }
267 #endif
268 
270 {
271  delete hiop_problem;
272 }
273 
275 {
276  problem = &prob;
278 
279  if (hiop_problem) { delete hiop_problem; }
280 
281 #ifdef MFEM_USE_MPI
283 #else
285 #endif
286 }
287 
288 void HiopNlpOptimizer::Mult(const Vector &xt, Vector &x) const
289 {
290  MFEM_ASSERT(hiop_problem != NULL,
291  "Unspecified OptimizationProblem that must be solved.");
292 
294 
295  hiop::hiopNlpDenseConstraints hiopInstance(*hiop_problem);
296 
297  hiopInstance.options->SetNumericValue("rel_tolerance", rel_tol);
298  hiopInstance.options->SetNumericValue("tolerance", abs_tol);
299  hiopInstance.options->SetIntegerValue("max_iter", max_iter);
300 
301  hiopInstance.options->SetStringValue("fixed_var", "relax");
302  hiopInstance.options->SetNumericValue("fixed_var_tolerance", 1e-20);
303  hiopInstance.options->SetNumericValue("fixed_var_perturb", 1e-9);
304 
305  // 0: no output; 3: not too much
306  hiopInstance.options->SetIntegerValue("verbosity_level", print_level);
307 
308  // Use the IPM solver.
309  hiop::hiopAlgFilterIPM solver(&hiopInstance);
310  const hiop::hiopSolveStatus status = solver.run();
311  final_norm = solver.getObjective();
312  final_iter = solver.getNumIterations();
313 
314  if (status != hiop::Solve_Success && status != hiop::Solve_Success_RelTol)
315  {
316  converged = false;
317  MFEM_WARNING("HIOP returned with a non-success status: " << status);
318  }
319  else { converged = true; }
320 
321  // Copy the final solution in x.
322  solver.getSolution(x.GetData());
323 }
324 
325 } // mfem namespace
326 #endif // MFEM_USE_HIOP
virtual void SetOptimizationProblem(const OptimizationProblem &prob)
Definition: hiop.cpp:274
void setStartingPoint(const Vector &x0)
Definition: hiop.hpp:86
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:199
virtual void Mult(const Vector &xt, Vector &x) const
Solves the optimization problem with xt as initial guess.
Definition: hiop.cpp:288
virtual ~HiopNlpOptimizer()
Definition: hiop.cpp:269
prob_type prob
Definition: ex25.cpp:153
int problem
Definition: ex15.cpp:62
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
Abstract solver for OptimizationProblems.
Definition: solvers.hpp:634
Vector data type.
Definition: vector.hpp:60
RefCoord s[3]
const OptimizationProblem * problem
Definition: solvers.hpp:637
HiopOptimizationProblem * hiop_problem
Definition: hiop.hpp:176
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
Internal class - adapts the OptimizationProblem class to HiOp&#39;s interface.
Definition: hiop.hpp:32