/* glpnet02.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 "glpnet.h"

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

*  NAME

*

*  mc13d - permutations to block triangular form

*

*  SYNOPSIS

*

*  #include "glpnet.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;

}

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

#if 0

#include "glplib.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 */