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

lemon-project-template-glpk

view deps/glpk/src/glpapi10.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 /* glpapi10.c (solution checking 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"

27 void _glp_check_kkt(glp_prob *P, int sol, int cond, double *_ae_max,

28 int *_ae_ind, double *_re_max, int *_re_ind)

29 { /* check feasibility and optimality conditions */

30 int m = P->m;

31 int n = P->n;

32 GLPROW *row;

33 GLPCOL *col;

34 GLPAIJ *aij;

35 int i, j, ae_ind, re_ind;

36 double e, sp, sn, t, ae_max, re_max;

37 if (!(sol == GLP_SOL || sol == GLP_IPT || sol == GLP_MIP))

38 xerror("glp_check_kkt: sol = %d; invalid solution indicator\n",

39 sol);

40 if (!(cond == GLP_KKT_PE || cond == GLP_KKT_PB ||

41 cond == GLP_KKT_DE || cond == GLP_KKT_DB ||

42 cond == GLP_KKT_CS))

43 xerror("glp_check_kkt: cond = %d; invalid condition indicator "

44 "\n", cond);

45 ae_max = re_max = 0.0;

46 ae_ind = re_ind = 0;

47 if (cond == GLP_KKT_PE)

48 { /* xR - A * xS = 0 */

49 for (i = 1; i <= m; i++)

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

51 sp = sn = 0.0;

52 /* t := xR[i] */

53 if (sol == GLP_SOL)

54 t = row->prim;

55 else if (sol == GLP_IPT)

56 t = row->pval;

57 else if (sol == GLP_MIP)

58 t = row->mipx;

59 else

60 xassert(sol != sol);

61 if (t >= 0.0) sp += t; else sn -= t;

62 for (aij = row->ptr; aij != NULL; aij = aij->r_next)

63 { col = aij->col;

64 /* t := - a[i,j] * xS[j] */

65 if (sol == GLP_SOL)

66 t = - aij->val * col->prim;

67 else if (sol == GLP_IPT)

68 t = - aij->val * col->pval;

69 else if (sol == GLP_MIP)

70 t = - aij->val * col->mipx;

71 else

72 xassert(sol != sol);

73 if (t >= 0.0) sp += t; else sn -= t;

74 }

75 /* absolute error */

76 e = fabs(sp - sn);

77 if (ae_max < e)

78 ae_max = e, ae_ind = i;

79 /* relative error */

80 e /= (1.0 + sp + sn);

81 if (re_max < e)

82 re_max = e, re_ind = i;

83 }

84 }

85 else if (cond == GLP_KKT_PB)

86 { /* lR <= xR <= uR */

87 for (i = 1; i <= m; i++)

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

89 /* t := xR[i] */

90 if (sol == GLP_SOL)

91 t = row->prim;

92 else if (sol == GLP_IPT)

93 t = row->pval;

94 else if (sol == GLP_MIP)

95 t = row->mipx;

96 else

97 xassert(sol != sol);

98 /* check lower bound */

99 if (row->type == GLP_LO || row->type == GLP_DB ||

100 row->type == GLP_FX)

101 { if (t < row->lb)

102 { /* absolute error */

103 e = row->lb - t;

104 if (ae_max < e)

105 ae_max = e, ae_ind = i;

106 /* relative error */

107 e /= (1.0 + fabs(row->lb));

108 if (re_max < e)

109 re_max = e, re_ind = i;

110 }

111 }

112 /* check upper bound */

113 if (row->type == GLP_UP || row->type == GLP_DB ||

114 row->type == GLP_FX)

115 { if (t > row->ub)

116 { /* absolute error */

117 e = t - row->ub;

118 if (ae_max < e)

119 ae_max = e, ae_ind = i;

120 /* relative error */

121 e /= (1.0 + fabs(row->ub));

122 if (re_max < e)

123 re_max = e, re_ind = i;

124 }

125 }

126 }

127 /* lS <= xS <= uS */

128 for (j = 1; j <= n; j++)

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

130 /* t := xS[j] */

131 if (sol == GLP_SOL)

132 t = col->prim;

133 else if (sol == GLP_IPT)

134 t = col->pval;

135 else if (sol == GLP_MIP)

136 t = col->mipx;

137 else

138 xassert(sol != sol);

139 /* check lower bound */

140 if (col->type == GLP_LO || col->type == GLP_DB ||

141 col->type == GLP_FX)

142 { if (t < col->lb)

143 { /* absolute error */

144 e = col->lb - t;

145 if (ae_max < e)

146 ae_max = e, ae_ind = m+j;

147 /* relative error */

148 e /= (1.0 + fabs(col->lb));

149 if (re_max < e)

150 re_max = e, re_ind = m+j;

151 }

152 }

153 /* check upper bound */

154 if (col->type == GLP_UP || col->type == GLP_DB ||

155 col->type == GLP_FX)

156 { if (t > col->ub)

157 { /* absolute error */

158 e = t - col->ub;

159 if (ae_max < e)

160 ae_max = e, ae_ind = m+j;

161 /* relative error */

162 e /= (1.0 + fabs(col->ub));

163 if (re_max < e)

164 re_max = e, re_ind = m+j;

165 }

166 }

167 }

168 }

169 else if (cond == GLP_KKT_DE)

170 { /* A' * (lambdaR - cR) + (lambdaS - cS) = 0 */

171 for (j = 1; j <= n; j++)

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

173 sp = sn = 0.0;

174 /* t := lambdaS[j] - cS[j] */

175 if (sol == GLP_SOL)

176 t = col->dual - col->coef;

177 else if (sol == GLP_IPT)

178 t = col->dval - col->coef;

179 else

180 xassert(sol != sol);

181 if (t >= 0.0) sp += t; else sn -= t;

182 for (aij = col->ptr; aij != NULL; aij = aij->c_next)

183 { row = aij->row;

184 /* t := a[i,j] * (lambdaR[i] - cR[i]) */

185 if (sol == GLP_SOL)

186 t = aij->val * row->dual;

187 else if (sol == GLP_IPT)

188 t = aij->val * row->dval;

189 else

190 xassert(sol != sol);

191 if (t >= 0.0) sp += t; else sn -= t;

192 }

193 /* absolute error */

194 e = fabs(sp - sn);

195 if (ae_max < e)

196 ae_max = e, ae_ind = m+j;

197 /* relative error */

198 e /= (1.0 + sp + sn);

199 if (re_max < e)

200 re_max = e, re_ind = m+j;

201 }

202 }

203 else if (cond == GLP_KKT_DB)

204 { /* check lambdaR */

205 for (i = 1; i <= m; i++)

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

207 /* t := lambdaR[i] */

208 if (sol == GLP_SOL)

209 t = row->dual;

210 else if (sol == GLP_IPT)

211 t = row->dval;

212 else

213 xassert(sol != sol);

214 /* correct sign */

215 if (P->dir == GLP_MIN)

216 t = + t;

217 else if (P->dir == GLP_MAX)

218 t = - t;

219 else

220 xassert(P != P);

221 /* check for positivity */

222 if (row->type == GLP_FR || row->type == GLP_LO)

223 { if (t < 0.0)

224 { e = - t;

225 if (ae_max < e)

226 ae_max = re_max = e, ae_ind = re_ind = i;

227 }

228 }

229 /* check for negativity */

230 if (row->type == GLP_FR || row->type == GLP_UP)

231 { if (t > 0.0)

232 { e = + t;

233 if (ae_max < e)

234 ae_max = re_max = e, ae_ind = re_ind = i;

235 }

236 }

237 }

238 /* check lambdaS */

239 for (j = 1; j <= n; j++)

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

241 /* t := lambdaS[j] */

242 if (sol == GLP_SOL)

243 t = col->dual;

244 else if (sol == GLP_IPT)

245 t = col->dval;

246 else

247 xassert(sol != sol);

248 /* correct sign */

249 if (P->dir == GLP_MIN)

250 t = + t;

251 else if (P->dir == GLP_MAX)

252 t = - t;

253 else

254 xassert(P != P);

255 /* check for positivity */

256 if (col->type == GLP_FR || col->type == GLP_LO)

257 { if (t < 0.0)

258 { e = - t;

259 if (ae_max < e)

260 ae_max = re_max = e, ae_ind = re_ind = m+j;

261 }

262 }

263 /* check for negativity */

264 if (col->type == GLP_FR || col->type == GLP_UP)

265 { if (t > 0.0)

266 { e = + t;

267 if (ae_max < e)

268 ae_max = re_max = e, ae_ind = re_ind = m+j;

269 }

270 }

271 }

272 }

273 else

274 xassert(cond != cond);

275 if (_ae_max != NULL) *_ae_max = ae_max;

276 if (_ae_ind != NULL) *_ae_ind = ae_ind;

277 if (_re_max != NULL) *_re_max = re_max;

278 if (_re_ind != NULL) *_re_ind = re_ind;

279 return;

280 }

281

282 /* eof */