SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pdormtr.f
Go to the documentation of this file.
1 SUBROUTINE pdormtr( SIDE, UPLO, TRANS, M, N, A, IA, JA, DESCA,
2 $ TAU, C, IC, JC, DESCC, WORK, LWORK, INFO )
3*
4* -- ScaLAPACK 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 SIDE, TRANS, UPLO
11 INTEGER IA, IC, INFO, JA, JC, LWORK, M, N
12* ..
13* .. Array Arguments ..
14 INTEGER DESCA( * ), DESCC( * )
15 DOUBLE PRECISION A( * ), C( * ), TAU( * ), WORK( * )
16* ..
17*
18* Purpose
19* =======
20*
21* PDORMTR overwrites the general real M-by-N distributed matrix
22* sub( C ) = C(IC:IC+M-1,JC:JC+N-1) with
23*
24* SIDE = 'L' SIDE = 'R'
25* TRANS = 'N': Q * sub( C ) sub( C ) * Q
26* TRANS = 'T': Q**T * sub( C ) sub( C ) * Q**T
27*
28* where Q is a real orthogonal distributed matrix of order nq, with
29* nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the
30* product of nq-1 elementary reflectors, as returned by PDSYTRD:
31*
32* if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);
33*
34* if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).
35*
36* Notes
37* =====
38*
39* Each global data object is described by an associated description
40* vector. This vector stores the information required to establish
41* the mapping between an object element and its corresponding process
42* and memory location.
43*
44* Let A be a generic term for any 2D block cyclicly distributed array.
45* Such a global array has an associated description vector DESCA.
46* In the following comments, the character _ should be read as
47* "of the global array".
48*
49* NOTATION STORED IN EXPLANATION
50* --------------- -------------- --------------------------------------
51* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
52* DTYPE_A = 1.
53* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
54* the BLACS process grid A is distribu-
55* ted over. The context itself is glo-
56* bal, but the handle (the integer
57* value) may vary.
58* M_A (global) DESCA( M_ ) The number of rows in the global
59* array A.
60* N_A (global) DESCA( N_ ) The number of columns in the global
61* array A.
62* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
63* the rows of the array.
64* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
65* the columns of the array.
66* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
67* row of the array A is distributed.
68* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
69* first column of the array A is
70* distributed.
71* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
72* array. LLD_A >= MAX(1,LOCr(M_A)).
73*
74* Let K be the number of rows or columns of a distributed matrix,
75* and assume that its process grid has dimension p x q.
76* LOCr( K ) denotes the number of elements of K that a process
77* would receive if K were distributed over the p processes of its
78* process column.
79* Similarly, LOCc( K ) denotes the number of elements of K that a
80* process would receive if K were distributed over the q processes of
81* its process row.
82* The values of LOCr() and LOCc() may be determined via a call to the
83* ScaLAPACK tool function, NUMROC:
84* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
85* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
86* An upper bound for these quantities may be computed by:
87* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
88* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
89*
90* Arguments
91* =========
92*
93* SIDE (global input) CHARACTER
94* = 'L': apply Q or Q**T from the Left;
95* = 'R': apply Q or Q**T from the Right.
96*
97* UPLO (global input) CHARACTER
98* = 'U': Upper triangle of A(IA:*,JA:*) contains elementary
99* reflectors from PDSYTRD;
100* = 'L': Lower triangle of A(IA:*,JA:*) contains elementary
101* reflectors from PDSYTRD.
102*
103* TRANS (global input) CHARACTER
104* = 'N': No transpose, apply Q;
105* = 'T': Transpose, apply Q**T.
106*
107* M (global input) INTEGER
108* The number of rows to be operated on i.e the number of rows
109* of the distributed submatrix sub( C ). M >= 0.
110*
111* N (global input) INTEGER
112* The number of columns to be operated on i.e the number of
113* columns of the distributed submatrix sub( C ). N >= 0.
114*
115* A (local input) DOUBLE PRECISION pointer into the local memory
116* to an array of dimension (LLD_A,LOCc(JA+M-1)) if SIDE='L',
117* or (LLD_A,LOCc(JA+N-1)) if SIDE = 'R'. The vectors which
118* define the elementary reflectors, as returned by PDSYTRD.
119* If SIDE = 'L', LLD_A >= max(1,LOCr(IA+M-1));
120* if SIDE = 'R', LLD_A >= max(1,LOCr(IA+N-1)).
121*
122* IA (global input) INTEGER
123* The row index in the global array A indicating the first
124* row of sub( A ).
125*
126* JA (global input) INTEGER
127* The column index in the global array A indicating the
128* first column of sub( A ).
129*
130* DESCA (global and local input) INTEGER array of dimension DLEN_.
131* The array descriptor for the distributed matrix A.
132*
133* TAU (local input) DOUBLE PRECISION array, dimension LTAU, where
134* if SIDE = 'L' and UPLO = 'U', LTAU = LOCc(M_A),
135* if SIDE = 'L' and UPLO = 'L', LTAU = LOCc(JA+M-2),
136* if SIDE = 'R' and UPLO = 'U', LTAU = LOCc(N_A),
137* if SIDE = 'R' and UPLO = 'L', LTAU = LOCc(JA+N-2).
138* TAU(i) must contain the scalar factor of the elementary
139* reflector H(i), as returned by PDSYTRD. TAU is tied to the
140* distributed matrix A.
141*
142* C (local input/local output) DOUBLE PRECISION pointer into the
143* local memory to an array of dimension (LLD_C,LOCc(JC+N-1)).
144* On entry, the local pieces of the distributed matrix sub(C).
145* On exit, sub( C ) is overwritten by Q*sub( C ) or Q'*sub( C )
146* or sub( C )*Q' or sub( C )*Q.
147*
148* IC (global input) INTEGER
149* The row index in the global array C indicating the first
150* row of sub( C ).
151*
152* JC (global input) INTEGER
153* The column index in the global array C indicating the
154* first column of sub( C ).
155*
156* DESCC (global and local input) INTEGER array of dimension DLEN_.
157* The array descriptor for the distributed matrix C.
158*
159* WORK (local workspace/local output) DOUBLE PRECISION array,
160* dimension (LWORK)
161* On exit, WORK(1) returns the minimal and optimal LWORK.
162*
163* LWORK (local or global input) INTEGER
164* The dimension of the array WORK.
165* LWORK is local input and must be at least
166*
167* If UPLO = 'U',
168* IAA = IA, JAA = JA+1, ICC = IC, JCC = JC;
169* else UPLO = 'L',
170* IAA = IA+1, JAA = JA;
171* if SIDE = 'L',
172* ICC = IC+1; JCC = JC;
173* else
174* ICC = IC; JCC = JC+1;
175* end if
176* end if
177*
178* If SIDE = 'L',
179* MI = M-1; NI = N;
180* LWORK >= MAX( (NB_A*(NB_A-1))/2, (NqC0 + MpC0)*NB_A ) +
181* NB_A * NB_A
182* else if SIDE = 'R',
183* MI = M; MI = N-1;
184* LWORK >= MAX( (NB_A*(NB_A-1))/2, ( NqC0 + MAX( NpA0 +
185* NUMROC( NUMROC( NI+ICOFFC, NB_A, 0, 0, NPCOL ),
186* NB_A, 0, 0, LCMQ ), MpC0 ) )*NB_A ) +
187* NB_A * NB_A
188* end if
189*
190* where LCMQ = LCM / NPCOL with LCM = ICLM( NPROW, NPCOL ),
191*
192* IROFFA = MOD( IAA-1, MB_A ), ICOFFA = MOD( JAA-1, NB_A ),
193* IAROW = INDXG2P( IAA, MB_A, MYROW, RSRC_A, NPROW ),
194* NpA0 = NUMROC( NI+IROFFA, MB_A, MYROW, IAROW, NPROW ),
195*
196* IROFFC = MOD( ICC-1, MB_C ), ICOFFC = MOD( JCC-1, NB_C ),
197* ICROW = INDXG2P( ICC, MB_C, MYROW, RSRC_C, NPROW ),
198* ICCOL = INDXG2P( JCC, NB_C, MYCOL, CSRC_C, NPCOL ),
199* MpC0 = NUMROC( MI+IROFFC, MB_C, MYROW, ICROW, NPROW ),
200* NqC0 = NUMROC( NI+ICOFFC, NB_C, MYCOL, ICCOL, NPCOL ),
201*
202* ILCM, INDXG2P and NUMROC are ScaLAPACK tool functions;
203* MYROW, MYCOL, NPROW and NPCOL can be determined by calling
204* the subroutine BLACS_GRIDINFO.
205*
206* If LWORK = -1, then LWORK is global input and a workspace
207* query is assumed; the routine only calculates the minimum
208* and optimal size for all work arrays. Each of these
209* values is returned in the first entry of the corresponding
210* work array, and no error message is issued by PXERBLA.
211*
212*
213* INFO (global output) INTEGER
214* = 0: successful exit
215* < 0: If the i-th argument is an array and the j-entry had
216* an illegal value, then INFO = -(i*100+j), if the i-th
217* argument is a scalar and had an illegal value, then
218* INFO = -i.
219*
220* Alignment requirements
221* ======================
222*
223* The distributed submatrices A(IA:*, JA:*) and C(IC:IC+M-1,JC:JC+N-1)
224* must verify some alignment properties, namely the following
225* expressions should be true:
226*
227* If SIDE = 'L',
228* ( MB_A.EQ.MB_C .AND. IROFFA.EQ.IROFFC .AND. IAROW.EQ.ICROW )
229* If SIDE = 'R',
230* ( MB_A.EQ.NB_C .AND. IROFFA.EQ.ICOFFC )
231*
232* =====================================================================
233*
234* .. Parameters ..
235 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
236 $ lld_, mb_, m_, nb_, n_, rsrc_
237 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
238 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
239 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
240* ..
241* .. Local Scalars ..
242 LOGICAL LEFT, LQUERY, NOTRAN, UPPER
243 INTEGER IAA, IAROW, ICC, ICCOL, ICOFFC, ICROW, ICTXT,
244 $ iinfo, iroffa, iroffc, jaa, jcc, lcm, lcmq,
245 $ lwmin, mi, mpc0, mycol, myrow, ni, npa0, npcol,
246 $ nprow, nq, nqc0
247* ..
248* .. Local Arrays ..
249 INTEGER IDUM1( 4 ), IDUM2( 4 )
250* ..
251* .. External Subroutines ..
252 EXTERNAL blacs_gridinfo, chk1mat, pchk2mat, pdormql,
254* ..
255* .. External Functions ..
256 LOGICAL LSAME
257 INTEGER ILCM, INDXG2P, NUMROC
258 EXTERNAL ilcm, indxg2p, lsame, numroc
259* ..
260* .. Intrinsic Functions ..
261 INTRINSIC dble, ichar, max, mod
262* ..
263* .. Executable Statements ..
264*
265* Get grid parameters
266*
267 ictxt = desca( ctxt_ )
268 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
269*
270* Test the input parameters
271*
272 info = 0
273 IF( nprow.EQ.-1 ) THEN
274 info = -(900+ctxt_)
275 ELSE
276 left = lsame( side, 'L' )
277 notran = lsame( trans, 'N' )
278 upper = lsame( uplo, 'U' )
279*
280 IF( upper ) THEN
281 iaa = ia
282 jaa = ja+1
283 icc = ic
284 jcc = jc
285 ELSE
286 iaa = ia+1
287 jaa = ja
288 IF( left ) THEN
289 icc = ic + 1
290 jcc = jc
291 ELSE
292 icc = ic
293 jcc = jc + 1
294 END IF
295 END IF
296*
297* NQ is the order of Q
298*
299 IF( left ) THEN
300 nq = m
301 mi = m - 1
302 ni = n
303 CALL chk1mat( mi, 4, nq-1, 4, iaa, jaa, desca, 9, info )
304 ELSE
305 nq = n
306 mi = m
307 ni = n - 1
308 CALL chk1mat( ni, 5, nq-1, 5, iaa, jaa, desca, 9, info )
309 END IF
310 CALL chk1mat( mi, 4, ni, 5, icc, jcc, descc, 14, info )
311 IF( info.EQ.0 ) THEN
312 iroffa = mod( iaa-1, desca( mb_ ) )
313 iroffc = mod( icc-1, descc( mb_ ) )
314 icoffc = mod( jcc-1, descc( nb_ ) )
315 iarow = indxg2p( iaa, desca( mb_ ), myrow, desca( rsrc_ ),
316 $ nprow )
317 icrow = indxg2p( icc, descc( mb_ ), myrow, descc( rsrc_ ),
318 $ nprow )
319 iccol = indxg2p( jcc, descc( nb_ ), mycol, descc( csrc_ ),
320 $ npcol )
321 mpc0 = numroc( mi+iroffc, descc( mb_ ), myrow, icrow,
322 $ nprow )
323 nqc0 = numroc( ni+icoffc, descc( nb_ ), mycol, iccol,
324 $ npcol )
325*
326 IF( left ) THEN
327 lwmin = max( ( desca( nb_ ) * ( desca( nb_ ) - 1 ) ) / 2,
328 $ ( mpc0 + nqc0 ) * desca( nb_ ) ) +
329 $ desca( nb_ ) * desca( nb_ )
330 ELSE
331 npa0 = numroc( ni+iroffa, desca( mb_ ), myrow, iarow,
332 $ nprow )
333 lcm = ilcm( nprow, npcol )
334 lcmq = lcm / npcol
335 lwmin = max( ( desca( nb_ ) * ( desca( nb_ ) - 1 ) )
336 $ / 2, ( nqc0 + max( npa0 + numroc( numroc(
337 $ ni+icoffc, desca( nb_ ), 0, 0, npcol ),
338 $ desca( nb_ ), 0, 0, lcmq ), mpc0 ) ) *
339 $ desca( nb_ ) ) + desca( nb_ ) * desca( nb_ )
340 END IF
341*
342 work( 1 ) = dble( lwmin )
343 lquery = ( lwork.EQ.-1 )
344 IF( .NOT.left .AND. .NOT.lsame( side, 'R' ) ) THEN
345 info = -1
346 ELSE IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
347 info = -2
348 ELSE IF( .NOT.lsame( trans, 'N' ) .AND.
349 $ .NOT.lsame( trans, 'T' ) ) THEN
350 info = -3
351 ELSE IF( .NOT.left .AND. desca( mb_ ).NE.descc( nb_ ) ) THEN
352 info = -(900+nb_)
353 ELSE IF( left .AND. iroffa.NE.iroffc ) THEN
354 info = -12
355 ELSE IF( left .AND. iarow.NE.icrow ) THEN
356 info = -12
357 ELSE IF( .NOT.left .AND. iroffa.NE.icoffc ) THEN
358 info = -13
359 ELSE IF( left .AND. desca( mb_ ).NE.descc( mb_ ) ) THEN
360 info = -(1400+mb_)
361 ELSE IF( ictxt.NE.descc( ctxt_ ) ) THEN
362 info = -(1400+ctxt_)
363 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
364 info = -16
365 END IF
366 END IF
367*
368 IF( left ) THEN
369 idum1( 1 ) = ichar( 'L' )
370 ELSE
371 idum1( 1 ) = ichar( 'R' )
372 END IF
373 idum2( 1 ) = 1
374 IF( upper ) THEN
375 idum1( 2 ) = ichar( 'U' )
376 ELSE
377 idum1( 2 ) = ichar( 'L' )
378 END IF
379 idum2( 2 ) = 2
380 IF( notran ) THEN
381 idum1( 3 ) = ichar( 'N' )
382 ELSE
383 idum1( 3 ) = ichar( 'T' )
384 END IF
385 idum2( 3 ) = 3
386 IF( lwork.EQ.-1 ) THEN
387 idum1( 4 ) = -1
388 ELSE
389 idum1( 4 ) = 1
390 END IF
391 idum2( 4 ) = 16
392 IF( left ) THEN
393 CALL pchk2mat( mi, 4, nq-1, 4, iaa, jaa, desca, 9, mi, 4,
394 $ ni, 5, icc, jcc, descc, 14, 4, idum1, idum2,
395 $ info )
396 ELSE
397 CALL pchk2mat( ni, 5, nq-1, 5, iaa, jaa, desca, 9, mi, 4,
398 $ ni, 5, icc, jcc, descc, 14, 4, idum1, idum2,
399 $ info )
400 END IF
401 END IF
402*
403 IF( info.NE.0 ) THEN
404 CALL pxerbla( ictxt, 'PDORMTR', -info )
405 RETURN
406 ELSE IF( lquery ) THEN
407 RETURN
408 END IF
409*
410* Quick return if possible
411*
412 IF( m.EQ.0 .OR. n.EQ.0 .OR. nq.EQ.1 )
413 $ RETURN
414*
415 IF( upper ) THEN
416*
417* Q was determined by a call to PDSYTRD with UPLO = 'U'
418*
419 CALL pdormql( side, trans, mi, ni, nq-1, a, iaa, jaa, desca,
420 $ tau, c, icc, jcc, descc, work, lwork, iinfo )
421*
422 ELSE
423*
424* Q was determined by a call to PDSYTRD with UPLO = 'L'
425*
426 CALL pdormqr( side, trans, mi, ni, nq-1, a, iaa, jaa, desca,
427 $ tau, c, icc, jcc, descc, work, lwork, iinfo )
428*
429 END IF
430*
431 work( 1 ) = dble( lwmin )
432*
433 RETURN
434*
435* End of PDORMTR
436*
437 END
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
Definition chk1mat.f:3
#define max(A, B)
Definition pcgemr.c:180
subroutine pchk2mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, mb, mbpos0, nb, nbpos0, ib, jb, descb, descbpos0, nextra, ex, expos, info)
Definition pchkxmat.f:175
subroutine pdormql(side, trans, m, n, k, a, ia, ja, desca, tau, c, ic, jc, descc, work, lwork, info)
Definition pdormql.f:3
subroutine pdormqr(side, trans, m, n, k, a, ia, ja, desca, tau, c, ic, jc, descc, work, lwork, info)
Definition pdormqr.f:3
subroutine pdormtr(side, uplo, trans, m, n, a, ia, ja, desca, tau, c, ic, jc, descc, work, lwork, info)
Definition pdormtr.f:3
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2