EpetraExt Development
mmc13e.f
Go to the documentation of this file.
00001       subroutine   mmc13e   ( nrows , ncols , nhcols, nhrows, nscols,
00002      $                        sqindx, hrzcmp, rowstr, colind, colset, 
00003      $                        trycol, cbegin, lowlnk, prev  , colmrk, 
00004      $                        ccmstr, rcmstr, cnto  , rnto  , sqcmpn )
00005 
00006 c     ==================================================================
00007 c     ==================================================================
00008 c     ====  mmc13e -- lower block triangular form of square matrix  ====
00009 c     ==================================================================
00010 c     ==================================================================
00011 
00012 c     mmc13e :   modified from harwell mc13e by alex pothen and
00013 c                chin-ju fan
00014 c     bcs modifications, john lewis, sept. 1990
00015 
00016 c     finds the lower block triangular form of the square submatrix
00017 c     in the general block triangular form.  the square submatrix
00018 c     consists entirely of matched rows and columns.  therefore,
00019 c     with each row matched to its matching column, the submatrix
00020 c     has a nonzero diagonal, as required by duff's algorithm.
00021 c
00022 c     from a graph-theoretic standard, this is the same as considering
00023 c     the subgraph induced by sr and sc, if non-matching edges
00024 c     are directed from rows to columns, and matching edges are shrunk 
00025 c     into single vertices, the resulting directed graph has strongly 
00026 c     connected components.
00027 c
00028 c     mmc13e uses Tarjan's algorithm to find the strongly connected
00029 c     components by depth-first search. All the pairs have been visited
00030 c     will be labeled in the order they are visited, and associated a
00031 c     lowlink for each vertex, stored in the stack - lowlnk.
00032 c
00033 c     input variables :
00034 c
00035 c        nrows  -- number of rows in matrix
00036 c        ncols  -- number of columns in matrix
00037 c        nhcols -- number of columns in horizontal (underdetermined)
00038 c                  partition
00039 c        nhrows -- number of rows in horizontal (underdetermined)
00040 c                  partition
00041 c        nscols -- number of rows and columns in square partition
00042 c        sqindx -- index for SR and SC, for rows and columns
00043 c                  in the square partition
00044 c        hrzcmp -- number of components in vertical partition
00045 c        rowstr, colind
00046 c               -- the adjacency structure, stored by rows
00047 c        colset -- the row matched to a column (if any)
00048 c
00049 c     output variables :
00050 c
00051 c        sqcmpn -- number of components in the square partition
00052 c        ccmstr -- global component start vector
00053 c        rcmstr -- global component start vector
00054 c        cnto   -- new to old mapping for columns
00055 c        rnto   -- new to old mapping for rows
00056 c
00057 c     working variables  :
00058 c
00059 c        trycol -- pointer to next unsearched column for this row 
00060 c        cbegin -- is the beginning of the component.
00061 c        colmrk -- column mark vector.
00062 c                  on input, is negative for all columns
00063 c                            = sqindx for columns in sc
00064 c                  used temporarily as a stack to
00065 c                  store the depth-first numbering for each pair.
00066 c                  that is, is the position of pair i in the stack
00067 c                  if it is on it, is the permuted order of pair i for
00068 c                  those pairs whose final position has been found and
00069 c                  is otherwise zero for columns in sc and negative
00070 c                  for all other columns.
00071 c                  on output, is restored to original values
00072 c        lowlnk -- stores the lowlink for each pair.
00073 c                  is the smallest stack position of any pair to which
00074 c                  a path from pair i has been found. it is set to n+1
00075 c                  when pair i is removed from the stack.
00076 c      prev   -- is the pair at the end of the path when pair i was 
00077 c                  placed on the stack
00078 c        
00079 c     ==================================================================
00080 
00081 c     --------------
00082 c     ... parameters
00083 c     --------------
00084 
00085       integer          nrows , ncols , nhcols, nhrows, nscols, 
00086      $                 sqindx, hrzcmp, sqcmpn
00087 
00088       integer          rowstr (nrows+1), colind (*), colset (ncols)
00089 
00090       integer          trycol (nrows), cbegin (nscols), lowlnk (ncols),
00091      $                 prev   (ncols) 
00092 
00093       integer          colmrk (ncols), cnto (ncols), rnto (nrows),
00094      $                 ccmstr (*)    , rcmstr (*)
00095 
00096 c     -------------------
00097 c     ... local variables
00098 c     -------------------
00099 
00100       integer          cmpbeg, col   , compnt, count , passes, fcol  , 
00101      $                 fnlpos, frow  , rootcl, j     , pair  , scol  ,
00102      $                 stackp, xcol
00103 
00104 c     ==================================================================
00105 
00106 c     fnlpos  is the number of pairs whose positions in final ordering
00107 c             have been found.
00108 c     sqcmpn  is the number of components that have been found.
00109 c     count   is the number of pairs on the stack  (stack pointer)
00110 
00111 c     ------------------------------------------------------
00112 c     ... initialization for columns in the square partition
00113 c     ------------------------------------------------------
00114 
00115       fnlpos = 0
00116       sqcmpn = 0
00117 
00118       do 100 col = 1, ncols
00119          if  ( colmrk (col) .eq. sqindx )  then
00120             colmrk (col) = 0
00121          endif
00122  100  continue
00123 
00124       do 200 j = 1, nrows
00125          trycol (j) = rowstr (j)
00126  200  continue
00127 
00128 c     ----------------------------
00129 c     ... look for a starting pair
00130 c     ----------------------------
00131 
00132       do 700 rootcl = 1, ncols
00133 
00134          if  ( colmrk (rootcl) .eq. 0 )  then
00135 
00136 c           -------------------------------------
00137 c           ... put pair (rootcl, colset(rootcl))
00138 c               at beginning of stack
00139 c           -------------------------------------
00140 
00141             fcol            = rootcl
00142             count           = 1
00143             lowlnk (fcol)   = count
00144             colmrk (fcol)   = count
00145             cbegin (nscols) = fcol
00146 
00147 c           --------------------------------------------
00148 c           ... the body of this loop puts a new pair on
00149 c               the stack or backtracks
00150 c           --------------------------------------------
00151 
00152             do 600 passes = 1, 2*nscols - 1
00153 
00154                frow = colset (fcol)
00155                
00156 c              -------------------------------
00157 c              ... have all edges leaving pair  
00158 c                  (frow,fcol)  been searched?
00159 c              -------------------------------
00160 
00161                if  ( trycol (frow) .gt. 0 )  then
00162 
00163 c                 -----------------------------------------------
00164 c                 ... look at edges leaving from row "frow" until
00165 c                     we find a new column "scol" that has not
00166 c                     yet been encountered or until all edges are
00167 c                     exhausted.
00168 c                 -----------------------------------------------
00169 
00170                   do 300 xcol = trycol (frow), rowstr (frow+1)-1
00171                      
00172                      scol = colind (xcol)
00173                      if  ( colmrk (scol) .eq. 0 )  then
00174 
00175 c                       --------------------------------------
00176 c                       ... put new pair  (scol, colset(scol))
00177 c                           on the stack
00178 c                       --------------------------------------
00179 
00180                         trycol (frow)           = xcol + 1
00181                         prev (scol)             = fcol
00182                         fcol                    = scol
00183                         count                   = count + 1 
00184                         lowlnk (fcol)           = count 
00185                         colmrk (fcol)           = count 
00186                         cbegin (nscols+1-count) = fcol 
00187                         go to 600
00188 
00189                      else
00190      $               if  ( colmrk (scol) .gt. 0 )  then
00191 
00192 c                       -------------------------------------------
00193 c                       ... has scol been on stack already?  then
00194 c                           update value of low (fcol) if necessary
00195 c                       -------------------------------------------
00196 
00197                         if  (lowlnk (scol) .lt. lowlnk (fcol))  then
00198                            lowlnk (fcol) = lowlnk (scol)
00199                         endif
00200                      endif
00201  300              continue
00202 
00203 c                 ----------------------------------------
00204 c                 ... there are no more edges leaving frow
00205 c                 ----------------------------------------
00206 
00207                   trycol (frow) = -1
00208 
00209                endif
00210 
00211 c              ------------------------------
00212 c              
00213 c              ------------------------------
00214 
00215                if  ( lowlnk (fcol) .ge. colmrk (fcol) )  then
00216 c
00217 c                 -----------------------------------------------------
00218 c                 ... is  frow  the root of a block?  if so, we have 
00219 c                     found a component.  order the nodes in this
00220 c                     block by starting at the top of the stack and
00221 c                     working down to the root of the block
00222 c                 -----------------------------------------------------
00223 
00224                   sqcmpn = sqcmpn + 1
00225                   cmpbeg = fnlpos + 1
00226                   
00227                   do 400 stackp = nscols + 1 - count, nscols
00228                      pair          = cbegin (stackp)
00229                      fnlpos        = fnlpos + 1
00230                      colmrk (pair) = fnlpos
00231                      count         = count-1
00232                      lowlnk (pair) = nscols + 1
00233                      if  ( pair .eq. fcol )  go to 500
00234  400              continue
00235 
00236 c                 -------------------------------------------------------
00237 c                 ... record the starting position for the new component
00238 c                 -------------------------------------------------------
00239 
00240  500              cbegin (sqcmpn) = cmpbeg 
00241 
00242 c                 --------------------------------------------
00243 c                 ... are there any pairs left on the stack.
00244 c                     if so, backtrack.
00245 c                     if not, have all the pairs been ordered?
00246 c                 --------------------------------------------
00247 
00248                   if  ( count .eq. 0 )  then
00249                      if  ( fnlpos .lt. nscols )  then
00250                         go to 700
00251                      else
00252                         go to 800
00253                      endif
00254                   endif
00255                   
00256                endif
00257 c
00258 c              --------------------------------------
00259 c              ... backtrack to previous pair on path
00260 c              --------------------------------------
00261 
00262                scol = fcol
00263                fcol = prev (fcol)
00264                if  ( lowlnk (scol) .lt. lowlnk (fcol) )  then
00265                   lowlnk (fcol) = lowlnk (scol)
00266                endif
00267 
00268  600         continue
00269             
00270          endif
00271 
00272  700  continue
00273 
00274 c     ----------------------------------------
00275 c     ... put permutation in the required form
00276 c     ----------------------------------------
00277 
00278  800  do 900 compnt = 1, sqcmpn
00279          ccmstr (compnt + hrzcmp) = (cbegin (compnt) + nhcols)
00280          rcmstr (compnt + hrzcmp) = (cbegin (compnt) + nhcols) -
00281      $                              (nhcols - nhrows)
00282  900  continue
00283 
00284       ccmstr (hrzcmp + sqcmpn + 1)  = nhcols + nscols + 1
00285       rcmstr (hrzcmp + sqcmpn + 1)  = nhrows + nscols + 1
00286 
00287 c     ------------------------------------------------------
00288 c     ... note that columns not in the square partition have
00289 c         colmrk set negative.  diagonal entries in the
00290 c         square block all correspond to matching pairs.
00291 c     ------------------------------------------------------
00292 
00293       do 1000 col = 1, ncols
00294          j = colmrk (col)
00295          if  ( j .gt. 0 )  then
00296             cnto (nhcols + j) = col
00297             rnto (nhrows + j) = colset (col)
00298             colmrk (col)      = sqindx
00299          endif
00300  1000 continue
00301 
00302       return
00303 
00304       end
00305 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines