lemon-project-template-glpk

annotate 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
rev   line source
alpar@9 1 /* glpapi10.c (solution checking routines) */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #include "glpapi.h"
alpar@9 26
alpar@9 27 void _glp_check_kkt(glp_prob *P, int sol, int cond, double *_ae_max,
alpar@9 28 int *_ae_ind, double *_re_max, int *_re_ind)
alpar@9 29 { /* check feasibility and optimality conditions */
alpar@9 30 int m = P->m;
alpar@9 31 int n = P->n;
alpar@9 32 GLPROW *row;
alpar@9 33 GLPCOL *col;
alpar@9 34 GLPAIJ *aij;
alpar@9 35 int i, j, ae_ind, re_ind;
alpar@9 36 double e, sp, sn, t, ae_max, re_max;
alpar@9 37 if (!(sol == GLP_SOL || sol == GLP_IPT || sol == GLP_MIP))
alpar@9 38 xerror("glp_check_kkt: sol = %d; invalid solution indicator\n",
alpar@9 39 sol);
alpar@9 40 if (!(cond == GLP_KKT_PE || cond == GLP_KKT_PB ||
alpar@9 41 cond == GLP_KKT_DE || cond == GLP_KKT_DB ||
alpar@9 42 cond == GLP_KKT_CS))
alpar@9 43 xerror("glp_check_kkt: cond = %d; invalid condition indicator "
alpar@9 44 "\n", cond);
alpar@9 45 ae_max = re_max = 0.0;
alpar@9 46 ae_ind = re_ind = 0;
alpar@9 47 if (cond == GLP_KKT_PE)
alpar@9 48 { /* xR - A * xS = 0 */
alpar@9 49 for (i = 1; i <= m; i++)
alpar@9 50 { row = P->row[i];
alpar@9 51 sp = sn = 0.0;
alpar@9 52 /* t := xR[i] */
alpar@9 53 if (sol == GLP_SOL)
alpar@9 54 t = row->prim;
alpar@9 55 else if (sol == GLP_IPT)
alpar@9 56 t = row->pval;
alpar@9 57 else if (sol == GLP_MIP)
alpar@9 58 t = row->mipx;
alpar@9 59 else
alpar@9 60 xassert(sol != sol);
alpar@9 61 if (t >= 0.0) sp += t; else sn -= t;
alpar@9 62 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
alpar@9 63 { col = aij->col;
alpar@9 64 /* t := - a[i,j] * xS[j] */
alpar@9 65 if (sol == GLP_SOL)
alpar@9 66 t = - aij->val * col->prim;
alpar@9 67 else if (sol == GLP_IPT)
alpar@9 68 t = - aij->val * col->pval;
alpar@9 69 else if (sol == GLP_MIP)
alpar@9 70 t = - aij->val * col->mipx;
alpar@9 71 else
alpar@9 72 xassert(sol != sol);
alpar@9 73 if (t >= 0.0) sp += t; else sn -= t;
alpar@9 74 }
alpar@9 75 /* absolute error */
alpar@9 76 e = fabs(sp - sn);
alpar@9 77 if (ae_max < e)
alpar@9 78 ae_max = e, ae_ind = i;
alpar@9 79 /* relative error */
alpar@9 80 e /= (1.0 + sp + sn);
alpar@9 81 if (re_max < e)
alpar@9 82 re_max = e, re_ind = i;
alpar@9 83 }
alpar@9 84 }
alpar@9 85 else if (cond == GLP_KKT_PB)
alpar@9 86 { /* lR <= xR <= uR */
alpar@9 87 for (i = 1; i <= m; i++)
alpar@9 88 { row = P->row[i];
alpar@9 89 /* t := xR[i] */
alpar@9 90 if (sol == GLP_SOL)
alpar@9 91 t = row->prim;
alpar@9 92 else if (sol == GLP_IPT)
alpar@9 93 t = row->pval;
alpar@9 94 else if (sol == GLP_MIP)
alpar@9 95 t = row->mipx;
alpar@9 96 else
alpar@9 97 xassert(sol != sol);
alpar@9 98 /* check lower bound */
alpar@9 99 if (row->type == GLP_LO || row->type == GLP_DB ||
alpar@9 100 row->type == GLP_FX)
alpar@9 101 { if (t < row->lb)
alpar@9 102 { /* absolute error */
alpar@9 103 e = row->lb - t;
alpar@9 104 if (ae_max < e)
alpar@9 105 ae_max = e, ae_ind = i;
alpar@9 106 /* relative error */
alpar@9 107 e /= (1.0 + fabs(row->lb));
alpar@9 108 if (re_max < e)
alpar@9 109 re_max = e, re_ind = i;
alpar@9 110 }
alpar@9 111 }
alpar@9 112 /* check upper bound */
alpar@9 113 if (row->type == GLP_UP || row->type == GLP_DB ||
alpar@9 114 row->type == GLP_FX)
alpar@9 115 { if (t > row->ub)
alpar@9 116 { /* absolute error */
alpar@9 117 e = t - row->ub;
alpar@9 118 if (ae_max < e)
alpar@9 119 ae_max = e, ae_ind = i;
alpar@9 120 /* relative error */
alpar@9 121 e /= (1.0 + fabs(row->ub));
alpar@9 122 if (re_max < e)
alpar@9 123 re_max = e, re_ind = i;
alpar@9 124 }
alpar@9 125 }
alpar@9 126 }
alpar@9 127 /* lS <= xS <= uS */
alpar@9 128 for (j = 1; j <= n; j++)
alpar@9 129 { col = P->col[j];
alpar@9 130 /* t := xS[j] */
alpar@9 131 if (sol == GLP_SOL)
alpar@9 132 t = col->prim;
alpar@9 133 else if (sol == GLP_IPT)
alpar@9 134 t = col->pval;
alpar@9 135 else if (sol == GLP_MIP)
alpar@9 136 t = col->mipx;
alpar@9 137 else
alpar@9 138 xassert(sol != sol);
alpar@9 139 /* check lower bound */
alpar@9 140 if (col->type == GLP_LO || col->type == GLP_DB ||
alpar@9 141 col->type == GLP_FX)
alpar@9 142 { if (t < col->lb)
alpar@9 143 { /* absolute error */
alpar@9 144 e = col->lb - t;
alpar@9 145 if (ae_max < e)
alpar@9 146 ae_max = e, ae_ind = m+j;
alpar@9 147 /* relative error */
alpar@9 148 e /= (1.0 + fabs(col->lb));
alpar@9 149 if (re_max < e)
alpar@9 150 re_max = e, re_ind = m+j;
alpar@9 151 }
alpar@9 152 }
alpar@9 153 /* check upper bound */
alpar@9 154 if (col->type == GLP_UP || col->type == GLP_DB ||
alpar@9 155 col->type == GLP_FX)
alpar@9 156 { if (t > col->ub)
alpar@9 157 { /* absolute error */
alpar@9 158 e = t - col->ub;
alpar@9 159 if (ae_max < e)
alpar@9 160 ae_max = e, ae_ind = m+j;
alpar@9 161 /* relative error */
alpar@9 162 e /= (1.0 + fabs(col->ub));
alpar@9 163 if (re_max < e)
alpar@9 164 re_max = e, re_ind = m+j;
alpar@9 165 }
alpar@9 166 }
alpar@9 167 }
alpar@9 168 }
alpar@9 169 else if (cond == GLP_KKT_DE)
alpar@9 170 { /* A' * (lambdaR - cR) + (lambdaS - cS) = 0 */
alpar@9 171 for (j = 1; j <= n; j++)
alpar@9 172 { col = P->col[j];
alpar@9 173 sp = sn = 0.0;
alpar@9 174 /* t := lambdaS[j] - cS[j] */
alpar@9 175 if (sol == GLP_SOL)
alpar@9 176 t = col->dual - col->coef;
alpar@9 177 else if (sol == GLP_IPT)
alpar@9 178 t = col->dval - col->coef;
alpar@9 179 else
alpar@9 180 xassert(sol != sol);
alpar@9 181 if (t >= 0.0) sp += t; else sn -= t;
alpar@9 182 for (aij = col->ptr; aij != NULL; aij = aij->c_next)
alpar@9 183 { row = aij->row;
alpar@9 184 /* t := a[i,j] * (lambdaR[i] - cR[i]) */
alpar@9 185 if (sol == GLP_SOL)
alpar@9 186 t = aij->val * row->dual;
alpar@9 187 else if (sol == GLP_IPT)
alpar@9 188 t = aij->val * row->dval;
alpar@9 189 else
alpar@9 190 xassert(sol != sol);
alpar@9 191 if (t >= 0.0) sp += t; else sn -= t;
alpar@9 192 }
alpar@9 193 /* absolute error */
alpar@9 194 e = fabs(sp - sn);
alpar@9 195 if (ae_max < e)
alpar@9 196 ae_max = e, ae_ind = m+j;
alpar@9 197 /* relative error */
alpar@9 198 e /= (1.0 + sp + sn);
alpar@9 199 if (re_max < e)
alpar@9 200 re_max = e, re_ind = m+j;
alpar@9 201 }
alpar@9 202 }
alpar@9 203 else if (cond == GLP_KKT_DB)
alpar@9 204 { /* check lambdaR */
alpar@9 205 for (i = 1; i <= m; i++)
alpar@9 206 { row = P->row[i];
alpar@9 207 /* t := lambdaR[i] */
alpar@9 208 if (sol == GLP_SOL)
alpar@9 209 t = row->dual;
alpar@9 210 else if (sol == GLP_IPT)
alpar@9 211 t = row->dval;
alpar@9 212 else
alpar@9 213 xassert(sol != sol);
alpar@9 214 /* correct sign */
alpar@9 215 if (P->dir == GLP_MIN)
alpar@9 216 t = + t;
alpar@9 217 else if (P->dir == GLP_MAX)
alpar@9 218 t = - t;
alpar@9 219 else
alpar@9 220 xassert(P != P);
alpar@9 221 /* check for positivity */
alpar@9 222 if (row->type == GLP_FR || row->type == GLP_LO)
alpar@9 223 { if (t < 0.0)
alpar@9 224 { e = - t;
alpar@9 225 if (ae_max < e)
alpar@9 226 ae_max = re_max = e, ae_ind = re_ind = i;
alpar@9 227 }
alpar@9 228 }
alpar@9 229 /* check for negativity */
alpar@9 230 if (row->type == GLP_FR || row->type == GLP_UP)
alpar@9 231 { if (t > 0.0)
alpar@9 232 { e = + t;
alpar@9 233 if (ae_max < e)
alpar@9 234 ae_max = re_max = e, ae_ind = re_ind = i;
alpar@9 235 }
alpar@9 236 }
alpar@9 237 }
alpar@9 238 /* check lambdaS */
alpar@9 239 for (j = 1; j <= n; j++)
alpar@9 240 { col = P->col[j];
alpar@9 241 /* t := lambdaS[j] */
alpar@9 242 if (sol == GLP_SOL)
alpar@9 243 t = col->dual;
alpar@9 244 else if (sol == GLP_IPT)
alpar@9 245 t = col->dval;
alpar@9 246 else
alpar@9 247 xassert(sol != sol);
alpar@9 248 /* correct sign */
alpar@9 249 if (P->dir == GLP_MIN)
alpar@9 250 t = + t;
alpar@9 251 else if (P->dir == GLP_MAX)
alpar@9 252 t = - t;
alpar@9 253 else
alpar@9 254 xassert(P != P);
alpar@9 255 /* check for positivity */
alpar@9 256 if (col->type == GLP_FR || col->type == GLP_LO)
alpar@9 257 { if (t < 0.0)
alpar@9 258 { e = - t;
alpar@9 259 if (ae_max < e)
alpar@9 260 ae_max = re_max = e, ae_ind = re_ind = m+j;
alpar@9 261 }
alpar@9 262 }
alpar@9 263 /* check for negativity */
alpar@9 264 if (col->type == GLP_FR || col->type == GLP_UP)
alpar@9 265 { if (t > 0.0)
alpar@9 266 { e = + t;
alpar@9 267 if (ae_max < e)
alpar@9 268 ae_max = re_max = e, ae_ind = re_ind = m+j;
alpar@9 269 }
alpar@9 270 }
alpar@9 271 }
alpar@9 272 }
alpar@9 273 else
alpar@9 274 xassert(cond != cond);
alpar@9 275 if (_ae_max != NULL) *_ae_max = ae_max;
alpar@9 276 if (_ae_ind != NULL) *_ae_ind = ae_ind;
alpar@9 277 if (_re_max != NULL) *_re_max = re_max;
alpar@9 278 if (_re_ind != NULL) *_re_ind = re_ind;
alpar@9 279 return;
alpar@9 280 }
alpar@9 281
alpar@9 282 /* eof */