src/libs/ck-libs/metis/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