smbfactor.c
Go to the documentation of this file.
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 /*************************************************************************/
20 /*************************************************************************/
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 /*************************************************************************/
110 /*************************************************************************/
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 
Matrix3f m
void errexit(char *f_str,...)
Definition: error.c:54
idx_t * xadj
#define imalloc
Definition: gklib_rename.h:42
void ComputeFillIn(graph_t *graph, idx_t *perm, idx_t *iperm, size_t *r_maxlnz, size_t *r_opc)
Definition: smbfactor.c:21
#define ismalloc
Definition: gklib_rename.h:68
idx_t idx_t * xadj
idx_t nvtxs
idx_t idx_t idx_t idx_t idx_t idx_t * iperm
NonlinearFactorGraph graph
int32_t idx_t
idx_t smbfct(idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t *perm, idx_t *invp, idx_t *xlnz, idx_t *maxlnz, idx_t *xnzsub, idx_t *nzsub, idx_t *maxsub)
Definition: smbfactor.c:111
idx_t idx_t idx_t idx_t idx_t * perm
void gk_free(void **ptr1,...)
Definition: memory.c:202
idx_t * adjncy
idx_t idx_t idx_t * adjncy
std::ptrdiff_t j
#define LTERM
Definition: gk_defs.h:14


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:52