Routine: PDLARFB()  File: SRC\pdlarfb.f

 
 
# lines: 886
  # code: 886
  # comment: 0
  # blank:0
# Variables:89
# Callers:13
# Callings:0
# Words:505
# Keywords:286
 

 

 
..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PDLARFB applies a real block reflector Q or its transpose Q**T to a
  real distributed M-by-N matrix sub( C ) = C(IC:IC+M-1,JC:JC+N-1)
  from the left or the right.
  Notes
  =====
  Each global data object is described by an associated description
  vector.  This vector stores the information required to establish
  the mapping between an object element and its corresponding process
  and memory location.
  Let A be a generic term for any 2D block cyclicly distributed array.
  Such a global array has an associated description vector DESCA.
  In the following comments, the character _ should be read as
  "of the global array".
  NOTATION        STORED IN      EXPLANATION
  --------------- -------------- --------------------------------------
  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
                                 DTYPE_A = 1.
  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
                                 the BLACS process grid A is distribu-
                                 ted over. The context itself is glo-
                                 bal, but the handle (the integer
                                 value) may vary.
  M_A    (global) DESCA( M_ )    The number of rows in the global
                                 array A.
  N_A    (global) DESCA( N_ )    The number of columns in the global
                                 array A.
  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
                                 the rows of the array.
  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
                                 the columns of the array.
  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
                                 row of the array A is distributed.
  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
                                 first column of the array A is
                                 distributed.
  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
  Let K be the number of rows or columns of a distributed matrix,
  and assume that its process grid has dimension p x q.
  LOCr( K ) denotes the number of elements of K that a process
  would receive if K were distributed over the p processes of its
  process column.
  Similarly, LOCc( K ) denotes the number of elements of K that a
  process would receive if K were distributed over the q processes of
  its process row.
  The values of LOCr() and LOCc() may be determined via a call to the
  ScaLAPACK tool function, NUMROC:
          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
  An upper bound for these quantities may be computed by:
          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
  Arguments
  =========
  SIDE    (global input) CHARACTER
          = 'L': apply Q or Q**T from the Left;
          = 'R': apply Q or Q**T from the Right.
  TRANS   (global input) CHARACTER
          = 'N':  No transpose, apply Q;
          = 'T':  Transpose, apply Q**T.
  DIRECT  (global input) CHARACTER
          Indicates how Q is formed from a product of elementary
          reflectors
          = 'F': Q = H(1) H(2) . . . H(k) (Forward)
          = 'B': Q = H(k) . . . H(2) H(1) (Backward)
  STOREV  (global input) CHARACTER
          Indicates how the vectors which define the elementary
          reflectors are stored:
          = 'C': Columnwise
          = 'R': Rowwise
  M       (global input) INTEGER
          The number of rows to be operated on i.e the number of rows
          of the distributed submatrix sub( C ). M >= 0.
  N       (global input) INTEGER
          The number of columns to be operated on i.e the number of
          columns of the distributed submatrix sub( C ). N >= 0.
  K       (global input) INTEGER
          The order of the matrix T (= the number of elementary
          reflectors whose product defines the block reflector).
  V       (local input) DOUBLE PRECISION pointer into the local memory
          to an array of dimension ( LLD_V, LOCc(JV+K-1) ) if
          STOREV = 'C', ( LLD_V, LOCc(JV+M-1)) if STOREV = 'R' and
          SIDE = 'L', ( LLD_V, LOCc(JV+N-1) ) if STOREV = 'R' and
          SIDE = 'R'. It contains the local pieces of the distributed
          vectors V representing the Householder transformation.
          See further details.
          If STOREV = 'C' and SIDE = 'L', LLD_V >= MAX(1,LOCr(IV+M-1));
          if STOREV = 'C' and SIDE = 'R', LLD_V >= MAX(1,LOCr(IV+N-1));
          if STOREV = 'R', LLD_V >= LOCr(IV+K-1).
  IV      (global input) INTEGER
          The row index in the global array V indicating the first
          row of sub( V ).
  JV      (global input) INTEGER
          The column index in the global array V indicating the
          first column of sub( V ).
  DESCV   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix V.
  T       (local input) DOUBLE PRECISION array, dimension MB_V by MB_V
          if STOREV = 'R' and NB_V by NB_V if STOREV = 'C'. The trian-
          gular matrix T in the representation of the block reflector.
  C       (local input/local output) DOUBLE PRECISION pointer into the
          local memory to an array of dimension (LLD_C,LOCc(JC+N-1)).
          On entry, the M-by-N distributed matrix sub( C ). On exit,
          sub( C ) is overwritten by Q*sub( C ) or Q'*sub( C ) or
          sub( C )*Q or sub( C )*Q'.
  IC      (global input) INTEGER
          The row index in the global array C indicating the first
          row of sub( C ).
  JC      (global input) INTEGER
          The column index in the global array C indicating the
          first column of sub( C ).
  DESCC   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix C.
  WORK    (local workspace) DOUBLE PRECISION array, dimension (LWORK)
          If STOREV = 'C',
            if SIDE = 'L',
              LWORK >= ( NqC0 + MpC0 ) * K
            else if SIDE = 'R',
              LWORK >= ( NqC0 + MAX( NpV0 + NUMROC( NUMROC( N+ICOFFC,
                         NB_V, 0, 0, NPCOL ), NB_V, 0, 0, LCMQ ),
                         MpC0 ) ) * K
            end if
          else if STOREV = 'R',
            if SIDE = 'L',
              LWORK >= ( MpC0 + MAX( MqV0 + NUMROC( NUMROC( M+IROFFC,
                         MB_V, 0, 0, NPROW ), MB_V, 0, 0, LCMP ),
                         NqC0 ) ) * K
            else if SIDE = 'R',
              LWORK >= ( MpC0 + NqC0 ) * K
            end if
          end if
          where LCMQ = LCM / NPCOL with LCM = ICLM( NPROW, NPCOL ),
          IROFFV = MOD( IV-1, MB_V ), ICOFFV = MOD( JV-1, NB_V ),
          IVROW = INDXG2P( IV, MB_V, MYROW, RSRC_V, NPROW ),
          IVCOL = INDXG2P( JV, NB_V, MYCOL, CSRC_V, NPCOL ),
          MqV0 = NUMROC( M+ICOFFV, NB_V, MYCOL, IVCOL, NPCOL ),
          NpV0 = NUMROC( N+IROFFV, MB_V, MYROW, IVROW, NPROW ),
          IROFFC = MOD( IC-1, MB_C ), ICOFFC = MOD( JC-1, NB_C ),
          ICROW = INDXG2P( IC, MB_C, MYROW, RSRC_C, NPROW ),
          ICCOL = INDXG2P( JC, NB_C, MYCOL, CSRC_C, NPCOL ),
          MpC0 = NUMROC( M+IROFFC, MB_C, MYROW, ICROW, NPROW ),
          NpC0 = NUMROC( N+ICOFFC, MB_C, MYROW, ICROW, NPROW ),
          NqC0 = NUMROC( N+ICOFFC, NB_C, MYCOL, ICCOL, NPCOL ),
          ILCM, INDXG2P and NUMROC are ScaLAPACK tool functions;
          MYROW, MYCOL, NPROW and NPCOL can be determined by calling
          the subroutine BLACS_GRIDINFO.
  Alignment requirements
  ======================
  The distributed submatrices V(IV:*, JV:*) and C(IC:IC+M-1,JC:JC+N-1)
  must verify some alignment properties, namely the following
  expressions should be true:
  If STOREV = 'Columnwise'
    If SIDE = 'Left',
      ( MB_V.EQ.MB_C .AND. IROFFV.EQ.IROFFC .AND. IVROW.EQ.ICROW )
    If SIDE = 'Right',
      ( MB_V.EQ.NB_C .AND. IROFFV.EQ.ICOFFC )
  else if STOREV = 'Rowwise'
    If SIDE = 'Left',
      ( NB_V.EQ.MB_C .AND. ICOFFV.EQ.IROFFC )
    If SIDE = 'Right',
      ( NB_V.EQ.NB_C .AND. ICOFFV.EQ.ICOFFC .AND. IVCOL.EQ.ICCOL )
  end if
  =====================================================================
     .. Parameters ..
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE PDLARFB( SIDE , TRANS , DIRECT , STOREV , M , N , K , V , IV ,
002       $JV , DESCV , T , C , IC , JC , DESCC , WORK )
003  
004  *     -- ScaLAPACK auxiliary routine(version 1.7) --
005  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
006  *     and University of California , Berkeley.
007  *     May 1 , 1997
008  
009  *     .. Scalar Arguments ..
010        CHARACTER SIDE , TRANS , DIRECT , STOREV
011        INTEGER IC , IV , JC , JV , K , M , N
012        INTEGER BLOCK_CYCLIC_2D , CSRC_ , CTXT_ , DLEN_ , DTYPE_ ,
013       $LLD_ , MB_ , M_ , NB_ , N_ , RSRC_
014        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
015       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
016       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
017        DOUBLE PRECISION ONE , ZERO
018        PARAMETER( ONE = 1.0D + 0 , ZERO = 0.0D + 0 )
019  *     ..
020  *     .. Local Scalars ..
021        LOGICAL FORWARD
022        CHARACTER COLBTOP , ROWBTOP , TRANST , UPLO
023        INTEGER HEIGHT , IBASE , ICCOL , ICOFFC , ICOFFV , ICROW ,
024       $ICTXT , II , IIBEG , IIC , IIEND , IINXT , IIV ,
025       $ILASTCOL , ILASTROW , ILEFT , IOFF , IOFFC , IOFFV ,
026       $IPT , IPV , IPW , IPW1 , IRIGHT , IROFFC , IROFFV ,
027       $ITOP , IVCOL , IVROW , JJ , JJBEG , JJC , JJEND ,
028       $JJNXT , JJV , KP , KQ , LDC , LDV , LV , LW , MBV , MPC ,
029       $MPC0 , MQV , MQV0 , MYCOL , MYDIST , MYROW , NBV ,
030       $NPV , NPV0 , NPCOL , NPROW , NQC , NQC0 , WIDE
031  *     ..
032  *     .. External Subroutines ..
033        EXTERNAL BLACS_GRIDINFO , DGEBR2D , DGEBS2D , DGEMM ,
034       $DGSUM2D , DLACPY , DLASET , DTRBR2D ,
035       $DTRBS2D , DTRMM , INFOG1L , INFOG2L , PB_TOPGET ,
036       $PBDTRAN
037  *     ..
038  *     .. Intrinsic Functions ..
039        INTRINSIC MAX , MIN , MOD
040  *     ..
041  *     .. External Functions ..
042        LOGICAL LSAME
043        INTEGER ICEIL , NUMROC
044        EXTERNAL ICEIL , LSAME , NUMROC
045  *     ..
046  *     .. Executable Statements ..
047  
048  *     Quick return if possible
049  
050        IF( M.LE.0 .OR. N.LE.0 .OR. K.LE.0 )
051       $    RETURN
052  
053  *         Get grid parameters
054  
055            ICTXT = DESCC( CTXT_ )
056            CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
057  
058            IF( LSAME( TRANS , 'N' ) ) THEN
059                TRANST = 'T'
060            ELSE
061                TRANST = 'N'
062            END IF
063            FORWARD = LSAME( DIRECT , 'F' )
064            IF( FORWARD ) THEN
065                UPLO = 'U'
066            ELSE
067                UPLO = 'L'
068            END IF
069  
070            CALL INFOG2L( IV , JV , DESCV , NPROW , NPCOL , MYROW , MYCOL , IIV , JJV ,
071       $    IVROW , IVCOL )
072            CALL INFOG2L( IC , JC , DESCC , NPROW , NPCOL , MYROW , MYCOL , IIC , JJC ,
073       $    ICROW , ICCOL )
074            LDC = DESCC( LLD_ )
075            LDV = DESCV( LLD_ )
076            IIC = MIN( IIC , LDC )
077            IIV = MIN( IIV , LDV )
078            IROFFC = MOD( IC - 1 , DESCC( MB_ ) )
079            ICOFFC = MOD( JC - 1 , DESCC( NB_ ) )
080            MBV = DESCV( MB_ )
081            NBV = DESCV( NB_ )
082            IROFFV = MOD( IV - 1 , MBV )
083            ICOFFV = MOD( JV - 1 , NBV )
084            MPC = NUMROC( M + IROFFC , DESCC( MB_ ) , MYROW , ICROW , NPROW )
085            NQC = NUMROC( N + ICOFFC , DESCC( NB_ ) , MYCOL , ICCOL , NPCOL )
086            IF( MYCOL.EQ.ICCOL )
087       $        NQC = NQC - ICOFFC
088                IF( MYROW.EQ.ICROW )
089       $            MPC = MPC - IROFFC
090                    JJC = MIN( JJC , MAX( 1 , JJC + NQC - 1 ) )
091                    JJV = MIN( JJV , MAX( 1 , NUMROC( DESCV( N_ ) , NBV , MYCOL ,
092       $            DESCV( CSRC_ ) , NPCOL ) ) )
093                    IOFFC = IIC + ( JJC - 1 ) * LDC
094                    IOFFV = IIV + ( JJV - 1 ) * LDV
095  
096                    IF( LSAME( STOREV , 'C' ) ) THEN
097  
098  *                     V is stored columnwise
099  
100                        IF( LSAME( SIDE , 'L' ) ) THEN
101  
102  *                         Form Q*sub( C ) or Q'*sub( C )
103  
104  *                         Locally V( IOFFV ) is MPV x K , C( IOFFC ) is MPC x NQC
105  *                         WORK( IPV ) is MPC x K = V( IOFFV ) , MPC = MPV
106  *                         WORK( IPW ) is NQC x K = C( IOFFC )' * V( IOFFV )
107  
108                            IPV = 1
109                            IPW = IPV + MPC * K
110                            LV = MAX( 1 , MPC )
111                            LW = MAX( 1 , NQC )
112  
113  *                         Broadcast V to the other process columns.
114  
115                            CALL PB_TOPGET( ICTXT , 'Broadcast' , 'Rowwise' , ROWBTOP )
116                            IF( MYCOL.EQ.IVCOL ) THEN
117                                CALL DGEBS2D( ICTXT , 'Rowwise' , ROWBTOP , MPC , K ,
118       $                        V( IOFFV ) , LDV )
119                                IF( MYROW.EQ.IVROW )
120       $                            CALL DTRBS2D( ICTXT , 'Rowwise' , ROWBTOP , UPLO ,
121       $                            'Non unit' , K , K , T , NBV )
122                                    CALL DLACPY( 'All' , MPC , K , V( IOFFV ) , LDV , WORK( IPV ) ,
123       $                            LV )
124                                ELSE
125                                    CALL DGEBR2D( ICTXT , 'Rowwise' , ROWBTOP , MPC , K ,
126       $                            WORK( IPV ) , LV , MYROW , IVCOL )
127                                    IF( MYROW.EQ.IVROW )
128       $                                CALL DTRBR2D( ICTXT , 'Rowwise' , ROWBTOP , UPLO ,
129       $                                'Non unit' , K , K , T , NBV , MYROW , IVCOL )
130                                    END IF
131  
132                                    IF( FORWARD ) THEN
133  
134  *                                     WORK(IPV) =( V1 ) where V1 is unit lower triangular ,
135  *                                     ( V2 ) zeroes upper triangular part of V1
136  
137                                        MYDIST = MOD( MYROW - IVROW + NPROW , NPROW )
138                                        ITOP = MAX( 0 , MYDIST*MBV - IROFFV )
139                                        IIBEG = IIV
140                                        IIEND = IIBEG + MPC - 1
141                                        IINXT = MIN( ICEIL( IIBEG , MBV )*MBV , IIEND )
142  
143     10 CONTINUE
144        IF( K - ITOP .GT.0 ) THEN
145            CALL DLASET( 'Upper' , IINXT - IIBEG + 1 , K - ITOP , ZERO ,
146       $    ONE , WORK( IPV + IIBEG - IIV + ITOP*LV ) , LV )
147            MYDIST = MYDIST + NPROW
148            ITOP = MYDIST * MBV - IROFFV
149            IIBEG = IINXT + 1
150            IINXT = MIN( IINXT + MBV , IIEND )
151            GO TO 10
152        END IF
153  
154        ELSE
155  
156  *         WORK(IPV) =( V1 ) where V2 is unit upper triangular ,
157  *         ( V2 ) zeroes lower triangular part of V2
158  
159            JJ = JJV
160            IOFF = MOD( IV + M - K - 1 , MBV )
161            CALL INFOG1L( IV + M - K , MBV , NPROW , MYROW , DESCV( RSRC_ ) ,
162       $    II , ILASTROW )
163            KP = NUMROC( K + IOFF , MBV , MYROW , ILASTROW , NPROW )
164            IF( MYROW.EQ.ILASTROW )
165       $        KP = KP - IOFF
166                MYDIST = MOD( MYROW - ILASTROW + NPROW , NPROW )
167                ITOP = MYDIST * MBV - IOFF
168                IBASE = MIN( ITOP + MBV , K )
169                ITOP = MIN( MAX( 0 , ITOP ) , K )
170  
171     20 CONTINUE
172        IF( JJ.LE.( JJV + K - 1 ) ) THEN
173            HEIGHT = IBASE - ITOP
174            CALL DLASET( 'All' , KP , ITOP - JJ + JJV , ZERO , ZERO ,
175       $    WORK( IPV + II - IIV + (JJ - JJV)*LV ) , LV )
176            CALL DLASET( 'Lower' , KP , HEIGHT , ZERO , ONE ,
177       $    WORK( IPV + II - IIV + ITOP*LV ) , LV )
178            KP = MAX( 0 , KP - HEIGHT )
179            II = II + HEIGHT
180            JJ = JJV + IBASE
181            MYDIST = MYDIST + NPROW
182            ITOP = MYDIST * MBV - IOFF
183            IBASE = MIN( ITOP + MBV , K )
184            ITOP = MIN( ITOP , K )
185            GO TO 20
186        END IF
187  
188        END IF
189  
190  *     WORK( IPW ) = C( IOFFC )' * V(NQC x MPC x K) -> NQC x K
191  
192        IF( MPC.GT.0 ) THEN
193            CALL DGEMM( 'Transpose' , 'No transpose' , NQC , K , MPC ,
194       $    ONE , C( IOFFC ) , LDC , WORK( IPV ) , LV , ZERO ,
195       $    WORK( IPW ) , LW )
196        ELSE
197            CALL DLASET( 'All' , NQC , K , ZERO , ZERO , WORK( IPW ) , LW )
198        END IF
199  
200        CALL DGSUM2D( ICTXT , 'Columnwise' , ' ' , NQC , K , WORK( IPW ) ,
201       $LW , IVROW , MYCOL )
202  
203        IF( MYROW.EQ.IVROW ) THEN
204  
205  *         WORK( IPW ) = WORK( IPW ) * T' or WORK( IPW ) * T
206  
207            CALL DTRMM( 'Right' , UPLO , TRANST , 'Non unit' , NQC , K ,
208       $    ONE , T , NBV , WORK( IPW ) , LW )
209            CALL DGEBS2D( ICTXT , 'Columnwise' , ' ' , NQC , K ,
210       $    WORK( IPW ) , LW )
211        ELSE
212            CALL DGEBR2D( ICTXT , 'Columnwise' , ' ' , NQC , K ,
213       $    WORK( IPW ) , LW , IVROW , MYCOL )
214        END IF
215  
216  *     C C - V * W'
217  *     C( IOFFC ) = C( IOFFC ) - WORK( IPV ) * WORK( IPW )'
218  *     MPC x NQC MPC x K K x NQC
219  
220        CALL DGEMM( 'No transpose' , 'Transpose' , MPC , NQC , K , - ONE ,
221       $WORK( IPV ) , LV , WORK( IPW ) , LW , ONE ,
222       $C( IOFFC ) , LDC )
223  
224        ELSE
225  
226  *         Form sub( C )*Q or sub( C )*Q'
227  
228  *         ICOFFC = IROFFV is required by the current transposition
229  *         routine PBDTRAN
230  
231            NPV0 = NUMROC( N + IROFFV , MBV , MYROW , IVROW , NPROW )
232            IF( MYROW.EQ.IVROW ) THEN
233                NPV = NPV0 - IROFFV
234            ELSE
235                NPV = NPV0
236            END IF
237            IF( MYCOL.EQ.ICCOL ) THEN
238                NQC0 = NQC + ICOFFC
239            ELSE
240                NQC0 = NQC
241            END IF
242  
243  *         Locally V( IOFFV ) is NPV x K C( IOFFC ) is MPC x NQC
244  *         WORK( IPV ) is K x NQC0 =[ . V( IOFFV ) ]'
245  *         WORK( IPW ) is NPV0 x K =[ . V( IOFFV )' ]'
246  *         WORK( IPT ) is the workspace for PBDTRAN
247  
248            IPV = 1
249            IPW = IPV + K * NQC0
250            IPT = IPW + NPV0 * K
251            LV = MAX( 1 , K )
252            LW = MAX( 1 , NPV0 )
253  
254            IF( MYCOL.EQ.IVCOL ) THEN
255                IF( MYROW.EQ.IVROW ) THEN
256                    CALL DLASET( 'All' , IROFFV , K , ZERO , ZERO ,
257       $            WORK( IPW ) , LW )
258                    IPW1 = IPW + IROFFV
259                    CALL DLACPY( 'All' , NPV , K , V( IOFFV ) , LDV ,
260       $            WORK( IPW1 ) , LW )
261                ELSE
262                    IPW1 = IPW
263                    CALL DLACPY( 'All' , NPV , K , V( IOFFV ) , LDV ,
264       $            WORK( IPW1 ) , LW )
265                END IF
266  
267                IF( FORWARD ) THEN
268  
269  *                 WORK(IPW) =( . V1' V2' )' where V1 is unit lower
270  *                 triangular , zeroes upper triangular part of V1
271  
272                    MYDIST = MOD( MYROW - IVROW + NPROW , NPROW )
273                    ITOP = MAX( 0 , MYDIST*MBV - IROFFV )
274                    IIBEG = IIV
275                    IIEND = IIBEG + NPV - 1
276                    IINXT = MIN( ICEIL( IIBEG , MBV )*MBV , IIEND )
277  
278     30 CONTINUE
279        IF(( K - ITOP ).GT.0 ) THEN
280            CALL DLASET( 'Upper' , IINXT - IIBEG + 1 , K - ITOP , ZERO ,
281       $    ONE , WORK( IPW1 + IIBEG - IIV + ITOP*LW ) ,
282       $    LW )
283            MYDIST = MYDIST + NPROW
284            ITOP = MYDIST * MBV - IROFFV
285            IIBEG = IINXT + 1
286            IINXT = MIN( IINXT + MBV , IIEND )
287            GO TO 30
288        END IF
289  
290        ELSE
291  
292  *         WORK( IPW ) =( . V1' V2' )' where V2 is unit upper
293  *         triangular , zeroes lower triangular part of V2.
294  
295            JJ = JJV
296            CALL INFOG1L( IV + N - K , MBV , NPROW , MYROW ,
297       $    DESCV( RSRC_ ) , II , ILASTROW )
298            IOFF = MOD( IV + N - K - 1 , MBV )
299            KP = NUMROC( K + IOFF , MBV , MYROW , ILASTROW , NPROW )
300            IF( MYROW.EQ.ILASTROW )
301       $        KP = KP - IOFF
302                MYDIST = MOD( MYROW - ILASTROW + NPROW , NPROW )
303                ITOP = MYDIST * MBV - IOFF
304                IBASE = MIN( ITOP + MBV , K )
305                ITOP = MIN( MAX( 0 , ITOP ) , K )
306  
307     40 CONTINUE
308        IF( JJ.LE.( JJV + K - 1 ) ) THEN
309            HEIGHT = IBASE - ITOP
310            CALL DLASET( 'All' , KP , ITOP - JJ + JJV , ZERO , ZERO ,
311       $    WORK( IPW1 + II - IIV + (JJ - JJV)*LW ) , LW )
312            CALL DLASET( 'Lower' , KP , HEIGHT , ZERO , ONE ,
313       $    WORK( IPW1 + II - IIV + ITOP*LW ) , LW )
314            KP = MAX( 0 , KP - HEIGHT )
315            II = II + HEIGHT
316            JJ = JJV + IBASE
317            MYDIST = MYDIST + NPROW
318            ITOP = MYDIST * MBV - IOFF
319            IBASE = MIN( ITOP + MBV , K )
320            ITOP = MIN( ITOP , K )
321            GO TO 40
322        END IF
323        END IF
324        END IF
325  
326        CALL PBDTRAN( ICTXT , 'Columnwise' , 'Transpose' , N + IROFFV , K ,
327       $MBV , WORK( IPW ) , LW , ZERO , WORK( IPV ) , LV ,
328       $IVROW , IVCOL , - 1 , ICCOL , WORK( IPT ) )
329  
330  *     WORK( IPV ) =( . V' ) -> WORK( IPV ) = V' is K x NQC
331  
332        IF( MYCOL.EQ.ICCOL )
333       $    IPV = IPV + ICOFFC * LV
334  
335  *         WORK( IPW ) becomes MPC x K = C( IOFFC ) * V
336  *         WORK( IPW ) = C( IOFFC ) * V(MPC x NQC x K) -> MPC x K
337  
338            LW = MAX( 1 , MPC )
339  
340            IF( NQC.GT.0 ) THEN
341                CALL DGEMM( 'No transpose' , 'Transpose' , MPC , K , NQC ,
342       $        ONE , C( IOFFC ) , LDC , WORK( IPV ) , LV , ZERO ,
343       $        WORK( IPW ) , LW )
344            ELSE
345                CALL DLASET( 'All' , MPC , K , ZERO , ZERO , WORK( IPW ) , LW )
346            END IF
347  
348            CALL DGSUM2D( ICTXT , 'Rowwise' , ' ' , MPC , K , WORK( IPW ) ,
349       $    LW , MYROW , IVCOL )
350  
351  *         WORK( IPW ) = WORK( IPW ) * T' or WORK( IPW ) * T
352  
353            IF( MYCOL.EQ.IVCOL ) THEN
354                IF( MYROW.EQ.IVROW ) THEN
355  
356  *                 Broadcast the block reflector to the other rows.
357  
358                    CALL DTRBS2D( ICTXT , 'Columnwise' , ' ' , UPLO ,
359       $            'Non unit' , K , K , T , NBV )
360                ELSE
361                    CALL DTRBR2D( ICTXT , 'Columnwise' , ' ' , UPLO ,
362       $            'Non unit' , K , K , T , NBV , IVROW , MYCOL )
363                END IF
364                CALL DTRMM( 'Right' , UPLO , TRANS , 'Non unit' , MPC , K ,
365       $        ONE , T , NBV , WORK( IPW ) , LW )
366  
367                CALL DGEBS2D( ICTXT , 'Rowwise' , ' ' , MPC , K , WORK( IPW ) ,
368       $        LW )
369            ELSE
370                CALL DGEBR2D( ICTXT , 'Rowwise' , ' ' , MPC , K , WORK( IPW ) ,
371       $        LW , MYROW , IVCOL )
372            END IF
373  
374  *         C C - W * V'
375  *         C( IOFFC ) = C( IOFFC ) - WORK( IPW ) * WORK( IPV )
376  *         MPC x NQC MPC x K K x NQC
377  
378            CALL DGEMM( 'No transpose' , 'No transpose' , MPC , NQC , K ,
379       $    - ONE , WORK( IPW ) , LW , WORK( IPV ) , LV , ONE ,
380       $    C( IOFFC ) , LDC )
381        END IF
382  
383        ELSE
384  
385  *         V is stored rowwise
386  
387            IF( LSAME( SIDE , 'L' ) ) THEN
388  
389  *             Form Q*sub( C ) or Q'*sub( C )
390  
391  *             IROFFC = ICOFFV is required by the current transposition
392  *             routine PBDTRAN
393  
394                MQV0 = NUMROC( M + ICOFFV , NBV , MYCOL , IVCOL , NPCOL )
395                IF( MYCOL.EQ.IVCOL ) THEN
396                    MQV = MQV0 - ICOFFV
397                ELSE
398                    MQV = MQV0
399                END IF
400                IF( MYROW.EQ.ICROW ) THEN
401                    MPC0 = MPC + IROFFC
402                ELSE
403                    MPC0 = MPC
404                END IF
405  
406  *             Locally V( IOFFV ) is K x MQV , C( IOFFC ) is MPC x NQC
407  *             WORK( IPV ) is MPC0 x K =[ . V( IOFFV ) ]'
408  *             WORK( IPW ) is K x MQV0 =[ . V( IOFFV ) ]
409  *             WORK( IPT ) is the workspace for PBDTRAN
410  
411                IPV = 1
412                IPW = IPV + MPC0 * K
413                IPT = IPW + K * MQV0
414                LV = MAX( 1 , MPC0 )
415                LW = MAX( 1 , K )
416  
417                IF( MYROW.EQ.IVROW ) THEN
418                    IF( MYCOL.EQ.IVCOL ) THEN
419                        CALL DLASET( 'All' , K , ICOFFV , ZERO , ZERO ,
420       $                WORK( IPW ) , LW )
421                        IPW1 = IPW + ICOFFV * LW
422                        CALL DLACPY( 'All' , K , MQV , V( IOFFV ) , LDV ,
423       $                WORK( IPW1 ) , LW )
424                    ELSE
425                        IPW1 = IPW
426                        CALL DLACPY( 'All' , K , MQV , V( IOFFV ) , LDV ,
427       $                WORK( IPW1 ) , LW )
428                    END IF
429  
430                    IF( FORWARD ) THEN
431  
432  *                     WORK( IPW ) =( . V1 V2 ) where V1 is unit upper
433  *                     triangular , zeroes lower triangular part of V1
434  
435                        MYDIST = MOD( MYCOL - IVCOL + NPCOL , NPCOL )
436                        ILEFT = MAX( 0 , MYDIST * NBV - ICOFFV )
437                        JJBEG = JJV
438                        JJEND = JJV + MQV - 1
439                        JJNXT = MIN( ICEIL( JJBEG , NBV ) * NBV , JJEND )
440  
441     50 CONTINUE
442        IF(( K - ILEFT ).GT.0 ) THEN
443            CALL DLASET( 'Lower' , K - ILEFT , JJNXT - JJBEG + 1 , ZERO ,
444       $    ONE ,
445       $    WORK( IPW1 + ILEFT + (JJBEG - JJV)*LW ) ,
446       $    LW )
447            MYDIST = MYDIST + NPCOL
448            ILEFT = MYDIST * NBV - ICOFFV
449            JJBEG = JJNXT + 1
450            JJNXT = MIN( JJNXT + NBV , JJEND )
451            GO TO 50
452        END IF
453  
454        ELSE
455  
456  *         WORK( IPW ) =( . V1 V2 ) where V2 is unit lower
457  *         triangular , zeroes upper triangular part of V2.
458  
459            II = IIV
460            CALL INFOG1L( JV + M - K , NBV , NPCOL , MYCOL ,
461       $    DESCV( CSRC_ ) , JJ , ILASTCOL )
462            IOFF = MOD( JV + M - K - 1 , NBV )
463            KQ = NUMROC( K + IOFF , NBV , MYCOL , ILASTCOL , NPCOL )
464            IF( MYCOL.EQ.ILASTCOL )
465       $        KQ = KQ - IOFF
466                MYDIST = MOD( MYCOL - ILASTCOL + NPCOL , NPCOL )
467                ILEFT = MYDIST * NBV - IOFF
468                IRIGHT = MIN( ILEFT + NBV , K )
469                ILEFT = MIN( MAX( 0 , ILEFT ) , K )
470  
471     60 CONTINUE
472        IF( II.LE.( IIV + K - 1 ) ) THEN
473            WIDE = IRIGHT - ILEFT
474            CALL DLASET( 'All' , ILEFT - II + IIV , KQ , ZERO , ZERO ,
475       $    WORK( IPW1 + II - IIV + (JJ - JJV)*LW ) , LW )
476            CALL DLASET( 'Upper' , WIDE , KQ , ZERO , ONE ,
477       $    WORK( IPW1 + ILEFT + (JJ - JJV)*LW ) , LW )
478            KQ = MAX( 0 , KQ - WIDE )
479            II = IIV + IRIGHT
480            JJ = JJ + WIDE
481            MYDIST = MYDIST + NPCOL
482            ILEFT = MYDIST * NBV - IOFF
483            IRIGHT = MIN( ILEFT + NBV , K )
484            ILEFT = MIN( ILEFT , K )
485            GO TO 60
486        END IF
487        END IF
488        END IF
489  
490  *     WORK( IPV ) = WORK( IPW )'(replicated) is MPC0 x K
491  
492        CALL PBDTRAN( ICTXT , 'Rowwise' , 'Transpose' , K , M + ICOFFV ,
493       $NBV , WORK( IPW ) , LW , ZERO , WORK( IPV ) , LV ,
494       $IVROW , IVCOL , ICROW , - 1 , WORK( IPT ) )
495  
496  *     WORK( IPV ) =( . V )' -> WORK( IPV ) = V' is MPC x K
497  
498        IF( MYROW.EQ.ICROW )
499       $    IPV = IPV + IROFFC
500  
501  *         WORK( IPW ) becomes NQC x K = C( IOFFC )' * V'
502  *         WORK( IPW ) = C( IOFFC )' * V'(NQC x MPC x K) -> NQC x K
503  
504            LW = MAX( 1 , NQC )
505  
506            IF( MPC.GT.0 ) THEN
507                CALL DGEMM( 'Transpose' , 'No transpose' , NQC , K , MPC ,
508       $        ONE , C( IOFFC ) , LDC , WORK( IPV ) , LV , ZERO ,
509       $        WORK( IPW ) , LW )
510            ELSE
511                CALL DLASET( 'All' , NQC , K , ZERO , ZERO , WORK( IPW ) , LW )
512            END IF
513  
514            CALL DGSUM2D( ICTXT , 'Columnwise' , ' ' , NQC , K , WORK( IPW ) ,
515       $    LW , IVROW , MYCOL )
516  
517  *         WORK( IPW ) = WORK( IPW ) * T' or WORK( IPW ) * T
518  
519            IF( MYROW.EQ.IVROW ) THEN
520                IF( MYCOL.EQ.IVCOL ) THEN
521  
522  *                 Broadcast the block reflector to the other columns.
523  
524                    CALL DTRBS2D( ICTXT , 'Rowwise' , ' ' , UPLO , 'Non unit' ,
525       $            K , K , T , MBV )
526                ELSE
527                    CALL DTRBR2D( ICTXT , 'Rowwise' , ' ' , UPLO , 'Non unit' ,
528       $            K , K , T , MBV , MYROW , IVCOL )
529                END IF
530                CALL DTRMM( 'Right' , UPLO , TRANST , 'Non unit' , NQC , K ,
531       $        ONE , T , MBV , WORK( IPW ) , LW )
532  
533                CALL DGEBS2D( ICTXT , 'Columnwise' , ' ' , NQC , K ,
534       $        WORK( IPW ) , LW )
535            ELSE
536                CALL DGEBR2D( ICTXT , 'Columnwise' , ' ' , NQC , K ,
537       $        WORK( IPW ) , LW , IVROW , MYCOL )
538            END IF
539  
540  *         C C - V' * W'
541  *         C( IOFFC ) = C( IOFFC ) - WORK( IPV ) * WORK( IPW )'
542  *         MPC x NQC MPC x K K x NQC
543  
544            CALL DGEMM( 'No transpose' , 'Transpose' , MPC , NQC , K , - ONE ,
545       $    WORK( IPV ) , LV , WORK( IPW ) , LW , ONE ,
546       $    C( IOFFC ) , LDC )
547  
548        ELSE
549  
550  *         Form Q*sub( C ) or Q'*sub( C )
551  
552  *         Locally V( IOFFV ) is K x NQV , C( IOFFC ) is MPC x NQC
553  *         WORK( IPV ) is K x NQV = V( IOFFV ) , NQV = NQC
554  *         WORK( IPW ) is MPC x K = C( IOFFC ) * V( IOFFV )'
555  
556            IPV = 1
557            IPW = IPV + K * NQC
558            LV = MAX( 1 , K )
559            LW = MAX( 1 , MPC )
560  
561  *         Broadcast V to the other process rows.
562  
563            CALL PB_TOPGET( ICTXT , 'Broadcast' , 'Columnwise' , COLBTOP )
564            IF( MYROW.EQ.IVROW ) THEN
565                CALL DGEBS2D( ICTXT , 'Columnwise' , COLBTOP , K , NQC ,
566       $        V( IOFFV ) , LDV )
567                IF( MYCOL.EQ.IVCOL )
568       $            CALL DTRBS2D( ICTXT , 'Columnwise' , COLBTOP , UPLO ,
569       $            'Non unit' , K , K , T , MBV )
570                    CALL DLACPY( 'All' , K , NQC , V( IOFFV ) , LDV , WORK( IPV ) ,
571       $            LV )
572                ELSE
573                    CALL DGEBR2D( ICTXT , 'Columnwise' , COLBTOP , K , NQC ,
574       $            WORK( IPV ) , LV , IVROW , MYCOL )
575                    IF( MYCOL.EQ.IVCOL )
576       $                CALL DTRBR2D( ICTXT , 'Columnwise' , COLBTOP , UPLO ,
577       $                'Non unit' , K , K , T , MBV , IVROW , MYCOL )
578                    END IF
579  
580                    IF( FORWARD ) THEN
581  
582  *                     WORK(IPW) =( V1 V2 ) where V1 is unit upper
583  *                     triangular , zeroes lower triangular part of V1
584  
585                        MYDIST = MOD( MYCOL - IVCOL + NPCOL , NPCOL )
586                        ILEFT = MAX( 0 , MYDIST * NBV - ICOFFV )
587                        JJBEG = JJV
588                        JJEND = JJV + NQC - 1
589                        JJNXT = MIN( ICEIL( JJBEG , NBV ) * NBV , JJEND )
590  
591     70 CONTINUE
592        IF(( K - ILEFT ).GT.0 ) THEN
593            CALL DLASET( 'Lower' , K - ILEFT , JJNXT - JJBEG + 1 , ZERO ,
594       $    ONE , WORK( IPV + ILEFT + (JJBEG - JJV)*LV ) ,
595       $    LV )
596            MYDIST = MYDIST + NPCOL
597            ILEFT = MYDIST * NBV - ICOFFV
598            JJBEG = JJNXT + 1
599            JJNXT = MIN( JJNXT + NBV , JJEND )
600            GO TO 70
601        END IF
602  
603        ELSE
604  
605  *         WORK( IPW ) =( . V1 V2 ) where V2 is unit lower
606  *         triangular , zeroes upper triangular part of V2.
607  
608            II = IIV
609            CALL INFOG1L( JV + N - K , NBV , NPCOL , MYCOL , DESCV( CSRC_ ) ,
610       $    JJ , ILASTCOL )
611            IOFF = MOD( JV + N - K - 1 , NBV )
612            KQ = NUMROC( K + IOFF , NBV , MYCOL , ILASTCOL , NPCOL )
613            IF( MYCOL.EQ.ILASTCOL )
614       $        KQ = KQ - IOFF
615                MYDIST = MOD( MYCOL - ILASTCOL + NPCOL , NPCOL )
616                ILEFT = MYDIST * NBV - IOFF
617                IRIGHT = MIN( ILEFT + NBV , K )
618                ILEFT = MIN( MAX( 0 , ILEFT ) , K )
619  
620     80 CONTINUE
621        IF( II.LE.( IIV + K - 1 ) ) THEN
622            WIDE = IRIGHT - ILEFT
623            CALL DLASET( 'All' , ILEFT - II + IIV , KQ , ZERO , ZERO ,
624       $    WORK( IPV + II - IIV + (JJ - JJV)*LV ) , LV )
625            CALL DLASET( 'Upper' , WIDE , KQ , ZERO , ONE ,
626       $    WORK( IPV + ILEFT + (JJ - JJV)*LV ) , LV )
627            KQ = MAX( 0 , KQ - WIDE )
628            II = IIV + IRIGHT
629            JJ = JJ + WIDE
630            MYDIST = MYDIST + NPCOL
631            ILEFT = MYDIST * NBV - IOFF
632            IRIGHT = MIN( ILEFT + NBV , K )
633            ILEFT = MIN( ILEFT , K )
634            GO TO 80
635        END IF
636  
637        END IF
638  
639  *     WORK( IPV ) is K x NQC = V = V( IOFFV )
640  *     WORK( IPW ) = C( IOFFC ) * V'(MPC x NQC x K) -> MPC x K
641  
642        IF( NQC.GT.0 ) THEN
643            CALL DGEMM( 'No Transpose' , 'Transpose' , MPC , K , NQC ,
644       $    ONE , C( IOFFC ) , LDC , WORK( IPV ) , LV , ZERO ,
645       $    WORK( IPW ) , LW )
646        ELSE
647            CALL DLASET( 'All' , MPC , K , ZERO , ZERO , WORK( IPW ) , LW )
648        END IF
649  
650        CALL DGSUM2D( ICTXT , 'Rowwise' , ' ' , MPC , K , WORK( IPW ) ,
651       $LW , MYROW , IVCOL )
652  
653  *     WORK( IPW ) = WORK( IPW ) * T' or WORK( IPW ) * T
654  
655        IF( MYCOL.EQ.IVCOL ) THEN
656            CALL DTRMM( 'Right' , UPLO , TRANS , 'Non unit' , MPC , K ,
657       $    ONE , T , MBV , WORK( IPW ) , LW )
658            CALL DGEBS2D( ICTXT , 'Rowwise' , ' ' , MPC , K , WORK( IPW ) ,
659       $    LW )
660        ELSE
661            CALL DGEBR2D( ICTXT , 'Rowwise' , ' ' , MPC , K , WORK( IPW ) ,
662       $    LW , MYROW , IVCOL )
663        END IF
664  
665  *     C C - W * V
666  *     C( IOFFC ) = C( IOFFC ) - WORK( IPW ) * WORK( IPV )
667  *     MPC x NQC MPC x K K x NQC
668  
669        CALL DGEMM( 'No transpose' , 'No transpose' , MPC , NQC , K ,
670       $- ONE , WORK( IPW ) , LW , WORK( IPV ) , LV , ONE ,
671       $C( IOFFC ) , LDC )
672  
673        END IF
674  
675        END IF
676  
677        RETURN
678  
679  *     End of PDLARFB
680  
681        END