1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpnet02.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,314 @@
     1.4 +/* glpnet02.c (permutations to block triangular form) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  This code is the result of translation of the Fortran subroutines
    1.10 +*  MC13D and MC13E associated with the following paper:
    1.11 +*
    1.12 +*  I.S.Duff, J.K.Reid, Algorithm 529: Permutations to block triangular
    1.13 +*  form, ACM Trans. on Math. Softw. 4 (1978), 189-192.
    1.14 +*
    1.15 +*  Use of ACM Algorithms is subject to the ACM Software Copyright and
    1.16 +*  License Agreement. See <http://www.acm.org/publications/policies>.
    1.17 +*
    1.18 +*  The translation was made by Andrew Makhorin <mao@gnu.org>.
    1.19 +*
    1.20 +*  GLPK is free software: you can redistribute it and/or modify it
    1.21 +*  under the terms of the GNU General Public License as published by
    1.22 +*  the Free Software Foundation, either version 3 of the License, or
    1.23 +*  (at your option) any later version.
    1.24 +*
    1.25 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.26 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.27 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.28 +*  License for more details.
    1.29 +*
    1.30 +*  You should have received a copy of the GNU General Public License
    1.31 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.32 +***********************************************************************/
    1.33 +
    1.34 +#include "glpnet.h"
    1.35 +
    1.36 +/***********************************************************************
    1.37 +*  NAME
    1.38 +*
    1.39 +*  mc13d - permutations to block triangular form
    1.40 +*
    1.41 +*  SYNOPSIS
    1.42 +*
    1.43 +*  #include "glpnet.h"
    1.44 +*  int mc13d(int n, const int icn[], const int ip[], const int lenr[],
    1.45 +*     int ior[], int ib[], int lowl[], int numb[], int prev[]);
    1.46 +*
    1.47 +*  DESCRIPTION
    1.48 +*
    1.49 +*  Given the column numbers of the nonzeros in each row of the sparse
    1.50 +*  matrix, the routine mc13d finds a symmetric permutation that makes
    1.51 +*  the matrix block lower triangular.
    1.52 +*
    1.53 +*  INPUT PARAMETERS
    1.54 +*
    1.55 +*  n     order of the matrix.
    1.56 +*
    1.57 +*  icn   array containing the column indices of the non-zeros. Those
    1.58 +*        belonging to a single row must be contiguous but the ordering
    1.59 +*        of column indices within each row is unimportant and wasted
    1.60 +*        space between rows is permitted.
    1.61 +*
    1.62 +*  ip    ip[i], i = 1,2,...,n, is the position in array icn of the
    1.63 +*        first column index of a non-zero in row i.
    1.64 +*
    1.65 +*  lenr  lenr[i], i = 1,2,...,n, is the number of non-zeros in row i.
    1.66 +*
    1.67 +*  OUTPUT PARAMETERS
    1.68 +*
    1.69 +*  ior   ior[i], i = 1,2,...,n, gives the position on the original
    1.70 +*        ordering of the row or column which is in position i in the
    1.71 +*        permuted form.
    1.72 +*
    1.73 +*  ib    ib[i], i = 1,2,...,num, is the row number in the permuted
    1.74 +*        matrix of the beginning of block i, 1 <= num <= n.
    1.75 +*
    1.76 +*  WORKING ARRAYS
    1.77 +*
    1.78 +*  arp   working array of length [1+n], where arp[0] is not used.
    1.79 +*        arp[i] is one less than the number of unsearched edges leaving
    1.80 +*        node i. At the end of the algorithm it is set to a permutation
    1.81 +*        which puts the matrix in block lower triangular form.
    1.82 +*
    1.83 +*  ib    working array of length [1+n], where ib[0] is not used.
    1.84 +*        ib[i] is the position in the ordering of the start of the ith
    1.85 +*        block. ib[n+1-i] holds the node number of the ith node on the
    1.86 +*        stack.
    1.87 +*
    1.88 +*  lowl  working array of length [1+n], where lowl[0] is not used.
    1.89 +*        lowl[i] is the smallest stack position of any node to which a
    1.90 +*        path from node i has been found. It is set to n+1 when node i
    1.91 +*        is removed from the stack.
    1.92 +*
    1.93 +*  numb  working array of length [1+n], where numb[0] is not used.
    1.94 +*        numb[i] is the position of node i in the stack if it is on it,
    1.95 +*        is the permuted order of node i for those nodes whose final
    1.96 +*        position has been found and is otherwise zero.
    1.97 +*
    1.98 +*  prev  working array of length [1+n], where prev[0] is not used.
    1.99 +*        prev[i] is the node at the end of the path when node i was
   1.100 +*        placed on the stack.
   1.101 +*
   1.102 +*  RETURNS
   1.103 +*
   1.104 +*  The routine mc13d returns num, the number of blocks found. */
   1.105 +
   1.106 +int mc13d(int n, const int icn[], const int ip[], const int lenr[],
   1.107 +      int ior[], int ib[], int lowl[], int numb[], int prev[])
   1.108 +{     int *arp = ior;
   1.109 +      int dummy, i, i1, i2, icnt, ii, isn, ist, ist1, iv, iw, j, lcnt,
   1.110 +         nnm1, num, stp;
   1.111 +      /* icnt is the number of nodes whose positions in final ordering
   1.112 +         have been found. */
   1.113 +      icnt = 0;
   1.114 +      /* num is the number of blocks that have been found. */
   1.115 +      num = 0;
   1.116 +      nnm1 = n + n - 1;
   1.117 +      /* Initialization of arrays. */
   1.118 +      for (j = 1; j <= n; j++)
   1.119 +      {  numb[j] = 0;
   1.120 +         arp[j] = lenr[j] - 1;
   1.121 +      }
   1.122 +      for (isn = 1; isn <= n; isn++)
   1.123 +      {  /* Look for a starting node. */
   1.124 +         if (numb[isn] != 0) continue;
   1.125 +         iv = isn;
   1.126 +         /* ist is the number of nodes on the stack ... it is the stack
   1.127 +            pointer. */
   1.128 +         ist = 1;
   1.129 +         /* Put node iv at beginning of stack. */
   1.130 +         lowl[iv] = numb[iv] = 1;
   1.131 +         ib[n] = iv;
   1.132 +         /* The body of this loop puts a new node on the stack or
   1.133 +            backtracks. */
   1.134 +         for (dummy = 1; dummy <= nnm1; dummy++)
   1.135 +         {  i1 = arp[iv];
   1.136 +            /* Have all edges leaving node iv been searched? */
   1.137 +            if (i1 >= 0)
   1.138 +            {  i2 = ip[iv] + lenr[iv] - 1;
   1.139 +               i1 = i2 - i1;
   1.140 +               /* Look at edges leaving node iv until one enters a new
   1.141 +                  node or all edges are exhausted. */
   1.142 +               for (ii = i1; ii <= i2; ii++)
   1.143 +               {  iw = icn[ii];
   1.144 +                  /* Has node iw been on stack already? */
   1.145 +                  if (numb[iw] == 0) goto L70;
   1.146 +                  /* Update value of lowl[iv] if necessary. */
   1.147 +                  if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw];
   1.148 +               }
   1.149 +               /* There are no more edges leaving node iv. */
   1.150 +               arp[iv] = -1;
   1.151 +            }
   1.152 +            /* Is node iv the root of a block? */
   1.153 +            if (lowl[iv] < numb[iv]) goto L60;
   1.154 +            /* Order nodes in a block. */
   1.155 +            num++;
   1.156 +            ist1 = n + 1 - ist;
   1.157 +            lcnt = icnt + 1;
   1.158 +            /* Peel block off the top of the stack starting at the top
   1.159 +               and working down to the root of the block. */
   1.160 +            for (stp = ist1; stp <= n; stp++)
   1.161 +            {  iw = ib[stp];
   1.162 +               lowl[iw] = n + 1;
   1.163 +               numb[iw] = ++icnt;
   1.164 +               if (iw == iv) break;
   1.165 +            }
   1.166 +            ist = n - stp;
   1.167 +            ib[num] = lcnt;
   1.168 +            /* Are there any nodes left on the stack? */
   1.169 +            if (ist != 0) goto L60;
   1.170 +            /* Have all the nodes been ordered? */
   1.171 +            if (icnt < n) break;
   1.172 +            goto L100;
   1.173 +L60:        /* Backtrack to previous node on path. */
   1.174 +            iw = iv;
   1.175 +            iv = prev[iv];
   1.176 +            /* Update value of lowl[iv] if necessary. */
   1.177 +            if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw];
   1.178 +            continue;
   1.179 +L70:        /* Put new node on the stack. */
   1.180 +            arp[iv] = i2 - ii - 1;
   1.181 +            prev[iw] = iv;
   1.182 +            iv = iw;
   1.183 +            lowl[iv] = numb[iv] = ++ist;
   1.184 +            ib[n+1-ist] = iv;
   1.185 +         }
   1.186 +      }
   1.187 +L100: /* Put permutation in the required form. */
   1.188 +      for (i = 1; i <= n; i++)
   1.189 +         arp[numb[i]] = i;
   1.190 +      return num;
   1.191 +}
   1.192 +
   1.193 +/**********************************************************************/
   1.194 +
   1.195 +#if 0
   1.196 +#include "glplib.h"
   1.197 +
   1.198 +void test(int n, int ipp);
   1.199 +
   1.200 +int main(void)
   1.201 +{     /* test program for routine mc13d */
   1.202 +      test( 1,   0);
   1.203 +      test( 2,   1);
   1.204 +      test( 2,   2);
   1.205 +      test( 3,   3);
   1.206 +      test( 4,   4);
   1.207 +      test( 5,  10);
   1.208 +      test(10,  10);
   1.209 +      test(10,  20);
   1.210 +      test(20,  20);
   1.211 +      test(20,  50);
   1.212 +      test(50,  50);
   1.213 +      test(50, 200);
   1.214 +      return 0;
   1.215 +}
   1.216 +
   1.217 +void fa01bs(int max, int *nrand);
   1.218 +
   1.219 +void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[]);
   1.220 +
   1.221 +void test(int n, int ipp)
   1.222 +{     int ip[1+50], icn[1+1000], ior[1+50], ib[1+51], iw[1+150],
   1.223 +         lenr[1+50];
   1.224 +      char a[1+50][1+50], hold[1+100];
   1.225 +      int i, ii, iblock, ij, index, j, jblock, jj, k9, num;
   1.226 +      xprintf("\n\n\nMatrix is of order %d and has %d off-diagonal non-"
   1.227 +         "zeros\n", n, ipp);
   1.228 +      for (j = 1; j <= n; j++)
   1.229 +      {  for (i = 1; i <= n; i++)
   1.230 +            a[i][j] = 0;
   1.231 +         a[j][j] = 1;
   1.232 +      }
   1.233 +      for (k9 = 1; k9 <= ipp; k9++)
   1.234 +      {  /* these statements should be replaced by calls to your
   1.235 +            favorite random number generator to place two pseudo-random
   1.236 +            numbers between 1 and n in the variables i and j */
   1.237 +         for (;;)
   1.238 +         {  fa01bs(n, &i);
   1.239 +            fa01bs(n, &j);
   1.240 +            if (!a[i][j]) break;
   1.241 +         }
   1.242 +         a[i][j] = 1;
   1.243 +      }
   1.244 +      /* setup converts matrix a[i,j] to required sparsity-oriented
   1.245 +         storage format */
   1.246 +      setup(n, a, ip, icn, lenr);
   1.247 +      num = mc13d(n, icn, ip, lenr, ior, ib, &iw[0], &iw[n], &iw[n+n]);
   1.248 +      /* output reordered matrix with blocking to improve clarity */
   1.249 +      xprintf("\nThe reordered matrix which has %d block%s is of the fo"
   1.250 +         "rm\n", num, num == 1 ? "" : "s");
   1.251 +      ib[num+1] = n + 1;
   1.252 +      index = 100;
   1.253 +      iblock = 1;
   1.254 +      for (i = 1; i <= n; i++)
   1.255 +      {  for (ij = 1; ij <= index; ij++)
   1.256 +            hold[ij] = ' ';
   1.257 +         if (i == ib[iblock])
   1.258 +         {  xprintf("\n");
   1.259 +            iblock++;
   1.260 +         }
   1.261 +         jblock = 1;
   1.262 +         index = 0;
   1.263 +         for (j = 1; j <= n; j++)
   1.264 +         {  if (j == ib[jblock])
   1.265 +            {  hold[++index] = ' ';
   1.266 +               jblock++;
   1.267 +            }
   1.268 +            ii = ior[i];
   1.269 +            jj = ior[j];
   1.270 +            hold[++index] = (char)(a[ii][jj] ? 'X' : '0');
   1.271 +         }
   1.272 +         xprintf("%.*s\n", index, &hold[1]);
   1.273 +      }
   1.274 +      xprintf("\nThe starting point for each block is given by\n");
   1.275 +      for (i = 1; i <= num; i++)
   1.276 +      {  if ((i - 1) % 12 == 0) xprintf("\n");
   1.277 +         xprintf(" %4d", ib[i]);
   1.278 +      }
   1.279 +      xprintf("\n");
   1.280 +      return;
   1.281 +}
   1.282 +
   1.283 +void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[])
   1.284 +{     int i, j, ind;
   1.285 +      for (i = 1; i <= n; i++)
   1.286 +         lenr[i] = 0;
   1.287 +      ind = 1;
   1.288 +      for (i = 1; i <= n; i++)
   1.289 +      {  ip[i] = ind;
   1.290 +         for (j = 1; j <= n; j++)
   1.291 +         {  if (a[i][j])
   1.292 +            {  lenr[i]++;
   1.293 +               icn[ind++] = j;
   1.294 +            }
   1.295 +         }
   1.296 +      }
   1.297 +      return;
   1.298 +}
   1.299 +
   1.300 +double g = 1431655765.0;
   1.301 +
   1.302 +double fa01as(int i)
   1.303 +{     /* random number generator */
   1.304 +      g = fmod(g * 9228907.0, 4294967296.0);
   1.305 +      if (i >= 0)
   1.306 +         return g / 4294967296.0;
   1.307 +      else
   1.308 +         return 2.0 * g / 4294967296.0 - 1.0;
   1.309 +}
   1.310 +
   1.311 +void fa01bs(int max, int *nrand)
   1.312 +{     *nrand = (int)(fa01as(1) * (double)max) + 1;
   1.313 +      return;
   1.314 +}
   1.315 +#endif
   1.316 +
   1.317 +/* eof */