-
Notifications
You must be signed in to change notification settings - Fork 0
/
boxconstr_lbfgsb_interface.cpp
175 lines (159 loc) · 5.21 KB
/
boxconstr_lbfgsb_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
164
165
166
167
168
169
170
171
172
173
174
175
// Copyright (C) 2023 Yixuan Qiu <yixuan.qiu@cos.name>
// Under MIT license
#include "interface.h"
void boxconstr_lbfgsb_stat(CUTEstStat& stat, bool verbose)
{
using Vector = Eigen::Matrix<doublereal, Eigen::Dynamic, 1>;
using IntVector = Eigen::VectorXi;
// 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 constrints.";
return;
}
if(verbose)
std::cout << "nvar = " << CUTEst_nvar << std::endl;
// Reserve memory for variables, bounds, and multipliers,
// 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;
}
if(verbose)
{
std::cout << "x0 = " << x.transpose().head(10) << " ... " << x.transpose().tail(10) << std::endl;
std::cout << "lb = " << lb.transpose().head(10) << " ... " << lb.transpose().tail(10) << std::endl;
std::cout << "ub = " << ub.transpose().head(10) << " ... " << ub.transpose().tail(10) << std::endl << 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;
const integer param_maxit = 10000;
// const doublereal param_factr = 0.0;
const doublereal param_factr = 1e7;
const doublereal param_pgtol = 1e-5;
// Bound type indicators
IntVector nbd(CUTEst_nvar);
const doublereal inf = std::numeric_limits<double>::infinity();
for(int i = 0; i < CUTEst_nvar; i++)
{
if(lb[i] == -inf)
{
nbd[i] = (ub[i] == inf) ? 0 : 3;
} else {
nbd[i] = (ub[i] == inf) ? 1 : 2;
}
}
// Objective function
CUTEstProblem fun(CUTEst_nvar);
doublereal fx;
Vector grad(CUTEst_nvar);
// Printing options
// Do not print
const integer iprint = -1;
// Working space and flags
Vector wa(2 * param_m * CUTEst_nvar + 11 * param_m * param_m + 5 * CUTEst_nvar + 8 * param_m);
IntVector iwa(3 * CUTEst_nvar);
integer itask;
integer icsave;
integer lsave[4];
integer isave[44];
double dsave[29];
// Optimization process
itask = 2;
int i = 0;
while (i < param_maxit)
{
// Call L-BFGS-B routine
setulb_(&CUTEst_nvar, ¶m_m, x.data(), lb.data(), ub.data(), nbd.data(),
&fx, grad.data(), ¶m_factr, ¶m_pgtol,
wa.data(), iwa.data(), &itask, &iprint,
&icsave, lsave, isave, dsave);
// std::cout << "i = " << i << ", itask = " << itask << std::endl;
if (itask == 4 || itask == 20 || itask == 21)
{
// Compute objective function value and gradient
fx = fun(x, grad);
// std::cout << " x = " << x.transpose() << std::endl;
// std::cout << " grad = " << grad.transpose() << std::endl;
// std::cout << " fx = " << fx << std::endl;
} else if (itask >= 6 && itask <= 8) {
// Converged
break;
} else if (itask == 1) {
// New x, update iteration number
i = isave[29];
} else {
// Errors
stat.prob = std::string(prob_name);
stat.nvar = CUTEst_nvar;
stat.flag = 1;
stat.msg = std::string("Solver abnormal exit. itask = ") +
std::to_string(itask);
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 = i;
stat.nfun = calls[0];
stat.objval = fx;
stat.proj_grad = dsave[12];
stat.setup_time = time[0];
stat.solve_time = time[1];
CUTEST_uterminate(&status);
}