alpar@1: /* glpnet02.c (permutations to block triangular form) */ 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: * MC13D and MC13E associated with the following paper: alpar@1: * alpar@1: * I.S.Duff, J.K.Reid, Algorithm 529: Permutations to block triangular alpar@1: * form, ACM Trans. on Math. Softw. 4 (1978), 189-192. 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: * mc13d - permutations to block triangular form alpar@1: * alpar@1: * SYNOPSIS alpar@1: * alpar@1: * #include "glpnet.h" alpar@1: * int mc13d(int n, const int icn[], const int ip[], const int lenr[], alpar@1: * int ior[], int ib[], int lowl[], int numb[], int prev[]); alpar@1: * alpar@1: * DESCRIPTION alpar@1: * alpar@1: * Given the column numbers of the nonzeros in each row of the sparse alpar@1: * matrix, the routine mc13d finds a symmetric permutation that makes alpar@1: * the matrix block lower triangular. alpar@1: * alpar@1: * INPUT PARAMETERS alpar@1: * alpar@1: * n order of the 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 PARAMETERS alpar@1: * alpar@1: * ior ior[i], i = 1,2,...,n, gives the position on the original alpar@1: * ordering of the row or column which is in position i in the alpar@1: * permuted form. alpar@1: * alpar@1: * ib ib[i], i = 1,2,...,num, is the row number in the permuted alpar@1: * matrix of the beginning of block i, 1 <= num <= n. alpar@1: * alpar@1: * WORKING ARRAYS 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 unsearched edges leaving alpar@1: * node i. At the end of the algorithm it is set to a permutation alpar@1: * which puts the matrix in block lower triangular form. alpar@1: * alpar@1: * ib working array of length [1+n], where ib[0] is not used. alpar@1: * ib[i] is the position in the ordering of the start of the ith alpar@1: * block. ib[n+1-i] holds the node number of the ith node on the alpar@1: * stack. alpar@1: * alpar@1: * lowl working array of length [1+n], where lowl[0] is not used. alpar@1: * lowl[i] is the smallest stack position of any node to which a alpar@1: * path from node i has been found. It is set to n+1 when node i alpar@1: * is removed from the stack. alpar@1: * alpar@1: * numb working array of length [1+n], where numb[0] is not used. alpar@1: * numb[i] is the position of node i in the stack if it is on it, alpar@1: * is the permuted order of node i for those nodes whose final alpar@1: * position has been found and is otherwise zero. alpar@1: * alpar@1: * prev working array of length [1+n], where prev[0] is not used. alpar@1: * prev[i] is the node at the end of the path when node i was alpar@1: * placed on the stack. alpar@1: * alpar@1: * RETURNS alpar@1: * alpar@1: * The routine mc13d returns num, the number of blocks found. */ alpar@1: alpar@1: int mc13d(int n, const int icn[], const int ip[], const int lenr[], alpar@1: int ior[], int ib[], int lowl[], int numb[], int prev[]) alpar@1: { int *arp = ior; alpar@1: int dummy, i, i1, i2, icnt, ii, isn, ist, ist1, iv, iw, j, lcnt, alpar@1: nnm1, num, stp; alpar@1: /* icnt is the number of nodes whose positions in final ordering alpar@1: have been found. */ alpar@1: icnt = 0; alpar@1: /* num is the number of blocks that have been found. */ alpar@1: num = 0; alpar@1: nnm1 = n + n - 1; alpar@1: /* Initialization of arrays. */ alpar@1: for (j = 1; j <= n; j++) alpar@1: { numb[j] = 0; alpar@1: arp[j] = lenr[j] - 1; alpar@1: } alpar@1: for (isn = 1; isn <= n; isn++) alpar@1: { /* Look for a starting node. */ alpar@1: if (numb[isn] != 0) continue; alpar@1: iv = isn; alpar@1: /* ist is the number of nodes on the stack ... it is the stack alpar@1: pointer. */ alpar@1: ist = 1; alpar@1: /* Put node iv at beginning of stack. */ alpar@1: lowl[iv] = numb[iv] = 1; alpar@1: ib[n] = iv; alpar@1: /* The body of this loop puts a new node on the stack or alpar@1: backtracks. */ alpar@1: for (dummy = 1; dummy <= nnm1; dummy++) alpar@1: { i1 = arp[iv]; alpar@1: /* Have all edges leaving node iv been searched? */ alpar@1: if (i1 >= 0) alpar@1: { i2 = ip[iv] + lenr[iv] - 1; alpar@1: i1 = i2 - i1; alpar@1: /* Look at edges leaving node iv until one enters a new alpar@1: node or all edges are exhausted. */ alpar@1: for (ii = i1; ii <= i2; ii++) alpar@1: { iw = icn[ii]; alpar@1: /* Has node iw been on stack already? */ alpar@1: if (numb[iw] == 0) goto L70; alpar@1: /* Update value of lowl[iv] if necessary. */ alpar@1: if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw]; alpar@1: } alpar@1: /* There are no more edges leaving node iv. */ alpar@1: arp[iv] = -1; alpar@1: } alpar@1: /* Is node iv the root of a block? */ alpar@1: if (lowl[iv] < numb[iv]) goto L60; alpar@1: /* Order nodes in a block. */ alpar@1: num++; alpar@1: ist1 = n + 1 - ist; alpar@1: lcnt = icnt + 1; alpar@1: /* Peel block off the top of the stack starting at the top alpar@1: and working down to the root of the block. */ alpar@1: for (stp = ist1; stp <= n; stp++) alpar@1: { iw = ib[stp]; alpar@1: lowl[iw] = n + 1; alpar@1: numb[iw] = ++icnt; alpar@1: if (iw == iv) break; alpar@1: } alpar@1: ist = n - stp; alpar@1: ib[num] = lcnt; alpar@1: /* Are there any nodes left on the stack? */ alpar@1: if (ist != 0) goto L60; alpar@1: /* Have all the nodes been ordered? */ alpar@1: if (icnt < n) break; alpar@1: goto L100; alpar@1: L60: /* Backtrack to previous node on path. */ alpar@1: iw = iv; alpar@1: iv = prev[iv]; alpar@1: /* Update value of lowl[iv] if necessary. */ alpar@1: if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw]; alpar@1: continue; alpar@1: L70: /* Put new node on the stack. */ alpar@1: arp[iv] = i2 - ii - 1; alpar@1: prev[iw] = iv; alpar@1: iv = iw; alpar@1: lowl[iv] = numb[iv] = ++ist; alpar@1: ib[n+1-ist] = iv; alpar@1: } alpar@1: } alpar@1: L100: /* Put permutation in the required form. */ alpar@1: for (i = 1; i <= n; i++) alpar@1: arp[numb[i]] = i; alpar@1: return num; alpar@1: } alpar@1: alpar@1: /**********************************************************************/ alpar@1: alpar@1: #if 0 alpar@1: #include "glplib.h" alpar@1: alpar@1: void test(int n, int ipp); alpar@1: alpar@1: int main(void) alpar@1: { /* test program for routine mc13d */ alpar@1: test( 1, 0); alpar@1: test( 2, 1); alpar@1: test( 2, 2); alpar@1: test( 3, 3); alpar@1: test( 4, 4); alpar@1: test( 5, 10); alpar@1: test(10, 10); alpar@1: test(10, 20); alpar@1: test(20, 20); alpar@1: test(20, 50); alpar@1: test(50, 50); alpar@1: test(50, 200); alpar@1: return 0; alpar@1: } alpar@1: alpar@1: void fa01bs(int max, int *nrand); alpar@1: alpar@1: void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[]); alpar@1: alpar@1: void test(int n, int ipp) alpar@1: { int ip[1+50], icn[1+1000], ior[1+50], ib[1+51], iw[1+150], alpar@1: lenr[1+50]; alpar@1: char a[1+50][1+50], hold[1+100]; alpar@1: int i, ii, iblock, ij, index, j, jblock, jj, k9, num; alpar@1: xprintf("\n\n\nMatrix is of order %d and has %d off-diagonal non-" alpar@1: "zeros\n", n, ipp); alpar@1: for (j = 1; j <= n; j++) alpar@1: { for (i = 1; i <= n; i++) alpar@1: a[i][j] = 0; alpar@1: a[j][j] = 1; alpar@1: } alpar@1: for (k9 = 1; k9 <= ipp; k9++) alpar@1: { /* these statements should be replaced by calls to your alpar@1: favorite random number generator to place two pseudo-random alpar@1: numbers between 1 and n in the variables i and j */ alpar@1: for (;;) alpar@1: { fa01bs(n, &i); alpar@1: fa01bs(n, &j); alpar@1: if (!a[i][j]) break; alpar@1: } alpar@1: a[i][j] = 1; alpar@1: } alpar@1: /* setup converts matrix a[i,j] to required sparsity-oriented alpar@1: storage format */ alpar@1: setup(n, a, ip, icn, lenr); alpar@1: num = mc13d(n, icn, ip, lenr, ior, ib, &iw[0], &iw[n], &iw[n+n]); alpar@1: /* output reordered matrix with blocking to improve clarity */ alpar@1: xprintf("\nThe reordered matrix which has %d block%s is of the fo" alpar@1: "rm\n", num, num == 1 ? "" : "s"); alpar@1: ib[num+1] = n + 1; alpar@1: index = 100; alpar@1: iblock = 1; alpar@1: for (i = 1; i <= n; i++) alpar@1: { for (ij = 1; ij <= index; ij++) alpar@1: hold[ij] = ' '; alpar@1: if (i == ib[iblock]) alpar@1: { xprintf("\n"); alpar@1: iblock++; alpar@1: } alpar@1: jblock = 1; alpar@1: index = 0; alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (j == ib[jblock]) alpar@1: { hold[++index] = ' '; alpar@1: jblock++; alpar@1: } alpar@1: ii = ior[i]; alpar@1: jj = ior[j]; alpar@1: hold[++index] = (char)(a[ii][jj] ? 'X' : '0'); alpar@1: } alpar@1: xprintf("%.*s\n", index, &hold[1]); alpar@1: } alpar@1: xprintf("\nThe starting point for each block is given by\n"); alpar@1: for (i = 1; i <= num; i++) alpar@1: { if ((i - 1) % 12 == 0) xprintf("\n"); alpar@1: xprintf(" %4d", ib[i]); alpar@1: } alpar@1: xprintf("\n"); alpar@1: return; alpar@1: } alpar@1: alpar@1: void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[]) alpar@1: { int i, j, ind; alpar@1: for (i = 1; i <= n; i++) alpar@1: lenr[i] = 0; alpar@1: ind = 1; alpar@1: for (i = 1; i <= n; i++) alpar@1: { ip[i] = ind; alpar@1: for (j = 1; j <= n; j++) alpar@1: { if (a[i][j]) alpar@1: { lenr[i]++; alpar@1: icn[ind++] = j; alpar@1: } alpar@1: } alpar@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 */