lemon-project-template-glpk

comparison deps/glpk/src/glpapi10.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:bda8949183e4
1 /* glpapi10.c (solution checking routines) */
2
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 *
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
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 ***********************************************************************/
24
25 #include "glpapi.h"
26
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 */