ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdlapv2.f
Go to the documentation of this file.
1  SUBROUTINE pdlapv2( DIREC, ROWCOL, M, N, A, IA, JA, DESCA, IPIV,
2  $ IP, JP, DESCIP )
3 *
4 * -- ScaLAPACK auxiliary routine (version 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * May 1, 1997
8 *
9 * .. Scalar Arguments ..
10  CHARACTER DIREC, ROWCOL
11  INTEGER IA, IP, JA, JP, M, N
12 * ..
13 * .. Array Arguments ..
14  INTEGER DESCA( * ), DESCIP( * ), IPIV( * )
15  DOUBLE PRECISION A( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PDLAPV2 applies either P (permutation matrix indicated by IPIV)
22 * or inv( P ) to a M-by-N distributed matrix sub( A ) denoting
23 * A(IA:IA+M-1,JA:JA+N-1), resulting in row or column pivoting. The
24 * pivot vector should be aligned with the distributed matrix A. For
25 * pivoting the rows of sub( A ), IPIV should be distributed along a
26 * process column and replicated over all process rows. Similarly,
27 * IPIV should be distributed along a process row and replicated over
28 * all process columns for column pivoting.
29 *
30 * Notes
31 * =====
32 *
33 * Each global data object is described by an associated description
34 * vector. This vector stores the information required to establish
35 * the mapping between an object element and its corresponding process
36 * and memory location.
37 *
38 * Let A be a generic term for any 2D block cyclicly distributed array.
39 * Such a global array has an associated description vector DESCA.
40 * In the following comments, the character _ should be read as
41 * "of the global array".
42 *
43 * NOTATION STORED IN EXPLANATION
44 * --------------- -------------- --------------------------------------
45 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
46 * DTYPE_A = 1.
47 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
48 * the BLACS process grid A is distribu-
49 * ted over. The context itself is glo-
50 * bal, but the handle (the integer
51 * value) may vary.
52 * M_A (global) DESCA( M_ ) The number of rows in the global
53 * array A.
54 * N_A (global) DESCA( N_ ) The number of columns in the global
55 * array A.
56 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
57 * the rows of the array.
58 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
59 * the columns of the array.
60 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
61 * row of the array A is distributed.
62 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
63 * first column of the array A is
64 * distributed.
65 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
66 * array. LLD_A >= MAX(1,LOCr(M_A)).
67 *
68 * Let K be the number of rows or columns of a distributed matrix,
69 * and assume that its process grid has dimension p x q.
70 * LOCr( K ) denotes the number of elements of K that a process
71 * would receive if K were distributed over the p processes of its
72 * process column.
73 * Similarly, LOCc( K ) denotes the number of elements of K that a
74 * process would receive if K were distributed over the q processes of
75 * its process row.
76 * The values of LOCr() and LOCc() may be determined via a call to the
77 * ScaLAPACK tool function, NUMROC:
78 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
79 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
80 * An upper bound for these quantities may be computed by:
81 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
82 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
83 *
84 * Arguments
85 * =========
86 *
87 * DIREC (global input) CHARACTER
88 * Specifies in which order the permutation is applied:
89 * = 'F' (Forward) Applies pivots Forward from top of matrix.
90 * Computes P * sub( A );
91 * = 'B' (Backward) Applies pivots Backward from bottom of
92 * matrix. Computes inv( P ) * sub( A ).
93 *
94 * ROWCOL (global input) CHARACTER
95 * Specifies if the rows or columns are to be permuted:
96 * = 'R' Rows will be permuted,
97 * = 'C' Columns will be permuted.
98 *
99 * M (global input) INTEGER
100 * The number of rows to be operated on, i.e. the number of rows
101 * of the distributed submatrix sub( A ). M >= 0.
102 *
103 * N (global input) INTEGER
104 * The number of columns to be operated on, i.e. the number of
105 * columns of the distributed submatrix sub( A ). N >= 0.
106 *
107 * A (local input/local output) DOUBLE PRECISION pointer into the
108 * local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
109 * On entry, this local array contains the local pieces of the
110 * distributed matrix sub( A ) to which the row or columns
111 * interchanges will be applied. On exit, this array contains
112 * the local pieces of the permuted distributed matrix.
113 *
114 * IA (global input) INTEGER
115 * The row index in the global array A indicating the first
116 * row of sub( A ).
117 *
118 * JA (global input) INTEGER
119 * The column index in the global array A indicating the
120 * first column of sub( A ).
121 *
122 * DESCA (global and local input) INTEGER array of dimension DLEN_.
123 * The array descriptor for the distributed matrix A.
124 *
125 * IPIV (input) INTEGER array, dimension >= LOCr(M_A)+MB_A if
126 * ROWCOL = 'R', LOCc(N_A)+NB_A otherwise. It contains
127 * the pivoting information. IPIV(i) is the global row (column),
128 * local row (column) i was swapped with. The last piece of the
129 * array of size MB_A (resp. NB_A) is used as workspace. IPIV is
130 * tied to the distributed matrix A.
131 *
132 * IP (global input) INTEGER
133 * IPIV's global row index, which points to the beginning of the
134 * submatrix which is to be operated on.
135 *
136 * JP (global input) INTEGER
137 * IPIV's global column index, which points to the beginning of
138 * the submatrix which is to be operated on.
139 *
140 * DESCIP (global and local input) INTEGER array of dimension 8
141 * The array descriptor for the distributed matrix IPIV.
142 *
143 * =====================================================================
144 *
145 * .. Parameters ..
146  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
147  $ lld_, mb_, m_, nb_, n_, rsrc_
148  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
149  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
150  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
151 * ..
152 * .. Local Scalars ..
153  LOGICAL FORWRD, ROWPVT
154  INTEGER I, IB, ICTXT, ICURCOL, ICURROW, IIP, IP1, ITMP,
155  $ ipvwrk, j, jb, jjp, jp1, k, ma, mba, mycol,
156  $ myrow, nba, npcol, nprow
157 * ..
158 * .. External Subroutines ..
159  EXTERNAL blacs_gridinfo, igebs2d, igebr2d, infog2l,
160  $ pdswap
161 * ..
162 * .. External Functions ..
163  LOGICAL LSAME
164  INTEGER ICEIL, NUMROC
165  EXTERNAL iceil, lsame, numroc
166 * ..
167 * .. Intrinsic Functions ..
168  INTRINSIC min, mod
169 * ..
170 * .. Executable Statements ..
171 *
172  rowpvt = lsame( rowcol, 'R' )
173  IF( rowpvt ) THEN
174  IF( m.LE.1 .OR. n.LT.1 )
175  $ RETURN
176  ELSE
177  IF( m.LT.1 .OR. n.LE.1 )
178  $ RETURN
179  END IF
180  forwrd = lsame( direc, 'F' )
181 *
182 *
183 * Get grid and matrix parameters
184 *
185  ma = desca( m_ )
186  mba = desca( mb_ )
187  nba = desca( nb_ )
188  ictxt = desca( ctxt_ )
189  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
190 *
191 * If I'm applying pivots from beginning to end (e.g., repeating
192 * pivoting done earlier). Thus this section computes P * sub( A ).
193 *
194  IF( forwrd ) THEN
195  CALL infog2l( ip, jp, descip, nprow, npcol, myrow, mycol,
196  $ iip, jjp, icurrow, icurcol )
197 *
198 * If I'm pivoting the rows of sub( A )
199 *
200  IF( rowpvt ) THEN
201  ipvwrk = numroc( descip( m_ ), descip( mb_ ), myrow,
202  $ descip( rsrc_ ), nprow ) + 1 -
203  $ descip( mb_ )
204 *
205 * Loop over rows of sub( A )
206 *
207  i = ia
208  ib = min( m, iceil( ia, mba ) * mba - ia + 1 )
209  10 CONTINUE
210 *
211 * Find local pointer into IPIV, and broadcast this block's
212 * pivot information to everyone in process column
213 *
214  IF( myrow.EQ.icurrow ) THEN
215  CALL igebs2d( ictxt, 'Columnwise', ' ', ib, 1,
216  $ ipiv( iip ), ib )
217  itmp = iip
218  iip = iip + ib
219  ELSE
220  itmp = ipvwrk
221  CALL igebr2d( ictxt, 'Columnwise', ' ', ib, 1,
222  $ ipiv( itmp ), ib, icurrow, mycol )
223  END IF
224 *
225 * Pivot the block of rows
226 *
227  DO 20 k = i, i+ib-1
228  ip1 = ipiv( itmp ) - ip + ia
229  IF( ip1.NE.k )
230  $ CALL pdswap( n, a, k, ja, desca, ma, a, ip1, ja,
231  $ desca, ma )
232  itmp = itmp + 1
233  20 CONTINUE
234 *
235 * Go on to next row of processes, increment row counter,
236 * and figure number of rows to pivot next
237 *
238  icurrow = mod( icurrow+1, nprow )
239  i = i + ib
240  ib = min( mba, m-i+ia )
241  IF( ib .GT. 0 ) GOTO 10
242 *
243 * If I am pivoting the columns of sub( A )
244 *
245  ELSE
246  ipvwrk = numroc( descip( n_ ), descip( nb_ ), mycol,
247  $ descip( csrc_ ), npcol ) + 1 -
248  $ descip( nb_ )
249 *
250 * Loop over columns of sub( A )
251 *
252  j = ja
253  jb = min( n, iceil( ja, nba ) * nba - ja + 1 )
254  30 CONTINUE
255 *
256 * Find local pointer into IPIV, and broadcast this block's
257 * pivot information to everyone in process row
258 *
259  IF( mycol.EQ.icurcol ) THEN
260  CALL igebs2d( ictxt, 'Rowwise', ' ', jb, 1,
261  $ ipiv( jjp ), jb )
262  itmp = jjp
263  jjp = jjp + jb
264  ELSE
265  itmp = ipvwrk
266  CALL igebr2d( ictxt, 'Rowwise', ' ', jb, 1,
267  $ ipiv( itmp ), jb, myrow, icurcol )
268  END IF
269 *
270 * Pivot the block of columns
271 *
272  DO 40 k = j, j+jb-1
273  jp1 = ipiv( itmp ) - jp + ja
274  IF( jp1.NE.k )
275  $ CALL pdswap( m, a, ia, k, desca, 1, a, ia, jp1,
276  $ desca, 1 )
277  itmp = itmp + 1
278  40 CONTINUE
279 *
280 * Go on to next column of processes, increment column
281 * counter, and figure number of columns to pivot next
282 *
283  icurcol = mod( icurcol+1, npcol )
284  j = j + jb
285  jb = min( nba, n-j+ja )
286  IF( jb .GT. 0 ) GOTO 30
287  END IF
288 *
289 * If I want to apply pivots in reverse order, i.e. reversing
290 * pivoting done earlier. Thus this section computes
291 * inv( P ) * sub( A ).
292 *
293  ELSE
294 *
295 * If I'm pivoting the rows of sub( A )
296 *
297  IF( rowpvt ) THEN
298  CALL infog2l( ip+m-1, jp, descip, nprow, npcol, myrow,
299  $ mycol, iip, jjp, icurrow, icurcol )
300 *
301  ipvwrk = numroc( descip( m_ ), descip( mb_ ), myrow,
302  $ descip( rsrc_ ), nprow ) + 1 -
303  $ descip( mb_ )
304 *
305 * If I'm not in the current process row, my IIP points out
306 * past end of pivot vector (since I don't own a piece of the
307 * last row). Adjust IIP so it points at last pivot entry.
308 *
309  IF( myrow.NE.icurrow ) iip = iip - 1
310 *
311 * Loop over rows in reverse order, starting at last row
312 *
313  i = ia + m - 1
314  ib = mod( i, mba )
315  IF( ib .EQ. 0 ) ib = mba
316  ib = min( ib, m )
317  50 CONTINUE
318 *
319 * Find local pointer into IPIV, and broadcast this block's
320 * pivot information to everyone in process column
321 *
322  IF( myrow.EQ.icurrow ) THEN
323  itmp = iip
324  iip = iip - ib
325  CALL igebs2d( ictxt, 'Columnwise', ' ', ib, 1,
326  $ ipiv( iip+1 ), ib )
327  ELSE
328  CALL igebr2d( ictxt, 'Columnwise', ' ', ib, 1,
329  $ ipiv( ipvwrk ), ib, icurrow, mycol )
330  itmp = ipvwrk + ib - 1
331  END IF
332 *
333 * Pivot the block of rows
334 *
335  DO 60 k = i, i-ib+1, -1
336  ip1 = ipiv( itmp ) - ip + ia
337  IF( ip1.NE.k )
338  $ CALL pdswap( n, a, k, ja, desca, ma, a, ip1, ja,
339  $ desca, ma )
340  itmp = itmp - 1
341  60 CONTINUE
342 *
343 * Go to previous row of processes, decrement row counter,
344 * and figure number of rows to be pivoted next
345 *
346  icurrow = mod( nprow+icurrow-1, nprow )
347  i = i - ib
348  ib = min( mba, i-ia+1 )
349  IF( ib .GT. 0 ) GOTO 50
350 *
351 * Otherwise, I'm pivoting the columns of sub( A )
352 *
353  ELSE
354  CALL infog2l( ip, jp+n-1, descip, nprow, npcol, myrow,
355  $ mycol, iip, jjp, icurrow, icurcol )
356  ipvwrk = numroc( descip( n_ ), descip( nb_ ), mycol,
357  $ descip( csrc_ ), npcol ) + 1 -
358  $ descip( nb_ )
359 *
360 * If I'm not in the current process column, my JJP points out
361 * past end of pivot vector (since I don't own a piece of the
362 * last column). Adjust JJP so it points at last pivot entry.
363 *
364  IF( mycol.NE.icurcol ) jjp = jjp - 1
365 *
366 * Loop over columns in reverse order starting at last column
367 *
368  j = ja + n - 1
369  jb = mod( j, nba )
370  IF( jb .EQ. 0 ) jb = nba
371  jb = min( jb, n )
372  70 CONTINUE
373 *
374 * Find local pointer into IPIV, and broadcast this block's
375 * pivot information to everyone in process row
376 *
377  IF( mycol.EQ.icurcol ) THEN
378  itmp = jjp
379  jjp = jjp - jb
380  CALL igebs2d( ictxt, 'Rowwise', ' ', jb, 1,
381  $ ipiv( jjp+1 ), jb )
382  ELSE
383  CALL igebr2d( ictxt, 'Rowwise', ' ', jb, 1,
384  $ ipiv( ipvwrk ), jb, myrow, icurcol )
385  itmp = ipvwrk + jb - 1
386  END IF
387 *
388 * Pivot a block of columns
389 *
390  DO 80 k = j, j-jb+1, -1
391  jp1 = ipiv( itmp ) - jp + ja
392  IF( jp1.NE.k )
393  $ CALL pdswap( m, a, ia, k, desca, 1, a, ia, jp1,
394  $ desca, 1 )
395  itmp = itmp - 1
396  80 CONTINUE
397 *
398 * Go to previous row of processes, decrement row counter,
399 * and figure number of rows to be pivoted next
400 *
401  icurcol = mod( npcol+icurcol-1, npcol )
402  j = j - jb
403  jb = min( nba, j-ja+1 )
404  IF( jb .GT. 0 ) GOTO 70
405  END IF
406 *
407  END IF
408 *
409  RETURN
410 *
411 * End PDLAPV2
412 *
413  END
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pdlapv2
subroutine pdlapv2(DIREC, ROWCOL, M, N, A, IA, JA, DESCA, IPIV, IP, JP, DESCIP)
Definition: pdlapv2.f:3
min
#define min(A, B)
Definition: pcgemr.c:181