COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpapi10.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 14 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 9.5 KB
Line 
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
27void _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 */
Note: See TracBrowser for help on using the repository browser.