COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpini01.c @ 10:5545663ca997

subpack-glpk
Last change on this file since 10:5545663ca997 was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 13 years ago

Import GLPK 4.47

File size: 22.1 KB
Line 
1/* glpini01.c */
2
3/***********************************************************************
4*  This code is part of GLPK (GNU Linear Programming Kit).
5*
6*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
8*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9*  E-mail: <mao@gnu.org>.
10*
11*  GLPK is free software: you can redistribute it and/or modify it
12*  under the terms of the GNU General Public License as published by
13*  the Free Software Foundation, either version 3 of the License, or
14*  (at your option) any later version.
15*
16*  GLPK is distributed in the hope that it will be useful, but WITHOUT
17*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19*  License for more details.
20*
21*  You should have received a copy of the GNU General Public License
22*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23***********************************************************************/
24
25#include "glpapi.h"
26
27/*----------------------------------------------------------------------
28-- triang - find maximal triangular part of a rectangular matrix.
29--
30-- *Synopsis*
31--
32-- int triang(int m, int n,
33--    void *info, int (*mat)(void *info, int k, int ndx[]),
34--    int rn[], int cn[]);
35--
36-- *Description*
37--
38-- For a given rectangular (sparse) matrix A with m rows and n columns
39-- the routine triang tries to find such permutation matrices P and Q
40-- that the first rows and columns of the matrix B = P*A*Q form a lower
41-- triangular submatrix of as greatest size as possible:
42--
43--                   1                       n
44--                1  * . . . . . . x x x x x x
45--                   * * . . . . . x x x x x x
46--                   * * * . . . . x x x x x x
47--                   * * * * . . . x x x x x x
48--    B = P*A*Q =    * * * * * . . x x x x x x
49--                   * * * * * * . x x x x x x
50--                   * * * * * * * x x x x x x
51--                   x x x x x x x x x x x x x
52--                   x x x x x x x x x x x x x
53--                m  x x x x x x x x x x x x x
54--
55-- where: '*' - elements of the lower triangular part, '.' - structural
56-- zeros, 'x' - other (either non-zero or zero) elements.
57--
58-- The parameter info is a transit pointer passed to the formal routine
59-- mat (see below).
60--
61-- The formal routine mat specifies the given matrix A in both row- and
62-- column-wise formats. In order to obtain an i-th row of the matrix A
63-- the routine triang calls the routine mat with the parameter k = +i,
64-- 1 <= i <= m. In response the routine mat should store column indices
65-- of (non-zero) elements of the i-th row to the locations ndx[1], ...,
66-- ndx[len], where len is number of non-zeros in the i-th row returned
67-- on exit. Analogously, in order to obtain a j-th column of the matrix
68-- A, the routine mat is called with the parameter k = -j, 1 <= j <= n,
69-- and should return pattern of the j-th column in the same way as for
70-- row patterns. Note that the routine mat may be called more than once
71-- for the same rows and columns.
72--
73-- On exit the routine computes two resultant arrays rn and cn, which
74-- define the permutation matrices P and Q, respectively. The array rn
75-- should have at least 1+m locations, where rn[i] = i' (1 <= i <= m)
76-- means that i-th row of the original matrix A corresponds to i'-th row
77-- of the matrix B = P*A*Q. Similarly, the array cn should have at least
78-- 1+n locations, where cn[j] = j' (1 <= j <= n) means that j-th column
79-- of the matrix A corresponds to j'-th column of the matrix B.
80--
81-- *Returns*
82--
83-- The routine triang returns the size of the lower tringular part of
84-- the matrix B = P*A*Q (see the figure above).
85--
86-- *Complexity*
87--
88-- The time complexity of the routine triang is O(nnz), where nnz is
89-- number of non-zeros in the given matrix A.
90--
91-- *Algorithm*
92--
93-- The routine triang starts from the matrix B = P*Q*A, where P and Q
94-- are unity matrices, so initially B = A.
95--
96-- Before the next iteration B = (B1 | B2 | B3), where B1 is partially
97-- built a lower triangular submatrix, B2 is the active submatrix, and
98-- B3 is a submatrix that contains rejected columns. Thus, the current
99-- matrix B looks like follows (initially k1 = 1 and k2 = n):
100--
101--       1         k1         k2         n
102--    1  x . . . . . . . . . . . . . # # #
103--       x x . . . . . . . . . . . . # # #
104--       x x x . . . . . . . . . . # # # #
105--       x x x x . . . . . . . . . # # # #
106--       x x x x x . . . . . . . # # # # #
107--    k1 x x x x x * * * * * * * # # # # #
108--       x x x x x * * * * * * * # # # # #
109--       x x x x x * * * * * * * # # # # #
110--       x x x x x * * * * * * * # # # # #
111--    m  x x x x x * * * * * * * # # # # #
112--       <--B1---> <----B2-----> <---B3-->
113--
114-- On each iteartion the routine looks for a singleton row, i.e. some
115-- row that has the only non-zero in the active submatrix B2. If such
116-- row exists and the corresponding non-zero is b[i,j], where (by the
117-- definition) k1 <= i <= m and k1 <= j <= k2, the routine permutes
118-- k1-th and i-th rows and k1-th and j-th columns of the matrix B (in
119-- order to place the element in the position b[k1,k1]), removes the
120-- k1-th column from the active submatrix B2, and adds this column to
121-- the submatrix B1. If no row singletons exist, but B2 is not empty
122-- yet, the routine chooses a j-th column, which has maximal number of
123-- non-zeros among other columns of B2, removes this column from B2 and
124-- adds it to the submatrix B3 in the hope that new row singletons will
125-- appear in the active submatrix. */
126
127static int triang(int m, int n,
128      void *info, int (*mat)(void *info, int k, int ndx[]),
129      int rn[], int cn[])
130{     int *ndx; /* int ndx[1+max(m,n)]; */
131      /* this array is used for querying row and column patterns of the
132         given matrix A (the third parameter to the routine mat) */
133      int *rs_len; /* int rs_len[1+m]; */
134      /* rs_len[0] is not used;
135         rs_len[i], 1 <= i <= m, is number of non-zeros in the i-th row
136         of the matrix A, which (non-zeros) belong to the current active
137         submatrix */
138      int *rs_head; /* int rs_head[1+n]; */
139      /* rs_head[len], 0 <= len <= n, is the number i of the first row
140         of the matrix A, for which rs_len[i] = len */
141      int *rs_prev; /* int rs_prev[1+m]; */
142      /* rs_prev[0] is not used;
143         rs_prev[i], 1 <= i <= m, is a number i' of the previous row of
144         the matrix A, for which rs_len[i] = rs_len[i'] (zero marks the
145         end of this linked list) */
146      int *rs_next; /* int rs_next[1+m]; */
147      /* rs_next[0] is not used;
148         rs_next[i], 1 <= i <= m, is a number i' of the next row of the
149         matrix A, for which rs_len[i] = rs_len[i'] (zero marks the end
150         this linked list) */
151      int cs_head;
152      /* is a number j of the first column of the matrix A, which has
153         maximal number of non-zeros among other columns */
154      int *cs_prev; /* cs_prev[1+n]; */
155      /* cs_prev[0] is not used;
156         cs_prev[j], 1 <= j <= n, is a number of the previous column of
157         the matrix A with the same or greater number of non-zeros than
158         in the j-th column (zero marks the end of this linked list) */
159      int *cs_next; /* cs_next[1+n]; */
160      /* cs_next[0] is not used;
161         cs_next[j], 1 <= j <= n, is a number of the next column of
162         the matrix A with the same or lesser number of non-zeros than
163         in the j-th column (zero marks the end of this linked list) */
164      int i, j, ii, jj, k1, k2, len, t, size = 0;
165      int *head, *rn_inv, *cn_inv;
166      if (!(m > 0 && n > 0))
167         xerror("triang: m = %d; n = %d; invalid dimension\n", m, n);
168      /* allocate working arrays */
169      ndx = xcalloc(1+(m >= n ? m : n), sizeof(int));
170      rs_len = xcalloc(1+m, sizeof(int));
171      rs_head = xcalloc(1+n, sizeof(int));
172      rs_prev = xcalloc(1+m, sizeof(int));
173      rs_next = xcalloc(1+m, sizeof(int));
174      cs_prev = xcalloc(1+n, sizeof(int));
175      cs_next = xcalloc(1+n, sizeof(int));
176      /* build linked lists of columns of the matrix A with the same
177         number of non-zeros */
178      head = rs_len; /* currently rs_len is used as working array */
179      for (len = 0; len <= m; len ++) head[len] = 0;
180      for (j = 1; j <= n; j++)
181      {  /* obtain length of the j-th column */
182         len = mat(info, -j, ndx);
183         xassert(0 <= len && len <= m);
184         /* include the j-th column in the corresponding linked list */
185         cs_prev[j] = head[len];
186         head[len] = j;
187      }
188      /* merge all linked lists of columns in one linked list, where
189         columns are ordered by descending of their lengths */
190      cs_head = 0;
191      for (len = 0; len <= m; len++)
192      {  for (j = head[len]; j != 0; j = cs_prev[j])
193         {  cs_next[j] = cs_head;
194            cs_head = j;
195         }
196      }
197      jj = 0;
198      for (j = cs_head; j != 0; j = cs_next[j])
199      {  cs_prev[j] = jj;
200         jj = j;
201      }
202      /* build initial doubly linked lists of rows of the matrix A with
203         the same number of non-zeros */
204      for (len = 0; len <= n; len++) rs_head[len] = 0;
205      for (i = 1; i <= m; i++)
206      {  /* obtain length of the i-th row */
207         rs_len[i] = len = mat(info, +i, ndx);
208         xassert(0 <= len && len <= n);
209         /* include the i-th row in the correspondng linked list */
210         rs_prev[i] = 0;
211         rs_next[i] = rs_head[len];
212         if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
213         rs_head[len] = i;
214      }
215      /* initially all rows and columns of the matrix A are active */
216      for (i = 1; i <= m; i++) rn[i] = 0;
217      for (j = 1; j <= n; j++) cn[j] = 0;
218      /* set initial bounds of the active submatrix */
219      k1 = 1, k2 = n;
220      /* main loop starts here */
221      while (k1 <= k2)
222      {  i = rs_head[1];
223         if (i != 0)
224         {  /* the i-th row of the matrix A is a row singleton, since
225               it has the only non-zero in the active submatrix */
226            xassert(rs_len[i] == 1);
227            /* determine the number j of an active column of the matrix
228               A, in which this non-zero is placed */
229            j = 0;
230            t = mat(info, +i, ndx);
231            xassert(0 <= t && t <= n);
232            for (t = t; t >= 1; t--)
233            {  jj = ndx[t];
234               xassert(1 <= jj && jj <= n);
235               if (cn[jj] == 0)
236               {  xassert(j == 0);
237                  j = jj;
238               }
239            }
240            xassert(j != 0);
241            /* the singleton is a[i,j]; move a[i,j] to the position
242               b[k1,k1] of the matrix B */
243            rn[i] = cn[j] = k1;
244            /* shift the left bound of the active submatrix */
245            k1++;
246            /* increase the size of the lower triangular part */
247            size++;
248         }
249         else
250         {  /* the current active submatrix has no row singletons */
251            /* remove an active column with maximal number of non-zeros
252               from the active submatrix */
253            j = cs_head;
254            xassert(j != 0);
255            cn[j] = k2;
256            /* shift the right bound of the active submatrix */
257            k2--;
258         }
259         /* the j-th column of the matrix A has been removed from the
260            active submatrix */
261         /* remove the j-th column from the linked list */
262         if (cs_prev[j] == 0)
263            cs_head = cs_next[j];
264         else
265            cs_next[cs_prev[j]] = cs_next[j];
266         if (cs_next[j] == 0)
267            /* nop */;
268         else
269            cs_prev[cs_next[j]] = cs_prev[j];
270         /* go through non-zeros of the j-th columns and update active
271            lengths of the corresponding rows */
272         t = mat(info, -j, ndx);
273         xassert(0 <= t && t <= m);
274         for (t = t; t >= 1; t--)
275         {  i = ndx[t];
276            xassert(1 <= i && i <= m);
277            /* the non-zero a[i,j] has left the active submatrix */
278            len = rs_len[i];
279            xassert(len >= 1);
280            /* remove the i-th row from the linked list of rows with
281               active length len */
282            if (rs_prev[i] == 0)
283               rs_head[len] = rs_next[i];
284            else
285               rs_next[rs_prev[i]] = rs_next[i];
286            if (rs_next[i] == 0)
287               /* nop */;
288            else
289               rs_prev[rs_next[i]] = rs_prev[i];
290            /* decrease the active length of the i-th row */
291            rs_len[i] = --len;
292            /* return the i-th row to the corresponding linked list */
293            rs_prev[i] = 0;
294            rs_next[i] = rs_head[len];
295            if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
296            rs_head[len] = i;
297         }
298      }
299      /* other rows of the matrix A, which are still active, correspond
300         to rows k1, ..., m of the matrix B (in arbitrary order) */
301      for (i = 1; i <= m; i++) if (rn[i] == 0) rn[i] = k1++;
302      /* but for columns this is not needed, because now the submatrix
303         B2 has no columns */
304      for (j = 1; j <= n; j++) xassert(cn[j] != 0);
305      /* perform some optional checks */
306      /* make sure that rn is a permutation of {1, ..., m} and cn is a
307         permutation of {1, ..., n} */
308      rn_inv = rs_len; /* used as working array */
309      for (ii = 1; ii <= m; ii++) rn_inv[ii] = 0;
310      for (i = 1; i <= m; i++)
311      {  ii = rn[i];
312         xassert(1 <= ii && ii <= m);
313         xassert(rn_inv[ii] == 0);
314         rn_inv[ii] = i;
315      }
316      cn_inv = rs_head; /* used as working array */
317      for (jj = 1; jj <= n; jj++) cn_inv[jj] = 0;
318      for (j = 1; j <= n; j++)
319      {  jj = cn[j];
320         xassert(1 <= jj && jj <= n);
321         xassert(cn_inv[jj] == 0);
322         cn_inv[jj] = j;
323      }
324      /* make sure that the matrix B = P*A*Q really has the form, which
325         was declared */
326      for (ii = 1; ii <= size; ii++)
327      {  int diag = 0;
328         i = rn_inv[ii];
329         t = mat(info, +i, ndx);
330         xassert(0 <= t && t <= n);
331         for (t = t; t >= 1; t--)
332         {  j = ndx[t];
333            xassert(1 <= j && j <= n);
334            jj = cn[j];
335            if (jj <= size) xassert(jj <= ii);
336            if (jj == ii)
337            {  xassert(!diag);
338               diag = 1;
339            }
340         }
341         xassert(diag);
342      }
343      /* free working arrays */
344      xfree(ndx);
345      xfree(rs_len);
346      xfree(rs_head);
347      xfree(rs_prev);
348      xfree(rs_next);
349      xfree(cs_prev);
350      xfree(cs_next);
351      /* return to the calling program */
352      return size;
353}
354
355/*----------------------------------------------------------------------
356-- adv_basis - construct advanced initial LP basis.
357--
358-- *Synopsis*
359--
360-- #include "glpini.h"
361-- void adv_basis(glp_prob *lp);
362--
363-- *Description*
364--
365-- The routine adv_basis constructs an advanced initial basis for an LP
366-- problem object, which the parameter lp points to.
367--
368-- In order to build the initial basis the routine does the following:
369--
370-- 1) includes in the basis all non-fixed auxiliary variables;
371--
372-- 2) includes in the basis as many as possible non-fixed structural
373--    variables preserving triangular form of the basis matrix;
374--
375-- 3) includes in the basis appropriate (fixed) auxiliary variables
376--    in order to complete the basis.
377--
378-- As a result the initial basis has minimum of fixed variables and the
379-- corresponding basis matrix is triangular. */
380
381static int mat(void *info, int k, int ndx[])
382{     /* this auxiliary routine returns the pattern of a given row or
383         a given column of the augmented constraint matrix A~ = (I|-A),
384         in which columns of fixed variables are implicitly cleared */
385      LPX *lp = info;
386      int m = lpx_get_num_rows(lp);
387      int n = lpx_get_num_cols(lp);
388      int typx, i, j, lll, len = 0;
389      if (k > 0)
390      {  /* the pattern of the i-th row is required */
391         i = +k;
392         xassert(1 <= i && i <= m);
393#if 0 /* 22/XII-2003 */
394         /* if the auxiliary variable x[i] is non-fixed, include its
395            element (placed in the i-th column) in the pattern */
396         lpx_get_row_bnds(lp, i, &typx, NULL, NULL);
397         if (typx != LPX_FX) ndx[++len] = i;
398         /* include in the pattern elements placed in columns, which
399            correspond to non-fixed structural varables */
400         i_beg = aa_ptr[i];
401         i_end = i_beg + aa_len[i] - 1;
402         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
403         {  j = m + sv_ndx[i_ptr];
404            lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL);
405            if (typx != LPX_FX) ndx[++len] = j;
406         }
407#else
408         lll = lpx_get_mat_row(lp, i, ndx, NULL);
409         for (k = 1; k <= lll; k++)
410         {  lpx_get_col_bnds(lp, ndx[k], &typx, NULL, NULL);
411            if (typx != LPX_FX) ndx[++len] = m + ndx[k];
412         }
413         lpx_get_row_bnds(lp, i, &typx, NULL, NULL);
414         if (typx != LPX_FX) ndx[++len] = i;
415#endif
416      }
417      else
418      {  /* the pattern of the j-th column is required */
419         j = -k;
420         xassert(1 <= j && j <= m+n);
421         /* if the (auxiliary or structural) variable x[j] is fixed,
422            the pattern of its column is empty */
423         if (j <= m)
424            lpx_get_row_bnds(lp, j, &typx, NULL, NULL);
425         else
426            lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL);
427         if (typx != LPX_FX)
428         {  if (j <= m)
429            {  /* x[j] is non-fixed auxiliary variable */
430               ndx[++len] = j;
431            }
432            else
433            {  /* x[j] is non-fixed structural variables */
434#if 0 /* 22/XII-2003 */
435               j_beg = aa_ptr[j];
436               j_end = j_beg + aa_len[j] - 1;
437               for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
438                  ndx[++len] = sv_ndx[j_ptr];
439#else
440               len = lpx_get_mat_col(lp, j-m, ndx, NULL);
441#endif
442            }
443         }
444      }
445      /* return the length of the row/column pattern */
446      return len;
447}
448
449static void adv_basis(glp_prob *lp)
450{     int m = lpx_get_num_rows(lp);
451      int n = lpx_get_num_cols(lp);
452      int i, j, jj, k, size;
453      int *rn, *cn, *rn_inv, *cn_inv;
454      int typx, *tagx = xcalloc(1+m+n, sizeof(int));
455      double lb, ub;
456      xprintf("Constructing initial basis...\n");
457#if 0 /* 13/V-2009 */
458      if (m == 0)
459         xerror("glp_adv_basis: problem has no rows\n");
460      if (n == 0)
461         xerror("glp_adv_basis: problem has no columns\n");
462#else
463      if (m == 0 || n == 0)
464      {  glp_std_basis(lp);
465         return;
466      }
467#endif
468      /* use the routine triang (see above) to find maximal triangular
469         part of the augmented constraint matrix A~ = (I|-A); in order
470         to prevent columns of fixed variables to be included in the
471         triangular part, such columns are implictly removed from the
472         matrix A~ by the routine adv_mat */
473      rn = xcalloc(1+m, sizeof(int));
474      cn = xcalloc(1+m+n, sizeof(int));
475      size = triang(m, m+n, lp, mat, rn, cn);
476      if (lpx_get_int_parm(lp, LPX_K_MSGLEV) >= 3)
477         xprintf("Size of triangular part = %d\n", size);
478      /* the first size rows and columns of the matrix P*A~*Q (where
479         P and Q are permutation matrices defined by the arrays rn and
480         cn) form a lower triangular matrix; build the arrays (rn_inv
481         and cn_inv), which define the matrices inv(P) and inv(Q) */
482      rn_inv = xcalloc(1+m, sizeof(int));
483      cn_inv = xcalloc(1+m+n, sizeof(int));
484      for (i = 1; i <= m; i++) rn_inv[rn[i]] = i;
485      for (j = 1; j <= m+n; j++) cn_inv[cn[j]] = j;
486      /* include the columns of the matrix A~, which correspond to the
487         first size columns of the matrix P*A~*Q, in the basis */
488      for (k = 1; k <= m+n; k++) tagx[k] = -1;
489      for (jj = 1; jj <= size; jj++)
490      {  j = cn_inv[jj];
491         /* the j-th column of A~ is the jj-th column of P*A~*Q */
492         tagx[j] = LPX_BS;
493      }
494      /* if size < m, we need to add appropriate columns of auxiliary
495         variables to the basis */
496      for (jj = size + 1; jj <= m; jj++)
497      {  /* the jj-th column of P*A~*Q should be replaced by the column
498            of the auxiliary variable, for which the only unity element
499            is placed in the position [jj,jj] */
500         i = rn_inv[jj];
501         /* the jj-th row of P*A~*Q is the i-th row of A~, but in the
502            i-th row of A~ the unity element belongs to the i-th column
503            of A~; therefore the disired column corresponds to the i-th
504            auxiliary variable (note that this column doesn't belong to
505            the triangular part found by the routine triang) */
506         xassert(1 <= i && i <= m);
507         xassert(cn[i] > size);
508         tagx[i] = LPX_BS;
509      }
510      /* free working arrays */
511      xfree(rn);
512      xfree(cn);
513      xfree(rn_inv);
514      xfree(cn_inv);
515      /* build tags of non-basic variables */
516      for (k = 1; k <= m+n; k++)
517      {  if (tagx[k] != LPX_BS)
518         {  if (k <= m)
519               lpx_get_row_bnds(lp, k, &typx, &lb, &ub);
520            else
521               lpx_get_col_bnds(lp, k-m, &typx, &lb, &ub);
522            switch (typx)
523            {  case LPX_FR:
524                  tagx[k] = LPX_NF; break;
525               case LPX_LO:
526                  tagx[k] = LPX_NL; break;
527               case LPX_UP:
528                  tagx[k] = LPX_NU; break;
529               case LPX_DB:
530                  tagx[k] =
531                     (fabs(lb) <= fabs(ub) ? LPX_NL : LPX_NU);
532                  break;
533               case LPX_FX:
534                  tagx[k] = LPX_NS; break;
535               default:
536                  xassert(typx != typx);
537            }
538         }
539      }
540      for (k = 1; k <= m+n; k++)
541      {  if (k <= m)
542            lpx_set_row_stat(lp, k, tagx[k]);
543         else
544            lpx_set_col_stat(lp, k-m, tagx[k]);
545      }
546      xfree(tagx);
547      return;
548}
549
550/***********************************************************************
551*  NAME
552*
553*  glp_adv_basis - construct advanced initial LP basis
554*
555*  SYNOPSIS
556*
557*  void glp_adv_basis(glp_prob *lp, int flags);
558*
559*  DESCRIPTION
560*
561*  The routine glp_adv_basis constructs an advanced initial basis for
562*  the specified problem object.
563*
564*  The parameter flags is reserved for use in the future and must be
565*  specified as zero. */
566
567void glp_adv_basis(glp_prob *lp, int flags)
568{     if (flags != 0)
569         xerror("glp_adv_basis: flags = %d; invalid flags\n", flags);
570      if (lp->m == 0 || lp->n == 0)
571         glp_std_basis(lp);
572      else
573         adv_basis(lp);
574      return;
575}
576
577/* eof */
Note: See TracBrowser for help on using the repository browser.