alpar@1: /* glpnet01.c (permutations for zero-free diagonal) */ alpar@1: alpar@1: /*********************************************************************** alpar@1: * This code is part of GLPK (GNU Linear Programming Kit). alpar@1: * alpar@1: * This code is the result of translation of the Fortran subroutines alpar@1: * MC21A and MC21B associated with the following paper: alpar@1: * alpar@1: * I.S.Duff, Algorithm 575: Permutations for zero-free diagonal, ACM alpar@1: * Trans. on Math. Softw. 7 (1981), 387-390. alpar@1: * alpar@1: * Use of ACM Algorithms is subject to the ACM Software Copyright and alpar@1: * License Agreement. See . alpar@1: * alpar@1: * The translation was made by Andrew Makhorin . alpar@1: * alpar@1: * GLPK is free software: you can redistribute it and/or modify it alpar@1: * under the terms of the GNU General Public License as published by alpar@1: * the Free Software Foundation, either version 3 of the License, or alpar@1: * (at your option) any later version. alpar@1: * alpar@1: * GLPK is distributed in the hope that it will be useful, but WITHOUT alpar@1: * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY alpar@1: * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public alpar@1: * License for more details. alpar@1: * alpar@1: * You should have received a copy of the GNU General Public License alpar@1: * along with GLPK. If not, see . alpar@1: ***********************************************************************/ alpar@1: alpar@1: #include "glpnet.h" alpar@1: alpar@1: /*********************************************************************** alpar@1: * NAME alpar@1: * alpar@1: * mc21a - permutations for zero-free diagonal alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpnet.h" alpar@1: * int mc21a(int n, const int icn[], const int ip[], const int lenr[], alpar@1: * int iperm[], int pr[], int arp[], int cv[], int out[]); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * Given the pattern of nonzeros of a sparse matrix, the routine mc21a alpar@1: * attempts to find a permutation of its rows that makes the matrix have alpar@1: * no zeros on its diagonal. alpar@1: * alpar@1: * INPUT PARAMETERS alpar@1: * alpar@1: * n order of matrix. alpar@1: * alpar@1: * icn array containing the column indices of the non-zeros. Those alpar@1: * belonging to a single row must be contiguous but the ordering alpar@1: * of column indices within each row is unimportant and wasted alpar@1: * space between rows is permitted. alpar@1: * alpar@1: * ip ip[i], i = 1,2,...,n, is the position in array icn of the alpar@1: * first column index of a non-zero in row i. alpar@1: * alpar@1: * lenr lenr[i], i = 1,2,...,n, is the number of non-zeros in row i. alpar@1: * alpar@1: * OUTPUT PARAMETER alpar@1: * alpar@1: * iperm contains permutation to make diagonal have the smallest alpar@1: * number of zeros on it. Elements (iperm[i], i), i = 1,2,...,n, alpar@1: * are non-zero at the end of the algorithm unless the matrix is alpar@1: * structurally singular. In this case, (iperm[i], i) will be alpar@1: * zero for n - numnz entries. alpar@1: * alpar@1: * WORKING ARRAYS alpar@1: * alpar@1: * pr working array of length [1+n], where pr[0] is not used. alpar@1: * pr[i] is the previous row to i in the depth first search. alpar@1: * alpar@1: * arp working array of length [1+n], where arp[0] is not used. alpar@1: * arp[i] is one less than the number of non-zeros in row i which alpar@1: * have not been scanned when looking for a cheap assignment. alpar@1: * alpar@1: * cv working array of length [1+n], where cv[0] is not used. alpar@1: * cv[i] is the most recent row extension at which column i was alpar@1: * visited. alpar@1: * alpar@1: * out working array of length [1+n], where out[0] is not used. alpar@1: * out[i] is one less than the number of non-zeros in row i alpar@1: * which have not been scanned during one pass through the main alpar@1: * loop. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine mc21a returns numnz, the number of non-zeros on diagonal alpar@1: * of permuted matrix. */ alpar@1: alpar@1: int mc21a(int n, const int icn[], const int ip[], const int lenr[], alpar@1: int iperm[], int pr[], int arp[], int cv[], int out[]) alpar@1: { int i, ii, in1, in2, j, j1, jord, k, kk, numnz; alpar@1: /* Initialization of arrays. */ alpar@1: for (i = 1; i <= n; i++) alpar@1: { arp[i] = lenr[i] - 1; alpar@1: cv[i] = iperm[i] = 0; alpar@1: } alpar@1: numnz = 0; alpar@1: /* Main loop. */ alpar@1: /* Each pass round this loop either results in a new assignment alpar@1: or gives a row with no assignment. */ alpar@1: for (jord = 1; jord <= n; jord++) alpar@1: { j = jord; alpar@1: pr[j] = -1; alpar@1: for (k = 1; k <= jord; k++) alpar@1: { /* Look for a cheap assignment. */ alpar@1: in1 = arp[j]; alpar@1: if (in1 >= 0) alpar@1: { in2 = ip[j] + lenr[j] - 1; alpar@1: in1 = in2 - in1; alpar@1: for (ii = in1; ii <= in2; ii++) alpar@1: { i = icn[ii]; alpar@1: if (iperm[i] == 0) goto L110; alpar@1: } alpar@1: /* No cheap assignment in row. */ alpar@1: arp[j] = -1; alpar@1: } alpar@1: /* Begin looking for assignment chain starting with row j.*/ alpar@1: out[j] = lenr[j] - 1; alpar@1: /* Inner loop. Extends chain by one or backtracks. */ alpar@1: for (kk = 1; kk <= jord; kk++) alpar@1: { in1 = out[j]; alpar@1: if (in1 >= 0) alpar@1: { in2 = ip[j] + lenr[j] - 1; alpar@1: in1 = in2 - in1; alpar@1: /* Forward scan. */ alpar@1: for (ii = in1; ii <= in2; ii++) alpar@1: { i = icn[ii]; alpar@1: if (cv[i] != jord) alpar@1: { /* Column i has not yet been accessed during alpar@1: this pass. */ alpar@1: j1 = j; alpar@1: j = iperm[i]; alpar@1: cv[i] = jord; alpar@1: pr[j] = j1; alpar@1: out[j1] = in2 - ii - 1; alpar@1: goto L100; alpar@1: } alpar@1: } alpar@1: } alpar@1: /* Backtracking step. */ alpar@1: j = pr[j]; alpar@1: if (j == -1) goto L130; alpar@1: } alpar@1: L100: ; alpar@1: } alpar@1: L110: /* New assignment is made. */ alpar@1: iperm[i] = j; alpar@1: arp[j] = in2 - ii - 1; alpar@1: numnz++; alpar@1: for (k = 1; k <= jord; k++) alpar@1: { j = pr[j]; alpar@1: if (j == -1) break; alpar@1: ii = ip[j] + lenr[j] - out[j] - 2; alpar@1: i = icn[ii]; alpar@1: iperm[i] = j; alpar@1: } alpar@1: L130: ; alpar@1: } alpar@1: /* If matrix is structurally singular, we now complete the alpar@1: permutation iperm. */ alpar@1: if (numnz < n) alpar@1: { for (i = 1; i <= n; i++) alpar@1: arp[i] = 0; alpar@1: k = 0; alpar@1: for (i = 1; i <= n; i++) alpar@1: { if (iperm[i] == 0) alpar@1: out[++k] = i; alpar@1: else alpar@1: arp[iperm[i]] = i; alpar@1: } alpar@1: k = 0; alpar@1: for (i = 1; i <= n; i++) alpar@1: { if (arp[i] == 0) alpar@1: iperm[out[++k]] = i; alpar@1: } alpar@1: } alpar@1: return numnz; alpar@1: } alpar@1: alpar@1: /**********************************************************************/ alpar@1: alpar@1: #if 0 alpar@1: #include "glplib.h" alpar@1: alpar@1: int sing; alpar@1: alpar@1: void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum, alpar@1: int iw[]); alpar@1: alpar@1: void fa01bs(int max, int *nrand); alpar@1: alpar@1: int main(void) alpar@1: { /* test program for the routine mc21a */ alpar@1: /* these runs on random matrices cause all possible statements in alpar@1: mc21a to be executed */ alpar@1: int i, iold, j, j1, j2, jj, knum, l, licn, n, nov4, num, numnz; alpar@1: int ip[1+21], icn[1+1000], iperm[1+20], lenr[1+20], iw1[1+80]; alpar@1: licn = 1000; alpar@1: /* run on random matrices of orders 1 through 20 */ alpar@1: for (n = 1; n <= 20; n++) alpar@1: { nov4 = n / 4; alpar@1: if (nov4 < 1) nov4 = 1; alpar@1: L10: fa01bs(nov4, &l); alpar@1: knum = l * n; alpar@1: /* knum is requested number of non-zeros in random matrix */ alpar@1: if (knum > licn) goto L10; alpar@1: /* if sing is false, matrix is guaranteed structurally alpar@1: non-singular */ alpar@1: sing = ((n / 2) * 2 == n); alpar@1: /* call to subroutine to generate random matrix */ alpar@1: ranmat(n, n, icn, ip, n+1, &knum, iw1); alpar@1: /* knum is now actual number of non-zeros in random matrix */ alpar@1: if (knum > licn) goto L10; alpar@1: xprintf("n = %2d; nz = %4d; sing = %d\n", n, knum, sing); alpar@1: /* set up array of row lengths */ alpar@1: for (i = 1; i <= n; i++) alpar@1: lenr[i] = ip[i+1] - ip[i]; alpar@1: /* call to mc21a */ alpar@1: numnz = mc21a(n, icn, ip, lenr, iperm, &iw1[0], &iw1[n], alpar@1: &iw1[n+n], &iw1[n+n+n]); alpar@1: /* testing to see if there are numnz non-zeros on the diagonal alpar@1: of the permuted matrix. */ alpar@1: num = 0; alpar@1: for (i = 1; i <= n; i++) alpar@1: { iold = iperm[i]; alpar@1: j1 = ip[iold]; alpar@1: j2 = j1 + lenr[iold] - 1; alpar@1: if (j2 < j1) continue; alpar@1: for (jj = j1; jj <= j2; jj++) alpar@1: { j = icn[jj]; alpar@1: if (j == i) alpar@1: { num++; alpar@1: break; alpar@1: } alpar@1: } alpar@1: } alpar@1: if (num != numnz) alpar@1: xprintf("Failure in mc21a, numnz = %d instead of %d\n", alpar@1: numnz, num); alpar@1: } alpar@1: return 0; alpar@1: } alpar@1: alpar@1: void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum, alpar@1: int iw[]) alpar@1: { /* subroutine to generate random matrix */ alpar@1: int i, ii, inum, j, lrow, matnum; alpar@1: inum = (*knum / n) * 2; alpar@1: if (inum > n-1) inum = n-1; alpar@1: matnum = 1; alpar@1: /* each pass through this loop generates a row of the matrix */ alpar@1: for (j = 1; j <= m; j++) alpar@1: { iptr[j] = matnum; alpar@1: if (!(sing || j > n)) alpar@1: icn[matnum++] = j; alpar@1: if (n == 1) continue; alpar@1: for (i = 1; i <= n; i++) iw[i] = 0; alpar@1: if (!sing) iw[j] = 1; alpar@1: fa01bs(inum, &lrow); alpar@1: lrow--; alpar@1: if (lrow == 0) continue; alpar@1: /* lrow off-diagonal non-zeros in row j of the matrix */ alpar@1: for (ii = 1; ii <= lrow; ii++) alpar@1: { for (;;) alpar@1: { fa01bs(n, &i); alpar@1: if (iw[i] != 1) break; alpar@1: } alpar@1: iw[i] = 1; alpar@1: icn[matnum++] = i; alpar@1: } alpar@1: } alpar@1: for (i = m+1; i <= nnnp1; i++) alpar@1: iptr[i] = matnum; alpar@1: *knum = matnum - 1; alpar@1: return; alpar@1: } alpar@1: alpar@1: double g = 1431655765.0; alpar@1: alpar@1: double fa01as(int i) alpar@1: { /* random number generator */ alpar@1: g = fmod(g * 9228907.0, 4294967296.0); alpar@1: if (i >= 0) alpar@1: return g / 4294967296.0; alpar@1: else alpar@1: return 2.0 * g / 4294967296.0 - 1.0; alpar@1: } alpar@1: alpar@1: void fa01bs(int max, int *nrand) alpar@1: { *nrand = (int)(fa01as(1) * (double)max) + 1; alpar@1: return; alpar@1: } alpar@1: #endif alpar@1: alpar@1: /* eof */