1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpnet02.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,314 @@
1.4 +/* glpnet02.c (permutations to block triangular form) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* This code is the result of translation of the Fortran subroutines
1.10 +* MC13D and MC13E associated with the following paper:
1.11 +*
1.12 +* I.S.Duff, J.K.Reid, Algorithm 529: Permutations to block triangular
1.13 +* form, ACM Trans. on Math. Softw. 4 (1978), 189-192.
1.14 +*
1.15 +* Use of ACM Algorithms is subject to the ACM Software Copyright and
1.16 +* License Agreement. See <http://www.acm.org/publications/policies>.
1.17 +*
1.18 +* The translation was made by Andrew Makhorin <mao@gnu.org>.
1.19 +*
1.20 +* GLPK is free software: you can redistribute it and/or modify it
1.21 +* under the terms of the GNU General Public License as published by
1.22 +* the Free Software Foundation, either version 3 of the License, or
1.23 +* (at your option) any later version.
1.24 +*
1.25 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.26 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.27 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.28 +* License for more details.
1.29 +*
1.30 +* You should have received a copy of the GNU General Public License
1.31 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.32 +***********************************************************************/
1.33 +
1.34 +#include "glpnet.h"
1.35 +
1.36 +/***********************************************************************
1.37 +* NAME
1.38 +*
1.39 +* mc13d - permutations to block triangular form
1.40 +*
1.41 +* SYNOPSIS
1.42 +*
1.43 +* #include "glpnet.h"
1.44 +* int mc13d(int n, const int icn[], const int ip[], const int lenr[],
1.45 +* int ior[], int ib[], int lowl[], int numb[], int prev[]);
1.46 +*
1.47 +* DESCRIPTION
1.48 +*
1.49 +* Given the column numbers of the nonzeros in each row of the sparse
1.50 +* matrix, the routine mc13d finds a symmetric permutation that makes
1.51 +* the matrix block lower triangular.
1.52 +*
1.53 +* INPUT PARAMETERS
1.54 +*
1.55 +* n order of the matrix.
1.56 +*
1.57 +* icn array containing the column indices of the non-zeros. Those
1.58 +* belonging to a single row must be contiguous but the ordering
1.59 +* of column indices within each row is unimportant and wasted
1.60 +* space between rows is permitted.
1.61 +*
1.62 +* ip ip[i], i = 1,2,...,n, is the position in array icn of the
1.63 +* first column index of a non-zero in row i.
1.64 +*
1.65 +* lenr lenr[i], i = 1,2,...,n, is the number of non-zeros in row i.
1.66 +*
1.67 +* OUTPUT PARAMETERS
1.68 +*
1.69 +* ior ior[i], i = 1,2,...,n, gives the position on the original
1.70 +* ordering of the row or column which is in position i in the
1.71 +* permuted form.
1.72 +*
1.73 +* ib ib[i], i = 1,2,...,num, is the row number in the permuted
1.74 +* matrix of the beginning of block i, 1 <= num <= n.
1.75 +*
1.76 +* WORKING ARRAYS
1.77 +*
1.78 +* arp working array of length [1+n], where arp[0] is not used.
1.79 +* arp[i] is one less than the number of unsearched edges leaving
1.80 +* node i. At the end of the algorithm it is set to a permutation
1.81 +* which puts the matrix in block lower triangular form.
1.82 +*
1.83 +* ib working array of length [1+n], where ib[0] is not used.
1.84 +* ib[i] is the position in the ordering of the start of the ith
1.85 +* block. ib[n+1-i] holds the node number of the ith node on the
1.86 +* stack.
1.87 +*
1.88 +* lowl working array of length [1+n], where lowl[0] is not used.
1.89 +* lowl[i] is the smallest stack position of any node to which a
1.90 +* path from node i has been found. It is set to n+1 when node i
1.91 +* is removed from the stack.
1.92 +*
1.93 +* numb working array of length [1+n], where numb[0] is not used.
1.94 +* numb[i] is the position of node i in the stack if it is on it,
1.95 +* is the permuted order of node i for those nodes whose final
1.96 +* position has been found and is otherwise zero.
1.97 +*
1.98 +* prev working array of length [1+n], where prev[0] is not used.
1.99 +* prev[i] is the node at the end of the path when node i was
1.100 +* placed on the stack.
1.101 +*
1.102 +* RETURNS
1.103 +*
1.104 +* The routine mc13d returns num, the number of blocks found. */
1.105 +
1.106 +int mc13d(int n, const int icn[], const int ip[], const int lenr[],
1.107 + int ior[], int ib[], int lowl[], int numb[], int prev[])
1.108 +{ int *arp = ior;
1.109 + int dummy, i, i1, i2, icnt, ii, isn, ist, ist1, iv, iw, j, lcnt,
1.110 + nnm1, num, stp;
1.111 + /* icnt is the number of nodes whose positions in final ordering
1.112 + have been found. */
1.113 + icnt = 0;
1.114 + /* num is the number of blocks that have been found. */
1.115 + num = 0;
1.116 + nnm1 = n + n - 1;
1.117 + /* Initialization of arrays. */
1.118 + for (j = 1; j <= n; j++)
1.119 + { numb[j] = 0;
1.120 + arp[j] = lenr[j] - 1;
1.121 + }
1.122 + for (isn = 1; isn <= n; isn++)
1.123 + { /* Look for a starting node. */
1.124 + if (numb[isn] != 0) continue;
1.125 + iv = isn;
1.126 + /* ist is the number of nodes on the stack ... it is the stack
1.127 + pointer. */
1.128 + ist = 1;
1.129 + /* Put node iv at beginning of stack. */
1.130 + lowl[iv] = numb[iv] = 1;
1.131 + ib[n] = iv;
1.132 + /* The body of this loop puts a new node on the stack or
1.133 + backtracks. */
1.134 + for (dummy = 1; dummy <= nnm1; dummy++)
1.135 + { i1 = arp[iv];
1.136 + /* Have all edges leaving node iv been searched? */
1.137 + if (i1 >= 0)
1.138 + { i2 = ip[iv] + lenr[iv] - 1;
1.139 + i1 = i2 - i1;
1.140 + /* Look at edges leaving node iv until one enters a new
1.141 + node or all edges are exhausted. */
1.142 + for (ii = i1; ii <= i2; ii++)
1.143 + { iw = icn[ii];
1.144 + /* Has node iw been on stack already? */
1.145 + if (numb[iw] == 0) goto L70;
1.146 + /* Update value of lowl[iv] if necessary. */
1.147 + if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw];
1.148 + }
1.149 + /* There are no more edges leaving node iv. */
1.150 + arp[iv] = -1;
1.151 + }
1.152 + /* Is node iv the root of a block? */
1.153 + if (lowl[iv] < numb[iv]) goto L60;
1.154 + /* Order nodes in a block. */
1.155 + num++;
1.156 + ist1 = n + 1 - ist;
1.157 + lcnt = icnt + 1;
1.158 + /* Peel block off the top of the stack starting at the top
1.159 + and working down to the root of the block. */
1.160 + for (stp = ist1; stp <= n; stp++)
1.161 + { iw = ib[stp];
1.162 + lowl[iw] = n + 1;
1.163 + numb[iw] = ++icnt;
1.164 + if (iw == iv) break;
1.165 + }
1.166 + ist = n - stp;
1.167 + ib[num] = lcnt;
1.168 + /* Are there any nodes left on the stack? */
1.169 + if (ist != 0) goto L60;
1.170 + /* Have all the nodes been ordered? */
1.171 + if (icnt < n) break;
1.172 + goto L100;
1.173 +L60: /* Backtrack to previous node on path. */
1.174 + iw = iv;
1.175 + iv = prev[iv];
1.176 + /* Update value of lowl[iv] if necessary. */
1.177 + if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw];
1.178 + continue;
1.179 +L70: /* Put new node on the stack. */
1.180 + arp[iv] = i2 - ii - 1;
1.181 + prev[iw] = iv;
1.182 + iv = iw;
1.183 + lowl[iv] = numb[iv] = ++ist;
1.184 + ib[n+1-ist] = iv;
1.185 + }
1.186 + }
1.187 +L100: /* Put permutation in the required form. */
1.188 + for (i = 1; i <= n; i++)
1.189 + arp[numb[i]] = i;
1.190 + return num;
1.191 +}
1.192 +
1.193 +/**********************************************************************/
1.194 +
1.195 +#if 0
1.196 +#include "glplib.h"
1.197 +
1.198 +void test(int n, int ipp);
1.199 +
1.200 +int main(void)
1.201 +{ /* test program for routine mc13d */
1.202 + test( 1, 0);
1.203 + test( 2, 1);
1.204 + test( 2, 2);
1.205 + test( 3, 3);
1.206 + test( 4, 4);
1.207 + test( 5, 10);
1.208 + test(10, 10);
1.209 + test(10, 20);
1.210 + test(20, 20);
1.211 + test(20, 50);
1.212 + test(50, 50);
1.213 + test(50, 200);
1.214 + return 0;
1.215 +}
1.216 +
1.217 +void fa01bs(int max, int *nrand);
1.218 +
1.219 +void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[]);
1.220 +
1.221 +void test(int n, int ipp)
1.222 +{ int ip[1+50], icn[1+1000], ior[1+50], ib[1+51], iw[1+150],
1.223 + lenr[1+50];
1.224 + char a[1+50][1+50], hold[1+100];
1.225 + int i, ii, iblock, ij, index, j, jblock, jj, k9, num;
1.226 + xprintf("\n\n\nMatrix is of order %d and has %d off-diagonal non-"
1.227 + "zeros\n", n, ipp);
1.228 + for (j = 1; j <= n; j++)
1.229 + { for (i = 1; i <= n; i++)
1.230 + a[i][j] = 0;
1.231 + a[j][j] = 1;
1.232 + }
1.233 + for (k9 = 1; k9 <= ipp; k9++)
1.234 + { /* these statements should be replaced by calls to your
1.235 + favorite random number generator to place two pseudo-random
1.236 + numbers between 1 and n in the variables i and j */
1.237 + for (;;)
1.238 + { fa01bs(n, &i);
1.239 + fa01bs(n, &j);
1.240 + if (!a[i][j]) break;
1.241 + }
1.242 + a[i][j] = 1;
1.243 + }
1.244 + /* setup converts matrix a[i,j] to required sparsity-oriented
1.245 + storage format */
1.246 + setup(n, a, ip, icn, lenr);
1.247 + num = mc13d(n, icn, ip, lenr, ior, ib, &iw[0], &iw[n], &iw[n+n]);
1.248 + /* output reordered matrix with blocking to improve clarity */
1.249 + xprintf("\nThe reordered matrix which has %d block%s is of the fo"
1.250 + "rm\n", num, num == 1 ? "" : "s");
1.251 + ib[num+1] = n + 1;
1.252 + index = 100;
1.253 + iblock = 1;
1.254 + for (i = 1; i <= n; i++)
1.255 + { for (ij = 1; ij <= index; ij++)
1.256 + hold[ij] = ' ';
1.257 + if (i == ib[iblock])
1.258 + { xprintf("\n");
1.259 + iblock++;
1.260 + }
1.261 + jblock = 1;
1.262 + index = 0;
1.263 + for (j = 1; j <= n; j++)
1.264 + { if (j == ib[jblock])
1.265 + { hold[++index] = ' ';
1.266 + jblock++;
1.267 + }
1.268 + ii = ior[i];
1.269 + jj = ior[j];
1.270 + hold[++index] = (char)(a[ii][jj] ? 'X' : '0');
1.271 + }
1.272 + xprintf("%.*s\n", index, &hold[1]);
1.273 + }
1.274 + xprintf("\nThe starting point for each block is given by\n");
1.275 + for (i = 1; i <= num; i++)
1.276 + { if ((i - 1) % 12 == 0) xprintf("\n");
1.277 + xprintf(" %4d", ib[i]);
1.278 + }
1.279 + xprintf("\n");
1.280 + return;
1.281 +}
1.282 +
1.283 +void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[])
1.284 +{ int i, j, ind;
1.285 + for (i = 1; i <= n; i++)
1.286 + lenr[i] = 0;
1.287 + ind = 1;
1.288 + for (i = 1; i <= n; i++)
1.289 + { ip[i] = ind;
1.290 + for (j = 1; j <= n; j++)
1.291 + { if (a[i][j])
1.292 + { lenr[i]++;
1.293 + icn[ind++] = j;
1.294 + }
1.295 + }
1.296 + }
1.297 + return;
1.298 +}
1.299 +
1.300 +double g = 1431655765.0;
1.301 +
1.302 +double fa01as(int i)
1.303 +{ /* random number generator */
1.304 + g = fmod(g * 9228907.0, 4294967296.0);
1.305 + if (i >= 0)
1.306 + return g / 4294967296.0;
1.307 + else
1.308 + return 2.0 * g / 4294967296.0 - 1.0;
1.309 +}
1.310 +
1.311 +void fa01bs(int max, int *nrand)
1.312 +{ *nrand = (int)(fa01as(1) * (double)max) + 1;
1.313 + return;
1.314 +}
1.315 +#endif
1.316 +
1.317 +/* eof */