1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpini01.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,577 @@
     1.4 +/* glpini01.c */
     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, 2011 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 +-- triang - find maximal triangular part of a rectangular matrix.
    1.32 +--
    1.33 +-- *Synopsis*
    1.34 +--
    1.35 +-- int triang(int m, int n,
    1.36 +--    void *info, int (*mat)(void *info, int k, int ndx[]),
    1.37 +--    int rn[], int cn[]);
    1.38 +--
    1.39 +-- *Description*
    1.40 +--
    1.41 +-- For a given rectangular (sparse) matrix A with m rows and n columns
    1.42 +-- the routine triang tries to find such permutation matrices P and Q
    1.43 +-- that the first rows and columns of the matrix B = P*A*Q form a lower
    1.44 +-- triangular submatrix of as greatest size as possible:
    1.45 +--
    1.46 +--                   1                       n
    1.47 +--                1  * . . . . . . x x x x x x
    1.48 +--                   * * . . . . . x x x x x x
    1.49 +--                   * * * . . . . x x x x x x
    1.50 +--                   * * * * . . . x x x x x x
    1.51 +--    B = P*A*Q =    * * * * * . . x x x x x x
    1.52 +--                   * * * * * * . x x x x x x
    1.53 +--                   * * * * * * * x x x x x x
    1.54 +--                   x x x x x x x x x x x x x
    1.55 +--                   x x x x x x x x x x x x x
    1.56 +--                m  x x x x x x x x x x x x x
    1.57 +--
    1.58 +-- where: '*' - elements of the lower triangular part, '.' - structural
    1.59 +-- zeros, 'x' - other (either non-zero or zero) elements.
    1.60 +--
    1.61 +-- The parameter info is a transit pointer passed to the formal routine
    1.62 +-- mat (see below).
    1.63 +--
    1.64 +-- The formal routine mat specifies the given matrix A in both row- and
    1.65 +-- column-wise formats. In order to obtain an i-th row of the matrix A
    1.66 +-- the routine triang calls the routine mat with the parameter k = +i,
    1.67 +-- 1 <= i <= m. In response the routine mat should store column indices
    1.68 +-- of (non-zero) elements of the i-th row to the locations ndx[1], ...,
    1.69 +-- ndx[len], where len is number of non-zeros in the i-th row returned
    1.70 +-- on exit. Analogously, in order to obtain a j-th column of the matrix
    1.71 +-- A, the routine mat is called with the parameter k = -j, 1 <= j <= n,
    1.72 +-- and should return pattern of the j-th column in the same way as for
    1.73 +-- row patterns. Note that the routine mat may be called more than once
    1.74 +-- for the same rows and columns.
    1.75 +--
    1.76 +-- On exit the routine computes two resultant arrays rn and cn, which
    1.77 +-- define the permutation matrices P and Q, respectively. The array rn
    1.78 +-- should have at least 1+m locations, where rn[i] = i' (1 <= i <= m)
    1.79 +-- means that i-th row of the original matrix A corresponds to i'-th row
    1.80 +-- of the matrix B = P*A*Q. Similarly, the array cn should have at least
    1.81 +-- 1+n locations, where cn[j] = j' (1 <= j <= n) means that j-th column
    1.82 +-- of the matrix A corresponds to j'-th column of the matrix B.
    1.83 +--
    1.84 +-- *Returns*
    1.85 +--
    1.86 +-- The routine triang returns the size of the lower tringular part of
    1.87 +-- the matrix B = P*A*Q (see the figure above).
    1.88 +--
    1.89 +-- *Complexity*
    1.90 +--
    1.91 +-- The time complexity of the routine triang is O(nnz), where nnz is
    1.92 +-- number of non-zeros in the given matrix A.
    1.93 +--
    1.94 +-- *Algorithm*
    1.95 +--
    1.96 +-- The routine triang starts from the matrix B = P*Q*A, where P and Q
    1.97 +-- are unity matrices, so initially B = A.
    1.98 +--
    1.99 +-- Before the next iteration B = (B1 | B2 | B3), where B1 is partially
   1.100 +-- built a lower triangular submatrix, B2 is the active submatrix, and
   1.101 +-- B3 is a submatrix that contains rejected columns. Thus, the current
   1.102 +-- matrix B looks like follows (initially k1 = 1 and k2 = n):
   1.103 +--
   1.104 +--       1         k1         k2         n
   1.105 +--    1  x . . . . . . . . . . . . . # # #
   1.106 +--       x x . . . . . . . . . . . . # # #
   1.107 +--       x x x . . . . . . . . . . # # # #
   1.108 +--       x x x x . . . . . . . . . # # # #
   1.109 +--       x x x x x . . . . . . . # # # # #
   1.110 +--    k1 x x x x x * * * * * * * # # # # #
   1.111 +--       x x x x x * * * * * * * # # # # #
   1.112 +--       x x x x x * * * * * * * # # # # #
   1.113 +--       x x x x x * * * * * * * # # # # #
   1.114 +--    m  x x x x x * * * * * * * # # # # #
   1.115 +--       <--B1---> <----B2-----> <---B3-->
   1.116 +--
   1.117 +-- On each iteartion the routine looks for a singleton row, i.e. some
   1.118 +-- row that has the only non-zero in the active submatrix B2. If such
   1.119 +-- row exists and the corresponding non-zero is b[i,j], where (by the
   1.120 +-- definition) k1 <= i <= m and k1 <= j <= k2, the routine permutes
   1.121 +-- k1-th and i-th rows and k1-th and j-th columns of the matrix B (in
   1.122 +-- order to place the element in the position b[k1,k1]), removes the
   1.123 +-- k1-th column from the active submatrix B2, and adds this column to
   1.124 +-- the submatrix B1. If no row singletons exist, but B2 is not empty
   1.125 +-- yet, the routine chooses a j-th column, which has maximal number of
   1.126 +-- non-zeros among other columns of B2, removes this column from B2 and
   1.127 +-- adds it to the submatrix B3 in the hope that new row singletons will
   1.128 +-- appear in the active submatrix. */
   1.129 +
   1.130 +static int triang(int m, int n,
   1.131 +      void *info, int (*mat)(void *info, int k, int ndx[]),
   1.132 +      int rn[], int cn[])
   1.133 +{     int *ndx; /* int ndx[1+max(m,n)]; */
   1.134 +      /* this array is used for querying row and column patterns of the
   1.135 +         given matrix A (the third parameter to the routine mat) */
   1.136 +      int *rs_len; /* int rs_len[1+m]; */
   1.137 +      /* rs_len[0] is not used;
   1.138 +         rs_len[i], 1 <= i <= m, is number of non-zeros in the i-th row
   1.139 +         of the matrix A, which (non-zeros) belong to the current active
   1.140 +         submatrix */
   1.141 +      int *rs_head; /* int rs_head[1+n]; */
   1.142 +      /* rs_head[len], 0 <= len <= n, is the number i of the first row
   1.143 +         of the matrix A, for which rs_len[i] = len */
   1.144 +      int *rs_prev; /* int rs_prev[1+m]; */
   1.145 +      /* rs_prev[0] is not used;
   1.146 +         rs_prev[i], 1 <= i <= m, is a number i' of the previous row of
   1.147 +         the matrix A, for which rs_len[i] = rs_len[i'] (zero marks the
   1.148 +         end of this linked list) */
   1.149 +      int *rs_next; /* int rs_next[1+m]; */
   1.150 +      /* rs_next[0] is not used;
   1.151 +         rs_next[i], 1 <= i <= m, is a number i' of the next row of the
   1.152 +         matrix A, for which rs_len[i] = rs_len[i'] (zero marks the end
   1.153 +         this linked list) */
   1.154 +      int cs_head;
   1.155 +      /* is a number j of the first column of the matrix A, which has
   1.156 +         maximal number of non-zeros among other columns */
   1.157 +      int *cs_prev; /* cs_prev[1+n]; */
   1.158 +      /* cs_prev[0] is not used;
   1.159 +         cs_prev[j], 1 <= j <= n, is a number of the previous column of
   1.160 +         the matrix A with the same or greater number of non-zeros than
   1.161 +         in the j-th column (zero marks the end of this linked list) */
   1.162 +      int *cs_next; /* cs_next[1+n]; */
   1.163 +      /* cs_next[0] is not used;
   1.164 +         cs_next[j], 1 <= j <= n, is a number of the next column of
   1.165 +         the matrix A with the same or lesser number of non-zeros than
   1.166 +         in the j-th column (zero marks the end of this linked list) */
   1.167 +      int i, j, ii, jj, k1, k2, len, t, size = 0;
   1.168 +      int *head, *rn_inv, *cn_inv;
   1.169 +      if (!(m > 0 && n > 0))
   1.170 +         xerror("triang: m = %d; n = %d; invalid dimension\n", m, n);
   1.171 +      /* allocate working arrays */
   1.172 +      ndx = xcalloc(1+(m >= n ? m : n), sizeof(int));
   1.173 +      rs_len = xcalloc(1+m, sizeof(int));
   1.174 +      rs_head = xcalloc(1+n, sizeof(int));
   1.175 +      rs_prev = xcalloc(1+m, sizeof(int));
   1.176 +      rs_next = xcalloc(1+m, sizeof(int));
   1.177 +      cs_prev = xcalloc(1+n, sizeof(int));
   1.178 +      cs_next = xcalloc(1+n, sizeof(int));
   1.179 +      /* build linked lists of columns of the matrix A with the same
   1.180 +         number of non-zeros */
   1.181 +      head = rs_len; /* currently rs_len is used as working array */
   1.182 +      for (len = 0; len <= m; len ++) head[len] = 0;
   1.183 +      for (j = 1; j <= n; j++)
   1.184 +      {  /* obtain length of the j-th column */
   1.185 +         len = mat(info, -j, ndx);
   1.186 +         xassert(0 <= len && len <= m);
   1.187 +         /* include the j-th column in the corresponding linked list */
   1.188 +         cs_prev[j] = head[len];
   1.189 +         head[len] = j;
   1.190 +      }
   1.191 +      /* merge all linked lists of columns in one linked list, where
   1.192 +         columns are ordered by descending of their lengths */
   1.193 +      cs_head = 0;
   1.194 +      for (len = 0; len <= m; len++)
   1.195 +      {  for (j = head[len]; j != 0; j = cs_prev[j])
   1.196 +         {  cs_next[j] = cs_head;
   1.197 +            cs_head = j;
   1.198 +         }
   1.199 +      }
   1.200 +      jj = 0;
   1.201 +      for (j = cs_head; j != 0; j = cs_next[j])
   1.202 +      {  cs_prev[j] = jj;
   1.203 +         jj = j;
   1.204 +      }
   1.205 +      /* build initial doubly linked lists of rows of the matrix A with
   1.206 +         the same number of non-zeros */
   1.207 +      for (len = 0; len <= n; len++) rs_head[len] = 0;
   1.208 +      for (i = 1; i <= m; i++)
   1.209 +      {  /* obtain length of the i-th row */
   1.210 +         rs_len[i] = len = mat(info, +i, ndx);
   1.211 +         xassert(0 <= len && len <= n);
   1.212 +         /* include the i-th row in the correspondng linked list */
   1.213 +         rs_prev[i] = 0;
   1.214 +         rs_next[i] = rs_head[len];
   1.215 +         if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
   1.216 +         rs_head[len] = i;
   1.217 +      }
   1.218 +      /* initially all rows and columns of the matrix A are active */
   1.219 +      for (i = 1; i <= m; i++) rn[i] = 0;
   1.220 +      for (j = 1; j <= n; j++) cn[j] = 0;
   1.221 +      /* set initial bounds of the active submatrix */
   1.222 +      k1 = 1, k2 = n;
   1.223 +      /* main loop starts here */
   1.224 +      while (k1 <= k2)
   1.225 +      {  i = rs_head[1];
   1.226 +         if (i != 0)
   1.227 +         {  /* the i-th row of the matrix A is a row singleton, since
   1.228 +               it has the only non-zero in the active submatrix */
   1.229 +            xassert(rs_len[i] == 1);
   1.230 +            /* determine the number j of an active column of the matrix
   1.231 +               A, in which this non-zero is placed */
   1.232 +            j = 0;
   1.233 +            t = mat(info, +i, ndx);
   1.234 +            xassert(0 <= t && t <= n);
   1.235 +            for (t = t; t >= 1; t--)
   1.236 +            {  jj = ndx[t];
   1.237 +               xassert(1 <= jj && jj <= n);
   1.238 +               if (cn[jj] == 0)
   1.239 +               {  xassert(j == 0);
   1.240 +                  j = jj;
   1.241 +               }
   1.242 +            }
   1.243 +            xassert(j != 0);
   1.244 +            /* the singleton is a[i,j]; move a[i,j] to the position
   1.245 +               b[k1,k1] of the matrix B */
   1.246 +            rn[i] = cn[j] = k1;
   1.247 +            /* shift the left bound of the active submatrix */
   1.248 +            k1++;
   1.249 +            /* increase the size of the lower triangular part */
   1.250 +            size++;
   1.251 +         }
   1.252 +         else
   1.253 +         {  /* the current active submatrix has no row singletons */
   1.254 +            /* remove an active column with maximal number of non-zeros
   1.255 +               from the active submatrix */
   1.256 +            j = cs_head;
   1.257 +            xassert(j != 0);
   1.258 +            cn[j] = k2;
   1.259 +            /* shift the right bound of the active submatrix */
   1.260 +            k2--;
   1.261 +         }
   1.262 +         /* the j-th column of the matrix A has been removed from the
   1.263 +            active submatrix */
   1.264 +         /* remove the j-th column from the linked list */
   1.265 +         if (cs_prev[j] == 0)
   1.266 +            cs_head = cs_next[j];
   1.267 +         else
   1.268 +            cs_next[cs_prev[j]] = cs_next[j];
   1.269 +         if (cs_next[j] == 0)
   1.270 +            /* nop */;
   1.271 +         else
   1.272 +            cs_prev[cs_next[j]] = cs_prev[j];
   1.273 +         /* go through non-zeros of the j-th columns and update active
   1.274 +            lengths of the corresponding rows */
   1.275 +         t = mat(info, -j, ndx);
   1.276 +         xassert(0 <= t && t <= m);
   1.277 +         for (t = t; t >= 1; t--)
   1.278 +         {  i = ndx[t];
   1.279 +            xassert(1 <= i && i <= m);
   1.280 +            /* the non-zero a[i,j] has left the active submatrix */
   1.281 +            len = rs_len[i];
   1.282 +            xassert(len >= 1);
   1.283 +            /* remove the i-th row from the linked list of rows with
   1.284 +               active length len */
   1.285 +            if (rs_prev[i] == 0)
   1.286 +               rs_head[len] = rs_next[i];
   1.287 +            else
   1.288 +               rs_next[rs_prev[i]] = rs_next[i];
   1.289 +            if (rs_next[i] == 0)
   1.290 +               /* nop */;
   1.291 +            else
   1.292 +               rs_prev[rs_next[i]] = rs_prev[i];
   1.293 +            /* decrease the active length of the i-th row */
   1.294 +            rs_len[i] = --len;
   1.295 +            /* return the i-th row to the corresponding linked list */
   1.296 +            rs_prev[i] = 0;
   1.297 +            rs_next[i] = rs_head[len];
   1.298 +            if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
   1.299 +            rs_head[len] = i;
   1.300 +         }
   1.301 +      }
   1.302 +      /* other rows of the matrix A, which are still active, correspond
   1.303 +         to rows k1, ..., m of the matrix B (in arbitrary order) */
   1.304 +      for (i = 1; i <= m; i++) if (rn[i] == 0) rn[i] = k1++;
   1.305 +      /* but for columns this is not needed, because now the submatrix
   1.306 +         B2 has no columns */
   1.307 +      for (j = 1; j <= n; j++) xassert(cn[j] != 0);
   1.308 +      /* perform some optional checks */
   1.309 +      /* make sure that rn is a permutation of {1, ..., m} and cn is a
   1.310 +         permutation of {1, ..., n} */
   1.311 +      rn_inv = rs_len; /* used as working array */
   1.312 +      for (ii = 1; ii <= m; ii++) rn_inv[ii] = 0;
   1.313 +      for (i = 1; i <= m; i++)
   1.314 +      {  ii = rn[i];
   1.315 +         xassert(1 <= ii && ii <= m);
   1.316 +         xassert(rn_inv[ii] == 0);
   1.317 +         rn_inv[ii] = i;
   1.318 +      }
   1.319 +      cn_inv = rs_head; /* used as working array */
   1.320 +      for (jj = 1; jj <= n; jj++) cn_inv[jj] = 0;
   1.321 +      for (j = 1; j <= n; j++)
   1.322 +      {  jj = cn[j];
   1.323 +         xassert(1 <= jj && jj <= n);
   1.324 +         xassert(cn_inv[jj] == 0);
   1.325 +         cn_inv[jj] = j;
   1.326 +      }
   1.327 +      /* make sure that the matrix B = P*A*Q really has the form, which
   1.328 +         was declared */
   1.329 +      for (ii = 1; ii <= size; ii++)
   1.330 +      {  int diag = 0;
   1.331 +         i = rn_inv[ii];
   1.332 +         t = mat(info, +i, ndx);
   1.333 +         xassert(0 <= t && t <= n);
   1.334 +         for (t = t; t >= 1; t--)
   1.335 +         {  j = ndx[t];
   1.336 +            xassert(1 <= j && j <= n);
   1.337 +            jj = cn[j];
   1.338 +            if (jj <= size) xassert(jj <= ii);
   1.339 +            if (jj == ii)
   1.340 +            {  xassert(!diag);
   1.341 +               diag = 1;
   1.342 +            }
   1.343 +         }
   1.344 +         xassert(diag);
   1.345 +      }
   1.346 +      /* free working arrays */
   1.347 +      xfree(ndx);
   1.348 +      xfree(rs_len);
   1.349 +      xfree(rs_head);
   1.350 +      xfree(rs_prev);
   1.351 +      xfree(rs_next);
   1.352 +      xfree(cs_prev);
   1.353 +      xfree(cs_next);
   1.354 +      /* return to the calling program */
   1.355 +      return size;
   1.356 +}
   1.357 +
   1.358 +/*----------------------------------------------------------------------
   1.359 +-- adv_basis - construct advanced initial LP basis.
   1.360 +--
   1.361 +-- *Synopsis*
   1.362 +--
   1.363 +-- #include "glpini.h"
   1.364 +-- void adv_basis(glp_prob *lp);
   1.365 +--
   1.366 +-- *Description*
   1.367 +--
   1.368 +-- The routine adv_basis constructs an advanced initial basis for an LP
   1.369 +-- problem object, which the parameter lp points to.
   1.370 +--
   1.371 +-- In order to build the initial basis the routine does the following:
   1.372 +--
   1.373 +-- 1) includes in the basis all non-fixed auxiliary variables;
   1.374 +--
   1.375 +-- 2) includes in the basis as many as possible non-fixed structural
   1.376 +--    variables preserving triangular form of the basis matrix;
   1.377 +--
   1.378 +-- 3) includes in the basis appropriate (fixed) auxiliary variables
   1.379 +--    in order to complete the basis.
   1.380 +--
   1.381 +-- As a result the initial basis has minimum of fixed variables and the
   1.382 +-- corresponding basis matrix is triangular. */
   1.383 +
   1.384 +static int mat(void *info, int k, int ndx[])
   1.385 +{     /* this auxiliary routine returns the pattern of a given row or
   1.386 +         a given column of the augmented constraint matrix A~ = (I|-A),
   1.387 +         in which columns of fixed variables are implicitly cleared */
   1.388 +      LPX *lp = info;
   1.389 +      int m = lpx_get_num_rows(lp);
   1.390 +      int n = lpx_get_num_cols(lp);
   1.391 +      int typx, i, j, lll, len = 0;
   1.392 +      if (k > 0)
   1.393 +      {  /* the pattern of the i-th row is required */
   1.394 +         i = +k;
   1.395 +         xassert(1 <= i && i <= m);
   1.396 +#if 0 /* 22/XII-2003 */
   1.397 +         /* if the auxiliary variable x[i] is non-fixed, include its
   1.398 +            element (placed in the i-th column) in the pattern */
   1.399 +         lpx_get_row_bnds(lp, i, &typx, NULL, NULL);
   1.400 +         if (typx != LPX_FX) ndx[++len] = i;
   1.401 +         /* include in the pattern elements placed in columns, which
   1.402 +            correspond to non-fixed structural varables */
   1.403 +         i_beg = aa_ptr[i];
   1.404 +         i_end = i_beg + aa_len[i] - 1;
   1.405 +         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
   1.406 +         {  j = m + sv_ndx[i_ptr];
   1.407 +            lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL);
   1.408 +            if (typx != LPX_FX) ndx[++len] = j;
   1.409 +         }
   1.410 +#else
   1.411 +         lll = lpx_get_mat_row(lp, i, ndx, NULL);
   1.412 +         for (k = 1; k <= lll; k++)
   1.413 +         {  lpx_get_col_bnds(lp, ndx[k], &typx, NULL, NULL);
   1.414 +            if (typx != LPX_FX) ndx[++len] = m + ndx[k];
   1.415 +         }
   1.416 +         lpx_get_row_bnds(lp, i, &typx, NULL, NULL);
   1.417 +         if (typx != LPX_FX) ndx[++len] = i;
   1.418 +#endif
   1.419 +      }
   1.420 +      else
   1.421 +      {  /* the pattern of the j-th column is required */
   1.422 +         j = -k;
   1.423 +         xassert(1 <= j && j <= m+n);
   1.424 +         /* if the (auxiliary or structural) variable x[j] is fixed,
   1.425 +            the pattern of its column is empty */
   1.426 +         if (j <= m)
   1.427 +            lpx_get_row_bnds(lp, j, &typx, NULL, NULL);
   1.428 +         else
   1.429 +            lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL);
   1.430 +         if (typx != LPX_FX)
   1.431 +         {  if (j <= m)
   1.432 +            {  /* x[j] is non-fixed auxiliary variable */
   1.433 +               ndx[++len] = j;
   1.434 +            }
   1.435 +            else
   1.436 +            {  /* x[j] is non-fixed structural variables */
   1.437 +#if 0 /* 22/XII-2003 */
   1.438 +               j_beg = aa_ptr[j];
   1.439 +               j_end = j_beg + aa_len[j] - 1;
   1.440 +               for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
   1.441 +                  ndx[++len] = sv_ndx[j_ptr];
   1.442 +#else
   1.443 +               len = lpx_get_mat_col(lp, j-m, ndx, NULL);
   1.444 +#endif
   1.445 +            }
   1.446 +         }
   1.447 +      }
   1.448 +      /* return the length of the row/column pattern */
   1.449 +      return len;
   1.450 +}
   1.451 +
   1.452 +static void adv_basis(glp_prob *lp)
   1.453 +{     int m = lpx_get_num_rows(lp);
   1.454 +      int n = lpx_get_num_cols(lp);
   1.455 +      int i, j, jj, k, size;
   1.456 +      int *rn, *cn, *rn_inv, *cn_inv;
   1.457 +      int typx, *tagx = xcalloc(1+m+n, sizeof(int));
   1.458 +      double lb, ub;
   1.459 +      xprintf("Constructing initial basis...\n");
   1.460 +#if 0 /* 13/V-2009 */
   1.461 +      if (m == 0)
   1.462 +         xerror("glp_adv_basis: problem has no rows\n");
   1.463 +      if (n == 0)
   1.464 +         xerror("glp_adv_basis: problem has no columns\n");
   1.465 +#else
   1.466 +      if (m == 0 || n == 0)
   1.467 +      {  glp_std_basis(lp);
   1.468 +         return;
   1.469 +      }
   1.470 +#endif
   1.471 +      /* use the routine triang (see above) to find maximal triangular
   1.472 +         part of the augmented constraint matrix A~ = (I|-A); in order
   1.473 +         to prevent columns of fixed variables to be included in the
   1.474 +         triangular part, such columns are implictly removed from the
   1.475 +         matrix A~ by the routine adv_mat */
   1.476 +      rn = xcalloc(1+m, sizeof(int));
   1.477 +      cn = xcalloc(1+m+n, sizeof(int));
   1.478 +      size = triang(m, m+n, lp, mat, rn, cn);
   1.479 +      if (lpx_get_int_parm(lp, LPX_K_MSGLEV) >= 3)
   1.480 +         xprintf("Size of triangular part = %d\n", size);
   1.481 +      /* the first size rows and columns of the matrix P*A~*Q (where
   1.482 +         P and Q are permutation matrices defined by the arrays rn and
   1.483 +         cn) form a lower triangular matrix; build the arrays (rn_inv
   1.484 +         and cn_inv), which define the matrices inv(P) and inv(Q) */
   1.485 +      rn_inv = xcalloc(1+m, sizeof(int));
   1.486 +      cn_inv = xcalloc(1+m+n, sizeof(int));
   1.487 +      for (i = 1; i <= m; i++) rn_inv[rn[i]] = i;
   1.488 +      for (j = 1; j <= m+n; j++) cn_inv[cn[j]] = j;
   1.489 +      /* include the columns of the matrix A~, which correspond to the
   1.490 +         first size columns of the matrix P*A~*Q, in the basis */
   1.491 +      for (k = 1; k <= m+n; k++) tagx[k] = -1;
   1.492 +      for (jj = 1; jj <= size; jj++)
   1.493 +      {  j = cn_inv[jj];
   1.494 +         /* the j-th column of A~ is the jj-th column of P*A~*Q */
   1.495 +         tagx[j] = LPX_BS;
   1.496 +      }
   1.497 +      /* if size < m, we need to add appropriate columns of auxiliary
   1.498 +         variables to the basis */
   1.499 +      for (jj = size + 1; jj <= m; jj++)
   1.500 +      {  /* the jj-th column of P*A~*Q should be replaced by the column
   1.501 +            of the auxiliary variable, for which the only unity element
   1.502 +            is placed in the position [jj,jj] */
   1.503 +         i = rn_inv[jj];
   1.504 +         /* the jj-th row of P*A~*Q is the i-th row of A~, but in the
   1.505 +            i-th row of A~ the unity element belongs to the i-th column
   1.506 +            of A~; therefore the disired column corresponds to the i-th
   1.507 +            auxiliary variable (note that this column doesn't belong to
   1.508 +            the triangular part found by the routine triang) */
   1.509 +         xassert(1 <= i && i <= m);
   1.510 +         xassert(cn[i] > size);
   1.511 +         tagx[i] = LPX_BS;
   1.512 +      }
   1.513 +      /* free working arrays */
   1.514 +      xfree(rn);
   1.515 +      xfree(cn);
   1.516 +      xfree(rn_inv);
   1.517 +      xfree(cn_inv);
   1.518 +      /* build tags of non-basic variables */
   1.519 +      for (k = 1; k <= m+n; k++)
   1.520 +      {  if (tagx[k] != LPX_BS)
   1.521 +         {  if (k <= m)
   1.522 +               lpx_get_row_bnds(lp, k, &typx, &lb, &ub);
   1.523 +            else
   1.524 +               lpx_get_col_bnds(lp, k-m, &typx, &lb, &ub);
   1.525 +            switch (typx)
   1.526 +            {  case LPX_FR:
   1.527 +                  tagx[k] = LPX_NF; break;
   1.528 +               case LPX_LO:
   1.529 +                  tagx[k] = LPX_NL; break;
   1.530 +               case LPX_UP:
   1.531 +                  tagx[k] = LPX_NU; break;
   1.532 +               case LPX_DB:
   1.533 +                  tagx[k] =
   1.534 +                     (fabs(lb) <= fabs(ub) ? LPX_NL : LPX_NU);
   1.535 +                  break;
   1.536 +               case LPX_FX:
   1.537 +                  tagx[k] = LPX_NS; break;
   1.538 +               default:
   1.539 +                  xassert(typx != typx);
   1.540 +            }
   1.541 +         }
   1.542 +      }
   1.543 +      for (k = 1; k <= m+n; k++)
   1.544 +      {  if (k <= m)
   1.545 +            lpx_set_row_stat(lp, k, tagx[k]);
   1.546 +         else
   1.547 +            lpx_set_col_stat(lp, k-m, tagx[k]);
   1.548 +      }
   1.549 +      xfree(tagx);
   1.550 +      return;
   1.551 +}
   1.552 +
   1.553 +/***********************************************************************
   1.554 +*  NAME
   1.555 +*
   1.556 +*  glp_adv_basis - construct advanced initial LP basis
   1.557 +*
   1.558 +*  SYNOPSIS
   1.559 +*
   1.560 +*  void glp_adv_basis(glp_prob *lp, int flags);
   1.561 +*
   1.562 +*  DESCRIPTION
   1.563 +*
   1.564 +*  The routine glp_adv_basis constructs an advanced initial basis for
   1.565 +*  the specified problem object.
   1.566 +*
   1.567 +*  The parameter flags is reserved for use in the future and must be
   1.568 +*  specified as zero. */
   1.569 +
   1.570 +void glp_adv_basis(glp_prob *lp, int flags)
   1.571 +{     if (flags != 0)
   1.572 +         xerror("glp_adv_basis: flags = %d; invalid flags\n", flags);
   1.573 +      if (lp->m == 0 || lp->n == 0)
   1.574 +         glp_std_basis(lp);
   1.575 +      else
   1.576 +         adv_basis(lp);
   1.577 +      return;
   1.578 +}
   1.579 +
   1.580 +/* eof */