← Back to branch summary

~elmer-csc-ubuntu/elmercsc/devel

« back to all changes in this revision

Viewing changes to elmergrid/src/metis-5.1.0/programs/smbfactor.c

  • Committer: GitHub
  • Author(s): Peter R
  • Date: 2017-10-04 13:47:14 UTC
  • mfrom: (575.1.17)
  • Revision ID: git-v1:91780af69137b74ec539fedc4ad70c9ec52d732a
Merge pull request #115 from ElmerCSC/metis_update

Updated to new Metis API for devel branch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
elmergrid/src/metis-5.1.0/programs/smbfactor.c
 
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 10154 2011-06-09 21:27:35Z karypis $
 
12
 *
 
13
 */
 
14
 
 
15
#include "metisbin.h"
 
16
 
 
17
 
 
18
/*************************************************************************/
 
19
/*! This function sets up data structures for fill-in computations */
 
20
/*************************************************************************/
 
21
void ComputeFillIn(graph_t *graph, idx_t *perm, idx_t *iperm, 
 
22
         size_t *r_maxlnz, size_t *r_opc)
 
23
{
 
24
  idx_t i, j, k, nvtxs, maxlnz, maxsub;
 
25
  idx_t *xadj, *adjncy;
 
26
  idx_t *xlnz, *xnzsub, *nzsub;
 
27
  size_t opc;
 
28
 
 
29
/*
 
30
  printf("\nSymbolic factorization... --------------------------------------------\n");
 
31
*/
 
32
 
 
33
  nvtxs  = graph->nvtxs;
 
34
  xadj   = graph->xadj;
 
35
  adjncy = graph->adjncy;
 
36
 
 
37
  maxsub = 8*(nvtxs+xadj[nvtxs]);
 
38
 
 
39
  /* Relabel the vertices so that it starts from 1 */
 
40
  for (i=0; i<xadj[nvtxs]; i++)
 
41
    adjncy[i]++;
 
42
  for (i=0; i<nvtxs+1; i++)
 
43
    xadj[i]++;
 
44
  for (i=0; i<nvtxs; i++) {
 
45
    iperm[i]++;
 
46
    perm[i]++;
 
47
  }
 
48
 
 
49
  /* Allocate the required memory */
 
50
  xlnz   = imalloc(nvtxs+2, "ComputeFillIn: xlnz");
 
51
  xnzsub = imalloc(nvtxs+2, "ComputeFillIn: xnzsub");
 
52
  nzsub  = imalloc(maxsub+1, "ComputeFillIn: nzsub");
 
53
 
 
54
  
 
55
  /* Call sparspak's routine. */
 
56
  if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) {
 
57
    printf("Realocating nzsub...\n");
 
58
    gk_free((void **)&nzsub, LTERM);
 
59
 
 
60
    maxsub *= 2;
 
61
    nzsub  = imalloc(maxsub+1, "ComputeFillIn: nzsub");
 
62
    if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) 
 
63
      errexit("MAXSUB is too small!");
 
64
  }
 
65
 
 
66
  for (i=0; i<nvtxs; i++)
 
67
    xlnz[i]--;
 
68
  for (opc=0, i=0; i<nvtxs; i++)
 
69
    opc += (xlnz[i+1]-xlnz[i])*(xlnz[i+1]-xlnz[i]) - (xlnz[i+1]-xlnz[i]);
 
70
 
 
71
  *r_maxlnz = maxlnz;
 
72
  *r_opc    = opc;
 
73
 
 
74
  gk_free((void **)&xlnz, &xnzsub, &nzsub, LTERM);
 
75
 
 
76
  /* Relabel the vertices so that it starts from 0 */
 
77
  for (i=0; i<nvtxs; i++) {
 
78
    iperm[i]--;
 
79
    perm[i]--;
 
80
  }
 
81
  for (i=0; i<nvtxs+1; i++)
 
82
    xadj[i]--;
 
83
  for (i=0; i<xadj[nvtxs]; i++)
 
84
    adjncy[i]--;
 
85
 
 
86
}
 
87
 
 
88
 
 
89
 
 
90
/*************************************************************************/
 
91
/*!
 
92
  PURPOSE - THIS ROUTINE PERFORMS SYMBOLIC FACTORIZATION               
 
93
  ON A PERMUTED LINEAR SYSTEM AND IT ALSO SETS UP THE               
 
94
  COMPRESSED DATA STRUCTURE FOR THE SYSTEM.                         
 
95
 
 
96
  INPUT PARAMETERS -                                               
 
97
     NEQNS - NUMBER OF EQUATIONS.                                 
 
98
     (XADJ, ADJNCY) - THE ADJACENCY STRUCTURE.                   
 
99
     (PERM, INVP) - THE PERMUTATION VECTOR AND ITS INVERSE.     
 
100
 
 
101
  UPDATED PARAMETERS -                                         
 
102
     MAXSUB - SIZE OF THE SUBSCRIPT ARRAY NZSUB.  ON RETURN,  
 
103
            IT CONTAINS THE NUMBER OF SUBSCRIPTS USED        
 
104
 
 
105
  OUTPUT PARAMETERS -                                       
 
106
     XLNZ - INDEX INTO THE NONZERO STORAGE VECTOR LNZ.   
 
107
     (XNZSUB, NZSUB) - THE COMPRESSED SUBSCRIPT VECTORS. 
 
108
     MAXLNZ - THE NUMBER OF NONZEROS FOUND.             
 
109
*/
 
110
/*************************************************************************/
 
111
idx_t smbfct(idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t *perm, idx_t *invp, 
 
112
               idx_t *xlnz, idx_t *maxlnz, idx_t *xnzsub, idx_t *nzsub, 
 
113
               idx_t *maxsub)
 
114
{
 
115
  /* Local variables */
 
116
  idx_t node, rchm, mrgk, lmax, i, j, k, m, nabor, nzbeg, nzend;
 
117
  idx_t kxsub, jstop, jstrt, mrkflg, inz, knz, flag;
 
118
  idx_t *mrglnk, *marker, *rchlnk;
 
119
 
 
120
  rchlnk = ismalloc(neqns+1, 0, "smbfct: rchlnk");
 
121
  marker = ismalloc(neqns+1, 0, "smbfct: marker");
 
122
  mrglnk = ismalloc(neqns+1, 0, "smbfct: mgrlnk");
 
123
 
 
124
  /* Parameter adjustments */
 
125
  --marker;
 
126
  --mrglnk;
 
127
  --rchlnk;
 
128
  --nzsub;
 
129
  --xnzsub;
 
130
  --xlnz;
 
131
  --invp;
 
132
  --perm;
 
133
  --adjncy;
 
134
  --xadj;
 
135
 
 
136
  /* Function Body */
 
137
  flag    = 0;
 
138
  nzbeg   = 1;
 
139
  nzend   = 0;
 
140
  xlnz[1] = 1;
 
141
 
 
142
  /* FOR EACH COLUMN KNZ COUNTS THE NUMBER OF NONZEROS IN COLUMN K ACCUMULATED IN RCHLNK. */
 
143
  for (k=1; k<=neqns; k++) {
 
144
    xnzsub[k] = nzend;
 
145
    node      = perm[k];
 
146
    knz       = 0;
 
147
    mrgk      = mrglnk[k];
 
148
    mrkflg    = 0;
 
149
    marker[k] = k;
 
150
    if (mrgk != 0) {
 
151
      assert(mrgk > 0 && mrgk <= neqns);
 
152
      marker[k] = marker[mrgk];
 
153
    }
 
154
 
 
155
    if (xadj[node] >= xadj[node+1]) {
 
156
      xlnz[k+1] = xlnz[k];
 
157
      continue;
 
158
    }
 
159
 
 
160
    /* USE RCHLNK TO LINK THROUGH THE STRUCTURE OF A(*,K) BELOW DIAGONAL */
 
161
    assert(k <= neqns && k > 0);
 
162
    rchlnk[k] = neqns+1;
 
163
    for (j=xadj[node]; j<xadj[node+1]; j++) {
 
164
      nabor = invp[adjncy[j]];
 
165
      if (nabor <= k) 
 
166
        continue;
 
167
      rchm = k;
 
168
 
 
169
      do {
 
170
        m    = rchm;
 
171
        assert(m > 0 && m <= neqns);
 
172
        rchm = rchlnk[m];
 
173
      } while (rchm <= nabor); 
 
174
 
 
175
      knz++;
 
176
      assert(m > 0 && m <= neqns);
 
177
      rchlnk[m]     = nabor;
 
178
      assert(nabor > 0 && nabor <= neqns);
 
179
      rchlnk[nabor] = rchm;
 
180
      assert(k > 0 && k <= neqns);
 
181
      if (marker[nabor] != marker[k]) 
 
182
        mrkflg = 1;
 
183
    }
 
184
 
 
185
 
 
186
    /* TEST FOR MASS SYMBOLIC ELIMINATION */
 
187
    lmax = 0;
 
188
    assert(mrgk >= 0 && mrgk <= neqns);
 
189
    if (mrkflg != 0 || mrgk == 0 || mrglnk[mrgk] != 0) 
 
190
      goto L350;
 
191
    xnzsub[k] = xnzsub[mrgk] + 1;
 
192
    knz = xlnz[mrgk + 1] - (xlnz[mrgk] + 1);
 
193
    goto L1400;
 
194
 
 
195
 
 
196
L350:
 
197
    /* LINK THROUGH EACH COLUMN I THAT AFFECTS L(*,K) */
 
198
    i = k;
 
199
    assert(i > 0 && i <= neqns);
 
200
    while ((i = mrglnk[i]) != 0) {
 
201
      assert(i > 0 && i <= neqns);
 
202
      inz   = xlnz[i+1] - (xlnz[i]+1);
 
203
      jstrt = xnzsub[i] + 1;
 
204
      jstop = xnzsub[i] + inz;
 
205
 
 
206
      if (inz > lmax) { 
 
207
        lmax      = inz;
 
208
        xnzsub[k] = jstrt;
 
209
      }
 
210
 
 
211
      /* MERGE STRUCTURE OF L(*,I) IN NZSUB INTO RCHLNK. */ 
 
212
      rchm = k;
 
213
      for (j=jstrt; j<=jstop; j++) {
 
214
        nabor = nzsub[j];
 
215
        do {
 
216
          m    = rchm;
 
217
          assert(m > 0 && m <= neqns);
 
218
          rchm = rchlnk[m];
 
219
        } while (rchm < nabor);
 
220
 
 
221
        if (rchm != nabor) {
 
222
          knz++;
 
223
          assert(m > 0 && m <= neqns);
 
224
          rchlnk[m]     = nabor;
 
225
          assert(nabor > 0 && nabor <= neqns);
 
226
          rchlnk[nabor] = rchm;
 
227
          rchm = nabor;
 
228
        }
 
229
      }
 
230
    }
 
231
 
 
232
 
 
233
    /* CHECK IF SUBSCRIPTS DUPLICATE THOSE OF ANOTHER COLUMN */
 
234
    if (knz == lmax) 
 
235
      goto L1400;
 
236
 
 
237
    /* OR IF TAIL OF K-1ST COLUMN MATCHES HEAD OF KTH */
 
238
    if (nzbeg > nzend) 
 
239
      goto L1200;
 
240
 
 
241
    assert(k > 0 && k <= neqns);
 
242
    i = rchlnk[k];
 
243
    for (jstrt = nzbeg; jstrt <= nzend; ++jstrt) {
 
244
      if (nzsub[jstrt] < i) 
 
245
        continue;
 
246
 
 
247
      if (nzsub[jstrt] == i) 
 
248
        goto L1000;
 
249
      else 
 
250
        goto L1200;
 
251
    }
 
252
    goto L1200;
 
253
 
 
254
 
 
255
L1000:
 
256
    xnzsub[k] = jstrt;
 
257
    for (j = jstrt; j <= nzend; ++j) {
 
258
      if (nzsub[j] != i) 
 
259
        goto L1200;
 
260
      
 
261
      assert(i > 0 && i <= neqns);
 
262
      i = rchlnk[i];
 
263
      if (i > neqns) 
 
264
        goto L1400;
 
265
    }
 
266
    nzend = jstrt - 1;
 
267
 
 
268
 
 
269
    /* COPY THE STRUCTURE OF L(*,K) FROM RCHLNK TO THE DATA STRUCTURE (XNZSUB, NZSUB) */
 
270
L1200:
 
271
    nzbeg = nzend + 1;
 
272
    nzend += knz;
 
273
 
 
274
    if (nzend >= *maxsub) {
 
275
      flag = 1; /* Out of memory */
 
276
      break;
 
277
    }
 
278
 
 
279
    i = k;
 
280
    for (j=nzbeg; j<=nzend; j++) {
 
281
      assert(i > 0 && i <= neqns);
 
282
      i = rchlnk[i];
 
283
      nzsub[j]  = i;
 
284
      assert(i > 0 && i <= neqns);
 
285
      marker[i] = k;
 
286
    }
 
287
    xnzsub[k] = nzbeg;
 
288
    assert(k > 0 && k <= neqns);
 
289
    marker[k] = k;
 
290
 
 
291
    /*
 
292
     * UPDATE THE VECTOR MRGLNK.  NOTE COLUMN L(*,K) JUST FOUND   
 
293
     * IS REQUIRED TO DETERMINE COLUMN L(*,J), WHERE              
 
294
     * L(J,K) IS THE FIRST NONZERO IN L(*,K) BELOW DIAGONAL.      
 
295
     */
 
296
L1400:
 
297
    if (knz > 1) { 
 
298
      kxsub = xnzsub[k];
 
299
      i = nzsub[kxsub];
 
300
      assert(i > 0 && i <= neqns);
 
301
      assert(k > 0 && k <= neqns);
 
302
      mrglnk[k] = mrglnk[i];
 
303
      mrglnk[i] = k;
 
304
    }
 
305
 
 
306
    xlnz[k + 1] = xlnz[k] + knz;
 
307
  }
 
308
 
 
309
  if (flag == 0) {
 
310
    *maxlnz = xlnz[neqns] - 1;
 
311
    *maxsub = xnzsub[neqns];
 
312
    xnzsub[neqns + 1] = xnzsub[neqns];
 
313
  }
 
314
 
 
315
 
 
316
  marker++;
 
317
  mrglnk++;
 
318
  rchlnk++;
 
319
  nzsub++;
 
320
  xnzsub++;
 
321
  xlnz++;
 
322
  invp++;
 
323
  perm++;
 
324
  adjncy++;
 
325
  xadj++;
 
326
 
 
327
  gk_free((void **)&rchlnk, &mrglnk, &marker, LTERM);
 
328
 
 
329
  return flag;
 
330
  
 
331
} 
 
332