|
EpetraExt
Development
|
|
EpetraExt
Development
|
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
1.7.6.1
