|
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 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 */ |