EpetraExt Development
Loading...
Searching...
No Matches
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
6c ==================================================================
7c ==================================================================
8c ==== mmc13e -- lower block triangular form of square matrix ====
9c ==================================================================
10c ==================================================================
11
12c mmc13e : modified from harwell mc13e by alex pothen and
13c chin-ju fan
14c bcs modifications, john lewis, sept. 1990
15
16c finds the lower block triangular form of the square submatrix
17c in the general block triangular form. the square submatrix
18c consists entirely of matched rows and columns. therefore,
19c with each row matched to its matching column, the submatrix
20c has a nonzero diagonal, as required by duff's algorithm.
21c
22c from a graph-theoretic standard, this is the same as considering
23c the subgraph induced by sr and sc, if non-matching edges
24c are directed from rows to columns, and matching edges are shrunk
25c into single vertices, the resulting directed graph has strongly
26c connected components.
27c
28c mmc13e uses Tarjan's algorithm to find the strongly connected
29c components by depth-first search. All the pairs have been visited
30c will be labeled in the order they are visited, and associated a
31c lowlink for each vertex, stored in the stack - lowlnk.
32c
33c input variables :
34c
35c nrows -- number of rows in matrix
36c ncols -- number of columns in matrix
37c nhcols -- number of columns in horizontal (underdetermined)
38c partition
39c nhrows -- number of rows in horizontal (underdetermined)
40c partition
41c nscols -- number of rows and columns in square partition
42c sqindx -- index for SR and SC, for rows and columns
43c in the square partition
44c hrzcmp -- number of components in vertical partition
45c rowstr, colind
46c -- the adjacency structure, stored by rows
47c colset -- the row matched to a column (if any)
48c
49c output variables :
50c
51c sqcmpn -- number of components in the square partition
52c ccmstr -- global component start vector
53c rcmstr -- global component start vector
54c cnto -- new to old mapping for columns
55c rnto -- new to old mapping for rows
56c
57c working variables :
58c
59c trycol -- pointer to next unsearched column for this row
60c cbegin -- is the beginning of the component.
61c colmrk -- column mark vector.
62c on input, is negative for all columns
63c = sqindx for columns in sc
64c used temporarily as a stack to
65c store the depth-first numbering for each pair.
66c that is, is the position of pair i in the stack
67c if it is on it, is the permuted order of pair i for
68c those pairs whose final position has been found and
69c is otherwise zero for columns in sc and negative
70c for all other columns.
71c on output, is restored to original values
72c lowlnk -- stores the lowlink for each pair.
73c is the smallest stack position of any pair to which
74c a path from pair i has been found. it is set to n+1
75c when pair i is removed from the stack.
76c prev -- is the pair at the end of the path when pair i was
77c placed on the stack
78c
79c ==================================================================
80
81c --------------
82c ... parameters
83c --------------
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
96c -------------------
97c ... local variables
98c -------------------
99
100 integer cmpbeg, col , compnt, count , passes, fcol ,
101 $ fnlpos, frow , rootcl, j , pair , scol ,
102 $ stackp, xcol
103
104c ==================================================================
105
106c fnlpos is the number of pairs whose positions in final ordering
107c have been found.
108c sqcmpn is the number of components that have been found.
109c count is the number of pairs on the stack (stack pointer)
110
111c ------------------------------------------------------
112c ... initialization for columns in the square partition
113c ------------------------------------------------------
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
128c ----------------------------
129c ... look for a starting pair
130c ----------------------------
131
132 do 700 rootcl = 1, ncols
133
134 if ( colmrk(rootcl) .eq. 0 ) then
135
136c -------------------------------------
137c ... put pair (rootcl, colset(rootcl))
138c at beginning of stack
139c -------------------------------------
140
141 fcol = rootcl
142 count = 1
143 lowlnk(fcol) = count
144 colmrk(fcol) = count
145 cbegin(nscols) = fcol
146
147c --------------------------------------------
148c ... the body of this loop puts a new pair on
149c the stack or backtracks
150c --------------------------------------------
151
152 do 600 passes = 1, 2*nscols - 1
153
154 frow = colset(fcol)
155
156c -------------------------------
157c ... have all edges leaving pair
158c (frow,fcol) been searched?
159c -------------------------------
160
161 if ( trycol(frow) .gt. 0 ) then
162
163c -----------------------------------------------
164c ... look at edges leaving from row "frow" until
165c we find a new column "scol" that has not
166c yet been encountered or until all edges are
167c exhausted.
168c -----------------------------------------------
169
170 do 300 xcol = trycol(frow), rowstr(frow+1)-1
171
172 scol = colind(xcol)
173 if ( colmrk(scol) .eq. 0 ) then
174
175c --------------------------------------
176c ... put new pair (scol, colset(scol))
177c on the stack
178c --------------------------------------
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
192c -------------------------------------------
193c ... has scol been on stack already? then
194c update value of low (fcol) if necessary
195c -------------------------------------------
196
197 if (lowlnk(scol) .lt. lowlnk(fcol)) then
198 lowlnk(fcol) = lowlnk(scol)
199 endif
200 endif
201 300 continue
202
203c ----------------------------------------
204c ... there are no more edges leaving frow
205c ----------------------------------------
206
207 trycol(frow) = -1
208
209 endif
210
211c ------------------------------
212c
213c ------------------------------
214
215 if ( lowlnk(fcol) .ge. colmrk(fcol) ) then
216c
217c -----------------------------------------------------
218c ... is frow the root of a block? if so, we have
219c found a component. order the nodes in this
220c block by starting at the top of the stack and
221c working down to the root of the block
222c -----------------------------------------------------
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
236c -------------------------------------------------------
237c ... record the starting position for the new component
238c -------------------------------------------------------
239
240 500 cbegin(sqcmpn) = cmpbeg
241
242c --------------------------------------------
243c ... are there any pairs left on the stack.
244c if so, backtrack.
245c if not, have all the pairs been ordered?
246c --------------------------------------------
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
257c
258c --------------------------------------
259c ... backtrack to previous pair on path
260c --------------------------------------
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
274c ----------------------------------------
275c ... put permutation in the required form
276c ----------------------------------------
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
287c ------------------------------------------------------
288c ... note that columns not in the square partition have
289c colmrk set negative. diagonal entries in the
290c square block all correspond to matching pairs.
291c ------------------------------------------------------
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
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:5