1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpnet01.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,301 @@
1.4 +/* glpnet01.c (permutations for zero-free diagonal) */
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 +* MC21A and MC21B associated with the following paper:
1.11 +*
1.12 +* I.S.Duff, Algorithm 575: Permutations for zero-free diagonal, ACM
1.13 +* Trans. on Math. Softw. 7 (1981), 387-390.
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 +* mc21a - permutations for zero-free diagonal
1.40 +*
1.41 +* SYNOPSIS
1.42 +*
1.43 +* #include "glpnet.h"
1.44 +* int mc21a(int n, const int icn[], const int ip[], const int lenr[],
1.45 +* int iperm[], int pr[], int arp[], int cv[], int out[]);
1.46 +*
1.47 +* DESCRIPTION
1.48 +*
1.49 +* Given the pattern of nonzeros of a sparse matrix, the routine mc21a
1.50 +* attempts to find a permutation of its rows that makes the matrix have
1.51 +* no zeros on its diagonal.
1.52 +*
1.53 +* INPUT PARAMETERS
1.54 +*
1.55 +* n order of 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 PARAMETER
1.68 +*
1.69 +* iperm contains permutation to make diagonal have the smallest
1.70 +* number of zeros on it. Elements (iperm[i], i), i = 1,2,...,n,
1.71 +* are non-zero at the end of the algorithm unless the matrix is
1.72 +* structurally singular. In this case, (iperm[i], i) will be
1.73 +* zero for n - numnz entries.
1.74 +*
1.75 +* WORKING ARRAYS
1.76 +*
1.77 +* pr working array of length [1+n], where pr[0] is not used.
1.78 +* pr[i] is the previous row to i in the depth first search.
1.79 +*
1.80 +* arp working array of length [1+n], where arp[0] is not used.
1.81 +* arp[i] is one less than the number of non-zeros in row i which
1.82 +* have not been scanned when looking for a cheap assignment.
1.83 +*
1.84 +* cv working array of length [1+n], where cv[0] is not used.
1.85 +* cv[i] is the most recent row extension at which column i was
1.86 +* visited.
1.87 +*
1.88 +* out working array of length [1+n], where out[0] is not used.
1.89 +* out[i] is one less than the number of non-zeros in row i
1.90 +* which have not been scanned during one pass through the main
1.91 +* loop.
1.92 +*
1.93 +* RETURNS
1.94 +*
1.95 +* The routine mc21a returns numnz, the number of non-zeros on diagonal
1.96 +* of permuted matrix. */
1.97 +
1.98 +int mc21a(int n, const int icn[], const int ip[], const int lenr[],
1.99 + int iperm[], int pr[], int arp[], int cv[], int out[])
1.100 +{ int i, ii, in1, in2, j, j1, jord, k, kk, numnz;
1.101 + /* Initialization of arrays. */
1.102 + for (i = 1; i <= n; i++)
1.103 + { arp[i] = lenr[i] - 1;
1.104 + cv[i] = iperm[i] = 0;
1.105 + }
1.106 + numnz = 0;
1.107 + /* Main loop. */
1.108 + /* Each pass round this loop either results in a new assignment
1.109 + or gives a row with no assignment. */
1.110 + for (jord = 1; jord <= n; jord++)
1.111 + { j = jord;
1.112 + pr[j] = -1;
1.113 + for (k = 1; k <= jord; k++)
1.114 + { /* Look for a cheap assignment. */
1.115 + in1 = arp[j];
1.116 + if (in1 >= 0)
1.117 + { in2 = ip[j] + lenr[j] - 1;
1.118 + in1 = in2 - in1;
1.119 + for (ii = in1; ii <= in2; ii++)
1.120 + { i = icn[ii];
1.121 + if (iperm[i] == 0) goto L110;
1.122 + }
1.123 + /* No cheap assignment in row. */
1.124 + arp[j] = -1;
1.125 + }
1.126 + /* Begin looking for assignment chain starting with row j.*/
1.127 + out[j] = lenr[j] - 1;
1.128 + /* Inner loop. Extends chain by one or backtracks. */
1.129 + for (kk = 1; kk <= jord; kk++)
1.130 + { in1 = out[j];
1.131 + if (in1 >= 0)
1.132 + { in2 = ip[j] + lenr[j] - 1;
1.133 + in1 = in2 - in1;
1.134 + /* Forward scan. */
1.135 + for (ii = in1; ii <= in2; ii++)
1.136 + { i = icn[ii];
1.137 + if (cv[i] != jord)
1.138 + { /* Column i has not yet been accessed during
1.139 + this pass. */
1.140 + j1 = j;
1.141 + j = iperm[i];
1.142 + cv[i] = jord;
1.143 + pr[j] = j1;
1.144 + out[j1] = in2 - ii - 1;
1.145 + goto L100;
1.146 + }
1.147 + }
1.148 + }
1.149 + /* Backtracking step. */
1.150 + j = pr[j];
1.151 + if (j == -1) goto L130;
1.152 + }
1.153 +L100: ;
1.154 + }
1.155 +L110: /* New assignment is made. */
1.156 + iperm[i] = j;
1.157 + arp[j] = in2 - ii - 1;
1.158 + numnz++;
1.159 + for (k = 1; k <= jord; k++)
1.160 + { j = pr[j];
1.161 + if (j == -1) break;
1.162 + ii = ip[j] + lenr[j] - out[j] - 2;
1.163 + i = icn[ii];
1.164 + iperm[i] = j;
1.165 + }
1.166 +L130: ;
1.167 + }
1.168 + /* If matrix is structurally singular, we now complete the
1.169 + permutation iperm. */
1.170 + if (numnz < n)
1.171 + { for (i = 1; i <= n; i++)
1.172 + arp[i] = 0;
1.173 + k = 0;
1.174 + for (i = 1; i <= n; i++)
1.175 + { if (iperm[i] == 0)
1.176 + out[++k] = i;
1.177 + else
1.178 + arp[iperm[i]] = i;
1.179 + }
1.180 + k = 0;
1.181 + for (i = 1; i <= n; i++)
1.182 + { if (arp[i] == 0)
1.183 + iperm[out[++k]] = i;
1.184 + }
1.185 + }
1.186 + return numnz;
1.187 +}
1.188 +
1.189 +/**********************************************************************/
1.190 +
1.191 +#if 0
1.192 +#include "glplib.h"
1.193 +
1.194 +int sing;
1.195 +
1.196 +void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum,
1.197 + int iw[]);
1.198 +
1.199 +void fa01bs(int max, int *nrand);
1.200 +
1.201 +int main(void)
1.202 +{ /* test program for the routine mc21a */
1.203 + /* these runs on random matrices cause all possible statements in
1.204 + mc21a to be executed */
1.205 + int i, iold, j, j1, j2, jj, knum, l, licn, n, nov4, num, numnz;
1.206 + int ip[1+21], icn[1+1000], iperm[1+20], lenr[1+20], iw1[1+80];
1.207 + licn = 1000;
1.208 + /* run on random matrices of orders 1 through 20 */
1.209 + for (n = 1; n <= 20; n++)
1.210 + { nov4 = n / 4;
1.211 + if (nov4 < 1) nov4 = 1;
1.212 +L10: fa01bs(nov4, &l);
1.213 + knum = l * n;
1.214 + /* knum is requested number of non-zeros in random matrix */
1.215 + if (knum > licn) goto L10;
1.216 + /* if sing is false, matrix is guaranteed structurally
1.217 + non-singular */
1.218 + sing = ((n / 2) * 2 == n);
1.219 + /* call to subroutine to generate random matrix */
1.220 + ranmat(n, n, icn, ip, n+1, &knum, iw1);
1.221 + /* knum is now actual number of non-zeros in random matrix */
1.222 + if (knum > licn) goto L10;
1.223 + xprintf("n = %2d; nz = %4d; sing = %d\n", n, knum, sing);
1.224 + /* set up array of row lengths */
1.225 + for (i = 1; i <= n; i++)
1.226 + lenr[i] = ip[i+1] - ip[i];
1.227 + /* call to mc21a */
1.228 + numnz = mc21a(n, icn, ip, lenr, iperm, &iw1[0], &iw1[n],
1.229 + &iw1[n+n], &iw1[n+n+n]);
1.230 + /* testing to see if there are numnz non-zeros on the diagonal
1.231 + of the permuted matrix. */
1.232 + num = 0;
1.233 + for (i = 1; i <= n; i++)
1.234 + { iold = iperm[i];
1.235 + j1 = ip[iold];
1.236 + j2 = j1 + lenr[iold] - 1;
1.237 + if (j2 < j1) continue;
1.238 + for (jj = j1; jj <= j2; jj++)
1.239 + { j = icn[jj];
1.240 + if (j == i)
1.241 + { num++;
1.242 + break;
1.243 + }
1.244 + }
1.245 + }
1.246 + if (num != numnz)
1.247 + xprintf("Failure in mc21a, numnz = %d instead of %d\n",
1.248 + numnz, num);
1.249 + }
1.250 + return 0;
1.251 +}
1.252 +
1.253 +void ranmat(int m, int n, int icn[], int iptr[], int nnnp1, int *knum,
1.254 + int iw[])
1.255 +{ /* subroutine to generate random matrix */
1.256 + int i, ii, inum, j, lrow, matnum;
1.257 + inum = (*knum / n) * 2;
1.258 + if (inum > n-1) inum = n-1;
1.259 + matnum = 1;
1.260 + /* each pass through this loop generates a row of the matrix */
1.261 + for (j = 1; j <= m; j++)
1.262 + { iptr[j] = matnum;
1.263 + if (!(sing || j > n))
1.264 + icn[matnum++] = j;
1.265 + if (n == 1) continue;
1.266 + for (i = 1; i <= n; i++) iw[i] = 0;
1.267 + if (!sing) iw[j] = 1;
1.268 + fa01bs(inum, &lrow);
1.269 + lrow--;
1.270 + if (lrow == 0) continue;
1.271 + /* lrow off-diagonal non-zeros in row j of the matrix */
1.272 + for (ii = 1; ii <= lrow; ii++)
1.273 + { for (;;)
1.274 + { fa01bs(n, &i);
1.275 + if (iw[i] != 1) break;
1.276 + }
1.277 + iw[i] = 1;
1.278 + icn[matnum++] = i;
1.279 + }
1.280 + }
1.281 + for (i = m+1; i <= nnnp1; i++)
1.282 + iptr[i] = matnum;
1.283 + *knum = matnum - 1;
1.284 + return;
1.285 +}
1.286 +
1.287 +double g = 1431655765.0;
1.288 +
1.289 +double fa01as(int i)
1.290 +{ /* random number generator */
1.291 + g = fmod(g * 9228907.0, 4294967296.0);
1.292 + if (i >= 0)
1.293 + return g / 4294967296.0;
1.294 + else
1.295 + return 2.0 * g / 4294967296.0 - 1.0;
1.296 +}
1.297 +
1.298 +void fa01bs(int max, int *nrand)
1.299 +{ *nrand = (int)(fa01as(1) * (double)max) + 1;
1.300 + return;
1.301 +}
1.302 +#endif
1.303 +
1.304 +/* eof */