/* glpnet01.c (permutations for zero-free diagonal) */

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

 *  This code is part of GLPK (GNU Linear Programming Kit).

 *

 *  This code is the result of translation of the Fortran subroutines

 *  MC21A and MC21B associated with the following paper:

 *

 *  I.S.Duff, Algorithm 575: Permutations for zero-free diagonal, ACM

 *  Trans. on Math. Softw. 7 (1981), 387-390.

 *

 *  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

 *

 *  mc21a - permutations for zero-free diagonal

 *

 *  SYNOPSIS

 *

 *  #include "glpnet.h"

 *  int mc21a(int n, const int icn[], const int ip[], const int lenr[],

 *     int iperm[], int pr[], int arp[], int cv[], int out[]);

 *

 *  DESCRIPTION

 *

 *  Given the pattern of nonzeros of a sparse matrix, the routine mc21a

 *  attempts to find a permutation of its rows that makes the matrix have

 *  no zeros on its diagonal.

 *

 *  INPUT PARAMETERS

 *

 *  n     order of 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 PARAMETER

 *

 *  iperm contains permutation to make diagonal have the smallest

 *        number of zeros on it. Elements (iperm[i], i), i = 1,2,...,n,

 *        are non-zero at the end of the algorithm unless the matrix is

 *        structurally singular. In this case, (iperm[i], i) will be

 *        zero for n - numnz entries.

 *

 *  WORKING ARRAYS

 *

 *  pr    working array of length [1+n], where pr[0] is not used.

 *        pr[i] is the previous row to i in the depth first search.

 *

 *  arp   working array of length [1+n], where arp[0] is not used.

 *        arp[i] is one less than the number of non-zeros in row i which

 *        have not been scanned when looking for a cheap assignment.

 *

 *  cv    working array of length [1+n], where cv[0] is not used.

 *        cv[i] is the most recent row extension at which column i was

 *        visited.

 *

 *  out   working array of length [1+n], where out[0] is not used.

 *        out[i] is one less than the number of non-zeros in row i

 *        which have not been scanned during one pass through the main

 *        loop.

 *

 *  RETURNS

 *

 *  The routine mc21a returns numnz, the number of non-zeros on diagonal

 *  of permuted matrix. */

 int mc21a(int n, const int icn[], const int ip[], const int lenr[],

       int iperm[], int pr[], int arp[], int cv[], int out[])

 {     int i, ii, in1, in2, j, j1, jord, k, kk, numnz;

       /* Initialization of arrays. */

       for (i = 1; i <= n; i++)

       {  arp[i] = lenr[i] - 1;

          cv[i] = iperm[i] = 0;

       }

       numnz = 0;

       /* Main loop. */

       /* Each pass round this loop either results in a new assignment

          or gives a row with no assignment. */

       for (jord = 1; jord <= n; jord++)

       {  j = jord;

          pr[j] = -1;

          for (k = 1; k <= jord; k++)

          {  /* Look for a cheap assignment. */

             in1 = arp[j];

             if (in1 >= 0)

             {  in2 = ip[j] + lenr[j] - 1;

                in1 = in2 - in1;

                for (ii = in1; ii <= in2; ii++)

                {  i = icn[ii];

                   if (iperm[i] == 0) goto L110;

                }

                /* No cheap assignment in row. */

                arp[j] = -1;

             }

             /* Begin looking for assignment chain starting with row j.*/

             out[j] = lenr[j] - 1;

             /* Inner loop. Extends chain by one or backtracks. */

             for (kk = 1; kk <= jord; kk++)

             {  in1 = out[j];

                if (in1 >= 0)

                {  in2 = ip[j] + lenr[j] - 1;

                   in1 = in2 - in1;

                   /* Forward scan. */

                   for (ii = in1; ii <= in2; ii++)

                   {  i = icn[ii];

                      if (cv[i] != jord)

                      {  /* Column i has not yet been accessed during

                            this pass. */

                         j1 = j;

                         j = iperm[i];

                         cv[i] = jord;

                         pr[j] = j1;

                         out[j1] = in2 - ii - 1;

                         goto L100;

                      }

                   }

                }

                /* Backtracking step. */

                j = pr[j];

                if (j == -1) goto L130;

             }

 L100:       ;

          }

 L110:    /* New assignment is made. */

          iperm[i] = j;

          arp[j] = in2 - ii - 1;

          numnz++;

          for (k = 1; k <= jord; k++)

          {  j = pr[j];

             if (j == -1) break;

             ii = ip[j] + lenr[j] - out[j] - 2;

             i = icn[ii];

             iperm[i] = j;

          }

 L130:    ;

       }

       /* If matrix is structurally singular, we now complete the

          permutation iperm. */

       if (numnz < n)

       {  for (i = 1; i <= n; i++)

             arp[i] = 0;

          k = 0;

          for (i = 1; i <= n; i++)

          {  if (iperm[i] == 0)

                out[++k] = i;

             else

                arp[iperm[i]] = i;

          }

          k = 0;

          for (i = 1; i <= n; i++)

          {  if (arp[i] == 0)

                iperm[out[++k]] = i;

          }

       }

       return numnz;

 }

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

 #if 0

 #include "glplib.h"

 int sing;

 void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum,

       int iw[]);

 void fa01bs(int max, int *nrand);

 int main(void)

 {     /* test program for the routine mc21a */

       /* these runs on random matrices cause all possible statements in

          mc21a to be executed */

       int i, iold, j, j1, j2, jj, knum, l, licn, n, nov4, num, numnz;

       int ip[1+21], icn[1+1000], iperm[1+20], lenr[1+20], iw1[1+80];

       licn = 1000;

       /* run on random matrices of orders 1 through 20 */

       for (n = 1; n <= 20; n++)

       {  nov4 = n / 4;

          if (nov4 < 1) nov4 = 1;

 L10:     fa01bs(nov4, &l);

          knum = l * n;

          /* knum is requested number of non-zeros in random matrix */

          if (knum > licn) goto L10;

          /* if sing is false, matrix is guaranteed structurally

             non-singular */

          sing = ((n / 2) * 2 == n);

          /* call to subroutine to generate random matrix */

          ranmat(n, n, icn, ip, n+1, &knum, iw1);

          /* knum is now actual number of non-zeros in random matrix */

          if (knum > licn) goto L10;

          xprintf("n = %2d; nz = %4d; sing = %d\n", n, knum, sing);

          /* set up array of row lengths */

          for (i = 1; i <= n; i++)

             lenr[i] = ip[i+1] - ip[i];

          /* call to mc21a */

          numnz = mc21a(n, icn, ip, lenr, iperm, &iw1[0], &iw1[n],

             &iw1[n+n], &iw1[n+n+n]);

          /* testing to see if there are numnz non-zeros on the diagonal

             of the permuted matrix. */

          num = 0;

          for (i = 1; i <= n; i++)

          {  iold = iperm[i];

             j1 = ip[iold];

             j2 = j1 + lenr[iold] - 1;

             if (j2 < j1) continue;

             for (jj = j1; jj <= j2; jj++)

             {  j = icn[jj];

                if (j == i)

                {  num++;

                   break;

                }

             }

          }

          if (num != numnz)

             xprintf("Failure in mc21a, numnz = %d instead of %d\n",

                numnz, num);

       }

       return 0;

 }

 void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum,

       int iw[])

 {     /* subroutine to generate random matrix */

       int i, ii, inum, j, lrow, matnum;

       inum = (*knum / n) * 2;

       if (inum > n-1) inum = n-1;

       matnum = 1;

       /* each pass through this loop generates a row of the matrix */

       for (j = 1; j <= m; j++)

       {  iptr[j] = matnum;

          if (!(sing || j > n))

             icn[matnum++] = j;

          if (n == 1) continue;

          for (i = 1; i <= n; i++) iw[i] = 0;

          if (!sing) iw[j] = 1;

          fa01bs(inum, &lrow);

          lrow--;

          if (lrow == 0) continue;

          /* lrow off-diagonal non-zeros in row j of the matrix */

          for (ii = 1; ii <= lrow; ii++)

          {  for (;;)

             {  fa01bs(n, &i);

                if (iw[i] != 1) break;

             }

             iw[i] = 1;

             icn[matnum++] = i;

          }

       }

       for (i = m+1; i <= nnnp1; i++)

          iptr[i] = matnum;

       *knum = matnum - 1;

       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 */