|
EpetraExt
Development
|
|
EpetraExt
Development
|
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) 00016 c last modified oct. 4, 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
1.7.6.1
