MFEM  v4.1.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-2020, 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  Jac[c][j] = constr_grads(idx_cons[c], j);
161  }
162  }
163 
164  return true;
165 }
166 
167 bool HiopOptimizationProblem::get_vecdistrib_info(long long global_n,
168  long long *cols)
169 {
170 #ifdef MFEM_USE_MPI
171  int nranks;
172  MPI_Comm_size(comm_, &nranks);
173 
174  long long *sizes = new long long[nranks];
175  MPI_Allgather(&ntdofs_loc, 1, MPI_LONG_LONG_INT, sizes, 1,
176  MPI_LONG_LONG_INT, comm_);
177  cols[0] = 0;
178  for (int r = 1; r <= nranks; r++)
179  {
180  cols[r] = sizes[r-1] + cols[r-1];
181  }
182 
183  delete [] sizes;
184  return true;
185 #else
186  // Returning false means that Hiop runs in non-distributed mode.
187  return false;
188 #endif
189 }
190 
191 void HiopOptimizationProblem::UpdateConstrValsGrads(const Vector x)
192 {
193  if (constr_info_is_current) { return; }
194 
195  // If needed (e.g. for CG spaces), communication should be handled by the
196  // operators' Mult() and GetGradient() methods.
197 
198  int cheight = 0;
199  if (problem.GetC())
200  {
201  cheight = problem.GetC()->Height();
202 
203  // Values of C.
204  Vector vals_C(constr_vals.GetData(), cheight);
205  problem.GetC()->Mult(x, vals_C);
206 
207  // Gradients C.
208  const Operator &oper_C = problem.GetC()->GetGradient(x);
209  const DenseMatrix *grad_C = dynamic_cast<const DenseMatrix *>(&oper_C);
210  MFEM_VERIFY(grad_C, "Hiop expects DenseMatrices as operator gradients.");
211  MFEM_ASSERT(grad_C->Height() == cheight && grad_C->Width() == ntdofs_loc,
212  "Incorrect dimensions of the C constraint gradient.");
213  for (int i = 0; i < cheight; i++)
214  {
215  for (int j = 0; j < ntdofs_loc; j++)
216  {
217  constr_grads(i, j) = (*grad_C)(i, j);
218  }
219  }
220  }
221 
222  if (problem.GetD())
223  {
224  const int dheight = problem.GetD()->Height();
225 
226  // Values of D.
227  Vector vals_D(constr_vals.GetData() + cheight, dheight);
228  problem.GetD()->Mult(x, vals_D);
229 
230  // Gradients of D.
231  const Operator &oper_D = problem.GetD()->GetGradient(x);
232  const DenseMatrix *grad_D = dynamic_cast<const DenseMatrix *>(&oper_D);
233  MFEM_VERIFY(grad_D, "Hiop expects DenseMatrices as operator gradients.");
234  MFEM_ASSERT(grad_D->Height() == dheight && grad_D->Width() == ntdofs_loc,
235  "Incorrect dimensions of the D constraint gradient.");
236  for (int i = 0; i < dheight; i++)
237  {
238  for (int j = 0; j < ntdofs_loc; j++)
239  {
240  constr_grads(i + cheight, j) = (*grad_D)(i, j);
241  }
242  }
243  }
244 
245  constr_info_is_current = true;
246 }
247 
248 HiopNlpOptimizer::HiopNlpOptimizer() : OptimizationSolver(), hiop_problem(NULL)
249 {
250 #ifdef MFEM_USE_MPI
251  // Set in case a serial driver uses a parallel MFEM build.
252  comm_ = MPI_COMM_WORLD;
253  int initialized, nret = MPI_Initialized(&initialized);
254  MFEM_ASSERT(MPI_SUCCESS == nret, "Failure in calling MPI_Initialized!");
255  if (!initialized)
256  {
257  nret = MPI_Init(NULL, NULL);
258  MFEM_ASSERT(MPI_SUCCESS == nret, "Failure in calling MPI_Init!");
259  }
260 #endif
261 }
262 
263 #ifdef MFEM_USE_MPI
265  : OptimizationSolver(_comm), hiop_problem(NULL), comm_(_comm) { }
266 #endif
267 
269 {
270  delete hiop_problem;
271 }
272 
274 {
275  problem = &prob;
277 
278  if (hiop_problem) { delete hiop_problem; }
279 
280 #ifdef MFEM_USE_MPI
282 #else
284 #endif
285 }
286 
287 void HiopNlpOptimizer::Mult(const Vector &xt, Vector &x) const
288 {
289  MFEM_ASSERT(hiop_problem != NULL,
290  "Unspecified OptimizationProblem that must be solved.");
291 
293 
294  hiop::hiopNlpDenseConstraints hiopInstance(*hiop_problem);
295 
296  hiopInstance.options->SetNumericValue("rel_tolerance", rel_tol);
297  hiopInstance.options->SetNumericValue("tolerance", abs_tol);
298  hiopInstance.options->SetIntegerValue("max_iter", max_iter);
299 
300  hiopInstance.options->SetStringValue("fixed_var", "relax");
301  hiopInstance.options->SetNumericValue("fixed_var_tolerance", 1e-20);
302  hiopInstance.options->SetNumericValue("fixed_var_perturb", 1e-9);
303 
304  // 0: no output; 3: not too much
305  hiopInstance.options->SetIntegerValue("verbosity_level", print_level);
306 
307  // Use the IPM solver.
308  hiop::hiopAlgFilterIPM solver(&hiopInstance);
309  const hiop::hiopSolveStatus status = solver.run();
310  final_norm = solver.getObjective();
311  final_iter = solver.getNumIterations();
312 
313  if (status != hiop::Solve_Success && status != hiop::Solve_Success_RelTol)
314  {
315  converged = false;
316  MFEM_WARNING("HIOP returned with a non-success status: " << status);
317  }
318  else { converged = true; }
319 
320  // Copy the final solution in x.
321  solver.getSolution(x.GetData());
322 }
323 
324 } // mfem namespace
325 #endif // MFEM_USE_HIOP
virtual void SetOptimizationProblem(const OptimizationProblem &prob)
Definition: hiop.cpp:273
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:166
virtual void Mult(const Vector &xt, Vector &x) const
Solves the optimization problem with xt as initial guess.
Definition: hiop.cpp:287
virtual ~HiopNlpOptimizer()
Definition: hiop.cpp:268
int problem
Definition: ex15.cpp:54
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
Abstract solver for OptimizationProblems.
Definition: solvers.hpp:394
Vector data type.
Definition: vector.hpp:48
const OptimizationProblem * problem
Definition: solvers.hpp:397
HiopOptimizationProblem * hiop_problem
Definition: hiop.hpp:175
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