Charm++: Charm Source Code Documentation

00001 /*
00002  * Copyright 1997, Regents of the University of Minnesota
00003  *
00004  * smbfactor.c
00005  *
00006  * This file performs the symbolic factorization of a matrix
00007  *
00008  * Started 8/1/97
00009  * George
00010  *
00011  * $Id: smbfactor.c 10154 2011-06-09 21:27:35Z karypis $
00012  *
00013  */
00014 
00015 #include "metisbin.h"
00016 
00017 
00018 /*************************************************************************/
00020 /*************************************************************************/
00021 void ComputeFillIn(graph_t *graph, idx_t *perm, idx_t *iperm, 
00022          size_t *r_maxlnz, size_t *r_opc)
00023 {
00024   idx_t i, j, k, nvtxs, maxlnz, maxsub;
00025   idx_t *xadj, *adjncy;
00026   idx_t *xlnz, *xnzsub, *nzsub;
00027   size_t opc;
00028 
00029 /*
00030   printf("\nSymbolic factorization... --------------------------------------------\n");
00031 */
00032 
00033   nvtxs  = graph->nvtxs;
00034   xadj   = graph->xadj;
00035   adjncy = graph->adjncy;
00036 
00037   maxsub = 8*(nvtxs+xadj[nvtxs]);
00038 
00039   /* Relabel the vertices so that it starts from 1 */
00040   for (i=0; i<xadj[nvtxs]; i++)
00041     adjncy[i]++;
00042   for (i=0; i<nvtxs+1; i++)
00043     xadj[i]++;
00044   for (i=0; i<nvtxs; i++) {
00045     iperm[i]++;
00046     perm[i]++;
00047   }
00048 
00049   /* Allocate the required memory */
00050   xlnz   = imalloc(nvtxs+2, "ComputeFillIn: xlnz");
00051   xnzsub = imalloc(nvtxs+2, "ComputeFillIn: xnzsub");
00052   nzsub  = imalloc(maxsub+1, "ComputeFillIn: nzsub");
00053 
00054   
00055   /* Call sparspak's routine. */
00056   if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) {
00057     printf("Realocating nzsub...\n");
00058     gk_free((void **)&nzsub, LTERM);
00059 
00060     maxsub *= 2;
00061     nzsub  = imalloc(maxsub+1, "ComputeFillIn: nzsub");
00062     if (smbfct(nvtxs, xadj, adjncy, perm, iperm, xlnz, &maxlnz, xnzsub, nzsub, &maxsub)) 
00063       errexit("MAXSUB is too small!");
00064   }
00065 
00066   for (i=0; i<nvtxs; i++)
00067     xlnz[i]--;
00068   for (opc=0, i=0; i<nvtxs; i++)
00069     opc += (xlnz[i+1]-xlnz[i])*(xlnz[i+1]-xlnz[i]) - (xlnz[i+1]-xlnz[i]);
00070 
00071   *r_maxlnz = maxlnz;
00072   *r_opc    = opc;
00073 
00074   gk_free((void **)&xlnz, &xnzsub, &nzsub, LTERM);
00075 
00076   /* Relabel the vertices so that it starts from 0 */
00077   for (i=0; i<nvtxs; i++) {
00078     iperm[i]--;
00079     perm[i]--;
00080   }
00081   for (i=0; i<nvtxs+1; i++)
00082     xadj[i]--;
00083   for (i=0; i<xadj[nvtxs]; i++)
00084     adjncy[i]--;
00085 
00086 }
00087 
00088 
00089 
00090 /*************************************************************************/
00110 /*************************************************************************/
00111 idx_t smbfct(idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t *perm, idx_t *invp, 
00112            idx_t *xlnz, idx_t *maxlnz, idx_t *xnzsub, idx_t *nzsub, 
00113                idx_t *maxsub)
00114 {
00115   /* Local variables */
00116   idx_t node, rchm, mrgk, lmax, i, j, k, m, nabor, nzbeg, nzend;
00117   idx_t kxsub, jstop, jstrt, mrkflg, inz, knz, flag;
00118   idx_t *mrglnk, *marker, *rchlnk;
00119 
00120   rchlnk = ismalloc(neqns+1, 0, "smbfct: rchlnk");
00121   marker = ismalloc(neqns+1, 0, "smbfct: marker");
00122   mrglnk = ismalloc(neqns+1, 0, "smbfct: mgrlnk");
00123 
00124   /* Parameter adjustments */
00125   --marker;
00126   --mrglnk;
00127   --rchlnk;
00128   --nzsub;
00129   --xnzsub;
00130   --xlnz;
00131   --invp;
00132   --perm;
00133   --adjncy;
00134   --xadj;
00135 
00136   /* Function Body */
00137   flag    = 0;
00138   nzbeg   = 1;
00139   nzend   = 0;
00140   xlnz[1] = 1;
00141 
00142   /* FOR EACH COLUMN KNZ COUNTS THE NUMBER OF NONZEROS IN COLUMN K ACCUMULATED IN RCHLNK. */
00143   for (k=1; k<=neqns; k++) {
00144     xnzsub[k] = nzend;
00145     node      = perm[k];
00146     knz       = 0;
00147     mrgk      = mrglnk[k];
00148     mrkflg    = 0;
00149     marker[k] = k;
00150     if (mrgk != 0) {
00151       assert(mrgk > 0 && mrgk <= neqns);
00152       marker[k] = marker[mrgk];
00153     }
00154 
00155     if (xadj[node] >= xadj[node+1]) {
00156       xlnz[k+1] = xlnz[k];
00157       continue;
00158     }
00159 
00160     /* USE RCHLNK TO LINK THROUGH THE STRUCTURE OF A(*,K) BELOW DIAGONAL */
00161     assert(k <= neqns && k > 0);
00162     rchlnk[k] = neqns+1;
00163     for (j=xadj[node]; j<xadj[node+1]; j++) {
00164       nabor = invp[adjncy[j]];
00165       if (nabor <= k) 
00166         continue;
00167       rchm = k;
00168 
00169       do {
00170         m    = rchm;
00171         assert(m > 0 && m <= neqns);
00172         rchm = rchlnk[m];
00173       } while (rchm <= nabor); 
00174 
00175       knz++;
00176       assert(m > 0 && m <= neqns);
00177       rchlnk[m]     = nabor;
00178       assert(nabor > 0 && nabor <= neqns);
00179       rchlnk[nabor] = rchm;
00180       assert(k > 0 && k <= neqns);
00181       if (marker[nabor] != marker[k]) 
00182         mrkflg = 1;
00183     }
00184 
00185 
00186     /* TEST FOR MASS SYMBOLIC ELIMINATION */
00187     lmax = 0;
00188     assert(mrgk >= 0 && mrgk <= neqns);
00189     if (mrkflg != 0 || mrgk == 0 || mrglnk[mrgk] != 0) 
00190       goto L350;
00191     xnzsub[k] = xnzsub[mrgk] + 1;
00192     knz = xlnz[mrgk + 1] - (xlnz[mrgk] + 1);
00193     goto L1400;
00194 
00195 
00196 L350:
00197     /* LINK THROUGH EACH COLUMN I THAT AFFECTS L(*,K) */
00198     i = k;
00199     assert(i > 0 && i <= neqns);
00200     while ((i = mrglnk[i]) != 0) {
00201       assert(i > 0 && i <= neqns);
00202       inz   = xlnz[i+1] - (xlnz[i]+1);
00203       jstrt = xnzsub[i] + 1;
00204       jstop = xnzsub[i] + inz;
00205 
00206       if (inz > lmax) { 
00207         lmax      = inz;
00208         xnzsub[k] = jstrt;
00209       }
00210 
00211       /* MERGE STRUCTURE OF L(*,I) IN NZSUB INTO RCHLNK. */ 
00212       rchm = k;
00213       for (j=jstrt; j<=jstop; j++) {
00214         nabor = nzsub[j];
00215         do {
00216           m    = rchm;
00217           assert(m > 0 && m <= neqns);
00218           rchm = rchlnk[m];
00219         } while (rchm < nabor);
00220 
00221         if (rchm != nabor) {
00222           knz++;
00223           assert(m > 0 && m <= neqns);
00224           rchlnk[m]     = nabor;
00225           assert(nabor > 0 && nabor <= neqns);
00226           rchlnk[nabor] = rchm;
00227           rchm = nabor;
00228         }
00229       }
00230     }
00231 
00232 
00233     /* CHECK IF SUBSCRIPTS DUPLICATE THOSE OF ANOTHER COLUMN */
00234     if (knz == lmax) 
00235       goto L1400;
00236 
00237     /* OR IF TAIL OF K-1ST COLUMN MATCHES HEAD OF KTH */
00238     if (nzbeg > nzend) 
00239       goto L1200;
00240 
00241     assert(k > 0 && k <= neqns);
00242     i = rchlnk[k];
00243     for (jstrt = nzbeg; jstrt <= nzend; ++jstrt) {
00244       if (nzsub[jstrt] < i) 
00245         continue;
00246 
00247       if (nzsub[jstrt] == i) 
00248         goto L1000;
00249       else 
00250         goto L1200;
00251     }
00252     goto L1200;
00253 
00254 
00255 L1000:
00256     xnzsub[k] = jstrt;
00257     for (j = jstrt; j <= nzend; ++j) {
00258       if (nzsub[j] != i) 
00259         goto L1200;
00260       
00261       assert(i > 0 && i <= neqns);
00262       i = rchlnk[i];
00263       if (i > neqns) 
00264         goto L1400;
00265     }
00266     nzend = jstrt - 1;
00267 
00268 
00269     /* COPY THE STRUCTURE OF L(*,K) FROM RCHLNK TO THE DATA STRUCTURE (XNZSUB, NZSUB) */
00270 L1200:
00271     nzbeg = nzend + 1;
00272     nzend += knz;
00273 
00274     if (nzend >= *maxsub) {
00275       flag = 1; /* Out of memory */
00276       break;
00277     }
00278 
00279     i = k;
00280     for (j=nzbeg; j<=nzend; j++) {
00281       assert(i > 0 && i <= neqns);
00282       i = rchlnk[i];
00283       nzsub[j]  = i;
00284       assert(i > 0 && i <= neqns);
00285       marker[i] = k;
00286     }
00287     xnzsub[k] = nzbeg;
00288     assert(k > 0 && k <= neqns);
00289     marker[k] = k;
00290 
00291     /*
00292      * UPDATE THE VECTOR MRGLNK.  NOTE COLUMN L(*,K) JUST FOUND   
00293      * IS REQUIRED TO DETERMINE COLUMN L(*,J), WHERE              
00294      * L(J,K) IS THE FIRST NONZERO IN L(*,K) BELOW DIAGONAL.      
00295      */
00296 L1400:
00297     if (knz > 1) { 
00298       kxsub = xnzsub[k];
00299       i = nzsub[kxsub];
00300       assert(i > 0 && i <= neqns);
00301       assert(k > 0 && k <= neqns);
00302       mrglnk[k] = mrglnk[i];
00303       mrglnk[i] = k;
00304     }
00305 
00306     xlnz[k + 1] = xlnz[k] + knz;
00307   }
00308 
00309   if (flag == 0) {
00310     *maxlnz = xlnz[neqns] - 1;
00311     *maxsub = xnzsub[neqns];
00312     xnzsub[neqns + 1] = xnzsub[neqns];
00313   }
00314 
00315 
00316   marker++;
00317   mrglnk++;
00318   rchlnk++;
00319   nzsub++;
00320   xnzsub++;
00321   xlnz++;
00322   invp++;
00323   perm++;
00324   adjncy++;
00325   xadj++;
00326 
00327   gk_free((void **)&rchlnk, &mrglnk, &marker, LTERM);
00328 
00329   return flag;
00330   
00331 } 
00332
libs/ck-libs/metis/programs/smbfactor.c
