EpetraExt Development
genbtf.f
Go to the documentation of this file.
00001       subroutine   genbtf   ( nrows , ncols , colstr, rowidx, rowstr,
00002      \$                        colidx, w     , rnto  , cnto  , nhrows,
00003      \$                        nhcols, hrzcmp, nsrows, sqcmpn, nvrows,
00004      \$                        nvcols, vrtcmp, rcmstr, ccmstr, msglvl,
00005      \$                        output )
00006
00007 c     ==================================================================
00008 c     ==================================================================
00009 c     ====  genbtf -- find the block triangular form (dulmadge-     ====
00010 c     ====            mendelson decomposition) of a general         ====
00011 c     ====            rectangular sparse matrix                     ====
00012 c     ==================================================================
00013 c     ==================================================================
00014 c
00015 c     created        sept. 14, 1990 (jgl)
00017
00018 c     algorithm by alex pothen and chin-ju fan
00019 c     this code based on code from alex pothen, penn state university
00020
00021 c     ... input variables
00022 c     -------------------
00023 c
00024 c        nrows  -- number of rows in matrix
00025 c        ncols  -- number of columns in matrix
00026 c        colstr, rowidx
00027 c               -- adjacency structure of matrix, where each
00028 c                  column of matrix is stored contiguously
00029 c                  (column-wise representation)
00030 c        rowstr, colidx
00031 c               -- adjacency structure of matrix, where each
00032 c                  row of matrix is stored contiguously
00033 c                  (row-wise representation)
00034 c                  (yes, two copies of the matrix)
00035 c
00036 c     ... temporary storage
00037 c     ---------------------
00038 c
00039 c        w      -- integer array of length 5*nrows + 5*ncols
00040 c
00041 c     ... output variables:
00042 c     ---------------------
00043 c
00044 c        rnto   -- the new to old permutation array for the rows
00045 c        cotn   -- the old to new permutation array for the columns
00046 c        nhrows, nhcols, hrzcmp
00047 c               -- number of rows, columns and connected components
00048 c                  in the horizontal (underdetermined) block
00049 c        nsrows, sqcmpn
00050 c               -- number of rows (and columns) and strong components
00051 c                  in the square (exactly determined) block
00052 c        nvrows, nvcols, vrtcmp
00053 c               -- number of rows, columns and connected components
00054 c                  in the vertical (overdetermined) block
00055 c        rcmstr -- index of first row in a diagonal block
00056 c                  (component starting row)
00057 c                  where
00058 c                      (rcmstr(1), ..., rcmstr(hrzcmp)) give the
00059 c                           indices for the components in the
00060 c                            horizontal block
00061 c                      (rcmstr(hrzcmp+1), ..., rcmstr(hrzcmp+sqcmpn))
00062 c                           give the indices for the components in the
00063 c                           square block
00064 c                      (rcmstr(hrzcmp+sqcmpn+1), ...,
00065 c                       rcmstr(hrzcmp+sqcmpn+vrtcmp)) give the indices
00066 c                           for the components in the vertical block
00067 c                       rcmstr(hrzcmp+sqcmpn+vrtcmp+1) is equal to
00068 c                           nrows+1 for convenience
00069 c        ccmstr -- index of first column in a diagonal block
00070 c                  (component starting column)
00071 c                  where
00072 c                      (ccmstr(1), ..., ccmstr(hrzcmp)) give the
00073 c                           indices for the components in the
00074 c                            horizontal block
00075 c                      (ccmstr(hrzcmp+1), ..., ccmstr(hrzcmp+sqcmpn))
00076 c                           give the indices for the components in the
00077 c                           square block, making this block itself
00078 c                           in block lower triangular form
00079 c                      (ccmstr(hrzcmp+sqcmpn+1), ...,
00080 c                       ccmstr(hrzcmp+sqcmpn+vrtcmp)) give the indices
00081 c                           for the components in the vertical block
00082 c                       ccmstr(hrzcmp+sqcmpn+vrtcmp+1) is equal to
00083 c                           ncols+1 for convenience
00084 c
00085 c                  note -- if the matrix has entirely empty rows,
00086 c                          these rows will be placed in the vertical
00087 c                          block, each as a component with one row
00088 c                          and zero columns.  similarly, entirely
00089 c                          empty columns will appear in the horizontal
00090 c                          block, each as a component with no rows and
00091 c                          one column.
00092 c
00093 c        msglvl -- message level
00094 c                  = 0 -- no output
00095 c                  = 1 -- timing and summary output
00096 c                  = 2 -- adds final permutation vectors
00097 c                  = 3 -- adds intermediate copies of vectros as
00098 c                         debugging output
00099 c
00100 c        output -- fortran unit number for printable output
00101 c
00102 c     efficiency note:
00103 c     ----------------
00104
00105 c        although it is not required by this code that the number
00106 c        of rows be larger than the number of columns, the first
00107 c        phase (the matching) will be faster in this case.  thus,
00108 c        in cases where the number of columns is substantially
00109 c        larger than the number of rows, it will probably be more
00110 c        efficient to apply this algorithm to the transpose of the
00111 c        matrix.  since the matrix is required with both row and
00112 c        column representations, applying the algorithm to the
00113 c        transposed matrix can be achieved simply by interchanging
00114 c        appropriate parameters in the call to  genbtf.
00115 c
00116 c     ==================================================================
00117
00118 c     --------------
00119 c     ... parameters
00120 c     --------------
00121
00122       integer         nrows , ncols , nhrows, nhcols, hrzcmp, nsrows,
00123      \$                sqcmpn, nvrows, nvcols, vrtcmp, msglvl, output
00124
00125       integer         colstr (ncols+1), rowidx (*),
00126      \$                rowstr (nrows+1), colidx (*),
00127      \$                w      (*)      ,
00128      \$                cnto   (ncols)  , rnto   (nrows),
00129      \$                rcmstr (nrows+1), ccmstr (ncols+1)
00130
00131 c     -------------------
00132 c     ... local variables
00133 c     -------------------
00134
00135       integer         cmk   , cmbase, cnbase, cst   , cw1   , cw2   ,
00136      \$                cw3   , i     , hindex, ncompn, nscols, rmk   ,
00137      \$                rnbase, rst   , rw1   , rw2   , rw3   , sqindx,
00138      \$                vindex
00139
00140       real            timeh , timem , times , timev , tmstrt
00141
00142       real            tarray (2)
00143
00144       real            etime
00145
00146 c     ==================================================================
00147
00148 c     --------------
00149 c     ... initialize
00150 c     --------------
00151
00152       vindex = -1
00153       sqindx = -2
00154       hindex = -3
00155
00156       cmk = 1
00157       cst = cmk + ncols
00158       rmk = cst + ncols
00159       rst = rmk + nrows
00160       rw1 = rst + nrows
00161       rw2 = rw1 + nrows
00162       rw3 = rw2 + nrows
00163       cw1 = rw3 + nrows
00164       cw2 = cw1 + ncols
00165       cw3 = cw2 + ncols
00166       call izero ( cw3 + ncols - 1, w, 1 )
00167
00168 c     ---------------------------------------
00169 c     ... algorithm consists of three stages:
00170 c         1.  find a maximum matching
00171 c         2.  find a coarse decomposition
00172 c         3.  find a fine decomposition
00173 c     ---------------------------------------
00174
00175 c     -----------------------------
00176 c     ... find the maximum matching
00177 c     -----------------------------
00178
00179       if  ( msglvl .ge. 1 )  then
00180          tmstrt = etime ( tarray )
00181       endif
00182
00183       call maxmatch ( nrows , ncols , colstr, rowidx, w(cw1), w(cmk),
00184      \$                w(rw2), w(cw2), w(cw3), w(rst), w(cst) )
00185
00186       do 100 i = 1, nrows
00187          w (rmk + i - 1) = sqindx
00188  100  continue
00189
00190       do 200 i = 1, ncols
00191          w (cmk + i - 1) = sqindx
00192  200  continue
00193
00194       if  ( msglvl .ge. 1 )  then
00195          timem = etime ( tarray ) - tmstrt
00196          if  ( msglvl .ge. 3 )  then
00197             call prtivs ( 'rowset', nrows, w(rst), output )
00198             call prtivs ( 'colset', ncols, w(cst), output )
00199          endif
00200       endif
00201
00202 c     ------------------------------------------------------------
00203 c     ... coarse partitioning -- divide the graph into three parts
00204 c     ------------------------------------------------------------
00205
00206 c     --------------------------------------------------------
00207 c     ... depth-first search from unmatched columns to get the
00208 c         horizontal submatrix
00209 c     --------------------------------------------------------
00210
00211       if  ( msglvl .ge. 1 )  then
00212          tmstrt = etime ( tarray )
00213       endif
00214
00215       call rectblk ( nrows , ncols , hindex, sqindx, colstr, rowidx,
00216      \$               w(cst), w(rst), w(cw1), w(cw2), w(cmk), w(rmk),
00217      \$               nhrows, nhcols )
00218
00219       if  ( msglvl .ge. 1 )  then
00220          timeh = etime ( tarray ) - tmstrt
00221          if  ( msglvl .ge. 3 )  then
00222             write ( output, * ) '0nhrows, nhcols', nhrows, nhcols
00223          endif
00224       endif
00225
00226 c     -----------------------------------------------------
00227 c     ... depth-first search from unmatched rows to get the
00228 c         vertical submatrix
00229 c     -----------------------------------------------------
00230
00231       if  ( msglvl .ge. 1 )  then
00232          tmstrt = etime ( tarray )
00233       endif
00234
00235       tmstrt = etime ( tarray )
00236
00237       call rectblk ( ncols , nrows , vindex, sqindx, rowstr, colidx,
00238      \$               w(rst), w(cst), w(rw1), w(rw2), w(rmk), w(cmk),
00239      \$               nvcols, nvrows )
00240
00241       if  ( msglvl .ge. 1 )  then
00242          timev = etime ( tarray ) - tmstrt
00243          if  ( msglvl .ge. 3 )  then
00244             write ( output, * ) '0nvrows, nvcols', nvrows, nvcols
00245          endif
00246       endif
00247
00248 c     ----------------------------------------
00249 c     ... the square submatrix is what is left
00250 c     ----------------------------------------
00251
00252       nscols = ncols - nhcols - nvcols
00253       nsrows = nrows - nhrows - nvrows
00254
00255       if  ( msglvl .ge. 1 )  then
00256          call corsum ( timem , timeh , timev , nhrows, nhcols, nsrows,
00257      \$                 nscols, nvrows, nvcols, output )
00258       endif
00259
00260 c     ----------------------------------------------
00261 c     ... begin the fine partitioning and create the
00262 c         new to old permutation vectors
00263 c     ----------------------------------------------
00264
00265 c     ---------------------------------------------------------
00266 c     ... find connected components in the horizontal submatrix
00267 c     ---------------------------------------------------------
00268
00269       if  ( nhcols .gt. 0 )  then
00270
00271          if  ( msglvl .ge. 1 )  then
00272             tmstrt = etime ( tarray )
00273          endif
00274
00275          cmbase  = 0
00276          rnbase  = 0
00277          cnbase  = 0
00278
00279          call concmp ( cmbase, cnbase, rnbase, hindex, ncols , nrows ,
00280      \$                 nhcols, nhrows, colstr, rowidx, rowstr, colidx,
00281      \$                 w(rw1), w(cw1), w(cw2), w(rw2), w(rw3), w(cw3),
00282      \$                 w(rmk), w(cmk), rcmstr, ccmstr, rnto  , cnto  ,
00283      \$                 hrzcmp )
00284
00285          if  ( msglvl .ge. 1 )  then
00286             timeh = etime ( tarray ) - tmstrt
00287             if  ( msglvl .ge. 3 )  then
00288                write ( output, * ) '0hrzcmp', hrzcmp
00289                call prtivs ( 'rcmstr', hrzcmp + 1, rcmstr, output )
00290                call prtivs ( 'ccmstr', hrzcmp + 1, ccmstr, output )
00291                call prtivs ( 'rnto', nrows, rnto, output )
00292                call prtivs ( 'cnto', ncols, cnto, output )
00293             endif
00294          endif
00295
00296       else
00297
00298          hrzcmp = 0
00299          timeh  = 0.0
00300
00301       endif
00302
00303       if  ( nsrows .gt. 0 ) then
00304
00305          if  ( msglvl .ge. 1 )  then
00306             tmstrt = etime ( tarray )
00307          endif
00308
00309 c        -----------------------------------------------------------
00310 c        ... find strongly connected components in square submatrix,
00311 c            putting this block into block lower triangular form.
00312 c        -----------------------------------------------------------
00313
00314          call mmc13e ( nrows , ncols , nhcols, nhrows, nsrows, sqindx,
00315      \$                 hrzcmp, rowstr, colidx, w(cst), w(rw1), w(rw2),
00316      \$                 w(cw1), w(cw2), w(cmk), ccmstr, rcmstr, cnto  ,
00317      \$                 rnto  , sqcmpn )
00318
00319          if  ( msglvl .ge. 1 )  then
00320             call strchk ( nrows , ncols , colstr, rowidx, nhrows,
00321      \$                    nhcols, nsrows, rnto  , cnto  , w(cst),
00322      \$                    w(rst), output )
00323
00324          endif
00325
00326          if  ( msglvl .ge. 1 )  then
00327             times = etime ( tarray ) - tmstrt
00328             if  ( msglvl .ge. 3 )  then
00329                ncompn = hrzcmp + sqcmpn + 1
00330                write ( output, * ) '0sqcmpn', sqcmpn
00331                call prtivs ( 'rcmstr', ncompn, rcmstr, output )
00332                call prtivs ( 'ccmstr', ncompn, ccmstr, output )
00333                call prtivs ( 'rnto', nrows, rnto, output )
00334                call prtivs ( 'cnto', ncols, cnto, output )
00335             endif
00336          endif
00337
00338       else
00339
00340          sqcmpn = 0
00341          times  = 0.0
00342
00343       endif
00344
00345       if  ( nvrows .gt. 0 )  then
00346
00347          cmbase = hrzcmp + sqcmpn
00348          rnbase = nhrows + nscols
00349          cnbase = nhcols + nscols
00350
00351 c        ---------------------------------------------------
00352 c        ... find connected components in vertical submatrix
00353 c        ---------------------------------------------------
00354
00355          if  ( msglvl .ge. 1 )  then
00356             tmstrt = etime ( tarray )
00357          endif
00358
00359          call concmp ( cmbase, rnbase, cnbase, vindex, nrows , ncols ,
00360      \$                 nvrows, nvcols, rowstr, colidx, colstr, rowidx,
00361      \$                 w(cw1), w(rw1), w(rw2), w(cw2), w(cw3), w(rw3),
00362      \$                 w(cmk), w(rmk), ccmstr, rcmstr, cnto  , rnto  ,
00363      \$                 vrtcmp )
00364
00365          if  ( msglvl .ge. 1 )  then
00366
00367             timev = etime ( tarray ) - tmstrt
00368
00369             if  ( msglvl .ge. 2 )  then
00370                call prtivs ( 'rnto', nrows, rnto, output )
00371                call prtivs ( 'cnto', ncols, cnto, output )
00372
00373                if  ( msglvl .ge. 3 )  then
00374                   ncompn = hrzcmp + sqcmpn + vrtcmp + 1
00375                   write ( output, * ) '0vrtcmp', vrtcmp
00376                   call prtivs ( 'rcmstr', ncompn, rcmstr, output )
00377                   call prtivs ( 'ccmstr', ncompn, ccmstr, output )
00378                endif
00379
00380             endif
00381          endif
00382
00383       else
00384
00385          vrtcmp = 0
00386          timev  = 0.0
00387
00388       endif
00389
00390       if  ( msglvl .ge. 1 )  then
00391          call finsum ( timeh , times , timev , hrzcmp, sqcmpn,
00392      \$                 vrtcmp, ccmstr, rcmstr, output )
00393       endif
00394
00395       return
00396
00397       end