src/glpapi11.c
changeset 1 c445c931472f
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/glpapi11.c	Mon Dec 06 13:09:21 2010 +0100
     1.3 @@ -0,0 +1,1217 @@
     1.4 +/* glpapi11.c (utility routines) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glpapi.h"
    1.29 +
    1.30 +int glp_print_sol(glp_prob *P, const char *fname)
    1.31 +{     /* write basic solution in printable format */
    1.32 +      XFILE *fp;
    1.33 +      GLPROW *row;
    1.34 +      GLPCOL *col;
    1.35 +      int i, j, t, ae_ind, re_ind, ret;
    1.36 +      double ae_max, re_max;
    1.37 +      xprintf("Writing basic solution to `%s'...\n", fname);
    1.38 +      fp = xfopen(fname, "w");
    1.39 +      if (fp == NULL)
    1.40 +      {  xprintf("Unable to create `%s' - %s\n", fname, xerrmsg());
    1.41 +         ret = 1;
    1.42 +         goto done;
    1.43 +      }
    1.44 +      xfprintf(fp, "%-12s%s\n", "Problem:",
    1.45 +         P->name == NULL ? "" : P->name);
    1.46 +      xfprintf(fp, "%-12s%d\n", "Rows:", P->m);
    1.47 +      xfprintf(fp, "%-12s%d\n", "Columns:", P->n);
    1.48 +      xfprintf(fp, "%-12s%d\n", "Non-zeros:", P->nnz);
    1.49 +      t = glp_get_status(P);
    1.50 +      xfprintf(fp, "%-12s%s\n", "Status:",
    1.51 +         t == GLP_OPT    ? "OPTIMAL" :
    1.52 +         t == GLP_FEAS   ? "FEASIBLE" :
    1.53 +         t == GLP_INFEAS ? "INFEASIBLE (INTERMEDIATE)" :
    1.54 +         t == GLP_NOFEAS ? "INFEASIBLE (FINAL)" :
    1.55 +         t == GLP_UNBND  ? "UNBOUNDED" :
    1.56 +         t == GLP_UNDEF  ? "UNDEFINED" : "???");
    1.57 +      xfprintf(fp, "%-12s%s%s%.10g (%s)\n", "Objective:",
    1.58 +         P->obj == NULL ? "" : P->obj,
    1.59 +         P->obj == NULL ? "" : " = ", P->obj_val,
    1.60 +         P->dir == GLP_MIN ? "MINimum" :
    1.61 +         P->dir == GLP_MAX ? "MAXimum" : "???");
    1.62 +      xfprintf(fp, "\n");
    1.63 +      xfprintf(fp, "   No.   Row name   St   Activity     Lower bound  "
    1.64 +         " Upper bound    Marginal\n");
    1.65 +      xfprintf(fp, "------ ------------ -- ------------- ------------- "
    1.66 +         "------------- -------------\n");
    1.67 +      for (i = 1; i <= P->m; i++)
    1.68 +      {  row = P->row[i];
    1.69 +         xfprintf(fp, "%6d ", i);
    1.70 +         if (row->name == NULL || strlen(row->name) <= 12)
    1.71 +            xfprintf(fp, "%-12s ", row->name == NULL ? "" : row->name);
    1.72 +         else
    1.73 +            xfprintf(fp, "%s\n%20s", row->name, "");
    1.74 +         xfprintf(fp, "%s ",
    1.75 +            row->stat == GLP_BS ? "B " :
    1.76 +            row->stat == GLP_NL ? "NL" :
    1.77 +            row->stat == GLP_NU ? "NU" :
    1.78 +            row->stat == GLP_NF ? "NF" :
    1.79 +            row->stat == GLP_NS ? "NS" : "??");
    1.80 +         xfprintf(fp, "%13.6g ",
    1.81 +            fabs(row->prim) <= 1e-9 ? 0.0 : row->prim);
    1.82 +         if (row->type == GLP_LO || row->type == GLP_DB ||
    1.83 +             row->type == GLP_FX)
    1.84 +            xfprintf(fp, "%13.6g ", row->lb);
    1.85 +         else
    1.86 +            xfprintf(fp, "%13s ", "");
    1.87 +         if (row->type == GLP_UP || row->type == GLP_DB)
    1.88 +            xfprintf(fp, "%13.6g ", row->ub);
    1.89 +         else
    1.90 +            xfprintf(fp, "%13s ", row->type == GLP_FX ? "=" : "");
    1.91 +         if (row->stat != GLP_BS)
    1.92 +         {  if (fabs(row->dual) <= 1e-9)
    1.93 +               xfprintf(fp, "%13s", "< eps");
    1.94 +            else
    1.95 +               xfprintf(fp, "%13.6g ", row->dual);
    1.96 +         }
    1.97 +         xfprintf(fp, "\n");
    1.98 +      }
    1.99 +      xfprintf(fp, "\n");
   1.100 +      xfprintf(fp, "   No. Column name  St   Activity     Lower bound  "
   1.101 +         " Upper bound    Marginal\n");
   1.102 +      xfprintf(fp, "------ ------------ -- ------------- ------------- "
   1.103 +         "------------- -------------\n");
   1.104 +      for (j = 1; j <= P->n; j++)
   1.105 +      {  col = P->col[j];
   1.106 +         xfprintf(fp, "%6d ", j);
   1.107 +         if (col->name == NULL || strlen(col->name) <= 12)
   1.108 +            xfprintf(fp, "%-12s ", col->name == NULL ? "" : col->name);
   1.109 +         else
   1.110 +            xfprintf(fp, "%s\n%20s", col->name, "");
   1.111 +         xfprintf(fp, "%s ",
   1.112 +            col->stat == GLP_BS ? "B " :
   1.113 +            col->stat == GLP_NL ? "NL" :
   1.114 +            col->stat == GLP_NU ? "NU" :
   1.115 +            col->stat == GLP_NF ? "NF" :
   1.116 +            col->stat == GLP_NS ? "NS" : "??");
   1.117 +         xfprintf(fp, "%13.6g ",
   1.118 +            fabs(col->prim) <= 1e-9 ? 0.0 : col->prim);
   1.119 +         if (col->type == GLP_LO || col->type == GLP_DB ||
   1.120 +             col->type == GLP_FX)
   1.121 +            xfprintf(fp, "%13.6g ", col->lb);
   1.122 +         else
   1.123 +            xfprintf(fp, "%13s ", "");
   1.124 +         if (col->type == GLP_UP || col->type == GLP_DB)
   1.125 +            xfprintf(fp, "%13.6g ", col->ub);
   1.126 +         else
   1.127 +            xfprintf(fp, "%13s ", col->type == GLP_FX ? "=" : "");
   1.128 +         if (col->stat != GLP_BS)
   1.129 +         {  if (fabs(col->dual) <= 1e-9)
   1.130 +               xfprintf(fp, "%13s", "< eps");
   1.131 +            else
   1.132 +               xfprintf(fp, "%13.6g ", col->dual);
   1.133 +         }
   1.134 +         xfprintf(fp, "\n");
   1.135 +      }
   1.136 +      xfprintf(fp, "\n");
   1.137 +      xfprintf(fp, "Karush-Kuhn-Tucker optimality conditions:\n");
   1.138 +      xfprintf(fp, "\n");
   1.139 +      _glp_check_kkt(P, GLP_SOL, GLP_KKT_PE, &ae_max, &ae_ind, &re_max,
   1.140 +         &re_ind);
   1.141 +      xfprintf(fp, "KKT.PE: max.abs.err = %.2e on row %d\n",
   1.142 +         ae_max, ae_ind);
   1.143 +      xfprintf(fp, "        max.rel.err = %.2e on row %d\n",
   1.144 +         re_max, re_ind);
   1.145 +      xfprintf(fp, "%8s%s\n", "",
   1.146 +         re_max <= 1e-9 ? "High quality" :
   1.147 +         re_max <= 1e-6 ? "Medium quality" :
   1.148 +         re_max <= 1e-3 ? "Low quality" : "PRIMAL SOLUTION IS WRONG");
   1.149 +      xfprintf(fp, "\n");
   1.150 +      _glp_check_kkt(P, GLP_SOL, GLP_KKT_PB, &ae_max, &ae_ind, &re_max,
   1.151 +         &re_ind);
   1.152 +      xfprintf(fp, "KKT.PB: max.abs.err = %.2e on %s %d\n",
   1.153 +            ae_max, ae_ind <= P->m ? "row" : "column",
   1.154 +            ae_ind <= P->m ? ae_ind : ae_ind - P->m);
   1.155 +      xfprintf(fp, "        max.rel.err = %.2e on %s %d\n",
   1.156 +            re_max, re_ind <= P->m ? "row" : "column",
   1.157 +            re_ind <= P->m ? re_ind : re_ind - P->m);
   1.158 +      xfprintf(fp, "%8s%s\n", "",
   1.159 +         re_max <= 1e-9 ? "High quality" :
   1.160 +         re_max <= 1e-6 ? "Medium quality" :
   1.161 +         re_max <= 1e-3 ? "Low quality" : "PRIMAL SOLUTION IS INFEASIBL"
   1.162 +            "E");
   1.163 +      xfprintf(fp, "\n");
   1.164 +      _glp_check_kkt(P, GLP_SOL, GLP_KKT_DE, &ae_max, &ae_ind, &re_max,
   1.165 +         &re_ind);
   1.166 +      xfprintf(fp, "KKT.DE: max.abs.err = %.2e on column %d\n",
   1.167 +         ae_max, ae_ind == 0 ? 0 : ae_ind - P->m);
   1.168 +      xfprintf(fp, "        max.rel.err = %.2e on column %d\n",
   1.169 +         re_max, re_ind == 0 ? 0 : re_ind - P->m);
   1.170 +      xfprintf(fp, "%8s%s\n", "",
   1.171 +         re_max <= 1e-9 ? "High quality" :
   1.172 +         re_max <= 1e-6 ? "Medium quality" :
   1.173 +         re_max <= 1e-3 ? "Low quality" : "DUAL SOLUTION IS WRONG");
   1.174 +      xfprintf(fp, "\n");
   1.175 +      _glp_check_kkt(P, GLP_SOL, GLP_KKT_DB, &ae_max, &ae_ind, &re_max,
   1.176 +         &re_ind);
   1.177 +      xfprintf(fp, "KKT.DB: max.abs.err = %.2e on %s %d\n",
   1.178 +            ae_max, ae_ind <= P->m ? "row" : "column",
   1.179 +            ae_ind <= P->m ? ae_ind : ae_ind - P->m);
   1.180 +      xfprintf(fp, "        max.rel.err = %.2e on %s %d\n",
   1.181 +            re_max, re_ind <= P->m ? "row" : "column",
   1.182 +            re_ind <= P->m ? re_ind : re_ind - P->m);
   1.183 +      xfprintf(fp, "%8s%s\n", "",
   1.184 +         re_max <= 1e-9 ? "High quality" :
   1.185 +         re_max <= 1e-6 ? "Medium quality" :
   1.186 +         re_max <= 1e-3 ? "Low quality" : "DUAL SOLUTION IS INFEASIBLE")
   1.187 +            ;
   1.188 +      xfprintf(fp, "\n");
   1.189 +      xfprintf(fp, "End of output\n");
   1.190 +      xfflush(fp);
   1.191 +      if (xferror(fp))
   1.192 +      {  xprintf("Write error on `%s' - %s\n", fname, xerrmsg());
   1.193 +         ret = 1;
   1.194 +         goto done;
   1.195 +      }
   1.196 +      ret = 0;
   1.197 +done: if (fp != NULL) xfclose(fp);
   1.198 +      return ret;
   1.199 +}
   1.200 +
   1.201 +/***********************************************************************
   1.202 +*  NAME
   1.203 +*
   1.204 +*  glp_read_sol - read basic solution from text file
   1.205 +*
   1.206 +*  SYNOPSIS
   1.207 +*
   1.208 +*  int glp_read_sol(glp_prob *lp, const char *fname);
   1.209 +*
   1.210 +*  DESCRIPTION
   1.211 +*
   1.212 +*  The routine glp_read_sol reads basic solution from a text file whose
   1.213 +*  name is specified by the parameter fname into the problem object.
   1.214 +*
   1.215 +*  For the file format see description of the routine glp_write_sol.
   1.216 +*
   1.217 +*  RETURNS
   1.218 +*
   1.219 +*  On success the routine returns zero, otherwise non-zero. */
   1.220 +
   1.221 +int glp_read_sol(glp_prob *lp, const char *fname)
   1.222 +{     glp_data *data;
   1.223 +      jmp_buf jump;
   1.224 +      int i, j, k, ret = 0;
   1.225 +      xprintf("Reading basic solution from `%s'...\n", fname);
   1.226 +      data = glp_sdf_open_file(fname);
   1.227 +      if (data == NULL)
   1.228 +      {  ret = 1;
   1.229 +         goto done;
   1.230 +      }
   1.231 +      if (setjmp(jump))
   1.232 +      {  ret = 1;
   1.233 +         goto done;
   1.234 +      }
   1.235 +      glp_sdf_set_jump(data, jump);
   1.236 +      /* number of rows, number of columns */
   1.237 +      k = glp_sdf_read_int(data);
   1.238 +      if (k != lp->m)
   1.239 +         glp_sdf_error(data, "wrong number of rows\n");
   1.240 +      k = glp_sdf_read_int(data);
   1.241 +      if (k != lp->n)
   1.242 +         glp_sdf_error(data, "wrong number of columns\n");
   1.243 +      /* primal status, dual status, objective value */
   1.244 +      k = glp_sdf_read_int(data);
   1.245 +      if (!(k == GLP_UNDEF || k == GLP_FEAS || k == GLP_INFEAS ||
   1.246 +            k == GLP_NOFEAS))
   1.247 +         glp_sdf_error(data, "invalid primal status\n");
   1.248 +      lp->pbs_stat = k;
   1.249 +      k = glp_sdf_read_int(data);
   1.250 +      if (!(k == GLP_UNDEF || k == GLP_FEAS || k == GLP_INFEAS ||
   1.251 +            k == GLP_NOFEAS))
   1.252 +         glp_sdf_error(data, "invalid dual status\n");
   1.253 +      lp->dbs_stat = k;
   1.254 +      lp->obj_val = glp_sdf_read_num(data);
   1.255 +      /* rows (auxiliary variables) */
   1.256 +      for (i = 1; i <= lp->m; i++)
   1.257 +      {  GLPROW *row = lp->row[i];
   1.258 +         /* status, primal value, dual value */
   1.259 +         k = glp_sdf_read_int(data);
   1.260 +         if (!(k == GLP_BS || k == GLP_NL || k == GLP_NU ||
   1.261 +               k == GLP_NF || k == GLP_NS))
   1.262 +            glp_sdf_error(data, "invalid row status\n");
   1.263 +         glp_set_row_stat(lp, i, k);
   1.264 +         row->prim = glp_sdf_read_num(data);
   1.265 +         row->dual = glp_sdf_read_num(data);
   1.266 +      }
   1.267 +      /* columns (structural variables) */
   1.268 +      for (j = 1; j <= lp->n; j++)
   1.269 +      {  GLPCOL *col = lp->col[j];
   1.270 +         /* status, primal value, dual value */
   1.271 +         k = glp_sdf_read_int(data);
   1.272 +         if (!(k == GLP_BS || k == GLP_NL || k == GLP_NU ||
   1.273 +               k == GLP_NF || k == GLP_NS))
   1.274 +            glp_sdf_error(data, "invalid column status\n");
   1.275 +         glp_set_col_stat(lp, j, k);
   1.276 +         col->prim = glp_sdf_read_num(data);
   1.277 +         col->dual = glp_sdf_read_num(data);
   1.278 +      }
   1.279 +      xprintf("%d lines were read\n", glp_sdf_line(data));
   1.280 +done: if (ret) lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
   1.281 +      if (data != NULL) glp_sdf_close_file(data);
   1.282 +      return ret;
   1.283 +}
   1.284 +
   1.285 +/***********************************************************************
   1.286 +*  NAME
   1.287 +*
   1.288 +*  glp_write_sol - write basic solution to text file
   1.289 +*
   1.290 +*  SYNOPSIS
   1.291 +*
   1.292 +*  int glp_write_sol(glp_prob *lp, const char *fname);
   1.293 +*
   1.294 +*  DESCRIPTION
   1.295 +*
   1.296 +*  The routine glp_write_sol writes the current basic solution to a
   1.297 +*  text file whose name is specified by the parameter fname. This file
   1.298 +*  can be read back with the routine glp_read_sol.
   1.299 +*
   1.300 +*  RETURNS
   1.301 +*
   1.302 +*  On success the routine returns zero, otherwise non-zero.
   1.303 +*
   1.304 +*  FILE FORMAT
   1.305 +*
   1.306 +*  The file created by the routine glp_write_sol is a plain text file,
   1.307 +*  which contains the following information:
   1.308 +*
   1.309 +*     m n
   1.310 +*     p_stat d_stat obj_val
   1.311 +*     r_stat[1] r_prim[1] r_dual[1]
   1.312 +*     . . .
   1.313 +*     r_stat[m] r_prim[m] r_dual[m]
   1.314 +*     c_stat[1] c_prim[1] c_dual[1]
   1.315 +*     . . .
   1.316 +*     c_stat[n] c_prim[n] c_dual[n]
   1.317 +*
   1.318 +*  where:
   1.319 +*  m is the number of rows (auxiliary variables);
   1.320 +*  n is the number of columns (structural variables);
   1.321 +*  p_stat is the primal status of the basic solution (GLP_UNDEF = 1,
   1.322 +*     GLP_FEAS = 2, GLP_INFEAS = 3, or GLP_NOFEAS = 4);
   1.323 +*  d_stat is the dual status of the basic solution (GLP_UNDEF = 1,
   1.324 +*     GLP_FEAS = 2, GLP_INFEAS = 3, or GLP_NOFEAS = 4);
   1.325 +*  obj_val is the objective value;
   1.326 +*  r_stat[i], i = 1,...,m, is the status of i-th row (GLP_BS = 1,
   1.327 +*     GLP_NL = 2, GLP_NU = 3, GLP_NF = 4, or GLP_NS = 5);
   1.328 +*  r_prim[i], i = 1,...,m, is the primal value of i-th row;
   1.329 +*  r_dual[i], i = 1,...,m, is the dual value of i-th row;
   1.330 +*  c_stat[j], j = 1,...,n, is the status of j-th column (GLP_BS = 1,
   1.331 +*     GLP_NL = 2, GLP_NU = 3, GLP_NF = 4, or GLP_NS = 5);
   1.332 +*  c_prim[j], j = 1,...,n, is the primal value of j-th column;
   1.333 +*  c_dual[j], j = 1,...,n, is the dual value of j-th column. */
   1.334 +
   1.335 +int glp_write_sol(glp_prob *lp, const char *fname)
   1.336 +{     XFILE *fp;
   1.337 +      int i, j, ret = 0;
   1.338 +      xprintf("Writing basic solution to `%s'...\n", fname);
   1.339 +      fp = xfopen(fname, "w");
   1.340 +      if (fp == NULL)
   1.341 +      {  xprintf("Unable to create `%s' - %s\n", fname, xerrmsg());
   1.342 +         ret = 1;
   1.343 +         goto done;
   1.344 +      }
   1.345 +      /* number of rows, number of columns */
   1.346 +      xfprintf(fp, "%d %d\n", lp->m, lp->n);
   1.347 +      /* primal status, dual status, objective value */
   1.348 +      xfprintf(fp, "%d %d %.*g\n", lp->pbs_stat, lp->dbs_stat, DBL_DIG,
   1.349 +         lp->obj_val);
   1.350 +      /* rows (auxiliary variables) */
   1.351 +      for (i = 1; i <= lp->m; i++)
   1.352 +      {  GLPROW *row = lp->row[i];
   1.353 +         /* status, primal value, dual value */
   1.354 +         xfprintf(fp, "%d %.*g %.*g\n", row->stat, DBL_DIG, row->prim,
   1.355 +            DBL_DIG, row->dual);
   1.356 +      }
   1.357 +      /* columns (structural variables) */
   1.358 +      for (j = 1; j <= lp->n; j++)
   1.359 +      {  GLPCOL *col = lp->col[j];
   1.360 +         /* status, primal value, dual value */
   1.361 +         xfprintf(fp, "%d %.*g %.*g\n", col->stat, DBL_DIG, col->prim,
   1.362 +            DBL_DIG, col->dual);
   1.363 +      }
   1.364 +      xfflush(fp);
   1.365 +      if (xferror(fp))
   1.366 +      {  xprintf("Write error on `%s' - %s\n", fname, xerrmsg());
   1.367 +         ret = 1;
   1.368 +         goto done;
   1.369 +      }
   1.370 +      xprintf("%d lines were written\n", 2 + lp->m + lp->n);
   1.371 +done: if (fp != NULL) xfclose(fp);
   1.372 +      return ret;
   1.373 +}
   1.374 +
   1.375 +/**********************************************************************/
   1.376 +
   1.377 +static char *format(char buf[13+1], double x)
   1.378 +{     /* format floating-point number in MPS/360-like style */
   1.379 +      if (x == -DBL_MAX)
   1.380 +         strcpy(buf, "         -Inf");
   1.381 +      else if (x == +DBL_MAX)
   1.382 +         strcpy(buf, "         +Inf");
   1.383 +      else if (fabs(x) <= 999999.99998)
   1.384 +      {  sprintf(buf, "%13.5f", x);
   1.385 +#if 1
   1.386 +         if (strcmp(buf, "      0.00000") == 0 ||
   1.387 +             strcmp(buf, "     -0.00000") == 0)
   1.388 +            strcpy(buf, "       .     ");
   1.389 +         else if (memcmp(buf, "      0.", 8) == 0)
   1.390 +            memcpy(buf, "       .", 8);
   1.391 +         else if (memcmp(buf, "     -0.", 8) == 0)
   1.392 +            memcpy(buf, "      -.", 8);
   1.393 +#endif
   1.394 +      }
   1.395 +      else
   1.396 +         sprintf(buf, "%13.6g", x);
   1.397 +      return buf;
   1.398 +}
   1.399 +
   1.400 +int glp_print_ranges(glp_prob *P, int len, const int list[],
   1.401 +      int flags, const char *fname)
   1.402 +{     /* print sensitivity analysis report */
   1.403 +      XFILE *fp = NULL;
   1.404 +      GLPROW *row;
   1.405 +      GLPCOL *col;
   1.406 +      int m, n, pass, k, t, numb, type, stat, var1, var2, count, page,
   1.407 +         ret;
   1.408 +      double lb, ub, slack, coef, prim, dual, value1, value2, coef1,
   1.409 +         coef2, obj1, obj2;
   1.410 +      const char *name, *limit;
   1.411 +      char buf[13+1];
   1.412 +      /* sanity checks */
   1.413 +      if (P == NULL || P->magic != GLP_PROB_MAGIC)
   1.414 +         xerror("glp_print_ranges: P = %p; invalid problem object\n",
   1.415 +            P);
   1.416 +      m = P->m, n = P->n;
   1.417 +      if (len < 0)
   1.418 +         xerror("glp_print_ranges: len = %d; invalid list length\n",
   1.419 +            len);
   1.420 +      if (len > 0)
   1.421 +      {  if (list == NULL)
   1.422 +            xerror("glp_print_ranges: list = %p: invalid parameter\n",
   1.423 +               list);
   1.424 +         for (t = 1; t <= len; t++)
   1.425 +         {  k = list[t];
   1.426 +            if (!(1 <= k && k <= m+n))
   1.427 +               xerror("glp_print_ranges: list[%d] = %d; row/column numb"
   1.428 +                  "er out of range\n", t, k);
   1.429 +         }
   1.430 +      }
   1.431 +      if (flags != 0)
   1.432 +         xerror("glp_print_ranges: flags = %d; invalid parameter\n",
   1.433 +            flags);
   1.434 +      if (fname == NULL)
   1.435 +         xerror("glp_print_ranges: fname = %p; invalid parameter\n",
   1.436 +            fname);
   1.437 +      if (glp_get_status(P) != GLP_OPT)
   1.438 +      {  xprintf("glp_print_ranges: optimal basic solution required\n");
   1.439 +         ret = 1;
   1.440 +         goto done;
   1.441 +      }
   1.442 +      if (!glp_bf_exists(P))
   1.443 +      {  xprintf("glp_print_ranges: basis factorization required\n");
   1.444 +         ret = 2;
   1.445 +         goto done;
   1.446 +      }
   1.447 +      /* start reporting */
   1.448 +      xprintf("Write sensitivity analysis report to `%s'...\n", fname);
   1.449 +      fp = xfopen(fname, "w");
   1.450 +      if (fp == NULL)
   1.451 +      {  xprintf("Unable to create `%s' - %s\n", fname, xerrmsg());
   1.452 +         ret = 3;
   1.453 +         goto done;
   1.454 +      }
   1.455 +      page = count = 0;
   1.456 +      for (pass = 1; pass <= 2; pass++)
   1.457 +      for (t = 1; t <= (len == 0 ? m+n : len); t++)
   1.458 +      {  if (t == 1) count = 0;
   1.459 +         k = (len == 0 ? t : list[t]);
   1.460 +         if (pass == 1 && k > m || pass == 2 && k <= m)
   1.461 +            continue;
   1.462 +         if (count == 0)
   1.463 +         {  xfprintf(fp, "GLPK %-4s - SENSITIVITY ANALYSIS REPORT%73sPa"
   1.464 +               "ge%4d\n", glp_version(), "", ++page);
   1.465 +            xfprintf(fp, "\n");
   1.466 +            xfprintf(fp, "%-12s%s\n", "Problem:",
   1.467 +               P->name == NULL ? "" : P->name);
   1.468 +            xfprintf(fp, "%-12s%s%s%.10g (%s)\n", "Objective:",
   1.469 +               P->obj == NULL ? "" : P->obj,
   1.470 +               P->obj == NULL ? "" : " = ", P->obj_val,
   1.471 +               P->dir == GLP_MIN ? "MINimum" :
   1.472 +               P->dir == GLP_MAX ? "MAXimum" : "???");
   1.473 +            xfprintf(fp, "\n");
   1.474 +            xfprintf(fp, "%6s %-12s %2s %13s %13s %13s  %13s %13s %13s "
   1.475 +               "%s\n", "No.", pass == 1 ? "Row name" : "Column name",
   1.476 +               "St", "Activity", pass == 1 ? "Slack" : "Obj coef",
   1.477 +               "Lower bound", "Activity", "Obj coef", "Obj value at",
   1.478 +               "Limiting");
   1.479 +            xfprintf(fp, "%6s %-12s %2s %13s %13s %13s  %13s %13s %13s "
   1.480 +               "%s\n", "", "", "", "", "Marginal", "Upper bound",
   1.481 +               "range", "range", "break point", "variable");
   1.482 +            xfprintf(fp, "------ ------------ -- ------------- --------"
   1.483 +               "----- -------------  ------------- ------------- ------"
   1.484 +               "------- ------------\n");
   1.485 +         }
   1.486 +         if (pass == 1)
   1.487 +         {  numb = k;
   1.488 +            xassert(1 <= numb && numb <= m);
   1.489 +            row = P->row[numb];
   1.490 +            name = row->name;
   1.491 +            type = row->type;
   1.492 +            lb = glp_get_row_lb(P, numb);
   1.493 +            ub = glp_get_row_ub(P, numb);
   1.494 +            coef = 0.0;
   1.495 +            stat = row->stat;
   1.496 +            prim = row->prim;
   1.497 +            if (type == GLP_FR)
   1.498 +               slack = - prim;
   1.499 +            else if (type == GLP_LO)
   1.500 +               slack = lb - prim;
   1.501 +            else if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
   1.502 +               slack = ub - prim;
   1.503 +            dual = row->dual;
   1.504 +         }
   1.505 +         else
   1.506 +         {  numb = k - m;
   1.507 +            xassert(1 <= numb && numb <= n);
   1.508 +            col = P->col[numb];
   1.509 +            name = col->name;
   1.510 +            lb = glp_get_col_lb(P, numb);
   1.511 +            ub = glp_get_col_ub(P, numb);
   1.512 +            coef = col->coef;
   1.513 +            stat = col->stat;
   1.514 +            prim = col->prim;
   1.515 +            slack = 0.0;
   1.516 +            dual = col->dual;
   1.517 +         }
   1.518 +         if (stat != GLP_BS)
   1.519 +         {  glp_analyze_bound(P, k, &value1, &var1, &value2, &var2);
   1.520 +            if (stat == GLP_NF)
   1.521 +               coef1 = coef2 = coef;
   1.522 +            else if (stat == GLP_NS)
   1.523 +               coef1 = -DBL_MAX, coef2 = +DBL_MAX;
   1.524 +            else if (stat == GLP_NL && P->dir == GLP_MIN ||
   1.525 +                     stat == GLP_NU && P->dir == GLP_MAX)
   1.526 +               coef1 = coef - dual, coef2 = +DBL_MAX;
   1.527 +            else
   1.528 +               coef1 = -DBL_MAX, coef2 = coef - dual;
   1.529 +            if (value1 == -DBL_MAX)
   1.530 +            {  if (dual < -1e-9)
   1.531 +                  obj1 = +DBL_MAX;
   1.532 +               else if (dual > +1e-9)
   1.533 +                  obj1 = -DBL_MAX;
   1.534 +               else
   1.535 +                  obj1 = P->obj_val;
   1.536 +            }
   1.537 +            else
   1.538 +               obj1 = P->obj_val + dual * (value1 - prim);
   1.539 +            if (value2 == +DBL_MAX)
   1.540 +            {  if (dual < -1e-9)
   1.541 +                  obj2 = -DBL_MAX;
   1.542 +               else if (dual > +1e-9)
   1.543 +                  obj2 = +DBL_MAX;
   1.544 +               else
   1.545 +                  obj2 = P->obj_val;
   1.546 +            }
   1.547 +            else
   1.548 +               obj2 = P->obj_val + dual * (value2 - prim);
   1.549 +         }
   1.550 +         else
   1.551 +         {  glp_analyze_coef(P, k, &coef1, &var1, &value1, &coef2,
   1.552 +               &var2, &value2);
   1.553 +            if (coef1 == -DBL_MAX)
   1.554 +            {  if (prim < -1e-9)
   1.555 +                  obj1 = +DBL_MAX;
   1.556 +               else if (prim > +1e-9)
   1.557 +                  obj1 = -DBL_MAX;
   1.558 +               else
   1.559 +                  obj1 = P->obj_val;
   1.560 +            }
   1.561 +            else
   1.562 +               obj1 = P->obj_val + (coef1 - coef) * prim;
   1.563 +            if (coef2 == +DBL_MAX)
   1.564 +            {  if (prim < -1e-9)
   1.565 +                  obj2 = -DBL_MAX;
   1.566 +               else if (prim > +1e-9)
   1.567 +                  obj2 = +DBL_MAX;
   1.568 +               else
   1.569 +                  obj2 = P->obj_val;
   1.570 +            }
   1.571 +            else
   1.572 +               obj2 = P->obj_val + (coef2 - coef) * prim;
   1.573 +         }
   1.574 +         /*** first line ***/
   1.575 +         /* row/column number */
   1.576 +         xfprintf(fp, "%6d", numb);
   1.577 +         /* row/column name */
   1.578 +         xfprintf(fp, " %-12.12s", name == NULL ? "" : name);
   1.579 +         if (name != NULL && strlen(name) > 12)
   1.580 +            xfprintf(fp, "%s\n%6s %12s", name+12, "", "");
   1.581 +         /* row/column status */
   1.582 +         xfprintf(fp, " %2s",
   1.583 +            stat == GLP_BS ? "BS" : stat == GLP_NL ? "NL" :
   1.584 +            stat == GLP_NU ? "NU" : stat == GLP_NF ? "NF" :
   1.585 +            stat == GLP_NS ? "NS" : "??");
   1.586 +         /* row/column activity */
   1.587 +         xfprintf(fp, " %s", format(buf, prim));
   1.588 +         /* row slack, column objective coefficient */
   1.589 +         xfprintf(fp, " %s", format(buf, k <= m ? slack : coef));
   1.590 +         /* row/column lower bound */
   1.591 +         xfprintf(fp, " %s", format(buf, lb));
   1.592 +         /* row/column activity range */
   1.593 +         xfprintf(fp, "  %s", format(buf, value1));
   1.594 +         /* row/column objective coefficient range */
   1.595 +         xfprintf(fp, " %s", format(buf, coef1));
   1.596 +         /* objective value at break point */
   1.597 +         xfprintf(fp, " %s", format(buf, obj1));
   1.598 +         /* limiting variable name */
   1.599 +         if (var1 != 0)
   1.600 +         {  if (var1 <= m)
   1.601 +               limit = glp_get_row_name(P, var1);
   1.602 +            else
   1.603 +               limit = glp_get_col_name(P, var1 - m);
   1.604 +            if (limit != NULL)
   1.605 +               xfprintf(fp, " %s", limit);
   1.606 +         }
   1.607 +         xfprintf(fp, "\n");
   1.608 +         /*** second line ***/
   1.609 +         xfprintf(fp, "%6s %-12s %2s %13s", "", "", "", "");
   1.610 +         /* row/column reduced cost */
   1.611 +         xfprintf(fp, " %s", format(buf, dual));
   1.612 +         /* row/column upper bound */
   1.613 +         xfprintf(fp, " %s", format(buf, ub));
   1.614 +         /* row/column activity range */
   1.615 +         xfprintf(fp, "  %s", format(buf, value2));
   1.616 +         /* row/column objective coefficient range */
   1.617 +         xfprintf(fp, " %s", format(buf, coef2));
   1.618 +         /* objective value at break point */
   1.619 +         xfprintf(fp, " %s", format(buf, obj2));
   1.620 +         /* limiting variable name */
   1.621 +         if (var2 != 0)
   1.622 +         {  if (var2 <= m)
   1.623 +               limit = glp_get_row_name(P, var2);
   1.624 +            else
   1.625 +               limit = glp_get_col_name(P, var2 - m);
   1.626 +            if (limit != NULL)
   1.627 +               xfprintf(fp, " %s", limit);
   1.628 +         }
   1.629 +         xfprintf(fp, "\n");
   1.630 +         xfprintf(fp, "\n");
   1.631 +         /* print 10 items per page */
   1.632 +         count = (count + 1) % 10;
   1.633 +      }
   1.634 +      xfprintf(fp, "End of report\n");
   1.635 +      xfflush(fp);
   1.636 +      if (xferror(fp))
   1.637 +      {  xprintf("Write error on `%s' - %s\n", fname, xerrmsg());
   1.638 +         ret = 4;
   1.639 +         goto done;
   1.640 +      }
   1.641 +      ret = 0;
   1.642 +done: if (fp != NULL) xfclose(fp);
   1.643 +      return ret;
   1.644 +}
   1.645 +
   1.646 +/**********************************************************************/
   1.647 +
   1.648 +int glp_print_ipt(glp_prob *P, const char *fname)
   1.649 +{     /* write interior-point solution in printable format */
   1.650 +      XFILE *fp;
   1.651 +      GLPROW *row;
   1.652 +      GLPCOL *col;
   1.653 +      int i, j, t, ae_ind, re_ind, ret;
   1.654 +      double ae_max, re_max;
   1.655 +      xprintf("Writing interior-point solution to `%s'...\n", fname);
   1.656 +      fp = xfopen(fname, "w");
   1.657 +      if (fp == NULL)
   1.658 +      {  xprintf("Unable to create `%s' - %s\n", fname, xerrmsg());
   1.659 +         ret = 1;
   1.660 +         goto done;
   1.661 +      }
   1.662 +      xfprintf(fp, "%-12s%s\n", "Problem:",
   1.663 +         P->name == NULL ? "" : P->name);
   1.664 +      xfprintf(fp, "%-12s%d\n", "Rows:", P->m);
   1.665 +      xfprintf(fp, "%-12s%d\n", "Columns:", P->n);
   1.666 +      xfprintf(fp, "%-12s%d\n", "Non-zeros:", P->nnz);
   1.667 +      t = glp_ipt_status(P);
   1.668 +      xfprintf(fp, "%-12s%s\n", "Status:",
   1.669 +         t == GLP_OPT    ? "OPTIMAL" :
   1.670 +         t == GLP_UNDEF  ? "UNDEFINED" :
   1.671 +         t == GLP_INFEAS ? "INFEASIBLE (INTERMEDIATE)" :
   1.672 +         t == GLP_NOFEAS ? "INFEASIBLE (FINAL)" : "???");
   1.673 +      xfprintf(fp, "%-12s%s%s%.10g (%s)\n", "Objective:",
   1.674 +         P->obj == NULL ? "" : P->obj,
   1.675 +         P->obj == NULL ? "" : " = ", P->ipt_obj,
   1.676 +         P->dir == GLP_MIN ? "MINimum" :
   1.677 +         P->dir == GLP_MAX ? "MAXimum" : "???");
   1.678 +      xfprintf(fp, "\n");
   1.679 +      xfprintf(fp, "   No.   Row name        Activity     Lower bound  "
   1.680 +         " Upper bound    Marginal\n");
   1.681 +      xfprintf(fp, "------ ------------    ------------- ------------- "
   1.682 +         "------------- -------------\n");
   1.683 +      for (i = 1; i <= P->m; i++)
   1.684 +      {  row = P->row[i];
   1.685 +         xfprintf(fp, "%6d ", i);
   1.686 +         if (row->name == NULL || strlen(row->name) <= 12)
   1.687 +            xfprintf(fp, "%-12s ", row->name == NULL ? "" : row->name);
   1.688 +         else
   1.689 +            xfprintf(fp, "%s\n%20s", row->name, "");
   1.690 +         xfprintf(fp, "%3s", "");
   1.691 +         xfprintf(fp, "%13.6g ",
   1.692 +            fabs(row->pval) <= 1e-9 ? 0.0 : row->pval);
   1.693 +         if (row->type == GLP_LO || row->type == GLP_DB ||
   1.694 +             row->type == GLP_FX)
   1.695 +            xfprintf(fp, "%13.6g ", row->lb);
   1.696 +         else
   1.697 +            xfprintf(fp, "%13s ", "");
   1.698 +         if (row->type == GLP_UP || row->type == GLP_DB)
   1.699 +            xfprintf(fp, "%13.6g ", row->ub);
   1.700 +         else
   1.701 +            xfprintf(fp, "%13s ", row->type == GLP_FX ? "=" : "");
   1.702 +         if (fabs(row->dval) <= 1e-9)
   1.703 +            xfprintf(fp, "%13s", "< eps");
   1.704 +         else
   1.705 +            xfprintf(fp, "%13.6g ", row->dval);
   1.706 +         xfprintf(fp, "\n");
   1.707 +      }
   1.708 +      xfprintf(fp, "\n");
   1.709 +      xfprintf(fp, "   No. Column name       Activity     Lower bound  "
   1.710 +         " Upper bound    Marginal\n");
   1.711 +      xfprintf(fp, "------ ------------    ------------- ------------- "
   1.712 +         "------------- -------------\n");
   1.713 +      for (j = 1; j <= P->n; j++)
   1.714 +      {  col = P->col[j];
   1.715 +         xfprintf(fp, "%6d ", j);
   1.716 +         if (col->name == NULL || strlen(col->name) <= 12)
   1.717 +            xfprintf(fp, "%-12s ", col->name == NULL ? "" : col->name);
   1.718 +         else
   1.719 +            xfprintf(fp, "%s\n%20s", col->name, "");
   1.720 +         xfprintf(fp, "%3s", "");
   1.721 +         xfprintf(fp, "%13.6g ",
   1.722 +            fabs(col->pval) <= 1e-9 ? 0.0 : col->pval);
   1.723 +         if (col->type == GLP_LO || col->type == GLP_DB ||
   1.724 +             col->type == GLP_FX)
   1.725 +            xfprintf(fp, "%13.6g ", col->lb);
   1.726 +         else
   1.727 +            xfprintf(fp, "%13s ", "");
   1.728 +         if (col->type == GLP_UP || col->type == GLP_DB)
   1.729 +            xfprintf(fp, "%13.6g ", col->ub);
   1.730 +         else
   1.731 +            xfprintf(fp, "%13s ", col->type == GLP_FX ? "=" : "");
   1.732 +         if (fabs(col->dval) <= 1e-9)
   1.733 +            xfprintf(fp, "%13s", "< eps");
   1.734 +         else
   1.735 +            xfprintf(fp, "%13.6g ", col->dval);
   1.736 +         xfprintf(fp, "\n");
   1.737 +      }
   1.738 +      xfprintf(fp, "\n");
   1.739 +      xfprintf(fp, "Karush-Kuhn-Tucker optimality conditions:\n");
   1.740 +      xfprintf(fp, "\n");
   1.741 +      _glp_check_kkt(P, GLP_IPT, GLP_KKT_PE, &ae_max, &ae_ind, &re_max,
   1.742 +         &re_ind);
   1.743 +      xfprintf(fp, "KKT.PE: max.abs.err = %.2e on row %d\n",
   1.744 +         ae_max, ae_ind);
   1.745 +      xfprintf(fp, "        max.rel.err = %.2e on row %d\n",
   1.746 +         re_max, re_ind);
   1.747 +      xfprintf(fp, "%8s%s\n", "",
   1.748 +         re_max <= 1e-9 ? "High quality" :
   1.749 +         re_max <= 1e-6 ? "Medium quality" :
   1.750 +         re_max <= 1e-3 ? "Low quality" : "PRIMAL SOLUTION IS WRONG");
   1.751 +      xfprintf(fp, "\n");
   1.752 +      _glp_check_kkt(P, GLP_IPT, GLP_KKT_PB, &ae_max, &ae_ind, &re_max,
   1.753 +         &re_ind);
   1.754 +      xfprintf(fp, "KKT.PB: max.abs.err = %.2e on %s %d\n",
   1.755 +            ae_max, ae_ind <= P->m ? "row" : "column",
   1.756 +            ae_ind <= P->m ? ae_ind : ae_ind - P->m);
   1.757 +      xfprintf(fp, "        max.rel.err = %.2e on %s %d\n",
   1.758 +            re_max, re_ind <= P->m ? "row" : "column",
   1.759 +            re_ind <= P->m ? re_ind : re_ind - P->m);
   1.760 +      xfprintf(fp, "%8s%s\n", "",
   1.761 +         re_max <= 1e-9 ? "High quality" :
   1.762 +         re_max <= 1e-6 ? "Medium quality" :
   1.763 +         re_max <= 1e-3 ? "Low quality" : "PRIMAL SOLUTION IS INFEASIBL"
   1.764 +            "E");
   1.765 +      xfprintf(fp, "\n");
   1.766 +      _glp_check_kkt(P, GLP_IPT, GLP_KKT_DE, &ae_max, &ae_ind, &re_max,
   1.767 +         &re_ind);
   1.768 +      xfprintf(fp, "KKT.DE: max.abs.err = %.2e on column %d\n",
   1.769 +         ae_max, ae_ind == 0 ? 0 : ae_ind - P->m);
   1.770 +      xfprintf(fp, "        max.rel.err = %.2e on column %d\n",
   1.771 +         re_max, re_ind == 0 ? 0 : re_ind - P->m);
   1.772 +      xfprintf(fp, "%8s%s\n", "",
   1.773 +         re_max <= 1e-9 ? "High quality" :
   1.774 +         re_max <= 1e-6 ? "Medium quality" :
   1.775 +         re_max <= 1e-3 ? "Low quality" : "DUAL SOLUTION IS WRONG");
   1.776 +      xfprintf(fp, "\n");
   1.777 +      _glp_check_kkt(P, GLP_IPT, GLP_KKT_DB, &ae_max, &ae_ind, &re_max,
   1.778 +         &re_ind);
   1.779 +      xfprintf(fp, "KKT.DB: max.abs.err = %.2e on %s %d\n",
   1.780 +            ae_max, ae_ind <= P->m ? "row" : "column",
   1.781 +            ae_ind <= P->m ? ae_ind : ae_ind - P->m);
   1.782 +      xfprintf(fp, "        max.rel.err = %.2e on %s %d\n",
   1.783 +            re_max, re_ind <= P->m ? "row" : "column",
   1.784 +            re_ind <= P->m ? re_ind : re_ind - P->m);
   1.785 +      xfprintf(fp, "%8s%s\n", "",
   1.786 +         re_max <= 1e-9 ? "High quality" :
   1.787 +         re_max <= 1e-6 ? "Medium quality" :
   1.788 +         re_max <= 1e-3 ? "Low quality" : "DUAL SOLUTION IS INFEASIBLE")
   1.789 +            ;
   1.790 +      xfprintf(fp, "\n");
   1.791 +      xfprintf(fp, "End of output\n");
   1.792 +      xfflush(fp);
   1.793 +      if (xferror(fp))
   1.794 +      {  xprintf("Write error on `%s' - %s\n", fname, xerrmsg());
   1.795 +         ret = 1;
   1.796 +         goto done;
   1.797 +      }
   1.798 +      ret = 0;
   1.799 +done: if (fp != NULL) xfclose(fp);
   1.800 +      return ret;
   1.801 +}
   1.802 +
   1.803 +/***********************************************************************
   1.804 +*  NAME
   1.805 +*
   1.806 +*  glp_read_ipt - read interior-point solution from text file
   1.807 +*
   1.808 +*  SYNOPSIS
   1.809 +*
   1.810 +*  int glp_read_ipt(glp_prob *lp, const char *fname);
   1.811 +*
   1.812 +*  DESCRIPTION
   1.813 +*
   1.814 +*  The routine glp_read_ipt reads interior-point solution from a text
   1.815 +*  file whose name is specified by the parameter fname into the problem
   1.816 +*  object.
   1.817 +*
   1.818 +*  For the file format see description of the routine glp_write_ipt.
   1.819 +*
   1.820 +*  RETURNS
   1.821 +*
   1.822 +*  On success the routine returns zero, otherwise non-zero. */
   1.823 +
   1.824 +int glp_read_ipt(glp_prob *lp, const char *fname)
   1.825 +{     glp_data *data;
   1.826 +      jmp_buf jump;
   1.827 +      int i, j, k, ret = 0;
   1.828 +      xprintf("Reading interior-point solution from `%s'...\n", fname);
   1.829 +      data = glp_sdf_open_file(fname);
   1.830 +      if (data == NULL)
   1.831 +      {  ret = 1;
   1.832 +         goto done;
   1.833 +      }
   1.834 +      if (setjmp(jump))
   1.835 +      {  ret = 1;
   1.836 +         goto done;
   1.837 +      }
   1.838 +      glp_sdf_set_jump(data, jump);
   1.839 +      /* number of rows, number of columns */
   1.840 +      k = glp_sdf_read_int(data);
   1.841 +      if (k != lp->m)
   1.842 +         glp_sdf_error(data, "wrong number of rows\n");
   1.843 +      k = glp_sdf_read_int(data);
   1.844 +      if (k != lp->n)
   1.845 +         glp_sdf_error(data, "wrong number of columns\n");
   1.846 +      /* solution status, objective value */
   1.847 +      k = glp_sdf_read_int(data);
   1.848 +      if (!(k == GLP_UNDEF || k == GLP_OPT))
   1.849 +         glp_sdf_error(data, "invalid solution status\n");
   1.850 +      lp->ipt_stat = k;
   1.851 +      lp->ipt_obj = glp_sdf_read_num(data);
   1.852 +      /* rows (auxiliary variables) */
   1.853 +      for (i = 1; i <= lp->m; i++)
   1.854 +      {  GLPROW *row = lp->row[i];
   1.855 +         /* primal value, dual value */
   1.856 +         row->pval = glp_sdf_read_num(data);
   1.857 +         row->dval = glp_sdf_read_num(data);
   1.858 +      }
   1.859 +      /* columns (structural variables) */
   1.860 +      for (j = 1; j <= lp->n; j++)
   1.861 +      {  GLPCOL *col = lp->col[j];
   1.862 +         /* primal value, dual value */
   1.863 +         col->pval = glp_sdf_read_num(data);
   1.864 +         col->dval = glp_sdf_read_num(data);
   1.865 +      }
   1.866 +      xprintf("%d lines were read\n", glp_sdf_line(data));
   1.867 +done: if (ret) lp->ipt_stat = GLP_UNDEF;
   1.868 +      if (data != NULL) glp_sdf_close_file(data);
   1.869 +      return ret;
   1.870 +}
   1.871 +
   1.872 +/***********************************************************************
   1.873 +*  NAME
   1.874 +*
   1.875 +*  glp_write_ipt - write interior-point solution to text file
   1.876 +*
   1.877 +*  SYNOPSIS
   1.878 +*
   1.879 +*  int glp_write_ipt(glp_prob *lp, const char *fname);
   1.880 +*
   1.881 +*  DESCRIPTION
   1.882 +*
   1.883 +*  The routine glp_write_ipt writes the current interior-point solution
   1.884 +*  to a text file whose name is specified by the parameter fname. This
   1.885 +*  file can be read back with the routine glp_read_ipt.
   1.886 +*
   1.887 +*  RETURNS
   1.888 +*
   1.889 +*  On success the routine returns zero, otherwise non-zero.
   1.890 +*
   1.891 +*  FILE FORMAT
   1.892 +*
   1.893 +*  The file created by the routine glp_write_ipt is a plain text file,
   1.894 +*  which contains the following information:
   1.895 +*
   1.896 +*     m n
   1.897 +*     stat obj_val
   1.898 +*     r_prim[1] r_dual[1]
   1.899 +*     . . .
   1.900 +*     r_prim[m] r_dual[m]
   1.901 +*     c_prim[1] c_dual[1]
   1.902 +*     . . .
   1.903 +*     c_prim[n] c_dual[n]
   1.904 +*
   1.905 +*  where:
   1.906 +*  m is the number of rows (auxiliary variables);
   1.907 +*  n is the number of columns (structural variables);
   1.908 +*  stat is the solution status (GLP_UNDEF = 1 or GLP_OPT = 5);
   1.909 +*  obj_val is the objective value;
   1.910 +*  r_prim[i], i = 1,...,m, is the primal value of i-th row;
   1.911 +*  r_dual[i], i = 1,...,m, is the dual value of i-th row;
   1.912 +*  c_prim[j], j = 1,...,n, is the primal value of j-th column;
   1.913 +*  c_dual[j], j = 1,...,n, is the dual value of j-th column. */
   1.914 +
   1.915 +int glp_write_ipt(glp_prob *lp, const char *fname)
   1.916 +{     XFILE *fp;
   1.917 +      int i, j, ret = 0;
   1.918 +      xprintf("Writing interior-point solution to `%s'...\n", fname);
   1.919 +      fp = xfopen(fname, "w");
   1.920 +      if (fp == NULL)
   1.921 +      {  xprintf("Unable to create `%s' - %s\n", fname, xerrmsg());
   1.922 +         ret = 1;
   1.923 +         goto done;
   1.924 +      }
   1.925 +      /* number of rows, number of columns */
   1.926 +      xfprintf(fp, "%d %d\n", lp->m, lp->n);
   1.927 +      /* solution status, objective value */
   1.928 +      xfprintf(fp, "%d %.*g\n", lp->ipt_stat, DBL_DIG, lp->ipt_obj);
   1.929 +      /* rows (auxiliary variables) */
   1.930 +      for (i = 1; i <= lp->m; i++)
   1.931 +      {  GLPROW *row = lp->row[i];
   1.932 +         /* primal value, dual value */
   1.933 +         xfprintf(fp, "%.*g %.*g\n", DBL_DIG, row->pval, DBL_DIG,
   1.934 +            row->dval);
   1.935 +      }
   1.936 +      /* columns (structural variables) */
   1.937 +      for (j = 1; j <= lp->n; j++)
   1.938 +      {  GLPCOL *col = lp->col[j];
   1.939 +         /* primal value, dual value */
   1.940 +         xfprintf(fp, "%.*g %.*g\n", DBL_DIG, col->pval, DBL_DIG,
   1.941 +            col->dval);
   1.942 +      }
   1.943 +      xfflush(fp);
   1.944 +      if (xferror(fp))
   1.945 +      {  xprintf("Write error on `%s' - %s\n", fname, xerrmsg());
   1.946 +         ret = 1;
   1.947 +         goto done;
   1.948 +      }
   1.949 +      xprintf("%d lines were written\n", 2 + lp->m + lp->n);
   1.950 +done: if (fp != NULL) xfclose(fp);
   1.951 +      return ret;
   1.952 +}
   1.953 +
   1.954 +/**********************************************************************/
   1.955 +
   1.956 +int glp_print_mip(glp_prob *P, const char *fname)
   1.957 +{     /* write MIP solution in printable format */
   1.958 +      XFILE *fp;
   1.959 +      GLPROW *row;
   1.960 +      GLPCOL *col;
   1.961 +      int i, j, t, ae_ind, re_ind, ret;
   1.962 +      double ae_max, re_max;
   1.963 +      xprintf("Writing MIP solution to `%s'...\n", fname);
   1.964 +      fp = xfopen(fname, "w");
   1.965 +      if (fp == NULL)
   1.966 +      {  xprintf("Unable to create `%s' - %s\n", fname, xerrmsg());
   1.967 +         ret = 1;
   1.968 +         goto done;
   1.969 +      }
   1.970 +      xfprintf(fp, "%-12s%s\n", "Problem:",
   1.971 +         P->name == NULL ? "" : P->name);
   1.972 +      xfprintf(fp, "%-12s%d\n", "Rows:", P->m);
   1.973 +      xfprintf(fp, "%-12s%d (%d integer, %d binary)\n", "Columns:",
   1.974 +         P->n, glp_get_num_int(P), glp_get_num_bin(P));
   1.975 +      xfprintf(fp, "%-12s%d\n", "Non-zeros:", P->nnz);
   1.976 +      t = glp_mip_status(P);
   1.977 +      xfprintf(fp, "%-12s%s\n", "Status:",
   1.978 +         t == GLP_OPT    ? "INTEGER OPTIMAL" :
   1.979 +         t == GLP_FEAS   ? "INTEGER NON-OPTIMAL" :
   1.980 +         t == GLP_NOFEAS ? "INTEGER EMPTY" :
   1.981 +         t == GLP_UNDEF  ? "INTEGER UNDEFINED" : "???");
   1.982 +      xfprintf(fp, "%-12s%s%s%.10g (%s)\n", "Objective:",
   1.983 +         P->obj == NULL ? "" : P->obj,
   1.984 +         P->obj == NULL ? "" : " = ", P->mip_obj,
   1.985 +         P->dir == GLP_MIN ? "MINimum" :
   1.986 +         P->dir == GLP_MAX ? "MAXimum" : "???");
   1.987 +      xfprintf(fp, "\n");
   1.988 +      xfprintf(fp, "   No.   Row name        Activity     Lower bound  "
   1.989 +         " Upper bound\n");
   1.990 +      xfprintf(fp, "------ ------------    ------------- ------------- "
   1.991 +         "-------------\n");
   1.992 +      for (i = 1; i <= P->m; i++)
   1.993 +      {  row = P->row[i];
   1.994 +         xfprintf(fp, "%6d ", i);
   1.995 +         if (row->name == NULL || strlen(row->name) <= 12)
   1.996 +            xfprintf(fp, "%-12s ", row->name == NULL ? "" : row->name);
   1.997 +         else
   1.998 +            xfprintf(fp, "%s\n%20s", row->name, "");
   1.999 +         xfprintf(fp, "%3s", "");
  1.1000 +         xfprintf(fp, "%13.6g ",
  1.1001 +            fabs(row->mipx) <= 1e-9 ? 0.0 : row->mipx);
  1.1002 +         if (row->type == GLP_LO || row->type == GLP_DB ||
  1.1003 +             row->type == GLP_FX)
  1.1004 +            xfprintf(fp, "%13.6g ", row->lb);
  1.1005 +         else
  1.1006 +            xfprintf(fp, "%13s ", "");
  1.1007 +         if (row->type == GLP_UP || row->type == GLP_DB)
  1.1008 +            xfprintf(fp, "%13.6g ", row->ub);
  1.1009 +         else
  1.1010 +            xfprintf(fp, "%13s ", row->type == GLP_FX ? "=" : "");
  1.1011 +         xfprintf(fp, "\n");
  1.1012 +      }
  1.1013 +      xfprintf(fp, "\n");
  1.1014 +      xfprintf(fp, "   No. Column name       Activity     Lower bound  "
  1.1015 +         " Upper bound\n");
  1.1016 +      xfprintf(fp, "------ ------------    ------------- ------------- "
  1.1017 +         "-------------\n");
  1.1018 +      for (j = 1; j <= P->n; j++)
  1.1019 +      {  col = P->col[j];
  1.1020 +         xfprintf(fp, "%6d ", j);
  1.1021 +         if (col->name == NULL || strlen(col->name) <= 12)
  1.1022 +            xfprintf(fp, "%-12s ", col->name == NULL ? "" : col->name);
  1.1023 +         else
  1.1024 +            xfprintf(fp, "%s\n%20s", col->name, "");
  1.1025 +         xfprintf(fp, "%s  ",
  1.1026 +            col->kind == GLP_CV ? " " :
  1.1027 +            col->kind == GLP_IV ? "*" : "?");
  1.1028 +         xfprintf(fp, "%13.6g ",
  1.1029 +            fabs(col->mipx) <= 1e-9 ? 0.0 : col->mipx);
  1.1030 +         if (col->type == GLP_LO || col->type == GLP_DB ||
  1.1031 +             col->type == GLP_FX)
  1.1032 +            xfprintf(fp, "%13.6g ", col->lb);
  1.1033 +         else
  1.1034 +            xfprintf(fp, "%13s ", "");
  1.1035 +         if (col->type == GLP_UP || col->type == GLP_DB)
  1.1036 +            xfprintf(fp, "%13.6g ", col->ub);
  1.1037 +         else
  1.1038 +            xfprintf(fp, "%13s ", col->type == GLP_FX ? "=" : "");
  1.1039 +         xfprintf(fp, "\n");
  1.1040 +      }
  1.1041 +      xfprintf(fp, "\n");
  1.1042 +      xfprintf(fp, "Integer feasibility conditions:\n");
  1.1043 +      xfprintf(fp, "\n");
  1.1044 +      _glp_check_kkt(P, GLP_MIP, GLP_KKT_PE, &ae_max, &ae_ind, &re_max,
  1.1045 +         &re_ind);
  1.1046 +      xfprintf(fp, "KKT.PE: max.abs.err = %.2e on row %d\n",
  1.1047 +         ae_max, ae_ind);
  1.1048 +      xfprintf(fp, "        max.rel.err = %.2e on row %d\n",
  1.1049 +         re_max, re_ind);
  1.1050 +      xfprintf(fp, "%8s%s\n", "",
  1.1051 +         re_max <= 1e-9 ? "High quality" :
  1.1052 +         re_max <= 1e-6 ? "Medium quality" :
  1.1053 +         re_max <= 1e-3 ? "Low quality" : "SOLUTION IS WRONG");
  1.1054 +      xfprintf(fp, "\n");
  1.1055 +      _glp_check_kkt(P, GLP_MIP, GLP_KKT_PB, &ae_max, &ae_ind, &re_max,
  1.1056 +         &re_ind);
  1.1057 +      xfprintf(fp, "KKT.PB: max.abs.err = %.2e on %s %d\n",
  1.1058 +            ae_max, ae_ind <= P->m ? "row" : "column",
  1.1059 +            ae_ind <= P->m ? ae_ind : ae_ind - P->m);
  1.1060 +      xfprintf(fp, "        max.rel.err = %.2e on %s %d\n",
  1.1061 +            re_max, re_ind <= P->m ? "row" : "column",
  1.1062 +            re_ind <= P->m ? re_ind : re_ind - P->m);
  1.1063 +      xfprintf(fp, "%8s%s\n", "",
  1.1064 +         re_max <= 1e-9 ? "High quality" :
  1.1065 +         re_max <= 1e-6 ? "Medium quality" :
  1.1066 +         re_max <= 1e-3 ? "Low quality" : "SOLUTION IS INFEASIBLE");
  1.1067 +      xfprintf(fp, "\n");
  1.1068 +      xfprintf(fp, "End of output\n");
  1.1069 +      xfflush(fp);
  1.1070 +      if (xferror(fp))
  1.1071 +      {  xprintf("Write error on `%s' - %s\n", fname, xerrmsg());
  1.1072 +         ret = 1;
  1.1073 +         goto done;
  1.1074 +      }
  1.1075 +      ret = 0;
  1.1076 +done: if (fp != NULL) xfclose(fp);
  1.1077 +      return ret;
  1.1078 +}
  1.1079 +
  1.1080 +/***********************************************************************
  1.1081 +*  NAME
  1.1082 +*
  1.1083 +*  glp_read_mip - read MIP solution from text file
  1.1084 +*
  1.1085 +*  SYNOPSIS
  1.1086 +*
  1.1087 +*  int glp_read_mip(glp_prob *mip, const char *fname);
  1.1088 +*
  1.1089 +*  DESCRIPTION
  1.1090 +*
  1.1091 +*  The routine glp_read_mip reads MIP solution from a text file whose
  1.1092 +*  name is specified by the parameter fname into the problem object.
  1.1093 +*
  1.1094 +*  For the file format see description of the routine glp_write_mip.
  1.1095 +*
  1.1096 +*  RETURNS
  1.1097 +*
  1.1098 +*  On success the routine returns zero, otherwise non-zero. */
  1.1099 +
  1.1100 +int glp_read_mip(glp_prob *mip, const char *fname)
  1.1101 +{     glp_data *data;
  1.1102 +      jmp_buf jump;
  1.1103 +      int i, j, k, ret = 0;
  1.1104 +      xprintf("Reading MIP solution from `%s'...\n", fname);
  1.1105 +      data = glp_sdf_open_file(fname);
  1.1106 +      if (data == NULL)
  1.1107 +      {  ret = 1;
  1.1108 +         goto done;
  1.1109 +      }
  1.1110 +      if (setjmp(jump))
  1.1111 +      {  ret = 1;
  1.1112 +         goto done;
  1.1113 +      }
  1.1114 +      glp_sdf_set_jump(data, jump);
  1.1115 +      /* number of rows, number of columns */
  1.1116 +      k = glp_sdf_read_int(data);
  1.1117 +      if (k != mip->m)
  1.1118 +         glp_sdf_error(data, "wrong number of rows\n");
  1.1119 +      k = glp_sdf_read_int(data);
  1.1120 +      if (k != mip->n)
  1.1121 +         glp_sdf_error(data, "wrong number of columns\n");
  1.1122 +      /* solution status, objective value */
  1.1123 +      k = glp_sdf_read_int(data);
  1.1124 +      if (!(k == GLP_UNDEF || k == GLP_OPT || k == GLP_FEAS ||
  1.1125 +            k == GLP_NOFEAS))
  1.1126 +         glp_sdf_error(data, "invalid solution status\n");
  1.1127 +      mip->mip_stat = k;
  1.1128 +      mip->mip_obj = glp_sdf_read_num(data);
  1.1129 +      /* rows (auxiliary variables) */
  1.1130 +      for (i = 1; i <= mip->m; i++)
  1.1131 +      {  GLPROW *row = mip->row[i];
  1.1132 +         row->mipx = glp_sdf_read_num(data);
  1.1133 +      }
  1.1134 +      /* columns (structural variables) */
  1.1135 +      for (j = 1; j <= mip->n; j++)
  1.1136 +      {  GLPCOL *col = mip->col[j];
  1.1137 +         col->mipx = glp_sdf_read_num(data);
  1.1138 +         if (col->kind == GLP_IV && col->mipx != floor(col->mipx))
  1.1139 +            glp_sdf_error(data, "non-integer column value");
  1.1140 +      }
  1.1141 +      xprintf("%d lines were read\n", glp_sdf_line(data));
  1.1142 +done: if (ret) mip->mip_stat = GLP_UNDEF;
  1.1143 +      if (data != NULL) glp_sdf_close_file(data);
  1.1144 +      return ret;
  1.1145 +}
  1.1146 +
  1.1147 +/***********************************************************************
  1.1148 +*  NAME
  1.1149 +*
  1.1150 +*  glp_write_mip - write MIP solution to text file
  1.1151 +*
  1.1152 +*  SYNOPSIS
  1.1153 +*
  1.1154 +*  int glp_write_mip(glp_prob *mip, const char *fname);
  1.1155 +*
  1.1156 +*  DESCRIPTION
  1.1157 +*
  1.1158 +*  The routine glp_write_mip writes the current MIP solution to a text
  1.1159 +*  file whose name is specified by the parameter fname. This file can
  1.1160 +*  be read back with the routine glp_read_mip.
  1.1161 +*
  1.1162 +*  RETURNS
  1.1163 +*
  1.1164 +*  On success the routine returns zero, otherwise non-zero.
  1.1165 +*
  1.1166 +*  FILE FORMAT
  1.1167 +*
  1.1168 +*  The file created by the routine glp_write_sol is a plain text file,
  1.1169 +*  which contains the following information:
  1.1170 +*
  1.1171 +*     m n
  1.1172 +*     stat obj_val
  1.1173 +*     r_val[1]
  1.1174 +*     . . .
  1.1175 +*     r_val[m]
  1.1176 +*     c_val[1]
  1.1177 +*     . . .
  1.1178 +*     c_val[n]
  1.1179 +*
  1.1180 +*  where:
  1.1181 +*  m is the number of rows (auxiliary variables);
  1.1182 +*  n is the number of columns (structural variables);
  1.1183 +*  stat is the solution status (GLP_UNDEF = 1, GLP_FEAS = 2,
  1.1184 +*     GLP_NOFEAS = 4, or GLP_OPT = 5);
  1.1185 +*  obj_val is the objective value;
  1.1186 +*  r_val[i], i = 1,...,m, is the value of i-th row;
  1.1187 +*  c_val[j], j = 1,...,n, is the value of j-th column. */
  1.1188 +
  1.1189 +int glp_write_mip(glp_prob *mip, const char *fname)
  1.1190 +{     XFILE *fp;
  1.1191 +      int i, j, ret = 0;
  1.1192 +      xprintf("Writing MIP solution to `%s'...\n", fname);
  1.1193 +      fp = xfopen(fname, "w");
  1.1194 +      if (fp == NULL)
  1.1195 +      {  xprintf("Unable to create `%s' - %s\n", fname, xerrmsg());
  1.1196 +         ret = 1;
  1.1197 +         goto done;
  1.1198 +      }
  1.1199 +      /* number of rows, number of columns */
  1.1200 +      xfprintf(fp, "%d %d\n", mip->m, mip->n);
  1.1201 +      /* solution status, objective value */
  1.1202 +      xfprintf(fp, "%d %.*g\n", mip->mip_stat, DBL_DIG, mip->mip_obj);
  1.1203 +      /* rows (auxiliary variables) */
  1.1204 +      for (i = 1; i <= mip->m; i++)
  1.1205 +         xfprintf(fp, "%.*g\n", DBL_DIG, mip->row[i]->mipx);
  1.1206 +      /* columns (structural variables) */
  1.1207 +      for (j = 1; j <= mip->n; j++)
  1.1208 +         xfprintf(fp, "%.*g\n", DBL_DIG, mip->col[j]->mipx);
  1.1209 +      xfflush(fp);
  1.1210 +      if (xferror(fp))
  1.1211 +      {  xprintf("Write error on `%s' - %s\n", fname, xerrmsg());
  1.1212 +         ret = 1;
  1.1213 +         goto done;
  1.1214 +      }
  1.1215 +      xprintf("%d lines were written\n", 2 + mip->m + mip->n);
  1.1216 +done: if (fp != NULL) xfclose(fp);
  1.1217 +      return ret;
  1.1218 +}
  1.1219 +
  1.1220 +/* eof */