Charm++: Charm Source Code Documentation

00001 /*
00002  * mmd.c
00003  *
00004  * **************************************************************
00005  * The following C function was developed from a FORTRAN subroutine
00006  * in SPARSPAK written by Eleanor Chu, Alan George, Joseph Liu
00007  * and Esmond Ng.
00008  * 
00009  * The FORTRAN-to-C transformation and modifications such as dynamic
00010  * memory allocation and deallocation were performed by Chunguang
00011  * Sun.
00012  * ************************************************************** 
00013  *
00014  * Taken from SMMS, George 12/13/94
00015  *
00016  * The meaning of invperm, and perm vectors is different from that
00017  * in genqmd_ of SparsPak
00018  *
00019  * $Id: mmd.c 5993 2009-01-07 02:09:57Z karypis $
00020  */
00021 
00022 #include "metislib.h"
00023 
00024 
00025 /*************************************************************************
00026 *  genmmd  -- multiple minimum external degree
00027 *  purpose -- this routine implements the minimum degree
00028 *     algorithm. it makes use of the implicit representation
00029 *     of elimination graphs by quotient graphs, and the notion
00030 *     of indistinguishable nodes. It also implements the modifications
00031 *     by multiple elimination and minimum external degree.
00032 *     Caution -- the adjacency vector adjncy will be destroyed.
00033 *  Input parameters --
00034 *     neqns -- number of equations.
00035 *     (xadj, adjncy) -- the adjacency structure.
00036 *     delta  -- tolerance value for multiple elimination.
00037 *     maxint -- maximum machine representable (short) integer
00038 *               (any smaller estimate will do) for marking nodes.
00039 *  Output parameters --
00040 *     perm -- the minimum degree ordering.
00041 *     invp -- the inverse of perm.
00042 *     *ncsub -- an upper bound on the number of nonzero subscripts
00043 *               for the compressed storage scheme.
00044 *  Working parameters --
00045 *     head -- vector for head of degree lists.
00046 *     invp  -- used temporarily for degree forward link.
00047 *     perm  -- used temporarily for degree backward link.
00048 *     qsize -- vector for size of supernodes.
00049 *     list -- vector for temporary linked lists.
00050 *     marker -- a temporary marker vector.
00051 *  Subroutines used -- mmdelm, mmdint, mmdnum, mmdupd.
00052 **************************************************************************/
00053 void genmmd(idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t *invp, idx_t *perm,
00054      idx_t delta, idx_t *head, idx_t *qsize, idx_t *list, idx_t *marker,
00055      idx_t maxint, idx_t *ncsub)
00056 {
00057     idx_t  ehead, i, mdeg, mdlmt, mdeg_node, nextmd, num, tag;
00058 
00059     if (neqns <= 0)  
00060       return;
00061 
00062     /* Adjust from C to Fortran */
00063     xadj--; adjncy--; invp--; perm--; head--; qsize--; list--; marker--;
00064 
00065     /* initialization for the minimum degree algorithm. */
00066     *ncsub = 0;
00067     mmdint(neqns, xadj, adjncy, head, invp, perm, qsize, list, marker);
00068 
00069     /*  'num' counts the number of ordered nodes plus 1. */
00070     num = 1;
00071 
00072     /* eliminate all isolated nodes. */
00073     nextmd = head[1];
00074     while (nextmd > 0) {
00075       mdeg_node = nextmd;
00076       nextmd = invp[mdeg_node];
00077       marker[mdeg_node] = maxint;
00078       invp[mdeg_node] = -num;
00079       num = num + 1;
00080     }
00081 
00082     /* search for node of the minimum degree. 'mdeg' is the current */
00083     /* minimum degree; 'tag' is used to facilitate marking nodes.   */
00084     if (num > neqns) 
00085       goto n1000;
00086     tag = 1;
00087     head[1] = 0;
00088     mdeg = 2;
00089 
00090     /* infinite loop here ! */
00091     while (1) {
00092       while (head[mdeg] <= 0) 
00093         mdeg++;
00094 
00095       /* use value of 'delta' to set up 'mdlmt', which governs */
00096       /* when a degree update is to be performed.              */
00097       mdlmt = mdeg + delta;
00098       ehead = 0;
00099 
00100 n500:
00101       mdeg_node = head[mdeg];
00102       while (mdeg_node <= 0) {
00103         mdeg++;
00104 
00105         if (mdeg > mdlmt) 
00106           goto n900;
00107         mdeg_node = head[mdeg];
00108       };
00109 
00110       /*  remove 'mdeg_node' from the degree structure. */
00111       nextmd = invp[mdeg_node];
00112       head[mdeg] = nextmd;
00113       if (nextmd > 0)  
00114         perm[nextmd] = -mdeg;
00115       invp[mdeg_node] = -num;
00116       *ncsub += mdeg + qsize[mdeg_node] - 2;
00117       if ((num+qsize[mdeg_node]) > neqns)  
00118         goto n1000;
00119 
00120       /*  eliminate 'mdeg_node' and perform quotient graph */
00121       /*  transformation. reset 'tag' value if necessary.    */
00122       tag++;
00123       if (tag >= maxint) {
00124         tag = 1;
00125         for (i = 1; i <= neqns; i++)
00126           if (marker[i] < maxint)  
00127             marker[i] = 0;
00128       };
00129 
00130       mmdelm(mdeg_node, xadj, adjncy, head, invp, perm, qsize, list, marker, maxint, tag);
00131 
00132       num += qsize[mdeg_node];
00133       list[mdeg_node] = ehead;
00134       ehead = mdeg_node;
00135       if (delta >= 0) 
00136         goto n500;
00137 
00138  n900:
00139       /* update degrees of the nodes involved in the  */
00140       /* minimum degree nodes elimination.            */
00141       if (num > neqns)  
00142         goto n1000;
00143       mmdupd( ehead, neqns, xadj, adjncy, delta, &mdeg, head, invp, perm, qsize, list, marker, maxint, &tag);
00144     }; /* end of -- while ( 1 ) -- */
00145 
00146 n1000:
00147     mmdnum( neqns, perm, invp, qsize );
00148 
00149     /* Adjust from Fortran back to C*/
00150     xadj++; adjncy++; invp++; perm++; head++; qsize++; list++; marker++;
00151 }
00152 
00153 
00154 /**************************************************************************
00155 *           mmdelm ...... multiple minimum degree elimination
00156 * Purpose -- This routine eliminates the node mdeg_node of minimum degree
00157 *     from the adjacency structure, which is stored in the quotient
00158 *     graph format. It also transforms the quotient graph representation
00159 *     of the elimination graph.
00160 * Input parameters --
00161 *     mdeg_node -- node of minimum degree.
00162 *     maxint -- estimate of maximum representable (short) integer.
00163 *     tag    -- tag value.
00164 * Updated parameters --
00165 *     (xadj, adjncy) -- updated adjacency structure.
00166 *     (head, forward, backward) -- degree doubly linked structure.
00167 *     qsize -- size of supernode.
00168 *     marker -- marker vector.
00169 *     list -- temporary linked list of eliminated nabors.
00170 ***************************************************************************/
00171 void mmdelm(idx_t mdeg_node, idx_t *xadj, idx_t *adjncy, idx_t *head, idx_t *forward,
00172      idx_t *backward, idx_t *qsize, idx_t *list, idx_t *marker, idx_t maxint, idx_t tag)
00173 {
00174     idx_t   element, i,   istop, istart, j,
00175           jstop, jstart, link,
00176           nabor, node, npv, nqnbrs, nxnode,
00177           pvnode, rlmt, rloc, rnode, xqnbr;
00178 
00179     /* find the reachable set of 'mdeg_node' and */
00180     /* place it in the data structure.           */
00181     marker[mdeg_node] = tag;
00182     istart = xadj[mdeg_node];
00183     istop = xadj[mdeg_node+1] - 1;
00184 
00185     /* 'element' points to the beginning of the list of  */
00186     /* eliminated nabors of 'mdeg_node', and 'rloc' gives the */
00187     /* storage location for the next reachable node.   */
00188     element = 0;
00189     rloc = istart;
00190     rlmt = istop;
00191     for ( i = istart; i <= istop; i++ ) {
00192         nabor = adjncy[i];
00193         if ( nabor == 0 ) break;
00194         if ( marker[nabor] < tag ) {
00195            marker[nabor] = tag;
00196            if ( forward[nabor] < 0 )  {
00197               list[nabor] = element;
00198               element = nabor;
00199            } else {
00200               adjncy[rloc] = nabor;
00201               rloc++;
00202            };
00203         }; /* end of -- if -- */
00204     }; /* end of -- for -- */
00205 
00206   /* merge with reachable nodes from generalized elements. */
00207   while ( element > 0 ) {
00208       adjncy[rlmt] = -element;
00209       link = element;
00210 
00211 n400:
00212       jstart = xadj[link];
00213       jstop = xadj[link+1] - 1;
00214       for ( j = jstart; j <= jstop; j++ ) {
00215           node = adjncy[j];
00216           link = -node;
00217           if ( node < 0 )  goto n400;
00218           if ( node == 0 ) break;
00219           if ((marker[node]<tag)&&(forward[node]>=0)) {
00220              marker[node] = tag;
00221              /*use storage from eliminated nodes if necessary.*/
00222              while ( rloc >= rlmt ) {
00223                    link = -adjncy[rlmt];
00224                    rloc = xadj[link];
00225                    rlmt = xadj[link+1] - 1;
00226              };
00227              adjncy[rloc] = node;
00228              rloc++;
00229           };
00230       }; /* end of -- for ( j = jstart; -- */
00231       element = list[element];
00232     };  /* end of -- while ( element > 0 ) -- */
00233     if ( rloc <= rlmt ) adjncy[rloc] = 0;
00234     /* for each node in the reachable set, do the following. */
00235     link = mdeg_node;
00236 
00237 n1100:
00238     istart = xadj[link];
00239     istop = xadj[link+1] - 1;
00240     for ( i = istart; i <= istop; i++ ) {
00241         rnode = adjncy[i];
00242         link = -rnode;
00243         if ( rnode < 0 ) goto n1100;
00244         if ( rnode == 0 ) return;
00245 
00246         /* 'rnode' is in the degree list structure. */
00247         pvnode = backward[rnode];
00248         if (( pvnode != 0 ) && ( pvnode != (-maxint) )) {
00249            /* then remove 'rnode' from the structure. */
00250            nxnode = forward[rnode];
00251            if ( nxnode > 0 ) backward[nxnode] = pvnode;
00252            if ( pvnode > 0 ) forward[pvnode] = nxnode;
00253            npv = -pvnode;
00254            if ( pvnode < 0 ) head[npv] = nxnode;
00255         };
00256 
00257         /* purge inactive quotient nabors of 'rnode'. */
00258         jstart = xadj[rnode];
00259         jstop = xadj[rnode+1] - 1;
00260         xqnbr = jstart;
00261         for ( j = jstart; j <= jstop; j++ ) {
00262             nabor = adjncy[j];
00263             if ( nabor == 0 ) break;
00264             if ( marker[nabor] < tag ) {
00265                 adjncy[xqnbr] = nabor;
00266                 xqnbr++;
00267             };
00268         };
00269 
00270         /* no active nabor after the purging. */
00271         nqnbrs = xqnbr - jstart;
00272         if ( nqnbrs <= 0 ) {
00273            /* merge 'rnode' with 'mdeg_node'. */
00274            qsize[mdeg_node] += qsize[rnode];
00275            qsize[rnode] = 0;
00276            marker[rnode] = maxint;
00277            forward[rnode] = -mdeg_node;
00278            backward[rnode] = -maxint;
00279         } else {
00280            /* flag 'rnode' for degree update, and  */
00281            /* add 'mdeg_node' as a nabor of 'rnode'.      */
00282            forward[rnode] = nqnbrs + 1;
00283            backward[rnode] = 0;
00284            adjncy[xqnbr] = mdeg_node;
00285            xqnbr++;
00286            if ( xqnbr <= jstop )  adjncy[xqnbr] = 0;
00287         };
00288       }; /* end of -- for ( i = istart; -- */
00289       return;
00290  }
00291 
00292 /***************************************************************************
00293 *    mmdint ---- mult minimum degree initialization
00294 *    purpose -- this routine performs initialization for the
00295 *       multiple elimination version of the minimum degree algorithm.
00296 *    input parameters --
00297 *       neqns  -- number of equations.
00298 *       (xadj, adjncy) -- adjacency structure.
00299 *    output parameters --
00300 *       (head, dfrow, backward) -- degree doubly linked structure.
00301 *       qsize -- size of supernode ( initialized to one).
00302 *       list -- linked list.
00303 *       marker -- marker vector.
00304 ****************************************************************************/
00305 idx_t  mmdint(idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t *head, idx_t *forward,
00306      idx_t *backward, idx_t *qsize, idx_t *list, idx_t *marker)
00307 {
00308     idx_t  fnode, ndeg, node;
00309 
00310     for ( node = 1; node <= neqns; node++ ) {
00311         head[node] = 0;
00312         qsize[node] = 1;
00313         marker[node] = 0;
00314         list[node] = 0;
00315     };
00316 
00317     /* initialize the degree doubly linked lists. */
00318     for ( node = 1; node <= neqns; node++ ) {
00319         ndeg = xadj[node+1] - xadj[node]/* + 1*/;   /* george */
00320         if (ndeg == 0)
00321           ndeg = 1;
00322         fnode = head[ndeg];
00323         forward[node] = fnode;
00324         head[ndeg] = node;
00325         if ( fnode > 0 ) backward[fnode] = node;
00326         backward[node] = -ndeg;
00327     };
00328     return 0;
00329 }
00330 
00331 /****************************************************************************
00332 * mmdnum --- multi minimum degree numbering
00333 * purpose -- this routine performs the final step in producing
00334 *    the permutation and inverse permutation vectors in the
00335 *    multiple elimination version of the minimum degree
00336 *    ordering algorithm.
00337 * input parameters --
00338 *     neqns -- number of equations.
00339 *     qsize -- size of supernodes at elimination.
00340 * updated parameters --
00341 *     invp -- inverse permutation vector. on input,
00342 *             if qsize[node] = 0, then node has been merged
00343 *             into the node -invp[node]; otherwise,
00344 *            -invp[node] is its inverse labelling.
00345 * output parameters --
00346 *     perm -- the permutation vector.
00347 ****************************************************************************/
00348 void mmdnum(idx_t neqns, idx_t *perm, idx_t *invp, idx_t *qsize)
00349 {
00350   idx_t father, nextf, node, nqsize, num, root;
00351 
00352   for ( node = 1; node <= neqns; node++ ) {
00353       nqsize = qsize[node];
00354       if ( nqsize <= 0 ) perm[node] = invp[node];
00355       if ( nqsize > 0 )  perm[node] = -invp[node];
00356   };
00357 
00358   /* for each node which has been merged, do the following. */
00359   for ( node = 1; node <= neqns; node++ ) {
00360       if ( perm[node] <= 0 )  {
00361 
00362      /* trace the merged tree until one which has not */
00363          /* been merged, call it root.                    */
00364          father = node;
00365          while ( perm[father] <= 0 )
00366             father = - perm[father];
00367 
00368          /* number node after root. */
00369          root = father;
00370          num = perm[root] + 1;
00371          invp[node] = -num;
00372          perm[root] = num;
00373 
00374          /* shorten the merged tree. */
00375          father = node;
00376          nextf = - perm[father];
00377          while ( nextf > 0 ) {
00378             perm[father] = -root;
00379             father = nextf;
00380             nextf = -perm[father];
00381          };
00382       };  /* end of -- if ( perm[node] <= 0 ) -- */
00383   }; /* end of -- for ( node = 1; -- */
00384 
00385   /* ready to compute perm. */
00386   for ( node = 1; node <= neqns; node++ ) {
00387         num = -invp[node];
00388         invp[node] = num;
00389         perm[num] = node;
00390   };
00391   return;
00392 }
00393 
00394 /****************************************************************************
00395 * mmdupd ---- multiple minimum degree update
00396 * purpose -- this routine updates the degrees of nodes after a
00397 *            multiple elimination step.
00398 * input parameters --
00399 *    ehead -- the beginning of the list of eliminated nodes
00400 *             (i.e., newly formed elements).
00401 *    neqns -- number of equations.
00402 *    (xadj, adjncy) -- adjacency structure.
00403 *    delta -- tolerance value for multiple elimination.
00404 *    maxint -- maximum machine representable (short) integer.
00405 * updated parameters --
00406 *    mdeg -- new minimum degree after degree update.
00407 *    (head, forward, backward) -- degree doubly linked structure.
00408 *    qsize -- size of supernode.
00409 *    list -- marker vector for degree update.
00410 *    *tag   -- tag value.
00411 ****************************************************************************/
00412 void mmdupd(idx_t ehead, idx_t neqns, idx_t *xadj, idx_t *adjncy, idx_t delta, idx_t *mdeg,
00413      idx_t *head, idx_t *forward, idx_t *backward, idx_t *qsize, idx_t *list,
00414      idx_t *marker, idx_t maxint, idx_t *tag)
00415 {
00416  idx_t  deg, deg0, element, enode, fnode, i, iq2, istop,
00417       istart, j, jstop, jstart, link, mdeg0, mtag, nabor,
00418       node, q2head, qxhead;
00419 
00420       mdeg0 = *mdeg + delta;
00421       element = ehead;
00422 
00423 n100:
00424       if ( element <= 0 ) return;
00425 
00426       /* for each of the newly formed element, do the following. */
00427       /* reset tag value if necessary.                           */
00428       mtag = *tag + mdeg0;
00429       if ( mtag >= maxint ) {
00430          *tag = 1;
00431          for ( i = 1; i <= neqns; i++ )
00432              if ( marker[i] < maxint ) marker[i] = 0;
00433          mtag = *tag + mdeg0;
00434       };
00435 
00436       /* create two linked lists from nodes associated with 'element': */
00437       /* one with two nabors (q2head) in the adjacency structure, and the*/
00438       /* other with more than two nabors (qxhead). also compute 'deg0',*/
00439       /* number of nodes in this element.                              */
00440       q2head = 0;
00441       qxhead = 0;
00442       deg0 = 0;
00443       link =element;
00444 
00445 n400:
00446       istart = xadj[link];
00447       istop = xadj[link+1] - 1;
00448       for ( i = istart; i <= istop; i++ ) {
00449           enode = adjncy[i];
00450           link = -enode;
00451           if ( enode < 0 )  goto n400;
00452           if ( enode == 0 ) break;
00453           if ( qsize[enode] != 0 ) {
00454              deg0 += qsize[enode];
00455              marker[enode] = mtag;
00456 
00457              /*'enode' requires a degree update*/
00458              if ( backward[enode] == 0 ) {
00459                 /* place either in qxhead or q2head list. */
00460                 if ( forward[enode] != 2 ) {
00461                      list[enode] = qxhead;
00462                      qxhead = enode;
00463                 } else {
00464                      list[enode] = q2head;
00465                      q2head = enode;
00466                 };
00467              };
00468           }; /* enf of -- if ( qsize[enode] != 0 ) -- */
00469       }; /* end of -- for ( i = istart; -- */
00470 
00471       /* for each node in q2 list, do the following. */
00472       enode = q2head;
00473       iq2 = 1;
00474 
00475 n900:
00476       if ( enode <= 0 ) goto n1500;
00477       if ( backward[enode] != 0 ) goto n2200;
00478       (*tag)++;
00479       deg = deg0;
00480 
00481       /* identify the other adjacent element nabor. */
00482       istart = xadj[enode];
00483       nabor = adjncy[istart];
00484       if ( nabor == element ) nabor = adjncy[istart+1];
00485       link = nabor;
00486       if ( forward[nabor] >= 0 ) {
00487            /* nabor is uneliminated, increase degree count. */
00488            deg += qsize[nabor];
00489            goto n2100;
00490       };
00491 
00492        /* the nabor is eliminated. for each node in the 2nd element */
00493        /* do the following.                                         */
00494 n1000:
00495        istart = xadj[link];
00496        istop = xadj[link+1] - 1;
00497        for ( i = istart; i <= istop; i++ ) {
00498            node = adjncy[i];
00499            link = -node;
00500            if ( node != enode ) {
00501                 if ( node < 0 ) goto n1000;
00502                 if ( node == 0 )  goto n2100;
00503                 if ( qsize[node] != 0 ) {
00504                      if ( marker[node] < *tag ) {
00505                         /* 'node' is not yet considered. */
00506                         marker[node] = *tag;
00507                         deg += qsize[node];
00508                      } else {
00509                         if ( backward[node] == 0 ) {
00510                              if ( forward[node] == 2 ) {
00511                                 /* 'node' is indistinguishable from 'enode'.*/
00512                                 /* merge them into a new supernode.         */
00513                                 qsize[enode] += qsize[node];
00514                                 qsize[node] = 0;
00515                                 marker[node] = maxint;
00516                                 forward[node] = -enode;
00517                                 backward[node] = -maxint;
00518                              } else {
00519                                 /* 'node' is outmacthed by 'enode' */
00520                 if (backward[node]==0) backward[node] = -maxint;
00521                              };
00522                         }; /* end of -- if ( backward[node] == 0 ) -- */
00523                     }; /* end of -- if ( marker[node] < *tag ) -- */
00524                 }; /* end of -- if ( qsize[node] != 0 ) -- */
00525               }; /* end of -- if ( node != enode ) -- */
00526           }; /* end of -- for ( i = istart; -- */
00527           goto n2100;
00528 
00529 n1500:
00530           /* for each 'enode' in the 'qx' list, do the following. */
00531           enode = qxhead;
00532           iq2 = 0;
00533 
00534 n1600:    if ( enode <= 0 )  goto n2300;
00535           if ( backward[enode] != 0 )  goto n2200;
00536           (*tag)++;
00537           deg = deg0;
00538 
00539           /*for each unmarked nabor of 'enode', do the following.*/
00540           istart = xadj[enode];
00541           istop = xadj[enode+1] - 1;
00542           for ( i = istart; i <= istop; i++ ) {
00543                 nabor = adjncy[i];
00544                 if ( nabor == 0 ) break;
00545                 if ( marker[nabor] < *tag ) {
00546                      marker[nabor] = *tag;
00547                      link = nabor;
00548                      if ( forward[nabor] >= 0 ) 
00549                           /*if uneliminated, include it in deg count.*/
00550                           deg += qsize[nabor];
00551                      else {
00552 n1700:
00553                           /* if eliminated, include unmarked nodes in this*/
00554                           /* element into the degree count.             */
00555                           jstart = xadj[link];
00556                           jstop = xadj[link+1] - 1;
00557                           for ( j = jstart; j <= jstop; j++ ) {
00558                                 node = adjncy[j];
00559                                 link = -node;
00560                                 if ( node < 0 ) goto n1700;
00561                                 if ( node == 0 ) break;
00562                                 if ( marker[node] < *tag ) {
00563                                     marker[node] = *tag;
00564                                     deg += qsize[node];
00565                                 };
00566                           }; /* end of -- for ( j = jstart; -- */
00567                      }; /* end of -- if ( forward[nabor] >= 0 ) -- */
00568                   }; /* end of -- if ( marker[nabor] < *tag ) -- */
00569           }; /* end of -- for ( i = istart; -- */
00570 
00571 n2100:
00572           /* update external degree of 'enode' in degree structure, */
00573           /* and '*mdeg' if necessary.                     */
00574           deg = deg - qsize[enode] + 1;
00575           fnode = head[deg];
00576           forward[enode] = fnode;
00577           backward[enode] = -deg;
00578           if ( fnode > 0 ) backward[fnode] = enode;
00579           head[deg] = enode;
00580           if ( deg < *mdeg ) *mdeg = deg;
00581 
00582 n2200:
00583           /* get next enode in current element. */
00584           enode = list[enode];
00585           if ( iq2 == 1 ) goto n900;
00586           goto n1600;
00587 
00588 n2300:
00589           /* get next element in the list. */
00590           *tag = mtag;
00591           element = list[element];
00592           goto n100;
00593     }
libs/ck-libs/metis/libmetis/mmd.c
