ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pslarzb.f
Go to the documentation of this file.
1  SUBROUTINE pslarzb( 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 1.7) --
5 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6 * and University of California, Berkeley.
7 * March 14, 2000
8 *
9 * .. Scalar Arguments ..
10  CHARACTER DIRECT, SIDE, STOREV, TRANS
11  INTEGER IC, IV, JC, JV, K, L, M, N
12 * ..
13 * .. Array Arguments ..
14  INTEGER DESCC( * ), DESCV( * )
15  REAL C( * ), T( * ), V( * ), WORK( * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * PSLARZB applies a real block reflector Q or its transpose Q**T to
22 * a real distributed M-by-N matrix sub( C ) = C(IC:IC+M-1,JC:JC+N-1)
23 * from the left or the right.
24 *
25 * Q is a product of k elementary reflectors as returned by PSTZRZF.
26 *
27 * Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
28 *
29 * Notes
30 * =====
31 *
32 * Each global data object is described by an associated description
33 * vector. This vector stores the information required to establish
34 * the mapping between an object element and its corresponding process
35 * and memory location.
36 *
37 * Let A be a generic term for any 2D block cyclicly distributed array.
38 * Such a global array has an associated description vector DESCA.
39 * In the following comments, the character _ should be read as
40 * "of the global array".
41 *
42 * NOTATION STORED IN EXPLANATION
43 * --------------- -------------- --------------------------------------
44 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
45 * DTYPE_A = 1.
46 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
47 * the BLACS process grid A is distribu-
48 * ted over. The context itself is glo-
49 * bal, but the handle (the integer
50 * value) may vary.
51 * M_A (global) DESCA( M_ ) The number of rows in the global
52 * array A.
53 * N_A (global) DESCA( N_ ) The number of columns in the global
54 * array A.
55 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
56 * the rows of the array.
57 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
58 * the columns of the array.
59 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
60 * row of the array A is distributed.
61 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
62 * first column of the array A is
63 * distributed.
64 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
65 * array. LLD_A >= MAX(1,LOCr(M_A)).
66 *
67 * Let K be the number of rows or columns of a distributed matrix,
68 * and assume that its process grid has dimension p x q.
69 * LOCr( K ) denotes the number of elements of K that a process
70 * would receive if K were distributed over the p processes of its
71 * process column.
72 * Similarly, LOCc( K ) denotes the number of elements of K that a
73 * process would receive if K were distributed over the q processes of
74 * its process row.
75 * The values of LOCr() and LOCc() may be determined via a call to the
76 * ScaLAPACK tool function, NUMROC:
77 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
78 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
79 * An upper bound for these quantities may be computed by:
80 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
81 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
82 *
83 * Arguments
84 * =========
85 *
86 * SIDE (global input) CHARACTER
87 * = 'L': apply Q or Q**T from the Left;
88 * = 'R': apply Q or Q**T from the Right.
89 *
90 * TRANS (global input) CHARACTER
91 * = 'N': No transpose, apply Q;
92 * = 'T': Transpose, apply Q**T.
93 *
94 * DIRECT (global input) CHARACTER
95 * Indicates how H is formed from a product of elementary
96 * reflectors
97 * = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
98 * = 'B': H = H(k) . . . H(2) H(1) (Backward)
99 *
100 * STOREV (global input) CHARACTER
101 * Indicates how the vectors which define the elementary
102 * reflectors are stored:
103 * = 'C': Columnwise (not supported yet)
104 * = 'R': Rowwise
105 *
106 * M (global input) INTEGER
107 * The number of rows to be operated on i.e the number of rows
108 * of the distributed submatrix sub( C ). M >= 0.
109 *
110 * N (global input) INTEGER
111 * The number of columns to be operated on i.e the number of
112 * columns of the distributed submatrix sub( C ). N >= 0.
113 *
114 * K (global input) INTEGER
115 * The order of the matrix T (= the number of elementary
116 * reflectors whose product defines the block reflector).
117 *
118 * L (global input) INTEGER
119 * The columns of the distributed submatrix sub( A ) containing
120 * the meaningful part of the Householder reflectors.
121 * If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
122 *
123 * V (local input) REAL pointer into the local memory
124 * to an array of dimension (LLD_V, LOCc(JV+M-1)) if SIDE = 'L',
125 * (LLD_V, LOCc(JV+N-1)) if SIDE = 'R'. It contains the local
126 * pieces of the distributed vectors V representing the
127 * Householder transformation as returned by PSTZRZF.
128 * LLD_V >= LOCr(IV+K-1).
129 *
130 * IV (global input) INTEGER
131 * The row index in the global array V indicating the first
132 * row of sub( V ).
133 *
134 * JV (global input) INTEGER
135 * The column index in the global array V indicating the
136 * first column of sub( V ).
137 *
138 * DESCV (global and local input) INTEGER array of dimension DLEN_.
139 * The array descriptor for the distributed matrix V.
140 *
141 * T (local input) REAL array, dimension MB_V by MB_V
142 * The lower triangular matrix T in the representation of the
143 * block reflector.
144 *
145 * C (local input/local output) REAL pointer into the
146 * local memory to an array of dimension (LLD_C,LOCc(JC+N-1)).
147 * On entry, the M-by-N distributed matrix sub( C ). On exit,
148 * sub( C ) is overwritten by Q*sub( C ) or Q'*sub( C ) or
149 * sub( C )*Q or sub( C )*Q'.
150 *
151 * IC (global input) INTEGER
152 * The row index in the global array C indicating the first
153 * row of sub( C ).
154 *
155 * JC (global input) INTEGER
156 * The column index in the global array C indicating the
157 * first column of sub( C ).
158 *
159 * DESCC (global and local input) INTEGER array of dimension DLEN_.
160 * The array descriptor for the distributed matrix C.
161 *
162 * WORK (local workspace) REAL array, dimension (LWORK)
163 * If STOREV = 'C',
164 * if SIDE = 'L',
165 * LWORK >= ( NqC0 + MpC0 ) * K
166 * else if SIDE = 'R',
167 * LWORK >= ( NqC0 + MAX( NpV0 + NUMROC( NUMROC( N+ICOFFC,
168 * NB_V, 0, 0, NPCOL ), NB_V, 0, 0, LCMQ ),
169 * MpC0 ) ) * K
170 * end if
171 * else if STOREV = 'R',
172 * if SIDE = 'L',
173 * LWORK >= ( MpC0 + MAX( MqV0 + NUMROC( NUMROC( M+IROFFC,
174 * MB_V, 0, 0, NPROW ), MB_V, 0, 0, LCMP ),
175 * NqC0 ) ) * K
176 * else if SIDE = 'R',
177 * LWORK >= ( MpC0 + NqC0 ) * K
178 * end if
179 * end if
180 *
181 * where LCMQ = LCM / NPCOL with LCM = ICLM( NPROW, NPCOL ),
182 *
183 * IROFFV = MOD( IV-1, MB_V ), ICOFFV = MOD( JV-1, NB_V ),
184 * IVROW = INDXG2P( IV, MB_V, MYROW, RSRC_V, NPROW ),
185 * IVCOL = INDXG2P( JV, NB_V, MYCOL, CSRC_V, NPCOL ),
186 * MqV0 = NUMROC( M+ICOFFV, NB_V, MYCOL, IVCOL, NPCOL ),
187 * NpV0 = NUMROC( N+IROFFV, MB_V, MYROW, IVROW, NPROW ),
188 *
189 * IROFFC = MOD( IC-1, MB_C ), ICOFFC = MOD( JC-1, NB_C ),
190 * ICROW = INDXG2P( IC, MB_C, MYROW, RSRC_C, NPROW ),
191 * ICCOL = INDXG2P( JC, NB_C, MYCOL, CSRC_C, NPCOL ),
192 * MpC0 = NUMROC( M+IROFFC, MB_C, MYROW, ICROW, NPROW ),
193 * NpC0 = NUMROC( N+ICOFFC, MB_C, MYROW, ICROW, NPROW ),
194 * NqC0 = NUMROC( N+ICOFFC, NB_C, MYCOL, ICCOL, NPCOL ),
195 *
196 * ILCM, INDXG2P and NUMROC are ScaLAPACK tool functions;
197 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling
198 * the subroutine BLACS_GRIDINFO.
199 *
200 * Alignment requirements
201 * ======================
202 *
203 * The distributed submatrices V(IV:*, JV:*) and C(IC:IC+M-1,JC:JC+N-1)
204 * must verify some alignment properties, namely the following
205 * expressions should be true:
206 *
207 * If STOREV = 'Columnwise'
208 * If SIDE = 'Left',
209 * ( MB_V.EQ.MB_C .AND. IROFFV.EQ.IROFFC .AND. IVROW.EQ.ICROW )
210 * If SIDE = 'Right',
211 * ( MB_V.EQ.NB_C .AND. IROFFV.EQ.ICOFFC )
212 * else if STOREV = 'Rowwise'
213 * If SIDE = 'Left',
214 * ( NB_V.EQ.MB_C .AND. ICOFFV.EQ.IROFFC )
215 * If SIDE = 'Right',
216 * ( NB_V.EQ.NB_C .AND. ICOFFV.EQ.ICOFFC .AND. IVCOL.EQ.ICCOL )
217 * end if
218 *
219 * =====================================================================
220 *
221 * .. Parameters ..
222  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
223  $ lld_, mb_, m_, nb_, n_, rsrc_
224  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
225  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
226  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
227  REAL ONE, ZERO
228  parameter( one = 1.0e+0, zero = 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, 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, infog2l,
244  $ pbsmatadd, pbstran, pb_topget, pxerbla,
245  $ sgebr2d, sgebs2d, sgemm,
246  $ sgsum2d, slamov, slaset, strbr2d,
247  $ strbs2d, strmm
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, 'PSLARZB', -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 = 'T'
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 PBSTRAN
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 PBSTRAN
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 slamov( 'All', k, mqv, v( ioffv ), ldv,
384  $ work( ipw+icoffv*lw ), lw )
385  ELSE
386  CALL slamov( '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 pbstran( ictxt, 'Rowwise', 'Transpose', k, m+icoffv,
394  $ 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 sgemm( 'Transpose', 'No transpose', nqc2, k, mpc2,
410  $ one, c( ioffc2 ), ldc, work( ipv ), lv, zero,
411  $ work( ipw ), lw )
412  ELSE
413  CALL slaset( '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 pbsmatadd( 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 sgsum2d( 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 strbs2d( ictxt, 'Rowwise', ' ', 'Lower', 'Non unit',
449  $ k, k, t, mbv )
450  ELSE
451  CALL strbr2d( ictxt, 'Rowwise', ' ', 'Lower', 'Non unit',
452  $ k, k, t, mbv, myrow, ivcol )
453  END IF
454  CALL strmm( 'Right', 'Lower', transt, 'Non unit', nqc2, k,
455  $ one, t, mbv, work( ipw ), lw )
456 *
457  CALL sgebs2d( ictxt, 'Columnwise', ' ', nqc2, k,
458  $ work( ipw ), lw )
459  ELSE
460  CALL sgebr2d( 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 pbsmatadd( 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  CALL sgemm( 'No transpose', 'Transpose', mpc2, nqc2, k, -one,
491  $ work( ipv ), lv, work( ipw ), lw, one,
492  $ c( ioffc2 ), ldc )
493 *
494  ELSE
495 *
496 * Form sub( C ) * Q or sub( C ) * Q'
497 *
498 * Locally V( IOFFV ) is K x NQV, C( IOFFC2 ) is MPC2 x NQC2
499 * WORK( IPV ) is K x NQV = V( IOFFV ), NQV = NQC2
500 * WORK( IPW ) is MPC2 x K = C( IOFFC2 ) * V( IOFFV )'
501 *
502  ipv = 1
503  ipw = ipv + k * nqc2
504  lv = max( 1, k )
505  lw = max( 1, mpc2 )
506 *
507 * Broadcast V to the other process rows.
508 *
509  CALL pb_topget( ictxt, 'Broadcast', 'Columnwise', colbtop )
510  IF( myrow.EQ.ivrow ) THEN
511  CALL sgebs2d( ictxt, 'Columnwise', colbtop, k, nqc2,
512  $ v( ioffv ), ldv )
513  IF( mycol.EQ.ivcol )
514  $ CALL strbs2d( ictxt, 'Columnwise', colbtop, 'Lower',
515  $ 'Non unit', k, k, t, mbv )
516  CALL slamov( 'All', k, nqc2, v( ioffv ), ldv, work( ipv ),
517  $ lv )
518  ELSE
519  CALL sgebr2d( ictxt, 'Columnwise', colbtop, k, nqc2,
520  $ work( ipv ), lv, ivrow, mycol )
521  IF( mycol.EQ.ivcol )
522  $ CALL strbr2d( ictxt, 'Columnwise', colbtop, 'Lower',
523  $ 'Non unit', k, k, t, mbv, ivrow, mycol )
524  END IF
525 *
526 * WORK( IPV ) is K x NQC2 = V = V( IOFFV )
527 * WORK( IPW ) = C( IOFFC2 ) * V' (MPC2 x NQC2 x K) -> MPC2 x K
528 *
529  IF( nqc2.GT.0 ) THEN
530  CALL sgemm( 'No Transpose', 'Transpose', mpc2, k, nqc2,
531  $ one, c( ioffc2 ), ldc, work( ipv ), lv, zero,
532  $ work( ipw ), lw )
533  ELSE
534  CALL slaset( 'All', mpc2, k, zero, zero, work( ipw ), lw )
535  END IF
536 *
537 * WORK( IPW ) = WORK( IPW ) + C1 ( MPC1 = MPC2 )
538 *
539  IF( nqc1.GT.0 ) THEN
540  mydist = mod( mycol-iccol1+npcol, npcol )
541  ileft = max( 0, mydist * nbc - icoffc1 )
542  jjbeg = jjc1
543  jjend = jjc1 + nqc1 - 1
544  jjnxt = min( iceil( jjbeg, nbc ) * nbc, jjend )
545 *
546  30 CONTINUE
547  IF( jjbeg.LE.jjnxt ) THEN
548  CALL pbsmatadd( ictxt, 'No transpose', mpc2,
549  $ jjnxt-jjbeg+1, one,
550  $ c( iic1+(jjbeg-1)*ldc ), ldc, one,
551  $ work( ipw+ileft*lw ), lw )
552  mydist = mydist + npcol
553  ileft = mydist * nbc - icoffc1
554  jjbeg = jjnxt +1
555  jjnxt = min( jjnxt+nbc, jjend )
556  GO TO 30
557  END IF
558  END IF
559 *
560  CALL sgsum2d( ictxt, 'Rowwise', ' ', mpc2, k, work( ipw ),
561  $ lw, myrow, ivcol )
562 *
563 * WORK( IPW ) = WORK( IPW ) * T' or WORK( IPW ) * T
564 *
565  IF( mycol.EQ.ivcol ) THEN
566  CALL strmm( 'Right', 'Lower', trans, 'Non unit', mpc2, k,
567  $ one, t, mbv, work( ipw ), lw )
568  CALL sgebs2d( ictxt, 'Rowwise', ' ', mpc2, k, work( ipw ),
569  $ lw )
570  ELSE
571  CALL sgebr2d( ictxt, 'Rowwise', ' ', mpc2, k, work( ipw ),
572  $ lw, myrow, ivcol )
573  END IF
574 *
575 * C1 = C1 - WORK( IPW )
576 *
577  IF( nqc1.GT.0 ) THEN
578  mydist = mod( mycol-iccol1+npcol, npcol )
579  ileft = max( 0, mydist * nbc - icoffc1 )
580  jjbeg = jjc1
581  jjend = jjc1 + nqc1 - 1
582  jjnxt = min( iceil( jjbeg, nbc ) * nbc, jjend )
583 *
584  40 CONTINUE
585  IF( jjbeg.LE.jjnxt ) THEN
586  CALL pbsmatadd( ictxt, 'No transpose', mpc2,
587  $ jjnxt-jjbeg+1, -one,
588  $ work( ipw+ileft*lw ), lw, one,
589  $ c( iic1+(jjbeg-1)*ldc ), ldc )
590  mydist = mydist + npcol
591  ileft = mydist * nbc - icoffc1
592  jjbeg = jjnxt +1
593  jjnxt = min( jjnxt+nbc, jjend )
594  GO TO 40
595  END IF
596  END IF
597 *
598 * C2 C2 - W * V
599 * C( IOFFC ) = C( IOFFC ) - WORK( IPW ) * WORK( IPV )
600 * MPC2 x NQC2 MPC2 x K K x NQC2
601 *
602  IF( ioffc2.GT.0 )
603  $ CALL sgemm( 'No transpose', 'No transpose', mpc2, nqc2, k,
604  $ -one, work( ipw ), lw, work( ipv ), lv, one,
605  $ c( ioffc2 ), ldc )
606 *
607  END IF
608 *
609  RETURN
610 *
611 * End of PSLARZB
612 *
613  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
pslarzb
subroutine pslarzb(SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, IV, JV, DESCV, T, C, IC, JC, DESCC, WORK)
Definition: pslarzb.f:3
pbstran
subroutine pbstran(ICONTXT, ADIST, TRANS, M, N, NB, A, LDA, BETA, C, LDC, IAROW, IACOL, ICROW, ICCOL, WORK)
Definition: pbstran.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
pbsmatadd
subroutine pbsmatadd(ICONTXT, MODE, M, N, ALPHA, A, LDA, BETA, B, LDB)
Definition: pbsmatadd.f:3
min
#define min(A, B)
Definition: pcgemr.c:181