SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pzlapv2.f
Go to the documentation of this file.
1 SUBROUTINE pzlapv2( 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 COMPLEX*16 A( * )
16* ..
17*
18* Purpose
19* =======
20*
21* PZLAPV2 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) COMPLEX*16 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 $ pzswap
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 pzswap( 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 pzswap( 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 pzswap( 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 pzswap( 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 PZLAPV2
412*
413 END
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
min
#define min(A, B)
Definition pcgemr.c:181
pzlapv2
subroutine pzlapv2(direc, rowcol, m, n, a, ia, ja, desca, ipiv, ip, jp, descip)
Definition pzlapv2.f:3