EpetraExt  Development
genbtf.f
Go to the documentation of this file.
1  subroutine genbtf ( nrows , ncols , colstr, rowidx, rowstr,
2  $ colidx, w , rnto , cnto , nhrows,
3  $ nhcols, hrzcmp, nsrows, sqcmpn, nvrows,
4  $ nvcols, vrtcmp, rcmstr, ccmstr, msglvl,
5  $ output )
6 
7 c ==================================================================
8 c ==================================================================
9 c ==== genbtf -- find the block triangular form (dulmadge- ====
10 c ==== mendelson decomposition) of a general ====
11 c ==== rectangular sparse matrix ====
12 c ==================================================================
13 c ==================================================================
14 c
15 c created sept. 14, 1990 (jgl)
16 c last modified oct. 4, 1990 (jgl)
17 
18 c algorithm by alex pothen and chin-ju fan
19 c this code based on code from alex pothen, penn state university
20 
21 c ... input variables
22 c -------------------
23 c
24 c nrows -- number of rows in matrix
25 c ncols -- number of columns in matrix
26 c colstr, rowidx
27 c -- adjacency structure of matrix, where each
28 c column of matrix is stored contiguously
29 c (column-wise representation)
30 c rowstr, colidx
31 c -- adjacency structure of matrix, where each
32 c row of matrix is stored contiguously
33 c (row-wise representation)
34 c (yes, two copies of the matrix)
35 c
36 c ... temporary storage
37 c ---------------------
38 c
39 c w -- integer array of length 5*nrows + 5*ncols
40 c
41 c ... output variables:
42 c ---------------------
43 c
44 c rnto -- the new to old permutation array for the rows
45 c cotn -- the old to new permutation array for the columns
46 c nhrows, nhcols, hrzcmp
47 c -- number of rows, columns and connected components
48 c in the horizontal (underdetermined) block
49 c nsrows, sqcmpn
50 c -- number of rows (and columns) and strong components
51 c in the square (exactly determined) block
52 c nvrows, nvcols, vrtcmp
53 c -- number of rows, columns and connected components
54 c in the vertical (overdetermined) block
55 c rcmstr -- index of first row in a diagonal block
56 c (component starting row)
57 c where
58 c (rcmstr(1), ..., rcmstr(hrzcmp)) give the
59 c indices for the components in the
60 c horizontal block
61 c (rcmstr(hrzcmp+1), ..., rcmstr(hrzcmp+sqcmpn))
62 c give the indices for the components in the
63 c square block
64 c (rcmstr(hrzcmp+sqcmpn+1), ...,
65 c rcmstr(hrzcmp+sqcmpn+vrtcmp)) give the indices
66 c for the components in the vertical block
67 c rcmstr(hrzcmp+sqcmpn+vrtcmp+1) is equal to
68 c nrows+1 for convenience
69 c ccmstr -- index of first column in a diagonal block
70 c (component starting column)
71 c where
72 c (ccmstr(1), ..., ccmstr(hrzcmp)) give the
73 c indices for the components in the
74 c horizontal block
75 c (ccmstr(hrzcmp+1), ..., ccmstr(hrzcmp+sqcmpn))
76 c give the indices for the components in the
77 c square block, making this block itself
78 c in block lower triangular form
79 c (ccmstr(hrzcmp+sqcmpn+1), ...,
80 c ccmstr(hrzcmp+sqcmpn+vrtcmp)) give the indices
81 c for the components in the vertical block
82 c ccmstr(hrzcmp+sqcmpn+vrtcmp+1) is equal to
83 c ncols+1 for convenience
84 c
85 c note -- if the matrix has entirely empty rows,
86 c these rows will be placed in the vertical
87 c block, each as a component with one row
88 c and zero columns. similarly, entirely
89 c empty columns will appear in the horizontal
90 c block, each as a component with no rows and
91 c one column.
92 c
93 c msglvl -- message level
94 c = 0 -- no output
95 c = 1 -- timing and summary output
96 c = 2 -- adds final permutation vectors
97 c = 3 -- adds intermediate copies of vectros as
98 c debugging output
99 c
100 c output -- fortran unit number for printable output
101 c
102 c efficiency note:
103 c ----------------
104 
105 c although it is not required by this code that the number
106 c of rows be larger than the number of columns, the first
107 c phase (the matching) will be faster in this case. thus,
108 c in cases where the number of columns is substantially
109 c larger than the number of rows, it will probably be more
110 c efficient to apply this algorithm to the transpose of the
111 c matrix. since the matrix is required with both row and
112 c column representations, applying the algorithm to the
113 c transposed matrix can be achieved simply by interchanging
114 c appropriate parameters in the call to genbtf.
115 c
116 c ==================================================================
117 
118 c --------------
119 c ... parameters
120 c --------------
121 
122  integer nrows , ncols , nhrows, nhcols, hrzcmp, nsrows,
123  $ sqcmpn, nvrows, nvcols, vrtcmp, msglvl, output
124 
125  integer colstr (ncols+1), rowidx (*),
126  $ rowstr(nrows+1), colidx(*),
127  $ w(*) ,
128  $ cnto(ncols) , rnto(nrows),
129  $ rcmstr(nrows+1), ccmstr(ncols+1)
130 
131 c -------------------
132 c ... local variables
133 c -------------------
134 
135  integer cmk , cmbase, cnbase, cst , cw1 , cw2 ,
136  $ cw3 , i , hindex, ncompn, nscols, rmk ,
137  $ rnbase, rst , rw1 , rw2 , rw3 , sqindx,
138  $ vindex
139 
140  real timeh , timem , times , timev , tmstrt
141 
142  real tarray (2)
143 
144  real etime
145 
146 c ==================================================================
147 
148 c --------------
149 c ... initialize
150 c --------------
151 
152  vindex = -1
153  sqindx = -2
154  hindex = -3
155 
156  cmk = 1
157  cst = cmk + ncols
158  rmk = cst + ncols
159  rst = rmk + nrows
160  rw1 = rst + nrows
161  rw2 = rw1 + nrows
162  rw3 = rw2 + nrows
163  cw1 = rw3 + nrows
164  cw2 = cw1 + ncols
165  cw3 = cw2 + ncols
166  call izero ( cw3 + ncols - 1, w, 1 )
167 
168 c ---------------------------------------
169 c ... algorithm consists of three stages:
170 c 1. find a maximum matching
171 c 2. find a coarse decomposition
172 c 3. find a fine decomposition
173 c ---------------------------------------
174 
175 c -----------------------------
176 c ... find the maximum matching
177 c -----------------------------
178 
179  if ( msglvl .ge. 1 ) then
180  tmstrt = etime( tarray )
181  endif
182 
183  call maxmatch ( nrows , ncols , colstr, rowidx, w(cw1), w(cmk),
184  $ w(rw2), w(cw2), w(cw3), w(rst), w(cst) )
185 
186  do 100 i = 1, nrows
187  w(rmk + i - 1) = sqindx
188  100 continue
189 
190  do 200 i = 1, ncols
191  w(cmk + i - 1) = sqindx
192  200 continue
193 
194  if ( msglvl .ge. 1 ) then
195  timem = etime( tarray ) - tmstrt
196  if ( msglvl .ge. 3 ) then
197  call prtivs ( 'rowset', nrows, w(rst), output )
198  call prtivs ( 'colset', ncols, w(cst), output )
199  endif
200  endif
201 
202 c ------------------------------------------------------------
203 c ... coarse partitioning -- divide the graph into three parts
204 c ------------------------------------------------------------
205 
206 c --------------------------------------------------------
207 c ... depth-first search from unmatched columns to get the
208 c horizontal submatrix
209 c --------------------------------------------------------
210 
211  if ( msglvl .ge. 1 ) then
212  tmstrt = etime( tarray )
213  endif
214 
215  call rectblk ( nrows , ncols , hindex, sqindx, colstr, rowidx,
216  $ w(cst), w(rst), w(cw1), w(cw2), w(cmk), w(rmk),
217  $ nhrows, nhcols )
218 
219  if ( msglvl .ge. 1 ) then
220  timeh = etime( tarray ) - tmstrt
221  if ( msglvl .ge. 3 ) then
222  write ( output, * ) '0nhrows, nhcols', nhrows, nhcols
223  endif
224  endif
225 
226 c -----------------------------------------------------
227 c ... depth-first search from unmatched rows to get the
228 c vertical submatrix
229 c -----------------------------------------------------
230 
231  if ( msglvl .ge. 1 ) then
232  tmstrt = etime( tarray )
233  endif
234 
235  tmstrt = etime( tarray )
236 
237  call rectblk ( ncols , nrows , vindex, sqindx, rowstr, colidx,
238  $ w(rst), w(cst), w(rw1), w(rw2), w(rmk), w(cmk),
239  $ nvcols, nvrows )
240 
241  if ( msglvl .ge. 1 ) then
242  timev = etime( tarray ) - tmstrt
243  if ( msglvl .ge. 3 ) then
244  write ( output, * ) '0nvrows, nvcols', nvrows, nvcols
245  endif
246  endif
247 
248 c ----------------------------------------
249 c ... the square submatrix is what is left
250 c ----------------------------------------
251 
252  nscols = ncols - nhcols - nvcols
253  nsrows = nrows - nhrows - nvrows
254 
255  if ( msglvl .ge. 1 ) then
256  call corsum ( timem , timeh , timev , nhrows, nhcols, nsrows,
257  $ nscols, nvrows, nvcols, output )
258  endif
259 
260 c ----------------------------------------------
261 c ... begin the fine partitioning and create the
262 c new to old permutation vectors
263 c ----------------------------------------------
264 
265 c ---------------------------------------------------------
266 c ... find connected components in the horizontal submatrix
267 c ---------------------------------------------------------
268 
269  if ( nhcols .gt. 0 ) then
270 
271  if ( msglvl .ge. 1 ) then
272  tmstrt = etime( tarray )
273  endif
274 
275  cmbase = 0
276  rnbase = 0
277  cnbase = 0
278 
279  call concmp ( cmbase, cnbase, rnbase, hindex, ncols , nrows ,
280  $ nhcols, nhrows, colstr, rowidx, rowstr, colidx,
281  $ w(rw1), w(cw1), w(cw2), w(rw2), w(rw3), w(cw3),
282  $ w(rmk), w(cmk), rcmstr, ccmstr, rnto , cnto ,
283  $ hrzcmp )
284 
285  if ( msglvl .ge. 1 ) then
286  timeh = etime( tarray ) - tmstrt
287  if ( msglvl .ge. 3 ) then
288  write ( output, * ) '0hrzcmp', hrzcmp
289  call prtivs ( 'rcmstr', hrzcmp + 1, rcmstr, output )
290  call prtivs ( 'ccmstr', hrzcmp + 1, ccmstr, output )
291  call prtivs ( 'rnto', nrows, rnto, output )
292  call prtivs ( 'cnto', ncols, cnto, output )
293  endif
294  endif
295 
296  else
297 
298  hrzcmp = 0
299  timeh = 0.0
300 
301  endif
302 
303  if ( nsrows .gt. 0 ) then
304 
305  if ( msglvl .ge. 1 ) then
306  tmstrt = etime( tarray )
307  endif
308 
309 c -----------------------------------------------------------
310 c ... find strongly connected components in square submatrix,
311 c putting this block into block lower triangular form.
312 c -----------------------------------------------------------
313 
314  call mmc13e ( nrows , ncols , nhcols, nhrows, nsrows, sqindx,
315  $ hrzcmp, rowstr, colidx, w(cst), w(rw1), w(rw2),
316  $ w(cw1), w(cw2), w(cmk), ccmstr, rcmstr, cnto ,
317  $ rnto , sqcmpn )
318 
319  if ( msglvl .ge. 1 ) then
320  call strchk ( nrows , ncols , colstr, rowidx, nhrows,
321  $ nhcols, nsrows, rnto , cnto , w(cst),
322  $ w(rst), output )
323 
324  endif
325 
326  if ( msglvl .ge. 1 ) then
327  times = etime( tarray ) - tmstrt
328  if ( msglvl .ge. 3 ) then
329  ncompn = hrzcmp + sqcmpn + 1
330  write ( output, * ) '0sqcmpn', sqcmpn
331  call prtivs ( 'rcmstr', ncompn, rcmstr, output )
332  call prtivs ( 'ccmstr', ncompn, ccmstr, output )
333  call prtivs ( 'rnto', nrows, rnto, output )
334  call prtivs ( 'cnto', ncols, cnto, output )
335  endif
336  endif
337 
338  else
339 
340  sqcmpn = 0
341  times = 0.0
342 
343  endif
344 
345  if ( nvrows .gt. 0 ) then
346 
347  cmbase = hrzcmp + sqcmpn
348  rnbase = nhrows + nscols
349  cnbase = nhcols + nscols
350 
351 c ---------------------------------------------------
352 c ... find connected components in vertical submatrix
353 c ---------------------------------------------------
354 
355  if ( msglvl .ge. 1 ) then
356  tmstrt = etime( tarray )
357  endif
358 
359  call concmp ( cmbase, rnbase, cnbase, vindex, nrows , ncols ,
360  $ nvrows, nvcols, rowstr, colidx, colstr, rowidx,
361  $ w(cw1), w(rw1), w(rw2), w(cw2), w(cw3), w(rw3),
362  $ w(cmk), w(rmk), ccmstr, rcmstr, cnto , rnto ,
363  $ vrtcmp )
364 
365  if ( msglvl .ge. 1 ) then
366 
367  timev = etime( tarray ) - tmstrt
368 
369  if ( msglvl .ge. 2 ) then
370  call prtivs ( 'rnto', nrows, rnto, output )
371  call prtivs ( 'cnto', ncols, cnto, output )
372 
373  if ( msglvl .ge. 3 ) then
374  ncompn = hrzcmp + sqcmpn + vrtcmp + 1
375  write ( output, * ) '0vrtcmp', vrtcmp
376  call prtivs ( 'rcmstr', ncompn, rcmstr, output )
377  call prtivs ( 'ccmstr', ncompn, ccmstr, output )
378  endif
379 
380  endif
381  endif
382 
383  else
384 
385  vrtcmp = 0
386  timev = 0.0
387 
388  endif
389 
390  if ( msglvl .ge. 1 ) then
391  call finsum ( timeh , times , timev , hrzcmp, sqcmpn,
392  $ vrtcmp, ccmstr, rcmstr, output )
393  endif
394 
395  return
396 
397  end
concmp
subroutine concmp(cmbase, rnbase, cnbase, vindex, nrows, ncols, nvrows, nvcols, rowstr, colidx, colstr, rowidx, predrw, nextrw, predcl, nextcl, ctab, rtab, colmrk, rowmrk, cmclad, cmrwad, cnto, rnto, numcmp)
Definition: concmp.f:6
corsum
subroutine corsum(tmmtch, timhrz, timvrt, nhrows, nhcols, nsrows, nscols, nvrows, nvcols, output)
Definition: corsum.f:3
finsum
subroutine finsum(timhrz, timesq, timvrt, hrzcmp, sqcmpn, vrtcmp, ccmstr, rcmstr, output)
Definition: finsum.f:3
genbtf
subroutine genbtf(nrows, ncols, colstr, rowidx, rowstr, colidx, w, rnto, cnto, nhrows, nhcols, hrzcmp, nsrows, sqcmpn, nvrows, nvcols, vrtcmp, rcmstr, ccmstr, msglvl, output)
Definition: genbtf.f:6
izero
subroutine izero(n, x, incx)
Definition: izero.f:2
maxmatch
subroutine maxmatch(nrows, ncols, colstr, rowind, prevcl, prevrw, marker, tryrow, nxtchp, rowset, colset)
Definition: maxmatch.f:4
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:5
prtivs
subroutine prtivs(title, n, x, output)
Definition: prtivs.f:2
rectblk
subroutine rectblk(nrows, ncols, marked, unmrkd, colstr, rowidx, colset, rowset, prevcl, tryrow, colmrk, rowmrk, nhrows, nhcols)
Definition: rectblk.f:4
strchk
subroutine strchk(nrows, ncols, colstr, rowidx, nhrows, nhcols, nsrows, rnto, cnto, colset, rowset, output)
Definition: strchk.f:4
Generated by doxygen 1.9.1