lemon-project-template-glpk

annotate deps/glpk/src/glpnet01.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 /* glpnet01.c (permutations for zero-free diagonal) */
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 * MC21A and MC21B associated with the following paper:
alpar@9 8 *
alpar@9 9 * I.S.Duff, Algorithm 575: Permutations for zero-free diagonal, ACM
alpar@9 10 * Trans. on Math. Softw. 7 (1981), 387-390.
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 * mc21a - permutations for zero-free diagonal
alpar@9 37 *
alpar@9 38 * SYNOPSIS
alpar@9 39 *
alpar@9 40 * #include "glpnet.h"
alpar@9 41 * int mc21a(int n, const int icn[], const int ip[], const int lenr[],
alpar@9 42 * int iperm[], int pr[], int arp[], int cv[], int out[]);
alpar@9 43 *
alpar@9 44 * DESCRIPTION
alpar@9 45 *
alpar@9 46 * Given the pattern of nonzeros of a sparse matrix, the routine mc21a
alpar@9 47 * attempts to find a permutation of its rows that makes the matrix have
alpar@9 48 * no zeros on its diagonal.
alpar@9 49 *
alpar@9 50 * INPUT PARAMETERS
alpar@9 51 *
alpar@9 52 * n order of 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 PARAMETER
alpar@9 65 *
alpar@9 66 * iperm contains permutation to make diagonal have the smallest
alpar@9 67 * number of zeros on it. Elements (iperm[i], i), i = 1,2,...,n,
alpar@9 68 * are non-zero at the end of the algorithm unless the matrix is
alpar@9 69 * structurally singular. In this case, (iperm[i], i) will be
alpar@9 70 * zero for n - numnz entries.
alpar@9 71 *
alpar@9 72 * WORKING ARRAYS
alpar@9 73 *
alpar@9 74 * pr working array of length [1+n], where pr[0] is not used.
alpar@9 75 * pr[i] is the previous row to i in the depth first search.
alpar@9 76 *
alpar@9 77 * arp working array of length [1+n], where arp[0] is not used.
alpar@9 78 * arp[i] is one less than the number of non-zeros in row i which
alpar@9 79 * have not been scanned when looking for a cheap assignment.
alpar@9 80 *
alpar@9 81 * cv working array of length [1+n], where cv[0] is not used.
alpar@9 82 * cv[i] is the most recent row extension at which column i was
alpar@9 83 * visited.
alpar@9 84 *
alpar@9 85 * out working array of length [1+n], where out[0] is not used.
alpar@9 86 * out[i] is one less than the number of non-zeros in row i
alpar@9 87 * which have not been scanned during one pass through the main
alpar@9 88 * loop.
alpar@9 89 *
alpar@9 90 * RETURNS
alpar@9 91 *
alpar@9 92 * The routine mc21a returns numnz, the number of non-zeros on diagonal
alpar@9 93 * of permuted matrix. */
alpar@9 94
alpar@9 95 int mc21a(int n, const int icn[], const int ip[], const int lenr[],
alpar@9 96 int iperm[], int pr[], int arp[], int cv[], int out[])
alpar@9 97 { int i, ii, in1, in2, j, j1, jord, k, kk, numnz;
alpar@9 98 /* Initialization of arrays. */
alpar@9 99 for (i = 1; i <= n; i++)
alpar@9 100 { arp[i] = lenr[i] - 1;
alpar@9 101 cv[i] = iperm[i] = 0;
alpar@9 102 }
alpar@9 103 numnz = 0;
alpar@9 104 /* Main loop. */
alpar@9 105 /* Each pass round this loop either results in a new assignment
alpar@9 106 or gives a row with no assignment. */
alpar@9 107 for (jord = 1; jord <= n; jord++)
alpar@9 108 { j = jord;
alpar@9 109 pr[j] = -1;
alpar@9 110 for (k = 1; k <= jord; k++)
alpar@9 111 { /* Look for a cheap assignment. */
alpar@9 112 in1 = arp[j];
alpar@9 113 if (in1 >= 0)
alpar@9 114 { in2 = ip[j] + lenr[j] - 1;
alpar@9 115 in1 = in2 - in1;
alpar@9 116 for (ii = in1; ii <= in2; ii++)
alpar@9 117 { i = icn[ii];
alpar@9 118 if (iperm[i] == 0) goto L110;
alpar@9 119 }
alpar@9 120 /* No cheap assignment in row. */
alpar@9 121 arp[j] = -1;
alpar@9 122 }
alpar@9 123 /* Begin looking for assignment chain starting with row j.*/
alpar@9 124 out[j] = lenr[j] - 1;
alpar@9 125 /* Inner loop. Extends chain by one or backtracks. */
alpar@9 126 for (kk = 1; kk <= jord; kk++)
alpar@9 127 { in1 = out[j];
alpar@9 128 if (in1 >= 0)
alpar@9 129 { in2 = ip[j] + lenr[j] - 1;
alpar@9 130 in1 = in2 - in1;
alpar@9 131 /* Forward scan. */
alpar@9 132 for (ii = in1; ii <= in2; ii++)
alpar@9 133 { i = icn[ii];
alpar@9 134 if (cv[i] != jord)
alpar@9 135 { /* Column i has not yet been accessed during
alpar@9 136 this pass. */
alpar@9 137 j1 = j;
alpar@9 138 j = iperm[i];
alpar@9 139 cv[i] = jord;
alpar@9 140 pr[j] = j1;
alpar@9 141 out[j1] = in2 - ii - 1;
alpar@9 142 goto L100;
alpar@9 143 }
alpar@9 144 }
alpar@9 145 }
alpar@9 146 /* Backtracking step. */
alpar@9 147 j = pr[j];
alpar@9 148 if (j == -1) goto L130;
alpar@9 149 }
alpar@9 150 L100: ;
alpar@9 151 }
alpar@9 152 L110: /* New assignment is made. */
alpar@9 153 iperm[i] = j;
alpar@9 154 arp[j] = in2 - ii - 1;
alpar@9 155 numnz++;
alpar@9 156 for (k = 1; k <= jord; k++)
alpar@9 157 { j = pr[j];
alpar@9 158 if (j == -1) break;
alpar@9 159 ii = ip[j] + lenr[j] - out[j] - 2;
alpar@9 160 i = icn[ii];
alpar@9 161 iperm[i] = j;
alpar@9 162 }
alpar@9 163 L130: ;
alpar@9 164 }
alpar@9 165 /* If matrix is structurally singular, we now complete the
alpar@9 166 permutation iperm. */
alpar@9 167 if (numnz < n)
alpar@9 168 { for (i = 1; i <= n; i++)
alpar@9 169 arp[i] = 0;
alpar@9 170 k = 0;
alpar@9 171 for (i = 1; i <= n; i++)
alpar@9 172 { if (iperm[i] == 0)
alpar@9 173 out[++k] = i;
alpar@9 174 else
alpar@9 175 arp[iperm[i]] = i;
alpar@9 176 }
alpar@9 177 k = 0;
alpar@9 178 for (i = 1; i <= n; i++)
alpar@9 179 { if (arp[i] == 0)
alpar@9 180 iperm[out[++k]] = i;
alpar@9 181 }
alpar@9 182 }
alpar@9 183 return numnz;
alpar@9 184 }
alpar@9 185
alpar@9 186 /**********************************************************************/
alpar@9 187
alpar@9 188 #if 0
alpar@9 189 #include "glplib.h"
alpar@9 190
alpar@9 191 int sing;
alpar@9 192
alpar@9 193 void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum,
alpar@9 194 int iw[]);
alpar@9 195
alpar@9 196 void fa01bs(int max, int *nrand);
alpar@9 197
alpar@9 198 int main(void)
alpar@9 199 { /* test program for the routine mc21a */
alpar@9 200 /* these runs on random matrices cause all possible statements in
alpar@9 201 mc21a to be executed */
alpar@9 202 int i, iold, j, j1, j2, jj, knum, l, licn, n, nov4, num, numnz;
alpar@9 203 int ip[1+21], icn[1+1000], iperm[1+20], lenr[1+20], iw1[1+80];
alpar@9 204 licn = 1000;
alpar@9 205 /* run on random matrices of orders 1 through 20 */
alpar@9 206 for (n = 1; n <= 20; n++)
alpar@9 207 { nov4 = n / 4;
alpar@9 208 if (nov4 < 1) nov4 = 1;
alpar@9 209 L10: fa01bs(nov4, &l);
alpar@9 210 knum = l * n;
alpar@9 211 /* knum is requested number of non-zeros in random matrix */
alpar@9 212 if (knum > licn) goto L10;
alpar@9 213 /* if sing is false, matrix is guaranteed structurally
alpar@9 214 non-singular */
alpar@9 215 sing = ((n / 2) * 2 == n);
alpar@9 216 /* call to subroutine to generate random matrix */
alpar@9 217 ranmat(n, n, icn, ip, n+1, &knum, iw1);
alpar@9 218 /* knum is now actual number of non-zeros in random matrix */
alpar@9 219 if (knum > licn) goto L10;
alpar@9 220 xprintf("n = %2d; nz = %4d; sing = %d\n", n, knum, sing);
alpar@9 221 /* set up array of row lengths */
alpar@9 222 for (i = 1; i <= n; i++)
alpar@9 223 lenr[i] = ip[i+1] - ip[i];
alpar@9 224 /* call to mc21a */
alpar@9 225 numnz = mc21a(n, icn, ip, lenr, iperm, &iw1[0], &iw1[n],
alpar@9 226 &iw1[n+n], &iw1[n+n+n]);
alpar@9 227 /* testing to see if there are numnz non-zeros on the diagonal
alpar@9 228 of the permuted matrix. */
alpar@9 229 num = 0;
alpar@9 230 for (i = 1; i <= n; i++)
alpar@9 231 { iold = iperm[i];
alpar@9 232 j1 = ip[iold];
alpar@9 233 j2 = j1 + lenr[iold] - 1;
alpar@9 234 if (j2 < j1) continue;
alpar@9 235 for (jj = j1; jj <= j2; jj++)
alpar@9 236 { j = icn[jj];
alpar@9 237 if (j == i)
alpar@9 238 { num++;
alpar@9 239 break;
alpar@9 240 }
alpar@9 241 }
alpar@9 242 }
alpar@9 243 if (num != numnz)
alpar@9 244 xprintf("Failure in mc21a, numnz = %d instead of %d\n",
alpar@9 245 numnz, num);
alpar@9 246 }
alpar@9 247 return 0;
alpar@9 248 }
alpar@9 249
alpar@9 250 void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum,
alpar@9 251 int iw[])
alpar@9 252 { /* subroutine to generate random matrix */
alpar@9 253 int i, ii, inum, j, lrow, matnum;
alpar@9 254 inum = (*knum / n) * 2;
alpar@9 255 if (inum > n-1) inum = n-1;
alpar@9 256 matnum = 1;
alpar@9 257 /* each pass through this loop generates a row of the matrix */
alpar@9 258 for (j = 1; j <= m; j++)
alpar@9 259 { iptr[j] = matnum;
alpar@9 260 if (!(sing || j > n))
alpar@9 261 icn[matnum++] = j;
alpar@9 262 if (n == 1) continue;
alpar@9 263 for (i = 1; i <= n; i++) iw[i] = 0;
alpar@9 264 if (!sing) iw[j] = 1;
alpar@9 265 fa01bs(inum, &lrow);
alpar@9 266 lrow--;
alpar@9 267 if (lrow == 0) continue;
alpar@9 268 /* lrow off-diagonal non-zeros in row j of the matrix */
alpar@9 269 for (ii = 1; ii <= lrow; ii++)
alpar@9 270 { for (;;)
alpar@9 271 { fa01bs(n, &i);
alpar@9 272 if (iw[i] != 1) break;
alpar@9 273 }
alpar@9 274 iw[i] = 1;
alpar@9 275 icn[matnum++] = i;
alpar@9 276 }
alpar@9 277 }
alpar@9 278 for (i = m+1; i <= nnnp1; i++)
alpar@9 279 iptr[i] = matnum;
alpar@9 280 *knum = matnum - 1;
alpar@9 281 return;
alpar@9 282 }
alpar@9 283
alpar@9 284 double g = 1431655765.0;
alpar@9 285
alpar@9 286 double fa01as(int i)
alpar@9 287 { /* random number generator */
alpar@9 288 g = fmod(g * 9228907.0, 4294967296.0);
alpar@9 289 if (i >= 0)
alpar@9 290 return g / 4294967296.0;
alpar@9 291 else
alpar@9 292 return 2.0 * g / 4294967296.0 - 1.0;
alpar@9 293 }
alpar@9 294
alpar@9 295 void fa01bs(int max, int *nrand)
alpar@9 296 { *nrand = (int)(fa01as(1) * (double)max) + 1;
alpar@9 297 return;
alpar@9 298 }
alpar@9 299 #endif
alpar@9 300
alpar@9 301 /* eof */