EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
mmc13e.f
Go to the documentation of this file.
1  subroutine mmc13e ( nrows , ncols , nhcols, nhrows, nscols,
2  $ sqindx, hrzcmp, rowstr, colind, colset,
3  $ trycol, cbegin, lowlnk, prev , colmrk,
4  $ ccmstr, rcmstr, cnto , rnto , sqcmpn )
5 
6 c ==================================================================
7 c ==================================================================
8 c ==== mmc13e -- lower block triangular form of square matrix ====
9 c ==================================================================
10 c ==================================================================
11 
12 c mmc13e : modified from harwell mc13e by alex pothen and
13 c chin-ju fan
14 c bcs modifications, john lewis, sept. 1990
15 
16 c finds the lower block triangular form of the square submatrix
17 c in the general block triangular form. the square submatrix
18 c consists entirely of matched rows and columns. therefore,
19 c with each row matched to its matching column, the submatrix
20 c has a nonzero diagonal, as required by duff's algorithm.
21 c
22 c from a graph-theoretic standard, this is the same as considering
23 c the subgraph induced by sr and sc, if non-matching edges
24 c are directed from rows to columns, and matching edges are shrunk
25 c into single vertices, the resulting directed graph has strongly
26 c connected components.
27 c
28 c mmc13e uses Tarjan's algorithm to find the strongly connected
29 c components by depth-first search. All the pairs have been visited
30 c will be labeled in the order they are visited, and associated a
31 c lowlink for each vertex, stored in the stack - lowlnk.
32 c
33 c input variables :
34 c
35 c nrows -- number of rows in matrix
36 c ncols -- number of columns in matrix
37 c nhcols -- number of columns in horizontal (underdetermined)
38 c partition
39 c nhrows -- number of rows in horizontal (underdetermined)
40 c partition
41 c nscols -- number of rows and columns in square partition
42 c sqindx -- index for SR and SC, for rows and columns
43 c in the square partition
44 c hrzcmp -- number of components in vertical partition
45 c rowstr, colind
46 c -- the adjacency structure, stored by rows
47 c colset -- the row matched to a column (if any)
48 c
49 c output variables :
50 c
51 c sqcmpn -- number of components in the square partition
52 c ccmstr -- global component start vector
53 c rcmstr -- global component start vector
54 c cnto -- new to old mapping for columns
55 c rnto -- new to old mapping for rows
56 c
57 c working variables :
58 c
59 c trycol -- pointer to next unsearched column for this row
60 c cbegin -- is the beginning of the component.
61 c colmrk -- column mark vector.
62 c on input, is negative for all columns
63 c = sqindx for columns in sc
64 c used temporarily as a stack to
65 c store the depth-first numbering for each pair.
66 c that is, is the position of pair i in the stack
67 c if it is on it, is the permuted order of pair i for
68 c those pairs whose final position has been found and
69 c is otherwise zero for columns in sc and negative
70 c for all other columns.
71 c on output, is restored to original values
72 c lowlnk -- stores the lowlink for each pair.
73 c is the smallest stack position of any pair to which
74 c a path from pair i has been found. it is set to n+1
75 c when pair i is removed from the stack.
76 c prev -- is the pair at the end of the path when pair i was
77 c placed on the stack
78 c
79 c ==================================================================
80 
81 c --------------
82 c ... parameters
83 c --------------
84 
85  integer nrows , ncols , nhcols, nhrows, nscols,
86  $ sqindx, hrzcmp, sqcmpn
87 
88  integer rowstr (nrows+1), colind (*), colset (ncols)
89 
90  integer trycol (nrows), cbegin (nscols), lowlnk (ncols),
91  $ prev(ncols)
92 
93  integer colmrk (ncols), cnto (ncols), rnto (nrows),
94  $ ccmstr(*) , rcmstr(*)
95 
96 c -------------------
97 c ... local variables
98 c -------------------
99 
100  integer cmpbeg, col , compnt, count , passes, fcol ,
101  $ fnlpos, frow , rootcl, j , pair , scol ,
102  $ stackp, xcol
103 
104 c ==================================================================
105 
106 c fnlpos is the number of pairs whose positions in final ordering
107 c have been found.
108 c sqcmpn is the number of components that have been found.
109 c count is the number of pairs on the stack (stack pointer)
110 
111 c ------------------------------------------------------
112 c ... initialization for columns in the square partition
113 c ------------------------------------------------------
114 
115  fnlpos = 0
116  sqcmpn = 0
117 
118  do 100 col = 1, ncols
119  if ( colmrk(col) .eq. sqindx ) then
120  colmrk(col) = 0
121  endif
122  100 continue
123 
124  do 200 j = 1, nrows
125  trycol(j) = rowstr(j)
126  200 continue
127 
128 c ----------------------------
129 c ... look for a starting pair
130 c ----------------------------
131 
132  do 700 rootcl = 1, ncols
133 
134  if ( colmrk(rootcl) .eq. 0 ) then
135 
136 c -------------------------------------
137 c ... put pair (rootcl, colset(rootcl))
138 c at beginning of stack
139 c -------------------------------------
140 
141  fcol = rootcl
142  count = 1
143  lowlnk(fcol) = count
144  colmrk(fcol) = count
145  cbegin(nscols) = fcol
146 
147 c --------------------------------------------
148 c ... the body of this loop puts a new pair on
149 c the stack or backtracks
150 c --------------------------------------------
151 
152  do 600 passes = 1, 2*nscols - 1
153 
154  frow = colset(fcol)
155 
156 c -------------------------------
157 c ... have all edges leaving pair
158 c (frow,fcol) been searched?
159 c -------------------------------
160 
161  if ( trycol(frow) .gt. 0 ) then
162 
163 c -----------------------------------------------
164 c ... look at edges leaving from row "frow" until
165 c we find a new column "scol" that has not
166 c yet been encountered or until all edges are
167 c exhausted.
168 c -----------------------------------------------
169 
170  do 300 xcol = trycol(frow), rowstr(frow+1)-1
171 
172  scol = colind(xcol)
173  if ( colmrk(scol) .eq. 0 ) then
174 
175 c --------------------------------------
176 c ... put new pair (scol, colset(scol))
177 c on the stack
178 c --------------------------------------
179 
180  trycol(frow) = xcol + 1
181  prev(scol) = fcol
182  fcol = scol
183  count = count + 1
184  lowlnk(fcol) = count
185  colmrk(fcol) = count
186  cbegin(nscols+1-count) = fcol
187  go to 600
188 
189  else
190  $ if ( colmrk(scol) .gt. 0 ) then
191 
192 c -------------------------------------------
193 c ... has scol been on stack already? then
194 c update value of low (fcol) if necessary
195 c -------------------------------------------
196 
197  if (lowlnk(scol) .lt. lowlnk(fcol)) then
198  lowlnk(fcol) = lowlnk(scol)
199  endif
200  endif
201  300 continue
202 
203 c ----------------------------------------
204 c ... there are no more edges leaving frow
205 c ----------------------------------------
206 
207  trycol(frow) = -1
208 
209  endif
210 
211 c ------------------------------
212 c
213 c ------------------------------
214 
215  if ( lowlnk(fcol) .ge. colmrk(fcol) ) then
216 c
217 c -----------------------------------------------------
218 c ... is frow the root of a block? if so, we have
219 c found a component. order the nodes in this
220 c block by starting at the top of the stack and
221 c working down to the root of the block
222 c -----------------------------------------------------
223 
224  sqcmpn = sqcmpn + 1
225  cmpbeg = fnlpos + 1
226 
227  do 400 stackp = nscols + 1 - count, nscols
228  pair = cbegin(stackp)
229  fnlpos = fnlpos + 1
230  colmrk(pair) = fnlpos
231  count = count-1
232  lowlnk(pair) = nscols + 1
233  if ( pair .eq. fcol ) go to 500
234  400 continue
235 
236 c -------------------------------------------------------
237 c ... record the starting position for the new component
238 c -------------------------------------------------------
239 
240  500 cbegin(sqcmpn) = cmpbeg
241 
242 c --------------------------------------------
243 c ... are there any pairs left on the stack.
244 c if so, backtrack.
245 c if not, have all the pairs been ordered?
246 c --------------------------------------------
247 
248  if ( count .eq. 0 ) then
249  if ( fnlpos .lt. nscols ) then
250  go to 700
251  else
252  go to 800
253  endif
254  endif
255 
256  endif
257 c
258 c --------------------------------------
259 c ... backtrack to previous pair on path
260 c --------------------------------------
261 
262  scol = fcol
263  fcol = prev(fcol)
264  if ( lowlnk(scol) .lt. lowlnk(fcol) ) then
265  lowlnk(fcol) = lowlnk(scol)
266  endif
267 
268  600 continue
269 
270  endif
271 
272  700 continue
273 
274 c ----------------------------------------
275 c ... put permutation in the required form
276 c ----------------------------------------
277 
278  800 do 900 compnt = 1, sqcmpn
279  ccmstr(compnt + hrzcmp) = (cbegin(compnt) + nhcols)
280  rcmstr(compnt + hrzcmp) = (cbegin(compnt) + nhcols) -
281  $ (nhcols - nhrows)
282  900 continue
283 
284  ccmstr(hrzcmp + sqcmpn + 1) = nhcols + nscols + 1
285  rcmstr(hrzcmp + sqcmpn + 1) = nhrows + nscols + 1
286 
287 c ------------------------------------------------------
288 c ... note that columns not in the square partition have
289 c colmrk set negative. diagonal entries in the
290 c square block all correspond to matching pairs.
291 c ------------------------------------------------------
292 
293  do 1000 col = 1, ncols
294  j = colmrk(col)
295  if ( j .gt. 0 ) then
296  cnto(nhcols + j) = col
297  rnto(nhrows + j) = colset(col)
298  colmrk(col) = sqindx
299  endif
300  1000 continue
301 
302  return
303 
304  end
305 
mmc13e
subroutine mmc13e(nrows, ncols, nhcols, nhrows, nscols, sqindx, hrzcmp, rowstr, colind, colset, trycol, cbegin, lowlnk, prev, colmrk, ccmstr, rcmstr, cnto, rnto, sqcmpn)
Definition: mmc13e.f:1
Generated by   doxygen 1.8.5