COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpnet01.c @ 10:5545663ca997

subpack-glpk
Last change on this file since 10:5545663ca997 was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 13 years ago

Import GLPK 4.47

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