EpetraExt Development
Loading...
Searching...
No Matches
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
7c ==================================================================
8c ==================================================================
9c ==== genbtf -- find the block triangular form (dulmadge- ====
10c ==== mendelson decomposition) of a general ====
11c ==== rectangular sparse matrix ====
12c ==================================================================
13c ==================================================================
14c
15c created sept. 14, 1990 (jgl)
16c last modified oct. 4, 1990 (jgl)
17
18c algorithm by alex pothen and chin-ju fan
19c this code based on code from alex pothen, penn state university
20
21c ... input variables
22c -------------------
23c
24c nrows -- number of rows in matrix
25c ncols -- number of columns in matrix
26c colstr, rowidx
27c -- adjacency structure of matrix, where each
28c column of matrix is stored contiguously
29c (column-wise representation)
30c rowstr, colidx
31c -- adjacency structure of matrix, where each
32c row of matrix is stored contiguously
33c (row-wise representation)
34c (yes, two copies of the matrix)
35c
36c ... temporary storage
37c ---------------------
38c
39c w -- integer array of length 5*nrows + 5*ncols
40c
41c ... output variables:
42c ---------------------
43c
44c rnto -- the new to old permutation array for the rows
45c cotn -- the old to new permutation array for the columns
46c nhrows, nhcols, hrzcmp
47c -- number of rows, columns and connected components
48c in the horizontal (underdetermined) block
49c nsrows, sqcmpn
50c -- number of rows (and columns) and strong components
51c in the square (exactly determined) block
52c nvrows, nvcols, vrtcmp
53c -- number of rows, columns and connected components
54c in the vertical (overdetermined) block
55c rcmstr -- index of first row in a diagonal block
56c (component starting row)
57c where
58c (rcmstr(1), ..., rcmstr(hrzcmp)) give the
59c indices for the components in the
60c horizontal block
61c (rcmstr(hrzcmp+1), ..., rcmstr(hrzcmp+sqcmpn))
62c give the indices for the components in the
63c square block
64c (rcmstr(hrzcmp+sqcmpn+1), ...,
65c rcmstr(hrzcmp+sqcmpn+vrtcmp)) give the indices
66c for the components in the vertical block
67c rcmstr(hrzcmp+sqcmpn+vrtcmp+1) is equal to
68c nrows+1 for convenience
69c ccmstr -- index of first column in a diagonal block
70c (component starting column)
71c where
72c (ccmstr(1), ..., ccmstr(hrzcmp)) give the
73c indices for the components in the
74c horizontal block
75c (ccmstr(hrzcmp+1), ..., ccmstr(hrzcmp+sqcmpn))
76c give the indices for the components in the
77c square block, making this block itself
78c in block lower triangular form
79c (ccmstr(hrzcmp+sqcmpn+1), ...,
80c ccmstr(hrzcmp+sqcmpn+vrtcmp)) give the indices
81c for the components in the vertical block
82c ccmstr(hrzcmp+sqcmpn+vrtcmp+1) is equal to
83c ncols+1 for convenience
84c
85c note -- if the matrix has entirely empty rows,
86c these rows will be placed in the vertical
87c block, each as a component with one row
88c and zero columns. similarly, entirely
89c empty columns will appear in the horizontal
90c block, each as a component with no rows and
91c one column.
92c
93c msglvl -- message level
94c = 0 -- no output
95c = 1 -- timing and summary output
96c = 2 -- adds final permutation vectors
97c = 3 -- adds intermediate copies of vectros as
98c debugging output
99c
100c output -- fortran unit number for printable output
101c
102c efficiency note:
103c ----------------
104
105c although it is not required by this code that the number
106c of rows be larger than the number of columns, the first
107c phase (the matching) will be faster in this case. thus,
108c in cases where the number of columns is substantially
109c larger than the number of rows, it will probably be more
110c efficient to apply this algorithm to the transpose of the
111c matrix. since the matrix is required with both row and
112c column representations, applying the algorithm to the
113c transposed matrix can be achieved simply by interchanging
114c appropriate parameters in the call to genbtf.
115c
116c ==================================================================
117
118c --------------
119c ... parameters
120c --------------
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
131c -------------------
132c ... local variables
133c -------------------
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
146c ==================================================================
147
148c --------------
149c ... initialize
150c --------------
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
168c ---------------------------------------
169c ... algorithm consists of three stages:
170c 1. find a maximum matching
171c 2. find a coarse decomposition
172c 3. find a fine decomposition
173c ---------------------------------------
174
175c -----------------------------
176c ... find the maximum matching
177c -----------------------------
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
202c ------------------------------------------------------------
203c ... coarse partitioning -- divide the graph into three parts
204c ------------------------------------------------------------
205
206c --------------------------------------------------------
207c ... depth-first search from unmatched columns to get the
208c horizontal submatrix
209c --------------------------------------------------------
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
226c -----------------------------------------------------
227c ... depth-first search from unmatched rows to get the
228c vertical submatrix
229c -----------------------------------------------------
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
248c ----------------------------------------
249c ... the square submatrix is what is left
250c ----------------------------------------
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
260c ----------------------------------------------
261c ... begin the fine partitioning and create the
262c new to old permutation vectors
263c ----------------------------------------------
264
265c ---------------------------------------------------------
266c ... find connected components in the horizontal submatrix
267c ---------------------------------------------------------
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
309c -----------------------------------------------------------
310c ... find strongly connected components in square submatrix,
311c putting this block into block lower triangular form.
312c -----------------------------------------------------------
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
351c ---------------------------------------------------
352c ... find connected components in vertical submatrix
353c ---------------------------------------------------
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
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
subroutine corsum(tmmtch, timhrz, timvrt, nhrows, nhcols, nsrows, nscols, nvrows, nvcols, output)
Definition corsum.f:3
subroutine finsum(timhrz, timesq, timvrt, hrzcmp, sqcmpn, vrtcmp, ccmstr, rcmstr, output)
Definition finsum.f:3
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
subroutine izero(n, x, incx)
Definition izero.f:2
subroutine maxmatch(nrows, ncols, colstr, rowind, prevcl, prevrw, marker, tryrow, nxtchp, rowset, colset)
Definition maxmatch.f:4
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
subroutine prtivs(title, n, x, output)
Definition prtivs.f:2
subroutine rectblk(nrows, ncols, marked, unmrkd, colstr, rowidx, colset, rowset, prevcl, tryrow, colmrk, rowmrk, nhrows, nhcols)
Definition rectblk.f:4
subroutine strchk(nrows, ncols, colstr, rowidx, nhrows, nhcols, nsrows, rnto, cnto, colset, rowset, output)
Definition strchk.f:4