glpk-cmake: src/glpapi10.c@c445c931472f

     1 /* glpapi10.c (solution checking routines) */

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

    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 }

   282 /* eof */

author	Alpar Juttner <alpar@cs.elte.hu>
	Mon, 06 Dec 2010 13:09:21 +0100
changeset 1	c445c931472f
permissions	-rw-r--r--