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