source: inundation-numpy-branch/pymetis/metis-4.0/Programs/smbfactor.c @ 3236

Last change on this file since 3236 was 2051, checked in by jack, 19 years ago

Python interface to metis. Currently provides only the
METIS_PartMeshNodal function, since that is what is currently needed for partitioning.
Module name is metis.
File size: 9.3 KB
Line 
1/*
2 * Copyright 1997, Regents of the University of Minnesota
3 *
4 * smbfactor.c
5 *
6 * This file performs the symbolic factorization of a matrix
7 *
8 * Started 8/1/97
9 * George
10 *
11 * $Id: smbfactor.c,v 1.1 1998/11/27 17:59:40 karypis Exp $
12 *
13 */
14
15#include <metis.h>
16
17
18/*************************************************************************
19* This function sets up data structures for fill-in computations
20**************************************************************************/
21void ComputeFillIn(GraphType *graph, idxtype *iperm)
22{
23  int i, j, k, nvtxs, maxlnz, maxsub;
24  idxtype *xadj, *adjncy;
25  idxtype *perm, *xlnz, *xnzsub, *nzsub;
26  double opc;
27
28/*
29  printf("\nSymbolic factorization... --------------------------------------------\n");
30*/
31
32  nvtxs = graph->nvtxs;
33  xadj = graph->xadj;
34  adjncy = graph->adjncy;
35
36  maxsub = 4*xadj[nvtxs];
37
38  /* Relabel the vertices so that it starts from 1 */
39  k = xadj[nvtxs];
40  for (i=0; i<k; i++)
41    adjncy[i]++;
42  for (i=0; i<nvtxs+1; i++)
43    xadj[i]++;
44
45  /* Allocate the required memory */
46  perm = idxmalloc(nvtxs+1, "ComputeFillIn: perm");
47  xlnz = idxmalloc(nvtxs+1, "ComputeFillIn: xlnz");
48  xnzsub = idxmalloc(nvtxs+1, "ComputeFillIn: xnzsub");
49  nzsub = idxmalloc(maxsub, "ComputeFillIn: nzsub");
50
51  /* Construct perm from iperm and change the numbering of iperm */
52  for (i=0; i<nvtxs; i++)
53    perm[iperm[i]] = i;
54  for (i=0; i<nvtxs; i++) {
55    iperm[i]++;
56    perm[i]++;
57  }
58 
59  /*
60   * Call sparspak routine.
61   */
62  if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) {
63    free(nzsub);
64
65    maxsub = 4*maxsub; 
66    nzsub = idxmalloc(maxsub, "ComputeFillIn: nzsub");
67    if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) 
68      errexit("MAXSUB is too small!");
69  }
70
71  opc = 0;
72  for (i=0; i<nvtxs; i++)
73    xlnz[i]--;
74  for (i=0; i<nvtxs; i++)
75    opc += (xlnz[i+1]-xlnz[i])*(xlnz[i+1]-xlnz[i]) - (xlnz[i+1]-xlnz[i]);
76
77  printf("  Nonzeros: %d, \tOperation Count: %6.4le\n", maxlnz, opc);
78
79
80  GKfree(&perm, &xlnz, &xnzsub, &nzsub, LTERM);
81
82
83  /* Relabel the vertices so that it starts from 0 */
84  for (i=0; i<nvtxs; i++)
85    iperm[i]--;
86  for (i=0; i<nvtxs+1; i++)
87    xadj[i]--;
88  k = xadj[nvtxs];
89  for (i=0; i<k; i++)
90    adjncy[i]--;
91
92}
93
94
95
96/*************************************************************************
97* This function sets up data structures for fill-in computations
98**************************************************************************/
99idxtype ComputeFillIn2(GraphType *graph, idxtype *iperm)
100{
101  int i, j, k, nvtxs, maxlnz, maxsub;
102  idxtype *xadj, *adjncy;
103  idxtype *perm, *xlnz, *xnzsub, *nzsub;
104  double opc;
105
106  nvtxs = graph->nvtxs;
107  xadj = graph->xadj;
108  adjncy = graph->adjncy;
109
110  maxsub = 4*xadj[nvtxs];
111
112  /* Relabel the vertices so that it starts from 1 */
113  k = xadj[nvtxs];
114  for (i=0; i<k; i++)
115    adjncy[i]++;
116  for (i=0; i<nvtxs+1; i++)
117    xadj[i]++;
118
119  /* Allocate the required memory */
120  perm = idxmalloc(nvtxs+1, "ComputeFillIn: perm");
121  xlnz = idxmalloc(nvtxs+1, "ComputeFillIn: xlnz");
122  xnzsub = idxmalloc(nvtxs+1, "ComputeFillIn: xnzsub");
123  nzsub = idxmalloc(maxsub, "ComputeFillIn: nzsub");
124
125  /* Construct perm from iperm and change the numbering of iperm */
126  for (i=0; i<nvtxs; i++)
127    perm[iperm[i]] = i;
128  for (i=0; i<nvtxs; i++) {
129    iperm[i]++;
130    perm[i]++;
131  }
132 
133  /*
134   * Call sparspak routine.
135   */
136  if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) {
137    free(nzsub);
138
139    maxsub = 4*maxsub; 
140    nzsub = idxmalloc(maxsub, "ComputeFillIn: nzsub");
141    if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) 
142      errexit("MAXSUB is too small!");
143  }
144
145  opc = 0;
146  for (i=0; i<nvtxs; i++)
147    xlnz[i]--;
148  for (i=0; i<nvtxs; i++)
149    opc += (xlnz[i+1]-xlnz[i])*(xlnz[i+1]-xlnz[i]) - (xlnz[i+1]-xlnz[i]);
150
151
152  GKfree(&perm, &xlnz, &xnzsub, &nzsub, LTERM);
153
154
155  /* Relabel the vertices so that it starts from 0 */
156  for (i=0; i<nvtxs; i++)
157    iperm[i]--;
158  for (i=0; i<nvtxs+1; i++)
159    xadj[i]--;
160  k = xadj[nvtxs];
161  for (i=0; i<k; i++)
162    adjncy[i]--;
163
164  return maxlnz;
165
166}
167
168
169/*****************************************************************         
170**********     SMBFCT ..... SYMBOLIC FACTORIZATION       *********
171******************************************************************         
172*   PURPOSE - THIS ROUTINE PERFORMS SYMBOLIC FACTORIZATION               
173*   ON A PERMUTED LINEAR SYSTEM AND IT ALSO SETS UP THE               
174*   COMPRESSED DATA STRUCTURE FOR THE SYSTEM.                         
175*
176*   INPUT PARAMETERS -                                               
177*      NEQNS - NUMBER OF EQUATIONS.                                 
178*      (XADJ, ADJNCY) - THE ADJACENCY STRUCTURE.                   
179*      (PERM, INVP) - THE PERMUTATION VECTOR AND ITS INVERSE.     
180*
181*   UPDATED PARAMETERS -                                         
182*      MAXSUB - SIZE OF THE SUBSCRIPT ARRAY NZSUB.  ON RETURN, 
183*             IT CONTAINS THE NUMBER OF SUBSCRIPTS USED       
184*
185*   OUTPUT PARAMETERS -                                       
186*      XLNZ - INDEX INTO THE NONZERO STORAGE VECTOR LNZ.   
187*      (XNZSUB, NZSUB) - THE COMPRESSED SUBSCRIPT VECTORS.
188*      MAXLNZ - THE NUMBER OF NONZEROS FOUND.             
189*
190*******************************************************************/
191int smbfct(int neqns, idxtype *xadj, idxtype *adjncy, idxtype *perm, idxtype *invp, idxtype *xlnz, 
192           int *maxlnz, idxtype *xnzsub, idxtype *nzsub, int *maxsub)
193{
194  /* Local variables */
195  int node, rchm, mrgk, lmax, i, j, k, m, nabor, nzbeg, nzend;
196  int kxsub, jstop, jstrt, mrkflg, inz, knz, flag;
197  idxtype *mrglnk, *marker, *rchlnk;
198
199  rchlnk = idxmalloc(neqns+1, "smbfct: rchlnk");
200  marker = idxsmalloc(neqns+1, 0, "smbfct: marker");
201  mrglnk = idxsmalloc(neqns+1, 0, "smbfct: mgrlnk");
202
203  /* Parameter adjustments */
204  --marker;
205  --mrglnk;
206  --rchlnk;
207  --nzsub;
208  --xnzsub;
209  --xlnz;
210  --invp;
211  --perm;
212  --adjncy;
213  --xadj;
214
215  /* Function Body */
216  flag = 0;
217  nzbeg = 1;
218  nzend = 0;
219  xlnz[1] = 1;
220
221  /* FOR EACH COLUMN KNZ COUNTS THE NUMBER OF NONZEROS IN COLUMN K ACCUMULATED IN RCHLNK. */
222  for (k = 1; k <= neqns; ++k) {
223    knz = 0;
224    mrgk = mrglnk[k];
225    mrkflg = 0;
226    marker[k] = k;
227    if (mrgk != 0) 
228      marker[k] = marker[mrgk];
229    xnzsub[k] = nzend;
230    node = perm[k];
231
232    if (xadj[node] >= xadj[node+1]) {
233      xlnz[k+1] = xlnz[k];
234      continue;
235    }
236
237    /* USE RCHLNK TO LINK THROUGH THE STRUCTURE OF A(*,K) BELOW DIAGONAL */
238    rchlnk[k] = neqns+1;
239    for (j=xadj[node]; j<xadj[node+1]; j++) {
240      nabor = invp[adjncy[j]];
241      if (nabor <= k) 
242        continue;
243      rchm = k;
244
245      do {
246        m = rchm;
247        rchm = rchlnk[m];
248      } while (rchm <= nabor); 
249
250      knz++;
251      rchlnk[m] = nabor;
252      rchlnk[nabor] = rchm;
253      if (marker[nabor] != marker[k]) 
254        mrkflg = 1;
255    }
256
257    /* TEST FOR MASS SYMBOLIC ELIMINATION */
258    lmax = 0;
259    if (mrkflg != 0 || mrgk == 0 || mrglnk[mrgk] != 0) 
260      goto L350;
261    xnzsub[k] = xnzsub[mrgk] + 1;
262    knz = xlnz[mrgk + 1] - (xlnz[mrgk] + 1);
263    goto L1400;
264
265
266    /* LINK THROUGH EACH COLUMN I THAT AFFECTS L(*,K) */
267L350:
268    i = k;
269    while ((i = mrglnk[i]) != 0) {
270      inz = xlnz[i+1] - (xlnz[i]+1);
271      jstrt = xnzsub[i] + 1;
272      jstop = xnzsub[i] + inz;
273
274      if (inz > lmax) { 
275        lmax = inz;
276        xnzsub[k] = jstrt;
277      }
278
279      /* MERGE STRUCTURE OF L(*,I) IN NZSUB INTO RCHLNK. */ 
280      rchm = k;
281      for (j = jstrt; j <= jstop; ++j) {
282        nabor = nzsub[j];
283        do {
284          m = rchm;
285          rchm = rchlnk[m];
286        } while (rchm < nabor);
287
288        if (rchm != nabor) {
289          knz++;
290          rchlnk[m] = nabor;
291          rchlnk[nabor] = rchm;
292          rchm = nabor;
293        }
294      }
295    }
296
297    /* CHECK IF SUBSCRIPTS DUPLICATE THOSE OF ANOTHER COLUMN */
298    if (knz == lmax) 
299      goto L1400;
300
301    /* OR IF TAIL OF K-1ST COLUMN MATCHES HEAD OF KTH */
302    if (nzbeg > nzend) 
303      goto L1200;
304
305    i = rchlnk[k];
306    for (jstrt = nzbeg; jstrt <= nzend; ++jstrt) {
307      if (nzsub[jstrt] < i) 
308        continue;
309
310      if (nzsub[jstrt] == i) 
311        goto L1000;
312      else 
313        goto L1200;
314    }
315    goto L1200;
316
317L1000:
318    xnzsub[k] = jstrt;
319    for (j = jstrt; j <= nzend; ++j) {
320      if (nzsub[j] != i) 
321        goto L1200;
322     
323      i = rchlnk[i];
324      if (i > neqns) 
325        goto L1400;
326    }
327    nzend = jstrt - 1;
328
329    /* COPY THE STRUCTURE OF L(*,K) FROM RCHLNK TO THE DATA STRUCTURE (XNZSUB, NZSUB) */
330L1200:
331    nzbeg = nzend + 1;
332    nzend += knz;
333
334    if (nzend > *maxsub) {
335      flag = 1; /* Out of memory */
336      break;
337    }
338
339    i = k;
340    for (j=nzbeg; j<=nzend; ++j) {
341      i = rchlnk[i];
342      nzsub[j] = i;
343      marker[i] = k;
344    }
345    xnzsub[k] = nzbeg;
346    marker[k] = k;
347
348    /*
349     * UPDATE THE VECTOR MRGLNK.  NOTE COLUMN L(*,K) JUST FOUND   
350     * IS REQUIRED TO DETERMINE COLUMN L(*,J), WHERE             
351     * L(J,K) IS THE FIRST NONZERO IN L(*,K) BELOW DIAGONAL.     
352     */
353L1400:
354    if (knz > 1) { 
355      kxsub = xnzsub[k];
356      i = nzsub[kxsub];
357      mrglnk[k] = mrglnk[i];
358      mrglnk[i] = k;
359    }
360
361    xlnz[k + 1] = xlnz[k] + knz;
362  }
363
364  if (flag == 0) {
365    *maxlnz = xlnz[neqns] - 1;
366    *maxsub = xnzsub[neqns];
367    xnzsub[neqns + 1] = xnzsub[neqns];
368  }
369
370  marker++;
371  mrglnk++;
372  rchlnk++;
373  nzsub++;
374  xnzsub++;
375  xlnz++;
376  invp++;
377  perm++;
378  adjncy++;
379  xadj++;
380  GKfree(&rchlnk, &mrglnk, &marker, LTERM);
381
382  return flag;
383 
384} 
385
Note: See TracBrowser for help on using the repository browser.