SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
subroutine pbsmatadd(icontxt, mode, m, n, alpha, a, lda, beta, b, ldb)
Definition pbsmatadd.f:3
subroutine pbstran(icontxt, adist, trans, m, n, nb, a, lda, beta, c, ldc, iarow, iacol, icrow, iccol, work)
Definition pbstran.f:3
#define max(A, B)
Definition pcgemr.c:180
#define min(A, B)
Definition pcgemr.c:181
subroutine pslarzb(side, trans, direct, storev, m, n, k, l, v, iv, jv, descv, t, c, ic, jc, descc, work)
Definition pslarzb.f:3
subroutine pxerbla(ictxt, srname, info)
Definition pxerbla.f:2