SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ pzunglq()

subroutine pzunglq ( integer  m,
integer  n,
integer  k,
complex*16, dimension( * )  a,
integer  ia,
integer  ja,
integer, dimension( * )  desca,
complex*16, dimension( * )  tau,
complex*16, dimension( * )  work,
integer  lwork,
integer  info 
)

Definition at line 1 of file pzunglq.f.

3*
4* -- ScaLAPACK routine (version 1.7) --
5* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6* and University of California, Berkeley.
7* May 25, 2001
8*
9* .. Scalar Arguments ..
10 INTEGER IA, INFO, JA, K, LWORK, M, N
11* ..
12* .. Array Arguments ..
13 INTEGER DESCA( * )
14 COMPLEX*16 A( * ), TAU( * ), WORK( * )
15* ..
16*
17* Purpose
18* =======
19*
20* PZUNGLQ generates an M-by-N complex distributed matrix Q denoting
21* A(IA:IA+M-1,JA:JA+N-1) with orthonormal rows, which is defined as
22* the first M rows of a product of K elementary reflectors of order N
23*
24* Q = H(k)' . . . H(2)' H(1)'
25*
26* as returned by PZGELQF.
27*
28* Notes
29* =====
30*
31* Each global data object is described by an associated description
32* vector. This vector stores the information required to establish
33* the mapping between an object element and its corresponding process
34* and memory location.
35*
36* Let A be a generic term for any 2D block cyclicly distributed array.
37* Such a global array has an associated description vector DESCA.
38* In the following comments, the character _ should be read as
39* "of the global array".
40*
41* NOTATION STORED IN EXPLANATION
42* --------------- -------------- --------------------------------------
43* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
44* DTYPE_A = 1.
45* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
46* the BLACS process grid A is distribu-
47* ted over. The context itself is glo-
48* bal, but the handle (the integer
49* value) may vary.
50* M_A (global) DESCA( M_ ) The number of rows in the global
51* array A.
52* N_A (global) DESCA( N_ ) The number of columns in the global
53* array A.
54* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
55* the rows of the array.
56* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
57* the columns of the array.
58* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
59* row of the array A is distributed.
60* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
61* first column of the array A is
62* distributed.
63* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
64* array. LLD_A >= MAX(1,LOCr(M_A)).
65*
66* Let K be the number of rows or columns of a distributed matrix,
67* and assume that its process grid has dimension p x q.
68* LOCr( K ) denotes the number of elements of K that a process
69* would receive if K were distributed over the p processes of its
70* process column.
71* Similarly, LOCc( K ) denotes the number of elements of K that a
72* process would receive if K were distributed over the q processes of
73* its process row.
74* The values of LOCr() and LOCc() may be determined via a call to the
75* ScaLAPACK tool function, NUMROC:
76* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
77* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
78* An upper bound for these quantities may be computed by:
79* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
80* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
81*
82* Arguments
83* =========
84*
85* M (global input) INTEGER
86* The number of rows to be operated on i.e the number of rows
87* of the distributed submatrix Q. M >= 0.
88*
89* N (global input) INTEGER
90* The number of columns to be operated on i.e the number of
91* columns of the distributed submatrix Q. N >= M >= 0.
92*
93* K (global input) INTEGER
94* The number of elementary reflectors whose product defines the
95* matrix Q. M >= K >= 0.
96*
97* A (local input/local output) COMPLEX*16 pointer into the
98* local memory to an array of dimension (LLD_A,LOCc(JA+N-1)).
99* On entry, the i-th row must contain the vector which defines
100* the elementary reflector H(i), IA <= i <= IA+K-1, as
101* returned by PZGELQF in the K rows of its distributed matrix
102* argument A(IA:IA+K-1,JA:*). On exit, this array contains the
103* local pieces of the M-by-N distributed matrix Q.
104*
105* IA (global input) INTEGER
106* The row index in the global array A indicating the first
107* row of sub( A ).
108*
109* JA (global input) INTEGER
110* The column index in the global array A indicating the
111* first column of sub( A ).
112*
113* DESCA (global and local input) INTEGER array of dimension DLEN_.
114* The array descriptor for the distributed matrix A.
115*
116* TAU (local input) COMPLEX*16, array, dimension LOCr(IA+K-1).
117* This array contains the scalar factors TAU(i) of the
118* elementary reflectors H(i) as returned by PZGELQF.
119* TAU is tied to the distributed matrix A.
120*
121* WORK (local workspace/local output) COMPLEX*16 array,
122* dimension (LWORK)
123* On exit, WORK(1) returns the minimal and optimal LWORK.
124*
125* LWORK (local or global input) INTEGER
126* The dimension of the array WORK.
127* LWORK is local input and must be at least
128* LWORK >= MB_A * ( MpA0 + NqA0 + MB_A ), where
129*
130* IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ),
131* IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ),
132* IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ),
133* MpA0 = NUMROC( M+IROFFA, MB_A, MYROW, IAROW, NPROW ),
134* NqA0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ),
135*
136* INDXG2P and NUMROC are ScaLAPACK tool functions;
137* MYROW, MYCOL, NPROW and NPCOL can be determined by calling
138* the subroutine BLACS_GRIDINFO.
139*
140* If LWORK = -1, then LWORK is global input and a workspace
141* query is assumed; the routine only calculates the minimum
142* and optimal size for all work arrays. Each of these
143* values is returned in the first entry of the corresponding
144* work array, and no error message is issued by PXERBLA.
145*
146*
147* INFO (global output) INTEGER
148* = 0: successful exit
149* < 0: If the i-th argument is an array and the j-entry had
150* an illegal value, then INFO = -(i*100+j), if the i-th
151* argument is a scalar and had an illegal value, then
152* INFO = -i.
153*
154* =====================================================================
155*
156* .. Parameters ..
157 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
158 $ LLD_, MB_, M_, NB_, N_, RSRC_
159 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
160 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
161 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
162 COMPLEX*16 ZERO
163 parameter( zero = ( 0.0d+0, 0.0d+0 ) )
164* ..
165* .. Local Scalars ..
166 LOGICAL LQUERY
167 CHARACTER COLBTOP, ROWBTOP
168 INTEGER I, IACOL, IAROW, IB, ICTXT, IINFO, IL, IN, IPW,
169 $ J, LWMIN, MPA0, MYCOL, MYROW, NPCOL, NPROW,
170 $ NQA0
171* ..
172* .. Local Arrays ..
173 INTEGER IDUM1( 2 ), IDUM2( 2 )
174* ..
175* .. External Subroutines ..
176 EXTERNAL blacs_gridinfo, chk1mat, pchk1mat, pb_topget,
177 $ pb_topset, pxerbla, pzlarfb, pzlarft,
178 $ pzlaset, pzungl2
179* ..
180* .. External Functions ..
181 INTEGER ICEIL, INDXG2P, NUMROC
182 EXTERNAL iceil, indxg2p, numroc
183* ..
184* .. Intrinsic Functions ..
185 INTRINSIC dble, dcmplx, max, min, mod
186* ..
187* .. Executable Statements ..
188*
189* Get grid parameters
190*
191 ictxt = desca( ctxt_ )
192 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
193*
194* Test the input parameters
195*
196 info = 0
197 IF( nprow.EQ.-1 ) THEN
198 info = -(700+ctxt_)
199 ELSE
200 CALL chk1mat( m, 1, n, 2, ia, ja, desca, 7, info )
201 IF( info.EQ.0 ) THEN
202 iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
203 $ nprow )
204 iacol = indxg2p( ja, desca( nb_ ), mycol, desca( csrc_ ),
205 $ npcol )
206 mpa0 = numroc( m+mod( ia-1, desca( mb_ ) ), desca( mb_ ),
207 $ myrow, iarow, nprow )
208 nqa0 = numroc( n+mod( ja-1, desca( nb_ ) ), desca( nb_ ),
209 $ mycol, iacol, npcol )
210 lwmin = desca( mb_ ) * ( mpa0 + nqa0 + desca( mb_ ) )
211*
212 work( 1 ) = dcmplx( dble( lwmin ) )
213 lquery = ( lwork.EQ.-1 )
214 IF( n.LT.m ) THEN
215 info = -2
216 ELSE IF( k.LT.0 .OR. k.GT.m ) THEN
217 info = -3
218 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
219 info = -10
220 END IF
221 END IF
222 idum1( 1 ) = k
223 idum2( 1 ) = 3
224 IF( lwork.EQ.-1 ) THEN
225 idum1( 2 ) = -1
226 ELSE
227 idum1( 2 ) = 1
228 END IF
229 idum2( 2 ) = 10
230 CALL pchk1mat( m, 1, n, 2, ia, ja, desca, 7, 2, idum1, idum2,
231 $ info )
232 END IF
233*
234 IF( info.NE.0 ) THEN
235 CALL pxerbla( ictxt, 'PZUNGLQ', -info )
236 RETURN
237 ELSE IF( lquery ) THEN
238 RETURN
239 END IF
240*
241* Quick return if possible
242*
243 IF( m.LE.0 )
244 $ RETURN
245*
246 ipw = desca( mb_ ) * desca( mb_ ) + 1
247 in = min( iceil( ia, desca( mb_ ) ) * desca( mb_ ), ia+k-1 )
248 il = max( ( (ia+k-2) / desca( mb_ ) ) * desca( mb_ ) + 1, ia )
249 CALL pb_topget( ictxt, 'Broadcast', 'Rowwise', rowbtop )
250 CALL pb_topget( ictxt, 'Broadcast', 'Columnwise', colbtop )
251 CALL pb_topset( ictxt, 'Broadcast', 'Rowwise', ' ' )
252 CALL pb_topset( ictxt, 'Broadcast', 'Columnwise', 'D-ring' )
253*
254 CALL pzlaset( 'All', ia+m-il, il-ia, zero, zero, a, il, ja,
255 $ desca )
256*
257* Use unblocked code for the last or only block.
258*
259 CALL pzungl2( ia+m-il, n-il+ia, ia+k-il, a, il, ja+il-ia, desca,
260 $ tau, work, lwork, iinfo )
261*
262* Is there at least one block of rows to loop over ?
263*
264 IF( il.GT.in+1 ) THEN
265*
266* Use blocked code
267*
268 DO 10 i = il-desca( mb_ ), in+1, -desca( mb_ )
269 ib = min( desca( mb_ ), ia+m-i )
270 j = ja + i - ia
271*
272 IF( i+ib.LE.ia+m-1 ) THEN
273*
274* Form the triangular factor of the block reflector
275* H = H(i) H(i+1) . . . H(i+ib-1)
276*
277 CALL pzlarft( 'Forward', 'Rowwise', n-i+ia, ib, a, i, j,
278 $ desca, tau, work, work( ipw ) )
279*
280* Apply H' to A(i+ib:ia+m-1,j:ja+n-1) from the right
281*
282 CALL pzlarfb( 'Right', 'Conjugate transpose', 'Forward',
283 $ 'Rowwise', m-i-ib+ia, n-i+ia, ib, a, i, j,
284 $ desca, work, a, i+ib, j, desca,
285 $ work( ipw ) )
286 END IF
287*
288* Apply H' to columns j:ja+n-1 of current block
289*
290 CALL pzungl2( ib, n-i+ia, ib, a, i, j, desca, tau, work,
291 $ lwork, iinfo )
292*
293* Set columns ia:i-1 of current block to zero
294*
295 CALL pzlaset( 'All', ib, i-ia, zero, zero, a, i, ja, desca )
296 10 CONTINUE
297*
298 END IF
299*
300* Handle first block separately
301*
302 IF( il.GT.ia ) THEN
303*
304 ib = in - ia + 1
305*
306* Form the triangular factor of the block reflector
307* H = H(i) H(i+1) . . . H(i+ib-1)
308*
309 CALL pzlarft( 'Forward', 'Rowwise', n, ib, a, ia, ja, desca,
310 $ tau, work, work( ipw ) )
311*
312* Apply H' to A(ia+ib:ia+m-1,ja:ja+n-1) from the right
313*
314 CALL pzlarfb( 'Right', 'Conjugate transpose', 'Forward',
315 $ 'Rowwise', m-ib, n, ib, a, ia, ja, desca, work,
316 $ a, ia+ib, ja, desca, work( ipw ) )
317*
318* Apply H' to columns ja:ja+n-1 of current block
319*
320 CALL pzungl2( ib, n, ib, a, ia, ja, desca, tau, work, lwork,
321 $ iinfo )
322*
323 END IF
324*
325 CALL pb_topset( ictxt, 'Broadcast', 'Rowwise', rowbtop )
326 CALL pb_topset( ictxt, 'Broadcast', 'Columnwise', colbtop )
327*
328 work( 1 ) = dcmplx( dble( lwmin ) )
329*
330 RETURN
331*
332* End of PZUNGLQ
333*
chk1mat
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
Definition chk1mat.f:3
iceil
integer function iceil(inum, idenom)
Definition iceil.f:2
indxg2p
integer function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
Definition indxg2p.f:2
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
pchk1mat
subroutine pchk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, nextra, ex, expos, info)
Definition pchkxmat.f:3
pxerbla
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2
pzlaset
subroutine pzlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
Definition pzblastst.f:7509
pzlarfb
subroutine pzlarfb(side, trans, direct, storev, m, n, k, v, iv, jv, descv, t, c, ic, jc, descc, work)
Definition pzlarfb.f:3
pzlarft
subroutine pzlarft(direct, storev, n, k, v, iv, jv, descv, tau, t, work)
Definition pzlarft.f:3
pzungl2
subroutine pzungl2(m, n, k, a, ia, ja, desca, tau, work, lwork, info)
Definition pzungl2.f:3
Here is the call graph for this function: