/* mc13d.c (permutations to block triangular form) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  This code is the result of translation of the Fortran subroutines
*  MC13D and MC13E associated with the following paper:
*
*  I.S.Duff, J.K.Reid, Algorithm 529: Permutations to block triangular
*  form, ACM Trans. on Math. Softw. 4 (1978), 189-192.
*
*  Use of ACM Algorithms is subject to the ACM Software Copyright and
*  License Agreement. See <http://www.acm.org/publications/policies>.
*
*  The translation was made by Andrew Makhorin <mao@gnu.org>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "mc13d.h"

/***********************************************************************
*  NAME
*
*  mc13d - permutations to block triangular form
*
*  SYNOPSIS
*
*  #include "mc13d.h"
*  int mc13d(int n, const int icn[], const int ip[], const int lenr[],
*     int ior[], int ib[], int lowl[], int numb[], int prev[]);
*
*  DESCRIPTION
*
*  Given the column numbers of the nonzeros in each row of the sparse
*  matrix, the routine mc13d finds a symmetric permutation that makes
*  the matrix block lower triangular.
*
*  INPUT PARAMETERS
*
*  n     order of the matrix.
*
*  icn   array containing the column indices of the non-zeros. Those
*        belonging to a single row must be contiguous but the ordering
*        of column indices within each row is unimportant and wasted
*        space between rows is permitted.
*
*  ip    ip[i], i = 1,2,...,n, is the position in array icn of the
*        first column index of a non-zero in row i.
*
*  lenr  lenr[i], i = 1,2,...,n, is the number of non-zeros in row i.
*
*  OUTPUT PARAMETERS
*
*  ior   ior[i], i = 1,2,...,n, gives the position on the original
*        ordering of the row or column which is in position i in the
*        permuted form.
*
*  ib    ib[i], i = 1,2,...,num, is the row number in the permuted
*        matrix of the beginning of block i, 1 <= num <= n.
*
*  WORKING ARRAYS
*
*  arp   working array of length [1+n], where arp[0] is not used.
*        arp[i] is one less than the number of unsearched edges leaving
*        node i. At the end of the algorithm it is set to a permutation
*        which puts the matrix in block lower triangular form.
*
*  ib    working array of length [1+n], where ib[0] is not used.
*        ib[i] is the position in the ordering of the start of the ith
*        block. ib[n+1-i] holds the node number of the ith node on the
*        stack.
*
*  lowl  working array of length [1+n], where lowl[0] is not used.
*        lowl[i] is the smallest stack position of any node to which a
*        path from node i has been found. It is set to n+1 when node i
*        is removed from the stack.
*
*  numb  working array of length [1+n], where numb[0] is not used.
*        numb[i] is the position of node i in the stack if it is on it,
*        is the permuted order of node i for those nodes whose final
*        position has been found and is otherwise zero.
*
*  prev  working array of length [1+n], where prev[0] is not used.
*        prev[i] is the node at the end of the path when node i was
*        placed on the stack.
*
*  RETURNS
*
*  The routine mc13d returns num, the number of blocks found. */

int mc13d(int n, const int icn[], const int ip[], const int lenr[],
      int ior[], int ib[], int lowl[], int numb[], int prev[])
{     int *arp = ior;
      int dummy, i, i1, i2, icnt, ii, isn, ist, ist1, iv, iw, j, lcnt,
         nnm1, num, stp;
      /* icnt is the number of nodes whose positions in final ordering
       * have been found. */
      icnt = 0;
      /* num is the number of blocks that have been found. */
      num = 0;
      nnm1 = n + n - 1;
      /* Initialization of arrays. */
      for (j = 1; j <= n; j++)
      {  numb[j] = 0;
         arp[j] = lenr[j] - 1;
      }
      for (isn = 1; isn <= n; isn++)
      {  /* Look for a starting node. */
         if (numb[isn] != 0) continue;
         iv = isn;
         /* ist is the number of nodes on the stack ... it is the stack
          * pointer. */
         ist = 1;
         /* Put node iv at beginning of stack. */
         lowl[iv] = numb[iv] = 1;
         ib[n] = iv;
         /* The body of this loop puts a new node on the stack or
          * backtracks. */
         for (dummy = 1; dummy <= nnm1; dummy++)
         {  i1 = arp[iv];
            /* Have all edges leaving node iv been searched? */
            if (i1 >= 0)
            {  i2 = ip[iv] + lenr[iv] - 1;
               i1 = i2 - i1;
               /* Look at edges leaving node iv until one enters a new
                * node or all edges are exhausted. */
               for (ii = i1; ii <= i2; ii++)
               {  iw = icn[ii];
                  /* Has node iw been on stack already? */
                  if (numb[iw] == 0) goto L70;
                  /* Update value of lowl[iv] if necessary. */
                  if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw];
               }
               /* There are no more edges leaving node iv. */
               arp[iv] = -1;
            }
            /* Is node iv the root of a block? */
            if (lowl[iv] < numb[iv]) goto L60;
            /* Order nodes in a block. */
            num++;
            ist1 = n + 1 - ist;
            lcnt = icnt + 1;
            /* Peel block off the top of the stack starting at the top
             * and working down to the root of the block. */
            for (stp = ist1; stp <= n; stp++)
            {  iw = ib[stp];
               lowl[iw] = n + 1;
               numb[iw] = ++icnt;
               if (iw == iv) break;
            }
            ist = n - stp;
            ib[num] = lcnt;
            /* Are there any nodes left on the stack? */
            if (ist != 0) goto L60;
            /* Have all the nodes been ordered? */
            if (icnt < n) break;
            goto L100;
L60:        /* Backtrack to previous node on path. */
            iw = iv;
            iv = prev[iv];
            /* Update value of lowl[iv] if necessary. */
            if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw];
            continue;
L70:        /* Put new node on the stack. */
            arp[iv] = i2 - ii - 1;
            prev[iw] = iv;
            iv = iw;
            lowl[iv] = numb[iv] = ++ist;
            ib[n+1-ist] = iv;
         }
      }
L100: /* Put permutation in the required form. */
      for (i = 1; i <= n; i++)
         arp[numb[i]] = i;
      return num;
}

/**********************************************************************/

#ifdef GLP_TEST
#include "glpenv.h"

void test(int n, int ipp);

int main(void)
{     /* test program for routine mc13d */
      test( 1,   0);
      test( 2,   1);
      test( 2,   2);
      test( 3,   3);
      test( 4,   4);
      test( 5,  10);
      test(10,  10);
      test(10,  20);
      test(20,  20);
      test(20,  50);
      test(50,  50);
      test(50, 200);
      return 0;
}

void fa01bs(int max, int *nrand);

void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[]);

void test(int n, int ipp)
{     int ip[1+50], icn[1+1000], ior[1+50], ib[1+51], iw[1+150],
         lenr[1+50];
      char a[1+50][1+50], hold[1+100];
      int i, ii, iblock, ij, index, j, jblock, jj, k9, num;
      xprintf("\n\n\nMatrix is of order %d and has %d off-diagonal non-"
         "zeros\n", n, ipp);
      for (j = 1; j <= n; j++)
      {  for (i = 1; i <= n; i++)
            a[i][j] = 0;
         a[j][j] = 1;
      }
      for (k9 = 1; k9 <= ipp; k9++)
      {  /* these statements should be replaced by calls to your
          * favorite random number generator to place two pseudo-random
          * numbers between 1 and n in the variables i and j */
         for (;;)
         {  fa01bs(n, &i);
            fa01bs(n, &j);
            if (!a[i][j]) break;
         }
         a[i][j] = 1;
      }
      /* setup converts matrix a[i,j] to required sparsity-oriented
       * storage format */
      setup(n, a, ip, icn, lenr);
      num = mc13d(n, icn, ip, lenr, ior, ib, &iw[0], &iw[n], &iw[n+n]);
      /* output reordered matrix with blocking to improve clarity */
      xprintf("\nThe reordered matrix which has %d block%s is of the fo"
         "rm\n", num, num == 1 ? "" : "s");
      ib[num+1] = n + 1;
      index = 100;
      iblock = 1;
      for (i = 1; i <= n; i++)
      {  for (ij = 1; ij <= index; ij++)
            hold[ij] = ' ';
         if (i == ib[iblock])
         {  xprintf("\n");
            iblock++;
         }
         jblock = 1;
         index = 0;
         for (j = 1; j <= n; j++)
         {  if (j == ib[jblock])
            {  hold[++index] = ' ';
               jblock++;
            }
            ii = ior[i];
            jj = ior[j];
            hold[++index] = (char)(a[ii][jj] ? 'X' : '0');
         }
         xprintf("%.*s\n", index, &hold[1]);
      }
      xprintf("\nThe starting point for each block is given by\n");
      for (i = 1; i <= num; i++)
      {  if ((i - 1) % 12 == 0) xprintf("\n");
         xprintf(" %4d", ib[i]);
      }
      xprintf("\n");
      return;
}

void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[])
{     int i, j, ind;
      for (i = 1; i <= n; i++)
         lenr[i] = 0;
      ind = 1;
      for (i = 1; i <= n; i++)
      {  ip[i] = ind;
         for (j = 1; j <= n; j++)
         {  if (a[i][j])
            {  lenr[i]++;
               icn[ind++] = j;
            }
         }
      }
      return;
}

double g = 1431655765.0;

double fa01as(int i)
{     /* random number generator */
      g = fmod(g * 9228907.0, 4294967296.0);
      if (i >= 0)
         return g / 4294967296.0;
      else
         return 2.0 * g / 4294967296.0 - 1.0;
}

void fa01bs(int max, int *nrand)
{     *nrand = (int)(fa01as(1) * (double)max) + 1;
      return;
}
#endif

/* eof */