glpnet02.c 10 KB

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