lemon-project-template-glpk

annotate deps/glpk/src/glpnet02.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
rev   line source
alpar@9 1 /* glpnet02.c (permutations to block triangular form) */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * This code is the result of translation of the Fortran subroutines
alpar@9 7 * MC13D and MC13E associated with the following paper:
alpar@9 8 *
alpar@9 9 * I.S.Duff, J.K.Reid, Algorithm 529: Permutations to block triangular
alpar@9 10 * form, ACM Trans. on Math. Softw. 4 (1978), 189-192.
alpar@9 11 *
alpar@9 12 * Use of ACM Algorithms is subject to the ACM Software Copyright and
alpar@9 13 * License Agreement. See <http://www.acm.org/publications/policies>.
alpar@9 14 *
alpar@9 15 * The translation was made by Andrew Makhorin <mao@gnu.org>.
alpar@9 16 *
alpar@9 17 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 18 * under the terms of the GNU General Public License as published by
alpar@9 19 * the Free Software Foundation, either version 3 of the License, or
alpar@9 20 * (at your option) any later version.
alpar@9 21 *
alpar@9 22 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 23 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 24 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 25 * License for more details.
alpar@9 26 *
alpar@9 27 * You should have received a copy of the GNU General Public License
alpar@9 28 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 29 ***********************************************************************/
alpar@9 30
alpar@9 31 #include "glpnet.h"
alpar@9 32
alpar@9 33 /***********************************************************************
alpar@9 34 * NAME
alpar@9 35 *
alpar@9 36 * mc13d - permutations to block triangular form
alpar@9 37 *
alpar@9 38 * SYNOPSIS
alpar@9 39 *
alpar@9 40 * #include "glpnet.h"
alpar@9 41 * int mc13d(int n, const int icn[], const int ip[], const int lenr[],
alpar@9 42 * int ior[], int ib[], int lowl[], int numb[], int prev[]);
alpar@9 43 *
alpar@9 44 * DESCRIPTION
alpar@9 45 *
alpar@9 46 * Given the column numbers of the nonzeros in each row of the sparse
alpar@9 47 * matrix, the routine mc13d finds a symmetric permutation that makes
alpar@9 48 * the matrix block lower triangular.
alpar@9 49 *
alpar@9 50 * INPUT PARAMETERS
alpar@9 51 *
alpar@9 52 * n order of the matrix.
alpar@9 53 *
alpar@9 54 * icn array containing the column indices of the non-zeros. Those
alpar@9 55 * belonging to a single row must be contiguous but the ordering
alpar@9 56 * of column indices within each row is unimportant and wasted
alpar@9 57 * space between rows is permitted.
alpar@9 58 *
alpar@9 59 * ip ip[i], i = 1,2,...,n, is the position in array icn of the
alpar@9 60 * first column index of a non-zero in row i.
alpar@9 61 *
alpar@9 62 * lenr lenr[i], i = 1,2,...,n, is the number of non-zeros in row i.
alpar@9 63 *
alpar@9 64 * OUTPUT PARAMETERS
alpar@9 65 *
alpar@9 66 * ior ior[i], i = 1,2,...,n, gives the position on the original
alpar@9 67 * ordering of the row or column which is in position i in the
alpar@9 68 * permuted form.
alpar@9 69 *
alpar@9 70 * ib ib[i], i = 1,2,...,num, is the row number in the permuted
alpar@9 71 * matrix of the beginning of block i, 1 <= num <= n.
alpar@9 72 *
alpar@9 73 * WORKING ARRAYS
alpar@9 74 *
alpar@9 75 * arp working array of length [1+n], where arp[0] is not used.
alpar@9 76 * arp[i] is one less than the number of unsearched edges leaving
alpar@9 77 * node i. At the end of the algorithm it is set to a permutation
alpar@9 78 * which puts the matrix in block lower triangular form.
alpar@9 79 *
alpar@9 80 * ib working array of length [1+n], where ib[0] is not used.
alpar@9 81 * ib[i] is the position in the ordering of the start of the ith
alpar@9 82 * block. ib[n+1-i] holds the node number of the ith node on the
alpar@9 83 * stack.
alpar@9 84 *
alpar@9 85 * lowl working array of length [1+n], where lowl[0] is not used.
alpar@9 86 * lowl[i] is the smallest stack position of any node to which a
alpar@9 87 * path from node i has been found. It is set to n+1 when node i
alpar@9 88 * is removed from the stack.
alpar@9 89 *
alpar@9 90 * numb working array of length [1+n], where numb[0] is not used.
alpar@9 91 * numb[i] is the position of node i in the stack if it is on it,
alpar@9 92 * is the permuted order of node i for those nodes whose final
alpar@9 93 * position has been found and is otherwise zero.
alpar@9 94 *
alpar@9 95 * prev working array of length [1+n], where prev[0] is not used.
alpar@9 96 * prev[i] is the node at the end of the path when node i was
alpar@9 97 * placed on the stack.
alpar@9 98 *
alpar@9 99 * RETURNS
alpar@9 100 *
alpar@9 101 * The routine mc13d returns num, the number of blocks found. */
alpar@9 102
alpar@9 103 int mc13d(int n, const int icn[], const int ip[], const int lenr[],
alpar@9 104 int ior[], int ib[], int lowl[], int numb[], int prev[])
alpar@9 105 { int *arp = ior;
alpar@9 106 int dummy, i, i1, i2, icnt, ii, isn, ist, ist1, iv, iw, j, lcnt,
alpar@9 107 nnm1, num, stp;
alpar@9 108 /* icnt is the number of nodes whose positions in final ordering
alpar@9 109 have been found. */
alpar@9 110 icnt = 0;
alpar@9 111 /* num is the number of blocks that have been found. */
alpar@9 112 num = 0;
alpar@9 113 nnm1 = n + n - 1;
alpar@9 114 /* Initialization of arrays. */
alpar@9 115 for (j = 1; j <= n; j++)
alpar@9 116 { numb[j] = 0;
alpar@9 117 arp[j] = lenr[j] - 1;
alpar@9 118 }
alpar@9 119 for (isn = 1; isn <= n; isn++)
alpar@9 120 { /* Look for a starting node. */
alpar@9 121 if (numb[isn] != 0) continue;
alpar@9 122 iv = isn;
alpar@9 123 /* ist is the number of nodes on the stack ... it is the stack
alpar@9 124 pointer. */
alpar@9 125 ist = 1;
alpar@9 126 /* Put node iv at beginning of stack. */
alpar@9 127 lowl[iv] = numb[iv] = 1;
alpar@9 128 ib[n] = iv;
alpar@9 129 /* The body of this loop puts a new node on the stack or
alpar@9 130 backtracks. */
alpar@9 131 for (dummy = 1; dummy <= nnm1; dummy++)
alpar@9 132 { i1 = arp[iv];
alpar@9 133 /* Have all edges leaving node iv been searched? */
alpar@9 134 if (i1 >= 0)
alpar@9 135 { i2 = ip[iv] + lenr[iv] - 1;
alpar@9 136 i1 = i2 - i1;
alpar@9 137 /* Look at edges leaving node iv until one enters a new
alpar@9 138 node or all edges are exhausted. */
alpar@9 139 for (ii = i1; ii <= i2; ii++)
alpar@9 140 { iw = icn[ii];
alpar@9 141 /* Has node iw been on stack already? */
alpar@9 142 if (numb[iw] == 0) goto L70;
alpar@9 143 /* Update value of lowl[iv] if necessary. */
alpar@9 144 if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw];
alpar@9 145 }
alpar@9 146 /* There are no more edges leaving node iv. */
alpar@9 147 arp[iv] = -1;
alpar@9 148 }
alpar@9 149 /* Is node iv the root of a block? */
alpar@9 150 if (lowl[iv] < numb[iv]) goto L60;
alpar@9 151 /* Order nodes in a block. */
alpar@9 152 num++;
alpar@9 153 ist1 = n + 1 - ist;
alpar@9 154 lcnt = icnt + 1;
alpar@9 155 /* Peel block off the top of the stack starting at the top
alpar@9 156 and working down to the root of the block. */
alpar@9 157 for (stp = ist1; stp <= n; stp++)
alpar@9 158 { iw = ib[stp];
alpar@9 159 lowl[iw] = n + 1;
alpar@9 160 numb[iw] = ++icnt;
alpar@9 161 if (iw == iv) break;
alpar@9 162 }
alpar@9 163 ist = n - stp;
alpar@9 164 ib[num] = lcnt;
alpar@9 165 /* Are there any nodes left on the stack? */
alpar@9 166 if (ist != 0) goto L60;
alpar@9 167 /* Have all the nodes been ordered? */
alpar@9 168 if (icnt < n) break;
alpar@9 169 goto L100;
alpar@9 170 L60: /* Backtrack to previous node on path. */
alpar@9 171 iw = iv;
alpar@9 172 iv = prev[iv];
alpar@9 173 /* Update value of lowl[iv] if necessary. */
alpar@9 174 if (lowl[iw] < lowl[iv]) lowl[iv] = lowl[iw];
alpar@9 175 continue;
alpar@9 176 L70: /* Put new node on the stack. */
alpar@9 177 arp[iv] = i2 - ii - 1;
alpar@9 178 prev[iw] = iv;
alpar@9 179 iv = iw;
alpar@9 180 lowl[iv] = numb[iv] = ++ist;
alpar@9 181 ib[n+1-ist] = iv;
alpar@9 182 }
alpar@9 183 }
alpar@9 184 L100: /* Put permutation in the required form. */
alpar@9 185 for (i = 1; i <= n; i++)
alpar@9 186 arp[numb[i]] = i;
alpar@9 187 return num;
alpar@9 188 }
alpar@9 189
alpar@9 190 /**********************************************************************/
alpar@9 191
alpar@9 192 #if 0
alpar@9 193 #include "glplib.h"
alpar@9 194
alpar@9 195 void test(int n, int ipp);
alpar@9 196
alpar@9 197 int main(void)
alpar@9 198 { /* test program for routine mc13d */
alpar@9 199 test( 1, 0);
alpar@9 200 test( 2, 1);
alpar@9 201 test( 2, 2);
alpar@9 202 test( 3, 3);
alpar@9 203 test( 4, 4);
alpar@9 204 test( 5, 10);
alpar@9 205 test(10, 10);
alpar@9 206 test(10, 20);
alpar@9 207 test(20, 20);
alpar@9 208 test(20, 50);
alpar@9 209 test(50, 50);
alpar@9 210 test(50, 200);
alpar@9 211 return 0;
alpar@9 212 }
alpar@9 213
alpar@9 214 void fa01bs(int max, int *nrand);
alpar@9 215
alpar@9 216 void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[]);
alpar@9 217
alpar@9 218 void test(int n, int ipp)
alpar@9 219 { int ip[1+50], icn[1+1000], ior[1+50], ib[1+51], iw[1+150],
alpar@9 220 lenr[1+50];
alpar@9 221 char a[1+50][1+50], hold[1+100];
alpar@9 222 int i, ii, iblock, ij, index, j, jblock, jj, k9, num;
alpar@9 223 xprintf("\n\n\nMatrix is of order %d and has %d off-diagonal non-"
alpar@9 224 "zeros\n", n, ipp);
alpar@9 225 for (j = 1; j <= n; j++)
alpar@9 226 { for (i = 1; i <= n; i++)
alpar@9 227 a[i][j] = 0;
alpar@9 228 a[j][j] = 1;
alpar@9 229 }
alpar@9 230 for (k9 = 1; k9 <= ipp; k9++)
alpar@9 231 { /* these statements should be replaced by calls to your
alpar@9 232 favorite random number generator to place two pseudo-random
alpar@9 233 numbers between 1 and n in the variables i and j */
alpar@9 234 for (;;)
alpar@9 235 { fa01bs(n, &i);
alpar@9 236 fa01bs(n, &j);
alpar@9 237 if (!a[i][j]) break;
alpar@9 238 }
alpar@9 239 a[i][j] = 1;
alpar@9 240 }
alpar@9 241 /* setup converts matrix a[i,j] to required sparsity-oriented
alpar@9 242 storage format */
alpar@9 243 setup(n, a, ip, icn, lenr);
alpar@9 244 num = mc13d(n, icn, ip, lenr, ior, ib, &iw[0], &iw[n], &iw[n+n]);
alpar@9 245 /* output reordered matrix with blocking to improve clarity */
alpar@9 246 xprintf("\nThe reordered matrix which has %d block%s is of the fo"
alpar@9 247 "rm\n", num, num == 1 ? "" : "s");
alpar@9 248 ib[num+1] = n + 1;
alpar@9 249 index = 100;
alpar@9 250 iblock = 1;
alpar@9 251 for (i = 1; i <= n; i++)
alpar@9 252 { for (ij = 1; ij <= index; ij++)
alpar@9 253 hold[ij] = ' ';
alpar@9 254 if (i == ib[iblock])
alpar@9 255 { xprintf("\n");
alpar@9 256 iblock++;
alpar@9 257 }
alpar@9 258 jblock = 1;
alpar@9 259 index = 0;
alpar@9 260 for (j = 1; j <= n; j++)
alpar@9 261 { if (j == ib[jblock])
alpar@9 262 { hold[++index] = ' ';
alpar@9 263 jblock++;
alpar@9 264 }
alpar@9 265 ii = ior[i];
alpar@9 266 jj = ior[j];
alpar@9 267 hold[++index] = (char)(a[ii][jj] ? 'X' : '0');
alpar@9 268 }
alpar@9 269 xprintf("%.*s\n", index, &hold[1]);
alpar@9 270 }
alpar@9 271 xprintf("\nThe starting point for each block is given by\n");
alpar@9 272 for (i = 1; i <= num; i++)
alpar@9 273 { if ((i - 1) % 12 == 0) xprintf("\n");
alpar@9 274 xprintf(" %4d", ib[i]);
alpar@9 275 }
alpar@9 276 xprintf("\n");
alpar@9 277 return;
alpar@9 278 }
alpar@9 279
alpar@9 280 void setup(int n, char a[1+50][1+50], int ip[], int icn[], int lenr[])
alpar@9 281 { int i, j, ind;
alpar@9 282 for (i = 1; i <= n; i++)
alpar@9 283 lenr[i] = 0;
alpar@9 284 ind = 1;
alpar@9 285 for (i = 1; i <= n; i++)
alpar@9 286 { ip[i] = ind;
alpar@9 287 for (j = 1; j <= n; j++)
alpar@9 288 { if (a[i][j])
alpar@9 289 { lenr[i]++;
alpar@9 290 icn[ind++] = j;
alpar@9 291 }
alpar@9 292 }
alpar@9 293 }
alpar@9 294 return;
alpar@9 295 }
alpar@9 296
alpar@9 297 double g = 1431655765.0;
alpar@9 298
alpar@9 299 double fa01as(int i)
alpar@9 300 { /* random number generator */
alpar@9 301 g = fmod(g * 9228907.0, 4294967296.0);
alpar@9 302 if (i >= 0)
alpar@9 303 return g / 4294967296.0;
alpar@9 304 else
alpar@9 305 return 2.0 * g / 4294967296.0 - 1.0;
alpar@9 306 }
alpar@9 307
alpar@9 308 void fa01bs(int max, int *nrand)
alpar@9 309 { *nrand = (int)(fa01as(1) * (double)max) + 1;
alpar@9 310 return;
alpar@9 311 }
alpar@9 312 #endif
alpar@9 313
alpar@9 314 /* eof */