src/glpmps.c
changeset 1 c445c931472f
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/src/glpmps.c	Mon Dec 06 13:09:21 2010 +0100
     1.3 @@ -0,0 +1,1401 @@
     1.4 +/* glpmps.c (MPS format 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 +/***********************************************************************
    1.31 +*  NAME
    1.32 +*
    1.33 +*  glp_init_mpscp - initialize MPS format control parameters
    1.34 +*
    1.35 +*  SYNOPSIS
    1.36 +*
    1.37 +*  void glp_init_mpscp(glp_mpscp *parm);
    1.38 +*
    1.39 +*  DESCRIPTION
    1.40 +*
    1.41 +*  The routine glp_init_mpscp initializes control parameters, which are
    1.42 +*  used by the MPS input/output routines glp_read_mps and glp_write_mps,
    1.43 +*  with default values.
    1.44 +*
    1.45 +*  Default values of the control parameters are stored in the glp_mpscp
    1.46 +*  structure, which the parameter parm points to. */
    1.47 +
    1.48 +void glp_init_mpscp(glp_mpscp *parm)
    1.49 +{     parm->blank = '\0';
    1.50 +      parm->obj_name = NULL;
    1.51 +      parm->tol_mps = 1e-12;
    1.52 +      return;
    1.53 +}
    1.54 +
    1.55 +static void check_parm(const char *func, const glp_mpscp *parm)
    1.56 +{     /* check control parameters */
    1.57 +      if (!(0x00 <= parm->blank && parm->blank <= 0xFF) ||
    1.58 +          !(parm->blank == '\0' || isprint(parm->blank)))
    1.59 +         xerror("%s: blank = 0x%02X; invalid parameter\n",
    1.60 +            func, parm->blank);
    1.61 +      if (!(parm->obj_name == NULL || strlen(parm->obj_name) <= 255))
    1.62 +         xerror("%s: obj_name = \"%.12s...\"; parameter too long\n",
    1.63 +            func, parm->obj_name);
    1.64 +      if (!(0.0 <= parm->tol_mps && parm->tol_mps < 1.0))
    1.65 +         xerror("%s: tol_mps = %g; invalid parameter\n",
    1.66 +            func, parm->tol_mps);
    1.67 +      return;
    1.68 +}
    1.69 +
    1.70 +/***********************************************************************
    1.71 +*  NAME
    1.72 +*
    1.73 +*  glp_read_mps - read problem data in MPS format
    1.74 +*
    1.75 +*  SYNOPSIS
    1.76 +*
    1.77 +*  int glp_read_mps(glp_prob *P, int fmt, const glp_mpscp *parm,
    1.78 +*     const char *fname);
    1.79 +*
    1.80 +*  DESCRIPTION
    1.81 +*
    1.82 +*  The routine glp_read_mps reads problem data in MPS format from a
    1.83 +*  text file.
    1.84 +*
    1.85 +*  The parameter fmt specifies the version of MPS format:
    1.86 +*
    1.87 +*  GLP_MPS_DECK - fixed (ancient) MPS format;
    1.88 +*  GLP_MPS_FILE - free (modern) MPS format.
    1.89 +*
    1.90 +*  The parameter parm is a pointer to the structure glp_mpscp, which
    1.91 +*  specifies control parameters used by the routine. If parm is NULL,
    1.92 +*  the routine uses default settings.
    1.93 +*
    1.94 +*  The character string fname specifies a name of the text file to be
    1.95 +*  read.
    1.96 +*
    1.97 +*  Note that before reading data the current content of the problem
    1.98 +*  object is completely erased with the routine glp_erase_prob.
    1.99 +*
   1.100 +*  RETURNS
   1.101 +*
   1.102 +*  If the operation was successful, the routine glp_read_mps returns
   1.103 +*  zero. Otherwise, it prints an error message and returns non-zero. */
   1.104 +
   1.105 +struct csa
   1.106 +{     /* common storage area */
   1.107 +      glp_prob *P;
   1.108 +      /* pointer to problem object */
   1.109 +      int deck;
   1.110 +      /* MPS format (0 - free, 1 - fixed) */
   1.111 +      const glp_mpscp *parm;
   1.112 +      /* pointer to control parameters */
   1.113 +      const char *fname;
   1.114 +      /* name of input MPS file */
   1.115 +      XFILE *fp;
   1.116 +      /* stream assigned to input MPS file */
   1.117 +      jmp_buf jump;
   1.118 +      /* label for go to in case of error */
   1.119 +      int recno;
   1.120 +      /* current record (card) number */
   1.121 +      int recpos;
   1.122 +      /* current record (card) position */
   1.123 +      int c;
   1.124 +      /* current character */
   1.125 +      int fldno;
   1.126 +      /* current field number */
   1.127 +      char field[255+1];
   1.128 +      /* current field content */
   1.129 +      int w80;
   1.130 +      /* warning 'record must not be longer than 80 chars' issued */
   1.131 +      int wef;
   1.132 +      /* warning 'extra fields detected beyond field 6' issued */
   1.133 +      int obj_row;
   1.134 +      /* objective row number */
   1.135 +      void *work1, *work2, *work3;
   1.136 +      /* working arrays */
   1.137 +};
   1.138 +
   1.139 +static void error(struct csa *csa, const char *fmt, ...)
   1.140 +{     /* print error message and terminate processing */
   1.141 +      va_list arg;
   1.142 +      xprintf("%s:%d: ", csa->fname, csa->recno);
   1.143 +      va_start(arg, fmt);
   1.144 +      xvprintf(fmt, arg);
   1.145 +      va_end(arg);
   1.146 +      longjmp(csa->jump, 1);
   1.147 +      /* no return */
   1.148 +}
   1.149 +
   1.150 +static void warning(struct csa *csa, const char *fmt, ...)
   1.151 +{     /* print warning message and continue processing */
   1.152 +      va_list arg;
   1.153 +      xprintf("%s:%d: warning: ", csa->fname, csa->recno);
   1.154 +      va_start(arg, fmt);
   1.155 +      xvprintf(fmt, arg);
   1.156 +      va_end(arg);
   1.157 +      return;
   1.158 +}
   1.159 +
   1.160 +static void read_char(struct csa *csa)
   1.161 +{     /* read next character */
   1.162 +      int c;
   1.163 +      if (csa->c == '\n')
   1.164 +         csa->recno++, csa->recpos = 0;
   1.165 +      csa->recpos++;
   1.166 +read: c = xfgetc(csa->fp);
   1.167 +      if (c < 0)
   1.168 +      {  if (xferror(csa->fp))
   1.169 +            error(csa, "read error - %s\n", xerrmsg());
   1.170 +         else if (csa->c == '\n')
   1.171 +            error(csa, "unexpected end of file\n");
   1.172 +         else
   1.173 +         {  warning(csa, "missing final end of line\n");
   1.174 +            c = '\n';
   1.175 +         }
   1.176 +      }
   1.177 +      else if (c == '\n')
   1.178 +         ;
   1.179 +      else if (csa->c == '\r')
   1.180 +      {  c = '\r';
   1.181 +         goto badc;
   1.182 +      }
   1.183 +      else if (csa->deck && c == '\r')
   1.184 +      {  csa->c = '\r';
   1.185 +         goto read;
   1.186 +      }
   1.187 +      else if (c == ' ')
   1.188 +         ;
   1.189 +      else if (isspace(c))
   1.190 +      {  if (csa->deck)
   1.191 +badc:       error(csa, "in fixed MPS format white-space character 0x%02"
   1.192 +               "X is not allowed\n", c);
   1.193 +         c = ' ';
   1.194 +      }
   1.195 +      else if (iscntrl(c))
   1.196 +         error(csa, "invalid control character 0x%02X\n", c);
   1.197 +      if (csa->deck && csa->recpos == 81 && c != '\n' && csa->w80 < 1)
   1.198 +      {  warning(csa, "in fixed MPS format record must not be longer th"
   1.199 +            "an 80 characters\n");
   1.200 +         csa->w80++;
   1.201 +      }
   1.202 +      csa->c = c;
   1.203 +      return;
   1.204 +}
   1.205 +
   1.206 +static int indicator(struct csa *csa, int name)
   1.207 +{     /* skip comment records and read possible indicator record */
   1.208 +      int ret;
   1.209 +      /* reset current field number */
   1.210 +      csa->fldno = 0;
   1.211 +loop: /* read the very first character of the next record */
   1.212 +      xassert(csa->c == '\n');
   1.213 +      read_char(csa);
   1.214 +      if (csa->c == ' ' || csa->c == '\n')
   1.215 +      {  /* data record */
   1.216 +         ret = 0;
   1.217 +      }
   1.218 +      else if (csa->c == '*')
   1.219 +      {  /* comment record */
   1.220 +         while (csa->c != '\n')
   1.221 +            read_char(csa);
   1.222 +         goto loop;
   1.223 +      }
   1.224 +      else
   1.225 +      {  /* indicator record */
   1.226 +         int len = 0;
   1.227 +         while (csa->c != ' ' && csa->c != '\n' && len < 12)
   1.228 +         {  csa->field[len++] = (char)csa->c;
   1.229 +            read_char(csa);
   1.230 +         }
   1.231 +         csa->field[len] = '\0';
   1.232 +         if (!(strcmp(csa->field, "NAME")    == 0 ||
   1.233 +               strcmp(csa->field, "ROWS")    == 0 ||
   1.234 +               strcmp(csa->field, "COLUMNS") == 0 ||
   1.235 +               strcmp(csa->field, "RHS")     == 0 ||
   1.236 +               strcmp(csa->field, "RANGES")  == 0 ||
   1.237 +               strcmp(csa->field, "BOUNDS")  == 0 ||
   1.238 +               strcmp(csa->field, "ENDATA")  == 0))
   1.239 +            error(csa, "invalid indicator record\n");
   1.240 +         if (!name)
   1.241 +         {  while (csa->c != '\n')
   1.242 +               read_char(csa);
   1.243 +         }
   1.244 +         ret = 1;
   1.245 +      }
   1.246 +      return ret;
   1.247 +}
   1.248 +
   1.249 +static void read_field(struct csa *csa)
   1.250 +{     /* read next field of the current data record */
   1.251 +      csa->fldno++;
   1.252 +      if (csa->deck)
   1.253 +      {  /* fixed MPS format */
   1.254 +         int beg, end, pos;
   1.255 +         /* determine predefined field positions */
   1.256 +         if (csa->fldno == 1)
   1.257 +            beg = 2, end = 3;
   1.258 +         else if (csa->fldno == 2)
   1.259 +            beg = 5, end = 12;
   1.260 +         else if (csa->fldno == 3)
   1.261 +            beg = 15, end = 22;
   1.262 +         else if (csa->fldno == 4)
   1.263 +            beg = 25, end = 36;
   1.264 +         else if (csa->fldno == 5)
   1.265 +            beg = 40, end = 47;
   1.266 +         else if (csa->fldno == 6)
   1.267 +            beg = 50, end = 61;
   1.268 +         else
   1.269 +            xassert(csa != csa);
   1.270 +         /* skip blanks preceding the current field */
   1.271 +         if (csa->c != '\n')
   1.272 +         {  pos = csa->recpos;
   1.273 +            while (csa->recpos < beg)
   1.274 +            {  if (csa->c == ' ')
   1.275 +                  ;
   1.276 +               else if (csa->c == '\n')
   1.277 +                  break;
   1.278 +               else
   1.279 +                  error(csa, "in fixed MPS format positions %d-%d must "
   1.280 +                     "be blank\n", pos, beg-1);
   1.281 +               read_char(csa);
   1.282 +            }
   1.283 +         }
   1.284 +         /* skip possible comment beginning in the field 3 or 5 */
   1.285 +         if ((csa->fldno == 3 || csa->fldno == 5) && csa->c == '$')
   1.286 +         {  while (csa->c != '\n')
   1.287 +               read_char(csa);
   1.288 +         }
   1.289 +         /* read the current field */
   1.290 +         for (pos = beg; pos <= end; pos++)
   1.291 +         {  if (csa->c == '\n') break;
   1.292 +            csa->field[pos-beg] = (char)csa->c;
   1.293 +            read_char(csa);
   1.294 +         }
   1.295 +         csa->field[pos-beg] = '\0';
   1.296 +         strtrim(csa->field);
   1.297 +         /* skip blanks following the last field */
   1.298 +         if (csa->fldno == 6 && csa->c != '\n')
   1.299 +         {  while (csa->recpos <= 72)
   1.300 +            {  if (csa->c == ' ')
   1.301 +                  ;
   1.302 +               else if (csa->c == '\n')
   1.303 +                  break;
   1.304 +               else
   1.305 +                  error(csa, "in fixed MPS format positions 62-72 must "
   1.306 +                     "be blank\n");
   1.307 +               read_char(csa);
   1.308 +            }
   1.309 +            while (csa->c != '\n')
   1.310 +               read_char(csa);
   1.311 +         }
   1.312 +      }
   1.313 +      else
   1.314 +      {  /* free MPS format */
   1.315 +         int len;
   1.316 +         /* skip blanks preceding the current field */
   1.317 +         while (csa->c == ' ')
   1.318 +            read_char(csa);
   1.319 +         /* skip possible comment */
   1.320 +         if (csa->c == '$')
   1.321 +         {  while (csa->c != '\n')
   1.322 +               read_char(csa);
   1.323 +         }
   1.324 +         /* read the current field */
   1.325 +         len = 0;
   1.326 +         while (!(csa->c == ' ' || csa->c == '\n'))
   1.327 +         {  if (len == 255)
   1.328 +               error(csa, "length of field %d exceeds 255 characters\n",
   1.329 +                  csa->fldno++);
   1.330 +            csa->field[len++] = (char)csa->c;
   1.331 +            read_char(csa);
   1.332 +         }
   1.333 +         csa->field[len] = '\0';
   1.334 +         /* skip anything following the last field (any extra fields
   1.335 +            are considered to be comments) */
   1.336 +         if (csa->fldno == 6)
   1.337 +         {  while (csa->c == ' ')
   1.338 +               read_char(csa);
   1.339 +            if (csa->c != '$' && csa->c != '\n' && csa->wef < 1)
   1.340 +            {  warning(csa, "some extra field(s) detected beyond field "
   1.341 +                  "6; field(s) ignored\n");
   1.342 +               csa->wef++;
   1.343 +            }
   1.344 +            while (csa->c != '\n')
   1.345 +               read_char(csa);
   1.346 +         }
   1.347 +      }
   1.348 +      return;
   1.349 +}
   1.350 +
   1.351 +static void patch_name(struct csa *csa, char *name)
   1.352 +{     /* process embedded blanks in symbolic name */
   1.353 +      int blank = csa->parm->blank;
   1.354 +      if (blank == '\0')
   1.355 +      {  /* remove emedded blanks */
   1.356 +         strspx(name);
   1.357 +      }
   1.358 +      else
   1.359 +      {  /* replace embedded blanks by specified character */
   1.360 +         for (; *name != '\0'; name++)
   1.361 +            if (*name == ' ') *name = (char)blank;
   1.362 +      }
   1.363 +      return;
   1.364 +}
   1.365 +
   1.366 +static double read_number(struct csa *csa)
   1.367 +{     /* read next field and convert it to floating-point number */
   1.368 +      double x;
   1.369 +      char *s;
   1.370 +      /* read next field */
   1.371 +      read_field(csa);
   1.372 +      xassert(csa->fldno == 4 || csa->fldno == 6);
   1.373 +      if (csa->field[0] == '\0')
   1.374 +         error(csa, "missing numeric value in field %d\n", csa->fldno);
   1.375 +      /* skip initial spaces of the field */
   1.376 +      for (s = csa->field; *s == ' '; s++);
   1.377 +      /* perform conversion */
   1.378 +      if (str2num(s, &x) != 0)
   1.379 +         error(csa, "cannot convert `%s' to floating-point number\n",
   1.380 +            s);
   1.381 +      return x;
   1.382 +}
   1.383 +
   1.384 +static void skip_field(struct csa *csa)
   1.385 +{     /* read and skip next field (assumed to be blank) */
   1.386 +      read_field(csa);
   1.387 +      if (csa->field[0] != '\0')
   1.388 +         error(csa, "field %d must be blank\n", csa->fldno);
   1.389 +      return;
   1.390 +}
   1.391 +
   1.392 +static void read_name(struct csa *csa)
   1.393 +{     /* read NAME indicator record */
   1.394 +      if (!(indicator(csa, 1) && strcmp(csa->field, "NAME") == 0))
   1.395 +         error(csa, "missing NAME indicator record\n");
   1.396 +      /* this indicator record looks like a data record; simulate that
   1.397 +         fields 1 and 2 were read */
   1.398 +      csa->fldno = 2;
   1.399 +      /* field 3: model name */
   1.400 +      read_field(csa), patch_name(csa, csa->field);
   1.401 +      if (csa->field[0] == '\0')
   1.402 +         warning(csa, "missing model name in field 3\n");
   1.403 +      else
   1.404 +         glp_set_prob_name(csa->P, csa->field);
   1.405 +      /* skip anything following field 3 */
   1.406 +      while (csa->c != '\n')
   1.407 +         read_char(csa);
   1.408 +      return;
   1.409 +}
   1.410 +
   1.411 +static void read_rows(struct csa *csa)
   1.412 +{     /* read ROWS section */
   1.413 +      int i, type;
   1.414 +loop: if (indicator(csa, 0)) goto done;
   1.415 +      /* field 1: row type */
   1.416 +      read_field(csa), strspx(csa->field);
   1.417 +      if (strcmp(csa->field, "N") == 0)
   1.418 +         type = GLP_FR;
   1.419 +      else if (strcmp(csa->field, "G") == 0)
   1.420 +         type = GLP_LO;
   1.421 +      else if (strcmp(csa->field, "L") == 0)
   1.422 +         type = GLP_UP;
   1.423 +      else if (strcmp(csa->field, "E") == 0)
   1.424 +         type = GLP_FX;
   1.425 +      else if (csa->field[0] == '\0')
   1.426 +         error(csa, "missing row type in field 1\n");
   1.427 +      else
   1.428 +         error(csa, "invalid row type in field 1\n");
   1.429 +      /* field 2: row name */
   1.430 +      read_field(csa), patch_name(csa, csa->field);
   1.431 +      if (csa->field[0] == '\0')
   1.432 +         error(csa, "missing row name in field 2\n");
   1.433 +      if (glp_find_row(csa->P, csa->field) != 0)
   1.434 +         error(csa, "row `%s' multiply specified\n", csa->field);
   1.435 +      i = glp_add_rows(csa->P, 1);
   1.436 +      glp_set_row_name(csa->P, i, csa->field);
   1.437 +      glp_set_row_bnds(csa->P, i, type, 0.0, 0.0);
   1.438 +      /* fields 3, 4, 5, and 6 must be blank */
   1.439 +      skip_field(csa);
   1.440 +      skip_field(csa);
   1.441 +      skip_field(csa);
   1.442 +      skip_field(csa);
   1.443 +      goto loop;
   1.444 +done: return;
   1.445 +}
   1.446 +
   1.447 +static void read_columns(struct csa *csa)
   1.448 +{     /* read COLUMNS section */
   1.449 +      int i, j, f, len, kind = GLP_CV, *ind;
   1.450 +      double aij, *val;
   1.451 +      char name[255+1], *flag;
   1.452 +      /* allocate working arrays */
   1.453 +      csa->work1 = ind = xcalloc(1+csa->P->m, sizeof(int));
   1.454 +      csa->work2 = val = xcalloc(1+csa->P->m, sizeof(double));
   1.455 +      csa->work3 = flag = xcalloc(1+csa->P->m, sizeof(char));
   1.456 +      memset(&flag[1], 0, csa->P->m);
   1.457 +      /* no current column exists */
   1.458 +      j = 0, len = 0;
   1.459 +loop: if (indicator(csa, 0)) goto done;
   1.460 +      /* field 1 must be blank */
   1.461 +      if (csa->deck)
   1.462 +      {  read_field(csa);
   1.463 +         if (csa->field[0] != '\0')
   1.464 +            error(csa, "field 1 must be blank\n");
   1.465 +      }
   1.466 +      else
   1.467 +         csa->fldno++;
   1.468 +      /* field 2: column or kind name */
   1.469 +      read_field(csa), patch_name(csa, csa->field);
   1.470 +      strcpy(name, csa->field);
   1.471 +      /* field 3: row name or keyword 'MARKER' */
   1.472 +      read_field(csa), patch_name(csa, csa->field);
   1.473 +      if (strcmp(csa->field, "'MARKER'") == 0)
   1.474 +      {  /* process kind data record */
   1.475 +         /* field 4 must be blank */
   1.476 +         if (csa->deck)
   1.477 +         {  read_field(csa);
   1.478 +            if (csa->field[0] != '\0')
   1.479 +               error(csa, "field 4 must be blank\n");
   1.480 +         }
   1.481 +         else
   1.482 +            csa->fldno++;
   1.483 +         /* field 5: keyword 'INTORG' or 'INTEND' */
   1.484 +         read_field(csa), patch_name(csa, csa->field);
   1.485 +         if (strcmp(csa->field, "'INTORG'") == 0)
   1.486 +            kind = GLP_IV;
   1.487 +         else if (strcmp(csa->field, "'INTEND'") == 0)
   1.488 +            kind = GLP_CV;
   1.489 +         else if (csa->field[0] == '\0')
   1.490 +            error(csa, "missing keyword in field 5\n");
   1.491 +         else
   1.492 +            error(csa, "invalid keyword in field 5\n");
   1.493 +         /* field 6 must be blank */
   1.494 +         skip_field(csa);
   1.495 +         goto loop;
   1.496 +      }
   1.497 +      /* process column name specified in field 2 */
   1.498 +      if (name[0] == '\0')
   1.499 +      {  /* the same column as in previous data record */
   1.500 +         if (j == 0)
   1.501 +            error(csa, "missing column name in field 2\n");
   1.502 +      }
   1.503 +      else if (j != 0 && strcmp(name, csa->P->col[j]->name) == 0)
   1.504 +      {  /* the same column as in previous data record */
   1.505 +         xassert(j != 0);
   1.506 +      }
   1.507 +      else
   1.508 +      {  /* store the current column */
   1.509 +         if (j != 0)
   1.510 +         {  glp_set_mat_col(csa->P, j, len, ind, val);
   1.511 +            while (len > 0) flag[ind[len--]] = 0;
   1.512 +         }
   1.513 +         /* create new column */
   1.514 +         if (glp_find_col(csa->P, name) != 0)
   1.515 +            error(csa, "column `%s' multiply specified\n", name);
   1.516 +         j = glp_add_cols(csa->P, 1);
   1.517 +         glp_set_col_name(csa->P, j, name);
   1.518 +         glp_set_col_kind(csa->P, j, kind);
   1.519 +         if (kind == GLP_CV)
   1.520 +            glp_set_col_bnds(csa->P, j, GLP_LO, 0.0, 0.0);
   1.521 +         else if (kind == GLP_IV)
   1.522 +            glp_set_col_bnds(csa->P, j, GLP_DB, 0.0, 1.0);
   1.523 +         else
   1.524 +            xassert(kind != kind);
   1.525 +      }
   1.526 +      /* process fields 3-4 and 5-6 */
   1.527 +      for (f = 3; f <= 5; f += 2)
   1.528 +      {  /* field 3 or 5: row name */
   1.529 +         if (f == 3)
   1.530 +         {  if (csa->field[0] == '\0')
   1.531 +               error(csa, "missing row name in field 3\n");
   1.532 +         }
   1.533 +         else
   1.534 +         {  read_field(csa), patch_name(csa, csa->field);
   1.535 +            if (csa->field[0] == '\0')
   1.536 +            {  /* if field 5 is blank, field 6 also must be blank */
   1.537 +               skip_field(csa);
   1.538 +               continue;
   1.539 +            }
   1.540 +         }
   1.541 +         i = glp_find_row(csa->P, csa->field);
   1.542 +         if (i == 0)
   1.543 +            error(csa, "row `%s' not found\n", csa->field);
   1.544 +         if (flag[i])
   1.545 +            error(csa, "duplicate coefficient in row `%s'\n",
   1.546 +               csa->field);
   1.547 +         /* field 4 or 6: coefficient value */
   1.548 +         aij = read_number(csa);
   1.549 +         if (fabs(aij) < csa->parm->tol_mps) aij = 0.0;
   1.550 +         len++, ind[len] = i, val[len] = aij, flag[i] = 1;
   1.551 +      }
   1.552 +      goto loop;
   1.553 +done: /* store the last column */
   1.554 +      if (j != 0)
   1.555 +         glp_set_mat_col(csa->P, j, len, ind, val);
   1.556 +      /* free working arrays */
   1.557 +      xfree(ind);
   1.558 +      xfree(val);
   1.559 +      xfree(flag);
   1.560 +      csa->work1 = csa->work2 = csa->work3 = NULL;
   1.561 +      return;
   1.562 +}
   1.563 +
   1.564 +static void read_rhs(struct csa *csa)
   1.565 +{     /* read RHS section */
   1.566 +      int i, f, v, type;
   1.567 +      double rhs;
   1.568 +      char name[255+1], *flag;
   1.569 +      /* allocate working array */
   1.570 +      csa->work3 = flag = xcalloc(1+csa->P->m, sizeof(char));
   1.571 +      memset(&flag[1], 0, csa->P->m);
   1.572 +      /* no current RHS vector exists */
   1.573 +      v = 0;
   1.574 +loop: if (indicator(csa, 0)) goto done;
   1.575 +      /* field 1 must be blank */
   1.576 +      if (csa->deck)
   1.577 +      {  read_field(csa);
   1.578 +         if (csa->field[0] != '\0')
   1.579 +            error(csa, "field 1 must be blank\n");
   1.580 +      }
   1.581 +      else
   1.582 +         csa->fldno++;
   1.583 +      /* field 2: RHS vector name */
   1.584 +      read_field(csa), patch_name(csa, csa->field);
   1.585 +      if (csa->field[0] == '\0')
   1.586 +      {  /* the same RHS vector as in previous data record */
   1.587 +         if (v == 0)
   1.588 +         {  warning(csa, "missing RHS vector name in field 2\n");
   1.589 +            goto blnk;
   1.590 +         }
   1.591 +      }
   1.592 +      else if (v != 0 && strcmp(csa->field, name) == 0)
   1.593 +      {  /* the same RHS vector as in previous data record */
   1.594 +         xassert(v != 0);
   1.595 +      }
   1.596 +      else
   1.597 +blnk: {  /* new RHS vector */
   1.598 +         if (v != 0)
   1.599 +            error(csa, "multiple RHS vectors not supported\n");
   1.600 +         v++;
   1.601 +         strcpy(name, csa->field);
   1.602 +      }
   1.603 +      /* process fields 3-4 and 5-6 */
   1.604 +      for (f = 3; f <= 5; f += 2)
   1.605 +      {  /* field 3 or 5: row name */
   1.606 +         read_field(csa), patch_name(csa, csa->field);
   1.607 +         if (csa->field[0] == '\0')
   1.608 +         {  if (f == 3)
   1.609 +               error(csa, "missing row name in field 3\n");
   1.610 +            else
   1.611 +            {  /* if field 5 is blank, field 6 also must be blank */
   1.612 +               skip_field(csa);
   1.613 +               continue;
   1.614 +            }
   1.615 +         }
   1.616 +         i = glp_find_row(csa->P, csa->field);
   1.617 +         if (i == 0)
   1.618 +            error(csa, "row `%s' not found\n", csa->field);
   1.619 +         if (flag[i])
   1.620 +            error(csa, "duplicate right-hand side for row `%s'\n",
   1.621 +               csa->field);
   1.622 +         /* field 4 or 6: right-hand side value */
   1.623 +         rhs = read_number(csa);
   1.624 +         if (fabs(rhs) < csa->parm->tol_mps) rhs = 0.0;
   1.625 +         type = csa->P->row[i]->type;
   1.626 +         if (type == GLP_FR)
   1.627 +         {  if (i == csa->obj_row)
   1.628 +               glp_set_obj_coef(csa->P, 0, rhs);
   1.629 +            else if (rhs != 0.0)
   1.630 +               warning(csa, "non-zero right-hand side for free row `%s'"
   1.631 +                  " ignored\n", csa->P->row[i]->name);
   1.632 +         }
   1.633 +         else
   1.634 +            glp_set_row_bnds(csa->P, i, type, rhs, rhs);
   1.635 +         flag[i] = 1;
   1.636 +      }
   1.637 +      goto loop;
   1.638 +done: /* free working array */
   1.639 +      xfree(flag);
   1.640 +      csa->work3 = NULL;
   1.641 +      return;
   1.642 +}
   1.643 +
   1.644 +static void read_ranges(struct csa *csa)
   1.645 +{     /* read RANGES section */
   1.646 +      int i, f, v, type;
   1.647 +      double rhs, rng;
   1.648 +      char name[255+1], *flag;
   1.649 +      /* allocate working array */
   1.650 +      csa->work3 = flag = xcalloc(1+csa->P->m, sizeof(char));
   1.651 +      memset(&flag[1], 0, csa->P->m);
   1.652 +      /* no current RANGES vector exists */
   1.653 +      v = 0;
   1.654 +loop: if (indicator(csa, 0)) goto done;
   1.655 +      /* field 1 must be blank */
   1.656 +      if (csa->deck)
   1.657 +      {  read_field(csa);
   1.658 +         if (csa->field[0] != '\0')
   1.659 +            error(csa, "field 1 must be blank\n");
   1.660 +      }
   1.661 +      else
   1.662 +         csa->fldno++;
   1.663 +      /* field 2: RANGES vector name */
   1.664 +      read_field(csa), patch_name(csa, csa->field);
   1.665 +      if (csa->field[0] == '\0')
   1.666 +      {  /* the same RANGES vector as in previous data record */
   1.667 +         if (v == 0)
   1.668 +         {  warning(csa, "missing RANGES vector name in field 2\n");
   1.669 +            goto blnk;
   1.670 +         }
   1.671 +      }
   1.672 +      else if (v != 0 && strcmp(csa->field, name) == 0)
   1.673 +      {  /* the same RANGES vector as in previous data record */
   1.674 +         xassert(v != 0);
   1.675 +      }
   1.676 +      else
   1.677 +blnk: {  /* new RANGES vector */
   1.678 +         if (v != 0)
   1.679 +            error(csa, "multiple RANGES vectors not supported\n");
   1.680 +         v++;
   1.681 +         strcpy(name, csa->field);
   1.682 +      }
   1.683 +      /* process fields 3-4 and 5-6 */
   1.684 +      for (f = 3; f <= 5; f += 2)
   1.685 +      {  /* field 3 or 5: row name */
   1.686 +         read_field(csa), patch_name(csa, csa->field);
   1.687 +         if (csa->field[0] == '\0')
   1.688 +         {  if (f == 3)
   1.689 +               error(csa, "missing row name in field 3\n");
   1.690 +            else
   1.691 +            {  /* if field 5 is blank, field 6 also must be blank */
   1.692 +               skip_field(csa);
   1.693 +               continue;
   1.694 +            }
   1.695 +         }
   1.696 +         i = glp_find_row(csa->P, csa->field);
   1.697 +         if (i == 0)
   1.698 +            error(csa, "row `%s' not found\n", csa->field);
   1.699 +         if (flag[i])
   1.700 +            error(csa, "duplicate range for row `%s'\n", csa->field);
   1.701 +         /* field 4 or 6: range value */
   1.702 +         rng = read_number(csa);
   1.703 +         if (fabs(rng) < csa->parm->tol_mps) rng = 0.0;
   1.704 +         type = csa->P->row[i]->type;
   1.705 +         if (type == GLP_FR)
   1.706 +            warning(csa, "range for free row `%s' ignored\n",
   1.707 +               csa->P->row[i]->name);
   1.708 +         else if (type == GLP_LO)
   1.709 +         {  rhs = csa->P->row[i]->lb;
   1.710 +            glp_set_row_bnds(csa->P, i, rhs == 0.0 ? GLP_FX : GLP_DB,
   1.711 +               rhs, rhs + fabs(rng));
   1.712 +         }
   1.713 +         else if (type == GLP_UP)
   1.714 +         {  rhs = csa->P->row[i]->ub;
   1.715 +            glp_set_row_bnds(csa->P, i, rhs == 0.0 ? GLP_FX : GLP_DB,
   1.716 +               rhs - fabs(rng), rhs);
   1.717 +         }
   1.718 +         else if (type == GLP_FX)
   1.719 +         {  rhs = csa->P->row[i]->lb;
   1.720 +            if (rng > 0.0)
   1.721 +               glp_set_row_bnds(csa->P, i, GLP_DB, rhs, rhs + rng);
   1.722 +            else if (rng < 0.0)
   1.723 +               glp_set_row_bnds(csa->P, i, GLP_DB, rhs + rng, rhs);
   1.724 +         }
   1.725 +         else
   1.726 +            xassert(type != type);
   1.727 +         flag[i] = 1;
   1.728 +      }
   1.729 +      goto loop;
   1.730 +done: /* free working array */
   1.731 +      xfree(flag);
   1.732 +      csa->work3 = NULL;
   1.733 +      return;
   1.734 +}
   1.735 +
   1.736 +static void read_bounds(struct csa *csa)
   1.737 +{     /* read BOUNDS section */
   1.738 +      GLPCOL *col;
   1.739 +      int j, v, mask, data;
   1.740 +      double bnd, lb, ub;
   1.741 +      char type[2+1], name[255+1], *flag;
   1.742 +      /* allocate working array */
   1.743 +      csa->work3 = flag = xcalloc(1+csa->P->n, sizeof(char));
   1.744 +      memset(&flag[1], 0, csa->P->n);
   1.745 +      /* no current BOUNDS vector exists */
   1.746 +      v = 0;
   1.747 +loop: if (indicator(csa, 0)) goto done;
   1.748 +      /* field 1: bound type */
   1.749 +      read_field(csa);
   1.750 +      if (strcmp(csa->field, "LO") == 0)
   1.751 +         mask = 0x01, data = 1;
   1.752 +      else if (strcmp(csa->field, "UP") == 0)
   1.753 +         mask = 0x10, data = 1;
   1.754 +      else if (strcmp(csa->field, "FX") == 0)
   1.755 +         mask = 0x11, data = 1;
   1.756 +      else if (strcmp(csa->field, "FR") == 0)
   1.757 +         mask = 0x11, data = 0;
   1.758 +      else if (strcmp(csa->field, "MI") == 0)
   1.759 +         mask = 0x01, data = 0;
   1.760 +      else if (strcmp(csa->field, "PL") == 0)
   1.761 +         mask = 0x10, data = 0;
   1.762 +      else if (strcmp(csa->field, "LI") == 0)
   1.763 +         mask = 0x01, data = 1;
   1.764 +      else if (strcmp(csa->field, "UI") == 0)
   1.765 +         mask = 0x10, data = 1;
   1.766 +      else if (strcmp(csa->field, "BV") == 0)
   1.767 +         mask = 0x11, data = 0;
   1.768 +      else if (csa->field[0] == '\0')
   1.769 +         error(csa, "missing bound type in field 1\n");
   1.770 +      else
   1.771 +         error(csa, "invalid bound type in field 1\n");
   1.772 +      strcpy(type, csa->field);
   1.773 +      /* field 2: BOUNDS vector name */
   1.774 +      read_field(csa), patch_name(csa, csa->field);
   1.775 +      if (csa->field[0] == '\0')
   1.776 +      {  /* the same BOUNDS vector as in previous data record */
   1.777 +         if (v == 0)
   1.778 +         {  warning(csa, "missing BOUNDS vector name in field 2\n");
   1.779 +            goto blnk;
   1.780 +         }
   1.781 +      }
   1.782 +      else if (v != 0 && strcmp(csa->field, name) == 0)
   1.783 +      {  /* the same BOUNDS vector as in previous data record */
   1.784 +         xassert(v != 0);
   1.785 +      }
   1.786 +      else
   1.787 +blnk: {  /* new BOUNDS vector */
   1.788 +         if (v != 0)
   1.789 +            error(csa, "multiple BOUNDS vectors not supported\n");
   1.790 +         v++;
   1.791 +         strcpy(name, csa->field);
   1.792 +      }
   1.793 +      /* field 3: column name */
   1.794 +      read_field(csa), patch_name(csa, csa->field);
   1.795 +      if (csa->field[0] == '\0')
   1.796 +         error(csa, "missing column name in field 3\n");
   1.797 +      j = glp_find_col(csa->P, csa->field);
   1.798 +      if (j == 0)
   1.799 +         error(csa, "column `%s' not found\n", csa->field);
   1.800 +      if ((flag[j] & mask) == 0x01)
   1.801 +         error(csa, "duplicate lower bound for column `%s'\n",
   1.802 +            csa->field);
   1.803 +      if ((flag[j] & mask) == 0x10)
   1.804 +         error(csa, "duplicate upper bound for column `%s'\n",
   1.805 +            csa->field);
   1.806 +      xassert((flag[j] & mask) == 0x00);
   1.807 +      /* field 4: bound value */
   1.808 +      if (data)
   1.809 +      {  bnd = read_number(csa);
   1.810 +         if (fabs(bnd) < csa->parm->tol_mps) bnd = 0.0;
   1.811 +      }
   1.812 +      else
   1.813 +         read_field(csa), bnd = 0.0;
   1.814 +      /* get current column bounds */
   1.815 +      col = csa->P->col[j];
   1.816 +      if (col->type == GLP_FR)
   1.817 +         lb = -DBL_MAX, ub = +DBL_MAX;
   1.818 +      else if (col->type == GLP_LO)
   1.819 +         lb = col->lb, ub = +DBL_MAX;
   1.820 +      else if (col->type == GLP_UP)
   1.821 +         lb = -DBL_MAX, ub = col->ub;
   1.822 +      else if (col->type == GLP_DB)
   1.823 +         lb = col->lb, ub = col->ub;
   1.824 +      else if (col->type == GLP_FX)
   1.825 +         lb = ub = col->lb;
   1.826 +      else
   1.827 +         xassert(col != col);
   1.828 +      /* change column bounds */
   1.829 +      if (strcmp(type, "LO") == 0)
   1.830 +         lb = bnd;
   1.831 +      else if (strcmp(type, "UP") == 0)
   1.832 +         ub = bnd;
   1.833 +      else if (strcmp(type, "FX") == 0)
   1.834 +         lb = ub = bnd;
   1.835 +      else if (strcmp(type, "FR") == 0)
   1.836 +         lb = -DBL_MAX, ub = +DBL_MAX;
   1.837 +      else if (strcmp(type, "MI") == 0)
   1.838 +         lb = -DBL_MAX;
   1.839 +      else if (strcmp(type, "PL") == 0)
   1.840 +         ub = +DBL_MAX;
   1.841 +      else if (strcmp(type, "LI") == 0)
   1.842 +      {  glp_set_col_kind(csa->P, j, GLP_IV);
   1.843 +         lb = ceil(bnd);
   1.844 +      }
   1.845 +      else if (strcmp(type, "UI") == 0)
   1.846 +      {  glp_set_col_kind(csa->P, j, GLP_IV);
   1.847 +         ub = floor(bnd);
   1.848 +      }
   1.849 +      else if (strcmp(type, "BV") == 0)
   1.850 +      {  glp_set_col_kind(csa->P, j, GLP_IV);
   1.851 +         lb = 0.0, ub = 1.0;
   1.852 +      }
   1.853 +      else
   1.854 +         xassert(type != type);
   1.855 +      /* set new column bounds */
   1.856 +      if (lb == -DBL_MAX && ub == +DBL_MAX)
   1.857 +         glp_set_col_bnds(csa->P, j, GLP_FR, lb, ub);
   1.858 +      else if (ub == +DBL_MAX)
   1.859 +         glp_set_col_bnds(csa->P, j, GLP_LO, lb, ub);
   1.860 +      else if (lb == -DBL_MAX)
   1.861 +         glp_set_col_bnds(csa->P, j, GLP_UP, lb, ub);
   1.862 +      else if (lb != ub)
   1.863 +         glp_set_col_bnds(csa->P, j, GLP_DB, lb, ub);
   1.864 +      else
   1.865 +         glp_set_col_bnds(csa->P, j, GLP_FX, lb, ub);
   1.866 +      flag[j] |= (char)mask;
   1.867 +      /* fields 5 and 6 must be blank */
   1.868 +      skip_field(csa);
   1.869 +      skip_field(csa);
   1.870 +      goto loop;
   1.871 +done: /* free working array */
   1.872 +      xfree(flag);
   1.873 +      csa->work3 = NULL;
   1.874 +      return;
   1.875 +}
   1.876 +
   1.877 +int glp_read_mps(glp_prob *P, int fmt, const glp_mpscp *parm,
   1.878 +      const char *fname)
   1.879 +{     /* read problem data in MPS format */
   1.880 +      glp_mpscp _parm;
   1.881 +      struct csa _csa, *csa = &_csa;
   1.882 +      int ret;
   1.883 +      xprintf("Reading problem data from `%s'...\n", fname);
   1.884 +      if (!(fmt == GLP_MPS_DECK || fmt == GLP_MPS_FILE))
   1.885 +         xerror("glp_read_mps: fmt = %d; invalid parameter\n", fmt);
   1.886 +      if (parm == NULL)
   1.887 +         glp_init_mpscp(&_parm), parm = &_parm;
   1.888 +      /* check control parameters */
   1.889 +      check_parm("glp_read_mps", parm);
   1.890 +      /* initialize common storage area */
   1.891 +      csa->P = P;
   1.892 +      csa->deck = (fmt == GLP_MPS_DECK);
   1.893 +      csa->parm = parm;
   1.894 +      csa->fname = fname;
   1.895 +      csa->fp = NULL;
   1.896 +      if (setjmp(csa->jump))
   1.897 +      {  ret = 1;
   1.898 +         goto done;
   1.899 +      }
   1.900 +      csa->recno = csa->recpos = 0;
   1.901 +      csa->c = '\n';
   1.902 +      csa->fldno = 0;
   1.903 +      csa->field[0] = '\0';
   1.904 +      csa->w80 = csa->wef = 0;
   1.905 +      csa->obj_row = 0;
   1.906 +      csa->work1 = csa->work2 = csa->work3 = NULL;
   1.907 +      /* erase problem object */
   1.908 +      glp_erase_prob(P);
   1.909 +      glp_create_index(P);
   1.910 +      /* open input MPS file */
   1.911 +      csa->fp = xfopen(fname, "r");
   1.912 +      if (csa->fp == NULL)
   1.913 +      {  xprintf("Unable to open `%s' - %s\n", fname, xerrmsg());
   1.914 +         ret = 1;
   1.915 +         goto done;
   1.916 +      }
   1.917 +      /* read NAME indicator record */
   1.918 +      read_name(csa);
   1.919 +      if (P->name != NULL)
   1.920 +         xprintf("Problem: %s\n", P->name);
   1.921 +      /* read ROWS section */
   1.922 +      if (!(indicator(csa, 0) && strcmp(csa->field, "ROWS") == 0))
   1.923 +         error(csa, "missing ROWS indicator record\n");
   1.924 +      read_rows(csa);
   1.925 +      /* determine objective row */
   1.926 +      if (parm->obj_name == NULL || parm->obj_name[0] == '\0')
   1.927 +      {  /* use the first row of N type */
   1.928 +         int i;
   1.929 +         for (i = 1; i <= P->m; i++)
   1.930 +         {  if (P->row[i]->type == GLP_FR)
   1.931 +            {  csa->obj_row = i;
   1.932 +               break;
   1.933 +            }
   1.934 +         }
   1.935 +         if (csa->obj_row == 0)
   1.936 +            warning(csa, "unable to determine objective row\n");
   1.937 +      }
   1.938 +      else
   1.939 +      {  /* use a row with specified name */
   1.940 +         int i;
   1.941 +         for (i = 1; i <= P->m; i++)
   1.942 +         {  xassert(P->row[i]->name != NULL);
   1.943 +            if (strcmp(parm->obj_name, P->row[i]->name) == 0)
   1.944 +            {  csa->obj_row = i;
   1.945 +               break;
   1.946 +            }
   1.947 +         }
   1.948 +         if (csa->obj_row == 0)
   1.949 +            error(csa, "objective row `%s' not found\n",
   1.950 +               parm->obj_name);
   1.951 +      }
   1.952 +      if (csa->obj_row != 0)
   1.953 +      {  glp_set_obj_name(P, P->row[csa->obj_row]->name);
   1.954 +         xprintf("Objective: %s\n", P->obj);
   1.955 +      }
   1.956 +      /* read COLUMNS section */
   1.957 +      if (strcmp(csa->field, "COLUMNS") != 0)
   1.958 +         error(csa, "missing COLUMNS indicator record\n");
   1.959 +      read_columns(csa);
   1.960 +      /* set objective coefficients */
   1.961 +      if (csa->obj_row != 0)
   1.962 +      {  GLPAIJ *aij;
   1.963 +         for (aij = P->row[csa->obj_row]->ptr; aij != NULL; aij =
   1.964 +            aij->r_next) glp_set_obj_coef(P, aij->col->j, aij->val);
   1.965 +      }
   1.966 +      /* read optional RHS section */
   1.967 +      if (strcmp(csa->field, "RHS") == 0)
   1.968 +         read_rhs(csa);
   1.969 +      /* read optional RANGES section */
   1.970 +      if (strcmp(csa->field, "RANGES") == 0)
   1.971 +         read_ranges(csa);
   1.972 +      /* read optional BOUNDS section */
   1.973 +      if (strcmp(csa->field, "BOUNDS") == 0)
   1.974 +         read_bounds(csa);
   1.975 +      /* read ENDATA indicator record */
   1.976 +      if (strcmp(csa->field, "ENDATA") != 0)
   1.977 +         error(csa, "invalid use of %s indicator record\n",
   1.978 +            csa->field);
   1.979 +      /* print some statistics */
   1.980 +      xprintf("%d row%s, %d column%s, %d non-zero%s\n",
   1.981 +         P->m, P->m == 1 ? "" : "s", P->n, P->n == 1 ? "" : "s",
   1.982 +         P->nnz, P->nnz == 1 ? "" : "s");
   1.983 +      if (glp_get_num_int(P) > 0)
   1.984 +      {  int ni = glp_get_num_int(P);
   1.985 +         int nb = glp_get_num_bin(P);
   1.986 +         if (ni == 1)
   1.987 +         {  if (nb == 0)
   1.988 +               xprintf("One variable is integer\n");
   1.989 +            else
   1.990 +               xprintf("One variable is binary\n");
   1.991 +         }
   1.992 +         else
   1.993 +         {  xprintf("%d integer variables, ", ni);
   1.994 +            if (nb == 0)
   1.995 +               xprintf("none");
   1.996 +            else if (nb == 1)
   1.997 +               xprintf("one");
   1.998 +            else if (nb == ni)
   1.999 +               xprintf("all");
  1.1000 +            else
  1.1001 +               xprintf("%d", nb);
  1.1002 +            xprintf(" of which %s binary\n", nb == 1 ? "is" : "are");
  1.1003 +         }
  1.1004 +      }
  1.1005 +      xprintf("%d records were read\n", csa->recno);
  1.1006 +      /* problem data has been successfully read */
  1.1007 +      glp_delete_index(P);
  1.1008 +      glp_sort_matrix(P);
  1.1009 +      ret = 0;
  1.1010 +done: if (csa->fp != NULL) xfclose(csa->fp);
  1.1011 +      if (csa->work1 != NULL) xfree(csa->work1);
  1.1012 +      if (csa->work2 != NULL) xfree(csa->work2);
  1.1013 +      if (csa->work3 != NULL) xfree(csa->work3);
  1.1014 +      if (ret != 0) glp_erase_prob(P);
  1.1015 +      return ret;
  1.1016 +}
  1.1017 +
  1.1018 +/***********************************************************************
  1.1019 +*  NAME
  1.1020 +*
  1.1021 +*  glp_write_mps - write problem data in MPS format
  1.1022 +*
  1.1023 +*  SYNOPSIS
  1.1024 +*
  1.1025 +*  int glp_write_mps(glp_prob *P, int fmt, const glp_mpscp *parm,
  1.1026 +*     const char *fname);
  1.1027 +*
  1.1028 +*  DESCRIPTION
  1.1029 +*
  1.1030 +*  The routine glp_write_mps writes problem data in MPS format to a
  1.1031 +*  text file.
  1.1032 +*
  1.1033 +*  The parameter fmt specifies the version of MPS format:
  1.1034 +*
  1.1035 +*  GLP_MPS_DECK - fixed (ancient) MPS format;
  1.1036 +*  GLP_MPS_FILE - free (modern) MPS format.
  1.1037 +*
  1.1038 +*  The parameter parm is a pointer to the structure glp_mpscp, which
  1.1039 +*  specifies control parameters used by the routine. If parm is NULL,
  1.1040 +*  the routine uses default settings.
  1.1041 +*
  1.1042 +*  The character string fname specifies a name of the text file to be
  1.1043 +*  written.
  1.1044 +*
  1.1045 +*  RETURNS
  1.1046 +*
  1.1047 +*  If the operation was successful, the routine glp_read_mps returns
  1.1048 +*  zero. Otherwise, it prints an error message and returns non-zero. */
  1.1049 +
  1.1050 +#define csa csa1
  1.1051 +
  1.1052 +struct csa
  1.1053 +{     /* common storage area */
  1.1054 +      glp_prob *P;
  1.1055 +      /* pointer to problem object */
  1.1056 +      int deck;
  1.1057 +      /* MPS format (0 - free, 1 - fixed) */
  1.1058 +      const glp_mpscp *parm;
  1.1059 +      /* pointer to control parameters */
  1.1060 +      char field[255+1];
  1.1061 +      /* field buffer */
  1.1062 +};
  1.1063 +
  1.1064 +static char *mps_name(struct csa *csa)
  1.1065 +{     /* make problem name */
  1.1066 +      char *f;
  1.1067 +      if (csa->P->name == NULL)
  1.1068 +         csa->field[0] = '\0';
  1.1069 +      else if (csa->deck)
  1.1070 +      {  strncpy(csa->field, csa->P->name, 8);
  1.1071 +         csa->field[8] = '\0';
  1.1072 +      }
  1.1073 +      else
  1.1074 +         strcpy(csa->field, csa->P->name);
  1.1075 +      for (f = csa->field; *f != '\0'; f++)
  1.1076 +         if (*f == ' ') *f = '_';
  1.1077 +      return csa->field;
  1.1078 +}
  1.1079 +
  1.1080 +static char *row_name(struct csa *csa, int i)
  1.1081 +{     /* make i-th row name */
  1.1082 +      char *f;
  1.1083 +      xassert(0 <= i && i <= csa->P->m);
  1.1084 +      if (i == 0 || csa->P->row[i]->name == NULL ||
  1.1085 +          csa->deck && strlen(csa->P->row[i]->name) > 8)
  1.1086 +         sprintf(csa->field, "R%07d", i);
  1.1087 +      else
  1.1088 +      {  strcpy(csa->field, csa->P->row[i]->name);
  1.1089 +         for (f = csa->field; *f != '\0'; f++)
  1.1090 +            if (*f == ' ') *f = '_';
  1.1091 +      }
  1.1092 +      return csa->field;
  1.1093 +}
  1.1094 +
  1.1095 +static char *col_name(struct csa *csa, int j)
  1.1096 +{     /* make j-th column name */
  1.1097 +      char *f;
  1.1098 +      xassert(1 <= j && j <= csa->P->n);
  1.1099 +      if (csa->P->col[j]->name == NULL ||
  1.1100 +          csa->deck && strlen(csa->P->col[j]->name) > 8)
  1.1101 +         sprintf(csa->field, "C%07d", j);
  1.1102 +      else
  1.1103 +      {  strcpy(csa->field, csa->P->col[j]->name);
  1.1104 +         for (f = csa->field; *f != '\0'; f++)
  1.1105 +            if (*f == ' ') *f = '_';
  1.1106 +      }
  1.1107 +      return csa->field;
  1.1108 +}
  1.1109 +
  1.1110 +static char *mps_numb(struct csa *csa, double val)
  1.1111 +{     /* format floating-point number */
  1.1112 +      int dig;
  1.1113 +      char *exp;
  1.1114 +      for (dig = 12; dig >= 6; dig--)
  1.1115 +      {  if (val != 0.0 && fabs(val) < 0.002)
  1.1116 +            sprintf(csa->field, "%.*E", dig-1, val);
  1.1117 +         else
  1.1118 +            sprintf(csa->field, "%.*G", dig, val);
  1.1119 +         exp = strchr(csa->field, 'E');
  1.1120 +         if (exp != NULL)
  1.1121 +            sprintf(exp+1, "%d", atoi(exp+1));
  1.1122 +         if (strlen(csa->field) <= 12) break;
  1.1123 +      }
  1.1124 +      xassert(strlen(csa->field) <= 12);
  1.1125 +      return csa->field;
  1.1126 +}
  1.1127 +
  1.1128 +int glp_write_mps(glp_prob *P, int fmt, const glp_mpscp *parm,
  1.1129 +      const char *fname)
  1.1130 +{     /* write problem data in MPS format */
  1.1131 +      glp_mpscp _parm;
  1.1132 +      struct csa _csa, *csa = &_csa;
  1.1133 +      XFILE *fp;
  1.1134 +      int out_obj, one_col = 0, empty = 0;
  1.1135 +      int i, j, recno, marker, count, gap, ret;
  1.1136 +      xprintf("Writing problem data to `%s'...\n", fname);
  1.1137 +      if (!(fmt == GLP_MPS_DECK || fmt == GLP_MPS_FILE))
  1.1138 +         xerror("glp_write_mps: fmt = %d; invalid parameter\n", fmt);
  1.1139 +      if (parm == NULL)
  1.1140 +         glp_init_mpscp(&_parm), parm = &_parm;
  1.1141 +      /* check control parameters */
  1.1142 +      check_parm("glp_write_mps", parm);
  1.1143 +      /* initialize common storage area */
  1.1144 +      csa->P = P;
  1.1145 +      csa->deck = (fmt == GLP_MPS_DECK);
  1.1146 +      csa->parm = parm;
  1.1147 +      /* create output MPS file */
  1.1148 +      fp = xfopen(fname, "w"), recno = 0;
  1.1149 +      if (fp == NULL)
  1.1150 +      {  xprintf("Unable to create `%s' - %s\n", fname, xerrmsg());
  1.1151 +         ret = 1;
  1.1152 +         goto done;
  1.1153 +      }
  1.1154 +      /* write comment records */
  1.1155 +      xfprintf(fp, "* %-*s%s\n", P->name == NULL ? 1 : 12, "Problem:",
  1.1156 +         P->name == NULL ? "" : P->name), recno++;
  1.1157 +      xfprintf(fp, "* %-12s%s\n", "Class:", glp_get_num_int(P) == 0 ?
  1.1158 +         "LP" : "MIP"), recno++;
  1.1159 +      xfprintf(fp, "* %-12s%d\n", "Rows:", P->m), recno++;
  1.1160 +      if (glp_get_num_int(P) == 0)
  1.1161 +         xfprintf(fp, "* %-12s%d\n", "Columns:", P->n), recno++;
  1.1162 +      else
  1.1163 +         xfprintf(fp, "* %-12s%d (%d integer, %d binary)\n",
  1.1164 +            "Columns:", P->n, glp_get_num_int(P), glp_get_num_bin(P)),
  1.1165 +            recno++;
  1.1166 +      xfprintf(fp, "* %-12s%d\n", "Non-zeros:", P->nnz), recno++;
  1.1167 +      xfprintf(fp, "* %-12s%s\n", "Format:", csa->deck ? "Fixed MPS" :
  1.1168 +         "Free MPS"), recno++;
  1.1169 +      xfprintf(fp, "*\n", recno++);
  1.1170 +      /* write NAME indicator record */
  1.1171 +      xfprintf(fp, "NAME%*s%s\n",
  1.1172 +         P->name == NULL ? 0 : csa->deck ? 10 : 1, "", mps_name(csa)),
  1.1173 +         recno++;
  1.1174 +#if 1
  1.1175 +      /* determine whether to write the objective row */
  1.1176 +      out_obj = 1;
  1.1177 +      for (i = 1; i <= P->m; i++)
  1.1178 +      {  if (P->row[i]->type == GLP_FR)
  1.1179 +         {  out_obj = 0;
  1.1180 +            break;
  1.1181 +         }
  1.1182 +      }
  1.1183 +#endif
  1.1184 +      /* write ROWS section */
  1.1185 +      xfprintf(fp, "ROWS\n"), recno++;
  1.1186 +      for (i = (out_obj ? 0 : 1); i <= P->m; i++)
  1.1187 +      {  int type;
  1.1188 +         type = (i == 0 ? GLP_FR : P->row[i]->type);
  1.1189 +         if (type == GLP_FR)
  1.1190 +            type = 'N';
  1.1191 +         else if (type == GLP_LO)
  1.1192 +            type = 'G';
  1.1193 +         else if (type == GLP_UP)
  1.1194 +            type = 'L';
  1.1195 +         else if (type == GLP_DB || type == GLP_FX)
  1.1196 +            type = 'E';
  1.1197 +         else
  1.1198 +            xassert(type != type);
  1.1199 +         xfprintf(fp, " %c%*s%s\n", type, csa->deck ? 2 : 1, "",
  1.1200 +            row_name(csa, i)), recno++;
  1.1201 +      }
  1.1202 +      /* write COLUMNS section */
  1.1203 +      xfprintf(fp, "COLUMNS\n"), recno++;
  1.1204 +      marker = 0;
  1.1205 +      for (j = 1; j <= P->n; j++)
  1.1206 +      {  GLPAIJ cj, *aij;
  1.1207 +         int kind;
  1.1208 +         kind = P->col[j]->kind;
  1.1209 +         if (kind == GLP_CV)
  1.1210 +         {  if (marker % 2 == 1)
  1.1211 +            {  /* close current integer block */
  1.1212 +               marker++;
  1.1213 +               xfprintf(fp, "%*sM%07d%*s'MARKER'%*s'INTEND'\n",
  1.1214 +                  csa->deck ? 4 : 1, "", marker, csa->deck ? 2 : 1, "",
  1.1215 +                  csa->deck ? 17 : 1, ""), recno++;
  1.1216 +            }
  1.1217 +         }
  1.1218 +         else if (kind == GLP_IV)
  1.1219 +         {  if (marker % 2 == 0)
  1.1220 +            {  /* open new integer block */
  1.1221 +               marker++;
  1.1222 +               xfprintf(fp, "%*sM%07d%*s'MARKER'%*s'INTORG'\n",
  1.1223 +                  csa->deck ? 4 : 1, "", marker, csa->deck ? 2 : 1, "",
  1.1224 +                  csa->deck ? 17 : 1, ""), recno++;
  1.1225 +            }
  1.1226 +         }
  1.1227 +         else
  1.1228 +            xassert(kind != kind);
  1.1229 +         if (out_obj && P->col[j]->coef != 0.0)
  1.1230 +         {  /* make fake objective coefficient */
  1.1231 +            aij = &cj;
  1.1232 +            aij->row = NULL;
  1.1233 +            aij->val = P->col[j]->coef;
  1.1234 +            aij->c_next = P->col[j]->ptr;
  1.1235 +         }
  1.1236 +         else
  1.1237 +            aij = P->col[j]->ptr;
  1.1238 +#if 1 /* FIXME */
  1.1239 +         if (aij == NULL)
  1.1240 +         {  /* empty column */
  1.1241 +            empty++;
  1.1242 +            xfprintf(fp, "%*s%-*s", csa->deck ? 4 : 1, "",
  1.1243 +               csa->deck ? 8 : 1, col_name(csa, j));
  1.1244 +            /* we need a row */
  1.1245 +            xassert(P->m > 0);
  1.1246 +            xfprintf(fp, "%*s%-*s",
  1.1247 +               csa->deck ? 2 : 1, "", csa->deck ? 8 : 1,
  1.1248 +               row_name(csa, 1));
  1.1249 +            xfprintf(fp, "%*s0%*s$ empty column\n",
  1.1250 +               csa->deck ? 13 : 1, "", csa->deck ? 3 : 1, ""), recno++;
  1.1251 +         }
  1.1252 +#endif
  1.1253 +         count = 0;
  1.1254 +         for (aij = aij; aij != NULL; aij = aij->c_next)
  1.1255 +         {  if (one_col || count % 2 == 0)
  1.1256 +               xfprintf(fp, "%*s%-*s", csa->deck ? 4 : 1, "",
  1.1257 +                  csa->deck ? 8 : 1, col_name(csa, j));
  1.1258 +            gap = (one_col || count % 2 == 0 ? 2 : 3);
  1.1259 +            xfprintf(fp, "%*s%-*s",
  1.1260 +               csa->deck ? gap : 1, "", csa->deck ? 8 : 1,
  1.1261 +               row_name(csa, aij->row == NULL ? 0 : aij->row->i));
  1.1262 +            xfprintf(fp, "%*s%*s",
  1.1263 +               csa->deck ? 2 : 1, "", csa->deck ? 12 : 1,
  1.1264 +               mps_numb(csa, aij->val)), count++;
  1.1265 +            if (one_col || count % 2 == 0)
  1.1266 +               xfprintf(fp, "\n"), recno++;
  1.1267 +         }
  1.1268 +         if (!(one_col || count % 2 == 0))
  1.1269 +            xfprintf(fp, "\n"), recno++;
  1.1270 +      }
  1.1271 +      if (marker % 2 == 1)
  1.1272 +      {  /* close last integer block */
  1.1273 +         marker++;
  1.1274 +         xfprintf(fp, "%*sM%07d%*s'MARKER'%*s'INTEND'\n",
  1.1275 +            csa->deck ? 4 : 1, "", marker, csa->deck ? 2 : 1, "",
  1.1276 +            csa->deck ? 17 : 1, ""), recno++;
  1.1277 +      }
  1.1278 +#if 1
  1.1279 +      if (empty > 0)
  1.1280 +         xprintf("Warning: problem has %d empty column(s)\n", empty);
  1.1281 +#endif
  1.1282 +      /* write RHS section */
  1.1283 +      xfprintf(fp, "RHS\n"), recno++;
  1.1284 +      count = 0;
  1.1285 +      for (i = (out_obj ? 0 : 1); i <= P->m; i++)
  1.1286 +      {  int type;
  1.1287 +         double rhs;
  1.1288 +         if (i == 0)
  1.1289 +            rhs = P->c0;
  1.1290 +         else
  1.1291 +         {  type = P->row[i]->type;
  1.1292 +            if (type == GLP_FR)
  1.1293 +               rhs = 0.0;
  1.1294 +            else if (type == GLP_LO)
  1.1295 +               rhs = P->row[i]->lb;
  1.1296 +            else if (type == GLP_UP)
  1.1297 +               rhs = P->row[i]->ub;
  1.1298 +            else if (type == GLP_DB || type == GLP_FX)
  1.1299 +               rhs = P->row[i]->lb;
  1.1300 +            else
  1.1301 +               xassert(type != type);
  1.1302 +         }
  1.1303 +         if (rhs != 0.0)
  1.1304 +         {  if (one_col || count % 2 == 0)
  1.1305 +               xfprintf(fp, "%*s%-*s", csa->deck ? 4 : 1, "",
  1.1306 +                  csa->deck ? 8 : 1, "RHS1");
  1.1307 +            gap = (one_col || count % 2 == 0 ? 2 : 3);
  1.1308 +            xfprintf(fp, "%*s%-*s",
  1.1309 +               csa->deck ? gap : 1, "", csa->deck ? 8 : 1,
  1.1310 +               row_name(csa, i));
  1.1311 +            xfprintf(fp, "%*s%*s",
  1.1312 +               csa->deck ? 2 : 1, "", csa->deck ? 12 : 1,
  1.1313 +               mps_numb(csa, rhs)), count++;
  1.1314 +            if (one_col || count % 2 == 0)
  1.1315 +               xfprintf(fp, "\n"), recno++;
  1.1316 +         }
  1.1317 +      }
  1.1318 +      if (!(one_col || count % 2 == 0))
  1.1319 +         xfprintf(fp, "\n"), recno++;
  1.1320 +      /* write RANGES section */
  1.1321 +      for (i = P->m; i >= 1; i--)
  1.1322 +         if (P->row[i]->type == GLP_DB) break;
  1.1323 +      if (i == 0) goto bnds;
  1.1324 +      xfprintf(fp, "RANGES\n"), recno++;
  1.1325 +      count = 0;
  1.1326 +      for (i = 1; i <= P->m; i++)
  1.1327 +      {  if (P->row[i]->type == GLP_DB)
  1.1328 +         {  if (one_col || count % 2 == 0)
  1.1329 +               xfprintf(fp, "%*s%-*s", csa->deck ? 4 : 1, "",
  1.1330 +                  csa->deck ? 8 : 1, "RNG1");
  1.1331 +            gap = (one_col || count % 2 == 0 ? 2 : 3);
  1.1332 +            xfprintf(fp, "%*s%-*s",
  1.1333 +               csa->deck ? gap : 1, "", csa->deck ? 8 : 1,
  1.1334 +               row_name(csa, i));
  1.1335 +            xfprintf(fp, "%*s%*s",
  1.1336 +               csa->deck ? 2 : 1, "", csa->deck ? 12 : 1,
  1.1337 +               mps_numb(csa, P->row[i]->ub - P->row[i]->lb)), count++;
  1.1338 +            if (one_col || count % 2 == 0)
  1.1339 +               xfprintf(fp, "\n"), recno++;
  1.1340 +         }
  1.1341 +      }
  1.1342 +      if (!(one_col || count % 2 == 0))
  1.1343 +         xfprintf(fp, "\n"), recno++;
  1.1344 +bnds: /* write BOUNDS section */
  1.1345 +      for (j = P->n; j >= 1; j--)
  1.1346 +         if (!(P->col[j]->type == GLP_LO && P->col[j]->lb == 0.0))
  1.1347 +            break;
  1.1348 +      if (j == 0) goto endt;
  1.1349 +      xfprintf(fp, "BOUNDS\n"), recno++;
  1.1350 +      for (j = 1; j <= P->n; j++)
  1.1351 +      {  int type, data[2];
  1.1352 +         double bnd[2];
  1.1353 +         char *spec[2];
  1.1354 +         spec[0] = spec[1] = NULL;
  1.1355 +         type = P->col[j]->type;
  1.1356 +         if (type == GLP_FR)
  1.1357 +            spec[0] = "FR", data[0] = 0;
  1.1358 +         else if (type == GLP_LO)
  1.1359 +         {  if (P->col[j]->lb != 0.0)
  1.1360 +               spec[0] = "LO", data[0] = 1, bnd[0] = P->col[j]->lb;
  1.1361 +            if (P->col[j]->kind == GLP_IV)
  1.1362 +               spec[1] = "PL", data[1] = 0;
  1.1363 +         }
  1.1364 +         else if (type == GLP_UP)
  1.1365 +         {  spec[0] = "MI", data[0] = 0;
  1.1366 +            spec[1] = "UP", data[1] = 1, bnd[1] = P->col[j]->ub;
  1.1367 +         }
  1.1368 +         else if (type == GLP_DB)
  1.1369 +         {  if (P->col[j]->lb != 0.0)
  1.1370 +               spec[0] = "LO", data[0] = 1, bnd[0] = P->col[j]->lb;
  1.1371 +            spec[1] = "UP", data[1] = 1, bnd[1] = P->col[j]->ub;
  1.1372 +         }
  1.1373 +         else if (type == GLP_FX)
  1.1374 +            spec[0] = "FX", data[0] = 1, bnd[0] = P->col[j]->lb;
  1.1375 +         else
  1.1376 +            xassert(type != type);
  1.1377 +         for (i = 0; i <= 1; i++)
  1.1378 +         {  if (spec[i] != NULL)
  1.1379 +            {  xfprintf(fp, " %s %-*s%*s%-*s", spec[i],
  1.1380 +                  csa->deck ? 8 : 1, "BND1", csa->deck ? 2 : 1, "",
  1.1381 +                  csa->deck ? 8 : 1, col_name(csa, j));
  1.1382 +               if (data[i])
  1.1383 +                  xfprintf(fp, "%*s%*s", csa->deck ? 2 : 1, "",
  1.1384 +                     csa->deck ? 12 : 1, mps_numb(csa, bnd[i]));
  1.1385 +               xfprintf(fp, "\n"), recno++;
  1.1386 +            }
  1.1387 +         }
  1.1388 +      }
  1.1389 +endt: /* write ENDATA indicator record */
  1.1390 +      xfprintf(fp, "ENDATA\n"), recno++;
  1.1391 +      xfflush(fp);
  1.1392 +      if (xferror(fp))
  1.1393 +      {  xprintf("Write error on `%s' - %s\n", fname, xerrmsg());
  1.1394 +         ret = 1;
  1.1395 +         goto done;
  1.1396 +      }
  1.1397 +      /* problem data has been successfully written */
  1.1398 +      xprintf("%d records were written\n", recno);
  1.1399 +      ret = 0;
  1.1400 +done: if (fp != NULL) xfclose(fp);
  1.1401 +      return ret;
  1.1402 +}
  1.1403 +
  1.1404 +/* eof */