lemon-project-template-glpk: 33de93886c88 deps/glpk/src/glpapi08.c

lemon-project-template-glpk

view deps/glpk/src/glpapi08.c @ 9:33de93886c88

Import GLPK 4.47

author	Alpar Juttner <alpar@cs.elte.hu>
date	Sun, 06 Nov 2011 20:59:10 +0100
parents
children

line source

1 /* glpapi08.c (interior-point method routines) */

3 /***********************************************************************

4 * This code is part of GLPK (GNU Linear Programming Kit).

5 *

7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,

9 * E-mail: <mao@gnu.org>.

10 *

11 * GLPK is free software: you can redistribute it and/or modify it

12 * under the terms of the GNU General Public License as published by

13 * the Free Software Foundation, either version 3 of the License, or

14 * (at your option) any later version.

15 *

16 * GLPK is distributed in the hope that it will be useful, but WITHOUT

17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY

18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public

19 * License for more details.

20 *

21 * You should have received a copy of the GNU General Public License

22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.

23 ***********************************************************************/

25 #include "glpapi.h"

26 #include "glpipm.h"

27 #include "glpnpp.h"

29 /***********************************************************************

30 * NAME

31 *

32 * glp_interior - solve LP problem with the interior-point method

33 *

34 * SYNOPSIS

35 *

36 * int glp_interior(glp_prob *P, const glp_iptcp *parm);

37 *

38 * The routine glp_interior is a driver to the LP solver based on the

39 * interior-point method.

40 *

41 * The interior-point solver has a set of control parameters. Values of

42 * the control parameters can be passed in a structure glp_iptcp, which

43 * the parameter parm points to.

44 *

45 * Currently this routine implements an easy variant of the primal-dual

46 * interior-point method based on Mehrotra's technique.

47 *

48 * This routine transforms the original LP problem to an equivalent LP

49 * problem in the standard formulation (all constraints are equalities,

50 * all variables are non-negative), calls the routine ipm_main to solve

51 * the transformed problem, and then transforms an obtained solution to

52 * the solution of the original problem.

53 *

54 * RETURNS

55 *

56 * 0 The LP problem instance has been successfully solved. This code

57 * does not necessarily mean that the solver has found optimal

58 * solution. It only means that the solution process was successful.

59 *

60 * GLP_EFAIL

61 * The problem has no rows/columns.

62 *

63 * GLP_ENOCVG

64 * Very slow convergence or divergence.

65 *

66 * GLP_EITLIM

67 * Iteration limit exceeded.

68 *

69 * GLP_EINSTAB

70 * Numerical instability on solving Newtonian system. */

72 static void transform(NPP *npp)

73 { /* transform LP to the standard formulation */

74 NPPROW *row, *prev_row;

75 NPPCOL *col, *prev_col;

76 for (row = npp->r_tail; row != NULL; row = prev_row)

77 { prev_row = row->prev;

78 if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)

79 npp_free_row(npp, row);

80 else if (row->lb == -DBL_MAX)

81 npp_leq_row(npp, row);

82 else if (row->ub == +DBL_MAX)

83 npp_geq_row(npp, row);

84 else if (row->lb != row->ub)

85 { if (fabs(row->lb) < fabs(row->ub))

86 npp_geq_row(npp, row);

87 else

88 npp_leq_row(npp, row);

89 }

90 }

91 for (col = npp->c_tail; col != NULL; col = prev_col)

92 { prev_col = col->prev;

93 if (col->lb == -DBL_MAX && col->ub == +DBL_MAX)

94 npp_free_col(npp, col);

95 else if (col->lb == -DBL_MAX)

96 npp_ubnd_col(npp, col);

97 else if (col->ub == +DBL_MAX)

98 { if (col->lb != 0.0)

99 npp_lbnd_col(npp, col);

100 }

101 else if (col->lb != col->ub)

102 { if (fabs(col->lb) < fabs(col->ub))

103 { if (col->lb != 0.0)

104 npp_lbnd_col(npp, col);

105 }

106 else

107 npp_ubnd_col(npp, col);

108 npp_dbnd_col(npp, col);

109 }

110 else

111 npp_fixed_col(npp, col);

112 }

113 for (row = npp->r_head; row != NULL; row = row->next)

114 xassert(row->lb == row->ub);

115 for (col = npp->c_head; col != NULL; col = col->next)

116 xassert(col->lb == 0.0 && col->ub == +DBL_MAX);

117 return;

118 }

119

120 int glp_interior(glp_prob *P, const glp_iptcp *parm)

121 { glp_iptcp _parm;

122 GLPROW *row;

123 GLPCOL *col;

124 NPP *npp = NULL;

125 glp_prob *prob = NULL;

126 int i, j, ret;

127 /* check control parameters */

128 if (parm == NULL)

129 glp_init_iptcp(&_parm), parm = &_parm;

130 if (!(parm->msg_lev == GLP_MSG_OFF ||

131 parm->msg_lev == GLP_MSG_ERR ||

132 parm->msg_lev == GLP_MSG_ON ||

133 parm->msg_lev == GLP_MSG_ALL))

134 xerror("glp_interior: msg_lev = %d; invalid parameter\n",

135 parm->msg_lev);

136 if (!(parm->ord_alg == GLP_ORD_NONE ||

137 parm->ord_alg == GLP_ORD_QMD ||

138 parm->ord_alg == GLP_ORD_AMD ||

139 parm->ord_alg == GLP_ORD_SYMAMD))

140 xerror("glp_interior: ord_alg = %d; invalid parameter\n",

141 parm->ord_alg);

142 /* interior-point solution is currently undefined */

143 P->ipt_stat = GLP_UNDEF;

144 P->ipt_obj = 0.0;

145 /* check bounds of double-bounded variables */

146 for (i = 1; i <= P->m; i++)

147 { row = P->row[i];

148 if (row->type == GLP_DB && row->lb >= row->ub)

149 { if (parm->msg_lev >= GLP_MSG_ERR)

150 xprintf("glp_interior: row %d: lb = %g, ub = %g; incorre"

151 "ct bounds\n", i, row->lb, row->ub);

152 ret = GLP_EBOUND;

153 goto done;

154 }

155 }

156 for (j = 1; j <= P->n; j++)

157 { col = P->col[j];

158 if (col->type == GLP_DB && col->lb >= col->ub)

159 { if (parm->msg_lev >= GLP_MSG_ERR)

160 xprintf("glp_interior: column %d: lb = %g, ub = %g; inco"

161 "rrect bounds\n", j, col->lb, col->ub);

162 ret = GLP_EBOUND;

163 goto done;

164 }

165 }

166 /* transform LP to the standard formulation */

167 if (parm->msg_lev >= GLP_MSG_ALL)

168 xprintf("Original LP has %d row(s), %d column(s), and %d non-z"

169 "ero(s)\n", P->m, P->n, P->nnz);

170 npp = npp_create_wksp();

171 npp_load_prob(npp, P, GLP_OFF, GLP_IPT, GLP_ON);

172 transform(npp);

173 prob = glp_create_prob();

174 npp_build_prob(npp, prob);

175 if (parm->msg_lev >= GLP_MSG_ALL)

176 xprintf("Working LP has %d row(s), %d column(s), and %d non-ze"

177 "ro(s)\n", prob->m, prob->n, prob->nnz);

178 #if 1

179 /* currently empty problem cannot be solved */

180 if (!(prob->m > 0 && prob->n > 0))

181 { if (parm->msg_lev >= GLP_MSG_ERR)

182 xprintf("glp_interior: unable to solve empty problem\n");

183 ret = GLP_EFAIL;

184 goto done;

185 }

186 #endif

187 /* scale the resultant LP */

188 { ENV *env = get_env_ptr();

189 int term_out = env->term_out;

190 env->term_out = GLP_OFF;

191 glp_scale_prob(prob, GLP_SF_EQ);

192 env->term_out = term_out;

193 }

194 /* warn about dense columns */

195 if (parm->msg_lev >= GLP_MSG_ON && prob->m >= 200)

196 { int len, cnt = 0;

197 for (j = 1; j <= prob->n; j++)

198 { len = glp_get_mat_col(prob, j, NULL, NULL);

199 if ((double)len >= 0.20 * (double)prob->m) cnt++;

200 }

201 if (cnt == 1)

202 xprintf("WARNING: PROBLEM HAS ONE DENSE COLUMN\n");

203 else if (cnt > 0)

204 xprintf("WARNING: PROBLEM HAS %d DENSE COLUMNS\n", cnt);

205 }

206 /* solve the transformed LP */

207 ret = ipm_solve(prob, parm);

208 /* postprocess solution from the transformed LP */

209 npp_postprocess(npp, prob);

210 /* and store solution to the original LP */

211 npp_unload_sol(npp, P);

212 done: /* free working program objects */

213 if (npp != NULL) npp_delete_wksp(npp);

214 if (prob != NULL) glp_delete_prob(prob);

215 /* return to the application program */

216 return ret;

217 }

218

219 /***********************************************************************

220 * NAME

221 *

222 * glp_init_iptcp - initialize interior-point solver control parameters

223 *

224 * SYNOPSIS

225 *

226 * void glp_init_iptcp(glp_iptcp *parm);

227 *

228 * DESCRIPTION

229 *

230 * The routine glp_init_iptcp initializes control parameters, which are

231 * used by the interior-point solver, with default values.

232 *

233 * Default values of the control parameters are stored in the glp_iptcp

234 * structure, which the parameter parm points to. */

235

236 void glp_init_iptcp(glp_iptcp *parm)

237 { parm->msg_lev = GLP_MSG_ALL;

238 parm->ord_alg = GLP_ORD_AMD;

239 return;

240 }

241

242 /***********************************************************************

243 * NAME

244 *

245 * glp_ipt_status - retrieve status of interior-point solution

246 *

247 * SYNOPSIS

248 *

249 * int glp_ipt_status(glp_prob *lp);

250 *

251 * RETURNS

252 *

253 * The routine glp_ipt_status reports the status of solution found by

254 * the interior-point solver as follows:

255 *

256 * GLP_UNDEF - interior-point solution is undefined;

257 * GLP_OPT - interior-point solution is optimal;

258 * GLP_INFEAS - interior-point solution is infeasible;

259 * GLP_NOFEAS - no feasible solution exists. */

260

261 int glp_ipt_status(glp_prob *lp)

262 { int ipt_stat = lp->ipt_stat;

263 return ipt_stat;

264 }

265

266 /***********************************************************************

267 * NAME

268 *

269 * glp_ipt_obj_val - retrieve objective value (interior point)

270 *

271 * SYNOPSIS

272 *

273 * double glp_ipt_obj_val(glp_prob *lp);

274 *

275 * RETURNS

276 *

277 * The routine glp_ipt_obj_val returns value of the objective function

278 * for interior-point solution. */

279

280 double glp_ipt_obj_val(glp_prob *lp)

281 { /*struct LPXCPS *cps = lp->cps;*/

282 double z;

283 z = lp->ipt_obj;

284 /*if (cps->round && fabs(z) < 1e-9) z = 0.0;*/

285 return z;

286 }

287

288 /***********************************************************************

289 * NAME

290 *

291 * glp_ipt_row_prim - retrieve row primal value (interior point)

292 *

293 * SYNOPSIS

294 *

295 * double glp_ipt_row_prim(glp_prob *lp, int i);

296 *

297 * RETURNS

298 *

299 * The routine glp_ipt_row_prim returns primal value of the auxiliary

300 * variable associated with i-th row. */

301

302 double glp_ipt_row_prim(glp_prob *lp, int i)

303 { /*struct LPXCPS *cps = lp->cps;*/

304 double pval;

305 if (!(1 <= i && i <= lp->m))

306 xerror("glp_ipt_row_prim: i = %d; row number out of range\n",

307 i);

308 pval = lp->row[i]->pval;

309 /*if (cps->round && fabs(pval) < 1e-9) pval = 0.0;*/

310 return pval;

311 }

312

313 /***********************************************************************

314 * NAME

315 *

316 * glp_ipt_row_dual - retrieve row dual value (interior point)

317 *

318 * SYNOPSIS

319 *

320 * double glp_ipt_row_dual(glp_prob *lp, int i);

321 *

322 * RETURNS

323 *

324 * The routine glp_ipt_row_dual returns dual value (i.e. reduced cost)

325 * of the auxiliary variable associated with i-th row. */

326

327 double glp_ipt_row_dual(glp_prob *lp, int i)

328 { /*struct LPXCPS *cps = lp->cps;*/

329 double dval;

330 if (!(1 <= i && i <= lp->m))

331 xerror("glp_ipt_row_dual: i = %d; row number out of range\n",

332 i);

333 dval = lp->row[i]->dval;

334 /*if (cps->round && fabs(dval) < 1e-9) dval = 0.0;*/

335 return dval;

336 }

337

338 /***********************************************************************

339 * NAME

340 *

341 * glp_ipt_col_prim - retrieve column primal value (interior point)

342 *

343 * SYNOPSIS

344 *

345 * double glp_ipt_col_prim(glp_prob *lp, int j);

346 *

347 * RETURNS

348 *

349 * The routine glp_ipt_col_prim returns primal value of the structural

350 * variable associated with j-th column. */

351

352 double glp_ipt_col_prim(glp_prob *lp, int j)

353 { /*struct LPXCPS *cps = lp->cps;*/

354 double pval;

355 if (!(1 <= j && j <= lp->n))

356 xerror("glp_ipt_col_prim: j = %d; column number out of range\n"

357 , j);

358 pval = lp->col[j]->pval;

359 /*if (cps->round && fabs(pval) < 1e-9) pval = 0.0;*/

360 return pval;

361 }

362

363 /***********************************************************************

364 * NAME

365 *

366 * glp_ipt_col_dual - retrieve column dual value (interior point)

367 *

368 * SYNOPSIS

369 *

370 * #include "glplpx.h"

371 * double glp_ipt_col_dual(glp_prob *lp, int j);

372 *

373 * RETURNS

374 *

375 * The routine glp_ipt_col_dual returns dual value (i.e. reduced cost)

376 * of the structural variable associated with j-th column. */

377

378 double glp_ipt_col_dual(glp_prob *lp, int j)

379 { /*struct LPXCPS *cps = lp->cps;*/

380 double dval;

381 if (!(1 <= j && j <= lp->n))

382 xerror("glp_ipt_col_dual: j = %d; column number out of range\n"

383 , j);

384 dval = lp->col[j]->dval;

385 /*if (cps->round && fabs(dval) < 1e-9) dval = 0.0;*/

386 return dval;

387 }

388

389 /* eof */