diff options
Diffstat (limited to 'test/monniaux/glpk-4.65/src/misc/mc13d.c')
| -rw-r--r-- | test/monniaux/glpk-4.65/src/misc/mc13d.c | 314 |
1 files changed, 314 insertions, 0 deletions
diff --git a/test/monniaux/glpk-4.65/src/misc/mc13d.c b/test/monniaux/glpk-4.65/src/misc/mc13d.c new file mode 100644 index 00000000..d8bab398 --- /dev/null +++ b/test/monniaux/glpk-4.65/src/misc/mc13d.c @@ -0,0 +1,314 @@ +/* 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 "env.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 */ |