PageRenderTime 50ms CodeModel.GetById 14ms app.highlight 26ms RepoModel.GetById 1ms app.codeStats 1ms

/farmR/src/glpnet02.c

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