-
Notifications
You must be signed in to change notification settings - Fork 0
/
unconstr_lbfgs_interface.cpp
163 lines (149 loc) · 4.77 KB
/
unconstr_lbfgs_interface.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
#include "interface.h"
#include <LBFGS.h>
using namespace LBFGSpp;
void unconstr_lbfgs_stat(CUTEstStat& stat, bool verbose)
{
using Vector = Eigen::Matrix<doublereal, Eigen::Dynamic, 1>;
// Open problem description file OUTSDIF.d
const char fname[] = "OUTSDIF.d";
// FORTRAN unit number for OUTSDIF.d
integer funit = 42;
// Exit flag from OPEN and CLOSE
integer ierr = 0;
FORTRAN_open(&funit, fname, &ierr);
if(ierr)
{
stat.flag = 2;
stat.msg = "Error opening file OUTSDIF.d.";
return;
}
// Determine problem size
// Exit flag from CUTEst tools
integer status;
// Number of variables
integer CUTEst_nvar;
// Number of general constraints
integer CUTEst_nconstr;
CUTEST_cdimen(&status, &funit, &CUTEst_nvar, &CUTEst_nconstr);
if(status)
{
stat.flag = 2;
stat.msg = "Error getting problem dimension.";
return;
}
if(CUTEst_nconstr > 0)
{
stat.flag = 2;
stat.msg = "Problem contains general constraints.";
return;
}
if(verbose)
std::cout << "nvar = " << CUTEst_nvar << std::endl;
// Reserve memory for variables and bounds,
// and call appropriate initialization routine for CUTEst
Vector x(CUTEst_nvar);
Vector lb(CUTEst_nvar);
Vector ub(CUTEst_nvar);
// FORTRAN unit number for error output
integer iout = 6;
// FORTRAN unit internal input/output
integer io_buffer = 11;
CUTEST_usetup(&status, &funit, &iout, &io_buffer,
&CUTEst_nvar, x.data(), lb.data(), ub.data());
if(status)
{
stat.flag = 2;
stat.msg = "Error setting up problem.";
CUTEST_uterminate(&status);
return;
}
// Even for unconstrained problems, lb and ub will be specified,
// but with a "fake" infinity value of +/- 1e20
// We need to make sure the problem is indeed unconstrained
const doublereal near_inf = 9.0e19;
if(lb.maxCoeff() > -near_inf || ub.minCoeff() < near_inf)
{
stat.flag = 2;
stat.msg = "Problem is not unconstrained.";
CUTEST_uterminate(&status);
return;
}
if(verbose)
std::cout << "x0 = " << x.transpose().head(10) << " ... " << x.transpose().tail(10) << std::endl;
// Problem name
char prob_name[16];
std::fill(prob_name, prob_name + 16, 0);
CUTEST_probname(&status, prob_name);
if(status)
{
stat.flag = 2;
stat.msg = "Error getting problem name.";
CUTEST_uterminate(&status);
return;
}
// Algorithm parameters
const integer param_m = 6;
// For very large problems, restrict to 1000 iterations
const integer param_maxit = (CUTEst_nvar < 50000) ? 10000 : 1000;
const doublereal param_eps = 1e-5;
// Machine precision
const doublereal param_xtol = std::numeric_limits<doublereal>::epsilon();
// Do not provide H0
const integer diagco = 0;
// Reuse the memory of lb
Vector& diag = lb;
// Objective function
CUTEstProblem fun(CUTEst_nvar);
doublereal fx;
Vector grad(CUTEst_nvar);
// Printing options
integer iprint[2];
iprint[0] = -1; // Do not print
iprint[1] = 0; // 0-3, larger value for more output
// Working space and flags
Vector work(CUTEst_nvar * (2 * param_m + 1) + 2 * param_m);
integer iflag = 0;
// Optimization process
integer i;
for (i = 0; i < param_maxit; i++)
{
// Compute objective function value and gradient
fx = fun(x, grad);
// Call L-BFGS routine
lbfgs_(&CUTEst_nvar, ¶m_m, x.data(), &fx, grad.data(),
&diagco, diag.data(), iprint, ¶m_eps, ¶m_xtol,
work.data(), &iflag);
// If iflag = 1, then continue iteration
if (iflag == 1)
{
continue;
} // If iflag = 0, then the solver finishes
else if (iflag == 0) {
break;
} // If iflag < 0, then some error occurs
else
{
// Errors
stat.prob = std::string(prob_name);
stat.nvar = CUTEst_nvar;
stat.flag = 1;
stat.msg = std::string("L-BFGS solver failed with code ") +
std::to_string(iflag);
CUTEST_uterminate(&status);
return;
}
}
doublereal calls[4], time[2];
CUTEST_ureport(&status, calls, time);
if(verbose)
std::cout << "x = " << x.transpose().head(5) << " ... " << x.transpose().tail(5) << std::endl << std::endl;
stat.prob = std::string(prob_name);
stat.nvar = CUTEst_nvar;
stat.niter = std::min(i + 1, param_maxit);
stat.nfun = calls[0];
stat.objval = fx;
stat.proj_grad = grad.norm();
stat.setup_time = time[0];
stat.solve_time = time[1];
CUTEST_uterminate(&status);
}