ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pclarzb.f
Go to the documentation of this file.
1  SUBROUTINE pclarzb( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V,
2  $ IV, JV, DESCV, T, C, IC, JC, DESCC, WORK )
3 *
4 * -- ScaLAPACK auxiliary routine (version 2.0.2) --
5 * Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
6 * May 1 2012
7 *
8 * .. Scalar Arguments ..
9  CHARACTER DIRECT, SIDE, STOREV, TRANS
10  INTEGER IC, IV, JC, JV, K, L, M, N
11 * ..
12 * .. Array Arguments ..
13  INTEGER DESCC( * ), DESCV( * )
14  COMPLEX C( * ), T( * ), V( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * PCLARZB applies a complex block reflector Q or its conjugate
21 * transpose Q**H to a complex M-by-N distributed matrix sub( C )
22 * denoting C(IC:IC+M-1,JC:JC+N-1), from the left or the right.
23 *
24 * Q is a product of k elementary reflectors as returned by PCTZRZF.
25 *
26 * Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
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 * SIDE (global input) CHARACTER
86 * = 'L': apply Q or Q**H from the Left;
87 * = 'R': apply Q or Q**H from the Right.
88 *
89 * TRANS (global input) CHARACTER
90 * = 'N': No transpose, apply Q;
91 * = 'C': Conjugate transpose, apply Q**H.
92 *
93 * DIRECT (global input) CHARACTER
94 * Indicates how H is formed from a product of elementary
95 * reflectors
96 * = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
97 * = 'B': H = H(k) . . . H(2) H(1) (Backward)
98 *
99 * STOREV (global input) CHARACTER
100 * Indicates how the vectors which define the elementary
101 * reflectors are stored:
102 * = 'C': Columnwise (not supported yet)
103 * = 'R': Rowwise
104 *
105 * M (global input) INTEGER
106 * The number of rows to be operated on i.e the number of rows
107 * of the distributed submatrix sub( C ). M >= 0.
108 *
109 * N (global input) INTEGER
110 * The number of columns to be operated on i.e the number of
111 * columns of the distributed submatrix sub( C ). N >= 0.
112 *
113 * K (global input) INTEGER
114 * The order of the matrix T (= the number of elementary
115 * reflectors whose product defines the block reflector).
116 *
117 * L (global input) INTEGER
118 * The columns of the distributed submatrix sub( A ) containing
119 * the meaningful part of the Householder reflectors.
120 * If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
121 *
122 * V (local input) COMPLEX pointer into the local memory
123 * to an array of dimension (LLD_V, LOCc(JV+M-1)) if SIDE = 'L',
124 * (LLD_V, LOCc(JV+N-1)) if SIDE = 'R'. It contains the local
125 * pieces of the distributed vectors V representing the
126 * Householder transformation as returned by PCTZRZF.
127 * LLD_V >= LOCr(IV+K-1).
128 *
129 * IV (global input) INTEGER
130 * The row index in the global array V indicating the first
131 * row of sub( V ).
132 *
133 * JV (global input) INTEGER
134 * The column index in the global array V indicating the
135 * first column of sub( V ).
136 *
137 * DESCV (global and local input) INTEGER array of dimension DLEN_.
138 * The array descriptor for the distributed matrix V.
139 *
140 * T (local input) COMPLEX array, dimension MB_V by MB_V
141 * The lower triangular matrix T in the representation of the
142 * block reflector.
143 *
144 * C (local input/local output) COMPLEX pointer into the
145 * local memory to an array of dimension (LLD_C,LOCc(JC+N-1)).
146 * On entry, the M-by-N distributed matrix sub( C ). On exit,
147 * sub( C ) is overwritten by Q*sub( C ) or Q'*sub( C ) or
148 * sub( C )*Q or sub( C )*Q'.
149 *
150 * IC (global input) INTEGER
151 * The row index in the global array C indicating the first
152 * row of sub( C ).
153 *
154 * JC (global input) INTEGER
155 * The column index in the global array C indicating the
156 * first column of sub( C ).
157 *
158 * DESCC (global and local input) INTEGER array of dimension DLEN_.
159 * The array descriptor for the distributed matrix C.
160 *
161 * WORK (local workspace) COMPLEX array, dimension (LWORK)
162 * If STOREV = 'C',
163 * if SIDE = 'L',
164 * LWORK >= ( NqC0 + MpC0 ) * K
165 * else if SIDE = 'R',
166 * LWORK >= ( NqC0 + MAX( NpV0 + NUMROC( NUMROC( N+ICOFFC,
167 * NB_V, 0, 0, NPCOL ), NB_V, 0, 0, LCMQ ),
168 * MpC0 ) ) * K
169 * end if
170 * else if STOREV = 'R',
171 * if SIDE = 'L',
172 * LWORK >= ( MpC0 + MAX( MqV0 + NUMROC( NUMROC( M+IROFFC,
173 * MB_V, 0, 0, NPROW ), MB_V, 0, 0, LCMP ),
174 * NqC0 ) ) * K
175 * else if SIDE = 'R',
176 * LWORK >= ( MpC0 + NqC0 ) * K
177 * end if
178 * end if
179 *
180 * where LCMQ = LCM / NPCOL with LCM = ICLM( NPROW, NPCOL ),
181 *
182 * IROFFV = MOD( IV-1, MB_V ), ICOFFV = MOD( JV-1, NB_V ),
183 * IVROW = INDXG2P( IV, MB_V, MYROW, RSRC_V, NPROW ),
184 * IVCOL = INDXG2P( JV, NB_V, MYCOL, CSRC_V, NPCOL ),
185 * MqV0 = NUMROC( M+ICOFFV, NB_V, MYCOL, IVCOL, NPCOL ),
186 * NpV0 = NUMROC( N+IROFFV, MB_V, MYROW, IVROW, NPROW ),
187 *
188 * IROFFC = MOD( IC-1, MB_C ), ICOFFC = MOD( JC-1, NB_C ),
189 * ICROW = INDXG2P( IC, MB_C, MYROW, RSRC_C, NPROW ),
190 * ICCOL = INDXG2P( JC, NB_C, MYCOL, CSRC_C, NPCOL ),
191 * MpC0 = NUMROC( M+IROFFC, MB_C, MYROW, ICROW, NPROW ),
192 * NpC0 = NUMROC( N+ICOFFC, MB_C, MYROW, ICROW, NPROW ),
193 * NqC0 = NUMROC( N+ICOFFC, NB_C, MYCOL, ICCOL, NPCOL ),
194 *
195 * ILCM, INDXG2P and NUMROC are ScaLAPACK tool functions;
196 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling
197 * the subroutine BLACS_GRIDINFO.
198 *
199 * Alignment requirements
200 * ======================
201 *
202 * The distributed submatrices V(IV:*, JV:*) and C(IC:IC+M-1,JC:JC+N-1)
203 * must verify some alignment properties, namely the following
204 * expressions should be true:
205 *
206 * If STOREV = 'Columnwise'
207 * If SIDE = 'Left',
208 * ( MB_V.EQ.MB_C .AND. IROFFV.EQ.IROFFC .AND. IVROW.EQ.ICROW )
209 * If SIDE = 'Right',
210 * ( MB_V.EQ.NB_C .AND. IROFFV.EQ.ICOFFC )
211 * else if STOREV = 'Rowwise'
212 * If SIDE = 'Left',
213 * ( NB_V.EQ.MB_C .AND. ICOFFV.EQ.IROFFC )
214 * If SIDE = 'Right',
215 * ( NB_V.EQ.NB_C .AND. ICOFFV.EQ.ICOFFC .AND. IVCOL.EQ.ICCOL )
216 * end if
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  COMPLEX ONE, ZERO
227  parameter( one = ( 1.0e+0, 0.0e+0 ),
228  $ zero = ( 0.0e+0, 0.0e+0 ) )
229 * ..
230 * .. Local Scalars ..
231  LOGICAL LEFT
232  CHARACTER COLBTOP, TRANST
233  INTEGER ICCOL1, ICCOL2, ICOFFC1, ICOFFC2, ICOFFV,
234  $ icrow1, icrow2, ictxt, iibeg, iic1, iic2,
235  $ iiend, iinxt, iiv, ileft, info, ioffc2, ioffv,
236  $ ipt, ipv, ipw, iroffc1, iroffc2, itop, ivcol,
237  $ ivrow, j, jjbeg, jjend, jjnxt, jjc1, jjc2, jjv,
238  $ ldc, ldv, lv, lw, mbc, mbv, mpc1, mpc2, mpc20,
239  $ mqv, mqv0, mycol, mydist, myrow, nbc, nbv,
240  $ npcol, nprow, nqc1, nqc2, nqcall, nqv
241 * ..
242 * .. External Subroutines ..
243  EXTERNAL blacs_abort, blacs_gridinfo, cgebr2d,
244  $ cgebs2d, cgemm, cgsum2d, clacgv,
245  $ clamov, claset, ctrbr2d, ctrbs2d,
246  $ ctrmm, infog2l, pbcmatadd, pbctran,
247  $ pb_topget, pxerbla
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC max, min, mod
251 * ..
252 * .. External Functions ..
253  LOGICAL LSAME
254  INTEGER ICEIL, NUMROC
255  EXTERNAL iceil, lsame, numroc
256 * ..
257 * .. Executable Statements ..
258 *
259 * Quick return if possible
260 *
261  IF( m.LE.0 .OR. n.LE.0 .OR. k.LE.0 )
262  $ RETURN
263 *
264 * Get grid parameters
265 *
266  ictxt = descc( ctxt_ )
267  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
268 *
269 * Check for currently supported options
270 *
271  info = 0
272  IF( .NOT.lsame( direct, 'B' ) ) THEN
273  info = -3
274  ELSE IF( .NOT.lsame( storev, 'R' ) ) THEN
275  info = -4
276  END IF
277  IF( info.NE.0 ) THEN
278  CALL pxerbla( ictxt, 'PCLARZB', -info )
279  CALL blacs_abort( ictxt, 1 )
280  RETURN
281  END IF
282 *
283  left = lsame( side, 'L' )
284  IF( lsame( trans, 'N' ) ) THEN
285  transt = 'C'
286  ELSE
287  transt = 'N'
288  END IF
289 *
290  CALL infog2l( iv, jv, descv, nprow, npcol, myrow, mycol, iiv, jjv,
291  $ ivrow, ivcol )
292  mbv = descv( mb_ )
293  nbv = descv( nb_ )
294  icoffv = mod( jv-1, nbv )
295  nqv = numroc( l+icoffv, nbv, mycol, ivcol, npcol )
296  IF( mycol.EQ.ivcol )
297  $ nqv = nqv - icoffv
298  ldv = descv( lld_ )
299  iiv = min( iiv, ldv )
300  jjv = min( jjv, max( 1, numroc( descv( n_ ), nbv, mycol,
301  $ descv( csrc_ ), npcol ) ) )
302  ioffv = iiv + ( jjv-1 ) * ldv
303  mbc = descc( mb_ )
304  nbc = descc( nb_ )
305  nqcall = numroc( descc( n_ ), nbc, mycol, descc( csrc_ ), npcol )
306  CALL infog2l( ic, jc, descc, nprow, npcol, myrow, mycol, iic1,
307  $ jjc1, icrow1, iccol1 )
308  ldc = descc( lld_ )
309  iic1 = min( iic1, ldc )
310  jjc1 = min( jjc1, max( 1, nqcall ) )
311 *
312  IF( left ) THEN
313  iroffc1 = mod( ic-1, mbc )
314  mpc1 = numroc( k+iroffc1, mbc, myrow, icrow1, nprow )
315  IF( myrow.EQ.icrow1 )
316  $ mpc1 = mpc1 - iroffc1
317  icoffc1 = mod( jc-1, nbc )
318  nqc1 = numroc( n+icoffc1, nbc, mycol, iccol1, npcol )
319  IF( mycol.EQ.iccol1 )
320  $ nqc1 = nqc1 - icoffc1
321  CALL infog2l( ic+m-l, jc, descc, nprow, npcol, myrow, mycol,
322  $ iic2, jjc2, icrow2, iccol2 )
323  iroffc2 = mod( ic+m-l-1, mbc )
324  mpc2 = numroc( l+iroffc2, mbc, myrow, icrow2, nprow )
325  IF( myrow.EQ.icrow2 )
326  $ mpc2 = mpc2 - iroffc2
327  icoffc2 = icoffc1
328  nqc2 = nqc1
329  ELSE
330  iroffc1 = mod( ic-1, mbc )
331  mpc1 = numroc( m+iroffc1, mbc, myrow, icrow1, nprow )
332  IF( myrow.EQ.icrow1 )
333  $ mpc1 = mpc1 - iroffc1
334  icoffc1 = mod( jc-1, nbc )
335  nqc1 = numroc( k+icoffc1, nbc, mycol, iccol1, npcol )
336  IF( mycol.EQ.iccol1 )
337  $ nqc1 = nqc1 - icoffc1
338  CALL infog2l( ic, jc+n-l, descc, nprow, npcol, myrow, mycol,
339  $ iic2, jjc2, icrow2, iccol2 )
340  iroffc2 = iroffc1
341  mpc2 = mpc1
342  icoffc2 = mod( jc+n-l-1, nbc )
343  nqc2 = numroc( l+icoffc2, nbc, mycol, iccol2, npcol )
344  IF( mycol.EQ.iccol2 )
345  $ nqc2 = nqc2 - icoffc2
346  END IF
347  iic2 = min( iic2, ldc )
348  jjc2 = min( jjc2, nqcall )
349  ioffc2 = iic2 + ( jjc2-1 ) * ldc
350 *
351  IF( lsame( side, 'L' ) ) THEN
352 *
353 * Form Q*sub( C ) or Q'*sub( C )
354 *
355 * IROFFC2 = ICOFFV is required by the current transposition
356 * routine PBCTRAN
357 *
358  mqv0 = numroc( m+icoffv, nbv, mycol, ivcol, npcol )
359  IF( mycol.EQ.ivcol ) THEN
360  mqv = mqv0 - icoffv
361  ELSE
362  mqv = mqv0
363  END IF
364  IF( myrow.EQ.icrow2 ) THEN
365  mpc20 = mpc2 + iroffc2
366  ELSE
367  mpc20 = mpc2
368  END IF
369 *
370 * Locally V( IOFFV ) is K x MQV, C( IOFFC2 ) is MPC2 x NQC2
371 * WORK( IPV ) is MPC20 x K = [ . V( IOFFV ) ]'
372 * WORK( IPW ) is K x MQV0 = [ . V( IOFFV ) ]
373 * WORK( IPT ) is the workspace for PBCTRAN
374 *
375  ipv = 1
376  ipw = ipv + mpc20 * k
377  ipt = ipw + k * mqv0
378  lv = max( 1, mpc20 )
379  lw = max( 1, k )
380 *
381  IF( myrow.EQ.ivrow ) THEN
382  IF( mycol.EQ.ivcol ) THEN
383  CALL clamov( 'All', k, mqv, v( ioffv ), ldv,
384  $ work( ipw+icoffv*lw ), lw )
385  ELSE
386  CALL clamov( 'All', k, mqv, v( ioffv ), ldv,
387  $ work( ipw ), lw )
388  END IF
389  END IF
390 *
391 * WORK( IPV ) = WORK( IPW )' (replicated) is MPC20 x K
392 *
393  CALL pbctran( ictxt, 'Rowwise', 'Conjugate transpose', k,
394  $ m+icoffv, descv( nb_ ), work( ipw ), lw, zero,
395  $ work( ipv ), lv, ivrow, ivcol, icrow2, -1,
396  $ work( ipt ) )
397 *
398 * WORK( IPV ) = ( . V )' -> WORK( IPV ) = V' is MPC2 x K
399 *
400  IF( myrow.EQ.icrow2 )
401  $ ipv = ipv + iroffc2
402 *
403 * WORK( IPW ) becomes NQC2 x K = C( IOFFC2 )' * V'
404 * WORK( IPW ) = C( IOFFC2 )' * V' (NQC2 x MPC2 x K) -> NQC2 x K
405 *
406  lw = max( 1, nqc2 )
407 *
408  IF( mpc2.GT.0 ) THEN
409  CALL cgemm( 'Transpose', 'No transpose', nqc2, k, mpc2,
410  $ one, c( ioffc2 ), ldc, work( ipv ), lv, zero,
411  $ work( ipw ), lw )
412  ELSE
413  CALL claset( 'All', nqc2, k, zero, zero, work( ipw ), lw )
414  END IF
415 *
416 * WORK( IPW ) = WORK( IPW ) + C1 ( NQC1 = NQC2 )
417 *
418  IF( mpc1.GT.0 ) THEN
419  mydist = mod( myrow-icrow1+nprow, nprow )
420  itop = max( 0, mydist * mbc - iroffc1 )
421  iibeg = iic1
422  iiend = iic1 + mpc1 - 1
423  iinxt = min( iceil( iibeg, mbc ) * mbc, iiend )
424 *
425  10 CONTINUE
426  IF( iibeg.LE.iinxt ) THEN
427  CALL pbcmatadd( ictxt, 'Transpose', nqc2, iinxt-iibeg+1,
428  $ one, c( iibeg+(jjc1-1)*ldc ), ldc, one,
429  $ work( ipw+itop ), lw )
430  mydist = mydist + nprow
431  itop = mydist * mbc - iroffc1
432  iibeg = iinxt +1
433  iinxt = min( iinxt+mbc, iiend )
434  GO TO 10
435  END IF
436  END IF
437 *
438  CALL cgsum2d( ictxt, 'Columnwise', ' ', nqc2, k, work( ipw ),
439  $ lw, ivrow, mycol )
440 *
441 * WORK( IPW ) = WORK( IPW ) * T' or WORK( IPW ) * T
442 *
443  IF( myrow.EQ.ivrow ) THEN
444  IF( mycol.EQ.ivcol ) THEN
445 *
446 * Broadcast the block reflector to the other columns.
447 *
448  CALL ctrbs2d( ictxt, 'Rowwise', ' ', 'Lower', 'Non unit',
449  $ k, k, t, mbv )
450  ELSE
451  CALL ctrbr2d( ictxt, 'Rowwise', ' ', 'Lower', 'Non unit',
452  $ k, k, t, mbv, myrow, ivcol )
453  END IF
454  CALL ctrmm( 'Right', 'Lower', transt, 'Non unit', nqc2, k,
455  $ one, t, mbv, work( ipw ), lw )
456 *
457  CALL cgebs2d( ictxt, 'Columnwise', ' ', nqc2, k,
458  $ work( ipw ), lw )
459  ELSE
460  CALL cgebr2d( ictxt, 'Columnwise', ' ', nqc2, k,
461  $ work( ipw ), lw, ivrow, mycol )
462  END IF
463 *
464 * C1 = C1 - WORK( IPW )
465 *
466  IF( mpc1.GT.0 ) THEN
467  mydist = mod( myrow-icrow1+nprow, nprow )
468  itop = max( 0, mydist * mbc - iroffc1 )
469  iibeg = iic1
470  iiend = iic1 + mpc1 - 1
471  iinxt = min( iceil( iibeg, mbc ) * mbc, iiend )
472 *
473  20 CONTINUE
474  IF( iibeg.LE.iinxt ) THEN
475  CALL pbcmatadd( ictxt, 'Transpose', iinxt-iibeg+1, nqc2,
476  $ -one, work( ipw+itop ), lw, one,
477  $ c( iibeg+(jjc1-1)*ldc ), ldc )
478  mydist = mydist + nprow
479  itop = mydist * mbc - iroffc1
480  iibeg = iinxt +1
481  iinxt = min( iinxt+mbc, iiend )
482  GO TO 20
483  END IF
484  END IF
485 *
486 * C2 C2 - V' * W'
487 * C( IOFFC2 ) = C( IOFFC2 ) - WORK( IPV ) * WORK( IPW )'
488 * MPC2 x NQC2 MPC2 x K K x NQC2
489 *
490  DO 30 j = 1, k
491  CALL clacgv( mpc2, work( ipv+(j-1)*lv ), 1 )
492  30 CONTINUE
493  CALL cgemm( 'No transpose', 'Transpose', mpc2, nqc2, k, -one,
494  $ work( ipv ), lv, work( ipw ), lw, one,
495  $ c( ioffc2 ), ldc )
496 *
497  ELSE
498 *
499 * Form sub( C ) * Q or sub( C ) * Q'
500 *
501 * Locally V( IOFFV ) is K x NQV, C( IOFFC2 ) is MPC2 x NQC2
502 * WORK( IPV ) is K x NQV = V( IOFFV ), NQV = NQC2
503 * WORK( IPW ) is MPC2 x K = C( IOFFC2 ) * V( IOFFV )'
504 *
505  ipv = 1
506  ipw = ipv + k * nqc2
507  lv = max( 1, k )
508  lw = max( 1, mpc2 )
509 *
510 * Broadcast V to the other process rows.
511 *
512  CALL pb_topget( ictxt, 'Broadcast', 'Columnwise', colbtop )
513  IF( myrow.EQ.ivrow ) THEN
514  CALL cgebs2d( ictxt, 'Columnwise', colbtop, k, nqc2,
515  $ v( ioffv ), ldv )
516  IF( mycol.EQ.ivcol )
517  $ CALL ctrbs2d( ictxt, 'Columnwise', colbtop, 'Lower',
518  $ 'Non unit', k, k, t, mbv )
519  CALL clamov( 'All', k, nqc2, v( ioffv ), ldv, work( ipv ),
520  $ lv )
521  ELSE
522  CALL cgebr2d( ictxt, 'Columnwise', colbtop, k, nqc2,
523  $ work( ipv ), lv, ivrow, mycol )
524  IF( mycol.EQ.ivcol )
525  $ CALL ctrbr2d( ictxt, 'Columnwise', colbtop, 'Lower',
526  $ 'Non unit', k, k, t, mbv, ivrow, mycol )
527  END IF
528 *
529 * WORK( IPV ) is K x NQC2 = V = V( IOFFV )
530 * WORK( IPW ) = C( IOFFC2 ) * V' (MPC2 x NQC2 x K) -> MPC2 x K
531 *
532  IF( nqc2.GT.0 ) THEN
533  CALL cgemm( 'No Transpose', 'Transpose', mpc2, k, nqc2,
534  $ one, c( ioffc2 ), ldc, work( ipv ), lv, zero,
535  $ work( ipw ), lw )
536  ELSE
537  CALL claset( 'All', mpc2, k, zero, zero, work( ipw ), lw )
538  END IF
539 *
540 * WORK( IPW ) = WORK( IPW ) + C1 ( MPC1 = MPC2 )
541 *
542  IF( nqc1.GT.0 ) THEN
543  mydist = mod( mycol-iccol1+npcol, npcol )
544  ileft = max( 0, mydist * nbc - icoffc1 )
545  jjbeg = jjc1
546  jjend = jjc1 + nqc1 - 1
547  jjnxt = min( iceil( jjbeg, nbc ) * nbc, jjend )
548 *
549  40 CONTINUE
550  IF( jjbeg.LE.jjnxt ) THEN
551  CALL pbcmatadd( ictxt, 'No transpose', mpc2,
552  $ jjnxt-jjbeg+1, one,
553  $ c( iic1+(jjbeg-1)*ldc ), ldc, one,
554  $ work( ipw+ileft*lw ), lw )
555  mydist = mydist + npcol
556  ileft = mydist * nbc - icoffc1
557  jjbeg = jjnxt +1
558  jjnxt = min( jjnxt+nbc, jjend )
559  GO TO 40
560  END IF
561  END IF
562 *
563  CALL cgsum2d( ictxt, 'Rowwise', ' ', mpc2, k, work( ipw ),
564  $ lw, myrow, ivcol )
565 *
566 * WORK( IPW ) = WORK( IPW ) * T' or WORK( IPW ) * T
567 *
568  IF( mycol.EQ.ivcol ) THEN
569  DO 50 j = 1, k
570  CALL clacgv( k-j+1, t( j+(j-1)*mbv ), 1 )
571  50 CONTINUE
572  CALL ctrmm( 'Right', 'Lower', trans, 'Non unit', mpc2, k,
573  $ one, t, mbv, work( ipw ), lw )
574  CALL cgebs2d( ictxt, 'Rowwise', ' ', mpc2, k, work( ipw ),
575  $ lw )
576  DO 60 j = 1, k
577  CALL clacgv( k-j+1, t( j+(j-1)*mbv ), 1 )
578  60 CONTINUE
579  ELSE
580  CALL cgebr2d( ictxt, 'Rowwise', ' ', mpc2, k, work( ipw ),
581  $ lw, myrow, ivcol )
582  END IF
583 *
584 * C1 = C1 - WORK( IPW )
585 *
586  IF( nqc1.GT.0 ) THEN
587  mydist = mod( mycol-iccol1+npcol, npcol )
588  ileft = max( 0, mydist * nbc - icoffc1 )
589  jjbeg = jjc1
590  jjend = jjc1 + nqc1 - 1
591  jjnxt = min( iceil( jjbeg, nbc ) * nbc, jjend )
592 *
593  70 CONTINUE
594  IF( jjbeg.LE.jjnxt ) THEN
595  CALL pbcmatadd( ictxt, 'No transpose', mpc2,
596  $ jjnxt-jjbeg+1, -one,
597  $ work( ipw+ileft*lw ), lw, one,
598  $ c( iic1+(jjbeg-1)*ldc ), ldc )
599  mydist = mydist + npcol
600  ileft = mydist * nbc - icoffc1
601  jjbeg = jjnxt +1
602  jjnxt = min( jjnxt+nbc, jjend )
603  GO TO 70
604  END IF
605  END IF
606 *
607 * C2 C2 - W * conjg( V )
608 * C( IOFFC ) = C( IOFFC ) - WORK( IPW ) * conjg( WORK( IPV ) )
609 * MPC2 x NQC2 MPC2 x K K x NQC2
610 *
611  DO 80 j = 1, nqc2
612  CALL clacgv( k, work( ipv+(j-1)*lv ), 1 )
613  80 CONTINUE
614  IF( ioffc2.GT.0 )
615  $ CALL cgemm( 'No transpose', 'No transpose', mpc2, nqc2, k,
616  $ -one, work( ipw ), lw, work( ipv ), lv, one,
617  $ c( ioffc2 ), ldc )
618 *
619  END IF
620 *
621  RETURN
622 *
623 * End of PCLARZB
624 *
625  END
max
#define max(A, B)
Definition: pcgemr.c:180
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pclarzb
subroutine pclarzb(SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, IV, JV, DESCV, T, C, IC, JC, DESCC, WORK)
Definition: pclarzb.f:3
pbcmatadd
subroutine pbcmatadd(ICONTXT, MODE, M, N, ALPHA, A, LDA, BETA, B, LDB)
Definition: pbcmatadd.f:3
pbctran
subroutine pbctran(ICONTXT, ADIST, TRANS, M, N, NB, A, LDA, BETA, C, LDC, IAROW, IACOL, ICROW, ICCOL, WORK)
Definition: pbctran.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
min
#define min(A, B)
Definition: pcgemr.c:181