ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
pzgebdrv.f
Go to the documentation of this file.
00001       SUBROUTINE PZGEBDRV( M, N, A, IA, JA, DESCA, D, E, TAUQ, TAUP,
00002      $                     WORK, INFO )
00003 *
00004 *  -- ScaLAPACK routine (version 1.7) --
00005 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00006 *     and University of California, Berkeley.
00007 *     May 1, 1997
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER              INFO, IA, JA, M, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER              DESCA( * )
00014       DOUBLE PRECISION   D( * ), E( * )
00015       COMPLEX*16         A( * ), TAUP( * ), TAUQ( * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  PZGEBDRV computes sub( A ) = A(IA:IA+M-1,JA:JA+N-1) from sub( A ),
00022 *  Q, P returned by PZGEBRD:
00023 *
00024 *                         sub( A ) := Q * B * P'.
00025 *
00026 *  Notes
00027 *  =====
00028 *
00029 *  Each global data object is described by an associated description
00030 *  vector.  This vector stores the information required to establish
00031 *  the mapping between an object element and its corresponding process
00032 *  and memory location.
00033 *
00034 *  Let A be a generic term for any 2D block cyclicly distributed array.
00035 *  Such a global array has an associated description vector DESCA.
00036 *  In the following comments, the character _ should be read as
00037 *  "of the global array".
00038 *
00039 *  NOTATION        STORED IN      EXPLANATION
00040 *  --------------- -------------- --------------------------------------
00041 *  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
00042 *                                 DTYPE_A = 1.
00043 *  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
00044 *                                 the BLACS process grid A is distribu-
00045 *                                 ted over. The context itself is glo-
00046 *                                 bal, but the handle (the integer
00047 *                                 value) may vary.
00048 *  M_A    (global) DESCA( M_ )    The number of rows in the global
00049 *                                 array A.
00050 *  N_A    (global) DESCA( N_ )    The number of columns in the global
00051 *                                 array A.
00052 *  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
00053 *                                 the rows of the array.
00054 *  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
00055 *                                 the columns of the array.
00056 *  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
00057 *                                 row of the array A is distributed.
00058 *  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
00059 *                                 first column of the array A is
00060 *                                 distributed.
00061 *  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
00062 *                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
00063 *
00064 *  Let K be the number of rows or columns of a distributed matrix,
00065 *  and assume that its process grid has dimension p x q.
00066 *  LOCr( K ) denotes the number of elements of K that a process
00067 *  would receive if K were distributed over the p processes of its
00068 *  process column.
00069 *  Similarly, LOCc( K ) denotes the number of elements of K that a
00070 *  process would receive if K were distributed over the q processes of
00071 *  its process row.
00072 *  The values of LOCr() and LOCc() may be determined via a call to the
00073 *  ScaLAPACK tool function, NUMROC:
00074 *          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
00075 *          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
00076 *  An upper bound for these quantities may be computed by:
00077 *          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
00078 *          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
00079 *
00080 *  Arguments
00081 *  =========
00082 *
00083 *  M       (global input) INTEGER
00084 *          The number of rows to be operated on, i.e. the number of rows
00085 *          of the distributed submatrix sub( A ). M >= 0.
00086 *
00087 *  N       (global input) INTEGER
00088 *          The number of columns to be operated on, i.e. the number of
00089 *          columns of the distributed submatrix sub( A ). N >= 0.
00090 *
00091 *  A       (local input/local output) COMPLEX*16 pointer into the
00092 *          local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
00093 *          On entry, this array contains the local pieces of sub( A )
00094 *          as returned by PZGEBRD. On exit, the original distribu-
00095 *          ted matrix sub( A ) is restored.
00096 *
00097 *  IA      (global input) INTEGER
00098 *          The row index in the global array A indicating the first
00099 *          row of sub( A ).
00100 *
00101 *  JA      (global input) INTEGER
00102 *          The column index in the global array A indicating the
00103 *          first column of sub( A ).
00104 *
00105 *  DESCA   (global and local input) INTEGER array of dimension DLEN_.
00106 *          The array descriptor for the distributed matrix A.
00107 *
00108 *  D       (local input) DOUBLE PRECISION array, dimension
00109 *          LOCc(JA+MIN(M,N)-1) if M >= N; LOCr(IA+MIN(M,N)-1) otherwise.
00110 *          The distributed diagonal elements of the bidiagonal matrix
00111 *          B: D(i) = A(i,i). D is tied to the distributed matrix A.
00112 *
00113 *  E       (local input) DOUBLE PRECISION array, dimension
00114 *          LOCr(IA+MIN(M,N)-1) if M >= N; LOCc(JA+MIN(M,N)-2) otherwise.
00115 *          The distributed off-diagonal elements of the bidiagonal
00116 *          distributed matrix B:
00117 *          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
00118 *          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
00119 *          E is tied to the distributed matrix A.
00120 *
00121 *  TAUQ    (local input) COMPLEX*16 array dimension
00122 *          LOCc(JA+MIN(M,N)-1). The scalar factors of the elementary
00123 *          reflectors which represent the unitary matrix Q. TAUQ is
00124 *          tied to the distributed matrix A. See Further Details.
00125 *
00126 *  TAUP    (local input) COMPLEX*16 array, dimension
00127 *          LOCr(IA+MIN(M,N)-1). The scalar factors of the elementary
00128 *          reflectors which represent the unitary matrix P. TAUP is
00129 *          tied to the distributed matrix A. See Further Details.
00130 *
00131 *  WORK    (local workspace) COMPLEX*16 array, dimension (LWORK)
00132 *          LWORK >= 2*NB*( MP + NQ + NB )
00133 *
00134 *          where NB = MB_A = NB_A,
00135 *          IROFFA = MOD( IA-1, NB ), ICOFFA = MOD( JA-1, NB ),
00136 *          IAROW = INDXG2P( IA, NB, MYROW, RSRC_A, NPROW ),
00137 *          IACOL = INDXG2P( JA, NB, MYCOL, CSRC_A, NPCOL ),
00138 *          MP = NUMROC( M+IROFFA, NB, MYROW, IAROW, NPROW ),
00139 *          NQ = NUMROC( N+ICOFFA, NB, MYCOL, IACOL, NPCOL ).
00140 *
00141 *          INDXG2P and NUMROC are ScaLAPACK tool functions;
00142 *          MYROW, MYCOL, NPROW and NPCOL can be determined by calling
00143 *          the subroutine BLACS_GRIDINFO.
00144 *
00145 *  INFO    (global output) INTEGER
00146 *          On exit, if INFO <> 0, a discrepancy has been found between
00147 *          the diagonal and off-diagonal elements of A and the copies
00148 *          contained in the arrays D and E.
00149 *
00150 *  =====================================================================
00151 *
00152 *     .. Parameters ..
00153       INTEGER            BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
00154      $                   LLD_, MB_, M_, NB_, N_, RSRC_
00155       PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
00156      $                     CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
00157      $                     RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
00158       DOUBLE PRECISION   REIGHT, RZERO
00159       PARAMETER          ( REIGHT = 8.0D+0, RZERO = 0.0D+0 )
00160       COMPLEX*16         ONE, ZERO
00161       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
00162      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
00163 *     ..
00164 *     .. Local Scalars ..
00165       INTEGER            I, IACOL, IAROW, ICTXT, IIA, IL, IPTP, IPTQ,
00166      $                   IPV, IPW, IPWK, IOFF, IV, J, JB, JJA, JL, JV,
00167      $                   K, MN, MP, MYCOL, MYROW, NB, NPCOL, NPROW, NQ
00168       DOUBLE PRECISION   ADDBND, D2, E2
00169       COMPLEX*16         D1, E1
00170 *     ..
00171 *     .. Local Arrays ..
00172       INTEGER            DESCD( DLEN_ ), DESCE( DLEN_ ), DESCV( DLEN_ ),
00173      $                   DESCW( DLEN_ )
00174 *     ..
00175 *     .. External Subroutines ..
00176       EXTERNAL           BLACS_GRIDINFO, DESCSET, IGSUM2D, INFOG2L,
00177      $                   PDELGET, PZLACPY, PZLARFB, PZLARFT,
00178      $                   PZLASET, PZELGET
00179 *     ..
00180 *     .. External Functions ..
00181       INTEGER            INDXG2P, NUMROC
00182       DOUBLE PRECISION   PDLAMCH
00183       EXTERNAL           INDXG2P, NUMROC, PDLAMCH
00184 *     ..
00185 *     .. Intrinsic Functions ..
00186       INTRINSIC          ABS, DCMPLX, MAX, MIN, MOD
00187 *     ..
00188 *     .. Executable Statements ..
00189 *
00190 *     Get grid parameters
00191 *
00192       ICTXT = DESCA( CTXT_ )
00193       CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
00194 *
00195       INFO = 0
00196       NB = DESCA( MB_ )
00197       IOFF = MOD( IA-1, NB )
00198       CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL, IIA, JJA,
00199      $              IAROW, IACOL )
00200       MP = NUMROC( M+IOFF, NB, MYROW, IAROW, NPROW )
00201       NQ = NUMROC( N+IOFF, NB, MYCOL, IACOL, NPCOL )
00202       IPV  = 1
00203       IPW  = IPV + MP*NB
00204       IPTP = IPW + NQ*NB
00205       IPTQ = IPTP + NB*NB
00206       IPWK = IPTQ + NB*NB
00207 *
00208       IV = 1
00209       JV = 1
00210       MN = MIN( M, N )
00211       IL = MAX( ( (IA+MN-2) / NB )*NB + 1, IA )
00212       JL = MAX( ( (JA+MN-2) / NB )*NB + 1, JA )
00213       IAROW = INDXG2P( IL, NB, MYROW, DESCA( RSRC_ ), NPROW )
00214       IACOL = INDXG2P( JL, NB, MYCOL, DESCA( CSRC_ ), NPCOL )
00215       CALL DESCSET( DESCV, IA+M-IL, NB, NB, NB, IAROW, IACOL, ICTXT,
00216      $              MAX( 1, MP ) )
00217       CALL DESCSET( DESCW, NB, JA+N-JL, NB, NB, IAROW, IACOL, ICTXT,
00218      $              NB )
00219 *
00220       ADDBND = REIGHT * PDLAMCH( ICTXT, 'eps' )
00221 *
00222 *     When A is an upper bidiagonal form
00223 *
00224       IF( M.GE.N ) THEN
00225 *
00226          CALL DESCSET( DESCD, 1, JA+MN-1, 1, DESCA( NB_ ), MYROW,
00227      $                 DESCA( CSRC_ ), DESCA( CTXT_ ), 1 )
00228          CALL DESCSET( DESCE, IA+MN-1, 1, DESCA( MB_ ), 1,
00229      $                 DESCA( RSRC_ ), MYCOL, DESCA( CTXT_ ),
00230      $                 DESCA( LLD_ ) )
00231 *
00232          DO 10 J = 0, MN-1
00233             D1 = ZERO
00234             E1 = ZERO
00235             D2 = RZERO
00236             E2 = RZERO
00237             CALL PDELGET( ' ', ' ', D2, D, 1, JA+J, DESCD )
00238             CALL PZELGET( 'Columnwise', ' ', D1, A, IA+J, JA+J, DESCA )
00239             IF( J.LT.(MN-1) ) THEN
00240                CALL PDELGET( ' ', ' ', E2, E, IA+J, 1, DESCE )
00241                CALL PZELGET( 'Rowwise', ' ', E1, A, IA+J, JA+J+1,
00242      $                       DESCA )
00243             END IF
00244 *
00245             IF( ( ABS( D1-DCMPLX( D2 ) ).GT.( ABS( D2 )*ADDBND ) ) .OR.
00246      $          ( ABS( E1-DCMPLX( E2 ) ).GT.( ABS( E2 )*ADDBND ) ) )
00247      $         INFO = INFO + 1
00248    10    CONTINUE
00249 *
00250          DO 20 J = JL, JA+NB-IOFF, -NB
00251             JB = MIN( JA+N-J, NB )
00252             I  = IA + J - JA
00253             K  = I - IA + 1
00254 *
00255 *           Compute upper triangular matrix TQ from TAUQ.
00256 *
00257             CALL PZLARFT( 'Forward', 'Columnwise', M-K+1, JB, A, I, J,
00258      $                    DESCA, TAUQ, WORK( IPTQ ), WORK( IPWK ) )
00259 *
00260 *           Copy Householder vectors into workspace.
00261 *
00262             CALL PZLACPY( 'Lower', M-K+1, JB, A, I, J, DESCA,
00263      $                    WORK( IPV ), IV, JV, DESCV )
00264             CALL PZLASET( 'Upper', M-K+1, JB, ZERO, ONE, WORK( IPV ),
00265      $                    IV, JV, DESCV )
00266 *
00267 *           Zero out the strict lower triangular part of A.
00268 *
00269             CALL PZLASET( 'Lower', M-K, JB, ZERO, ZERO, A, I+1, J,
00270      $                    DESCA )
00271 *
00272 *           Compute upper triangular matrix TP from TAUP.
00273 *
00274             CALL PZLARFT( 'Forward', 'Rowwise', N-K, JB, A, I, J+1,
00275      $                    DESCA, TAUP, WORK( IPTP ), WORK( IPWK ) )
00276 *
00277 *           Copy Householder vectors into workspace.
00278 *
00279             CALL PZLACPY( 'Upper', JB, N-K, A, I, J+1, DESCA,
00280      $                    WORK( IPW ), IV, JV+1, DESCW )
00281             CALL PZLASET( 'Lower', JB, N-K, ZERO, ONE, WORK( IPW ), IV,
00282      $                    JV+1, DESCW )
00283 *
00284 *           Zero out the strict+1 upper triangular part of A.
00285 *
00286             CALL PZLASET( 'Upper', JB, N-K-1, ZERO, ZERO, A, I, J+2,
00287      $                    DESCA )
00288 *
00289 *           Apply block Householder transformation from Left.
00290 *
00291             CALL PZLARFB( 'Left', 'No transpose', 'Forward',
00292      $                    'Columnwise', M-K+1, N-K+1, JB, WORK( IPV ),
00293      $                    IV, JV, DESCV, WORK( IPTQ ), A, I, J, DESCA,
00294      $                    WORK( IPWK ) )
00295 *
00296 *           Apply block Householder transformation from Right.
00297 *
00298             CALL PZLARFB( 'Right', 'Conjugate transpose', 'Forward',
00299      $                    'Rowwise', M-K+1, N-K, JB, WORK( IPW ), IV,
00300      $                    JV+1, DESCW, WORK( IPTP ), A, I, J+1, DESCA,
00301      $                    WORK( IPWK ) )
00302 *
00303             DESCV( M_ ) = DESCV( M_ ) + NB
00304             DESCV( RSRC_ ) = MOD( DESCV( RSRC_ ) + NPROW - 1, NPROW )
00305             DESCV( CSRC_ ) = MOD( DESCV( CSRC_ ) + NPCOL - 1, NPCOL )
00306             DESCW( N_ ) = DESCW( N_ ) + NB
00307             DESCW( RSRC_ ) = DESCV( RSRC_ )
00308             DESCW( CSRC_ ) = DESCV( CSRC_ )
00309 *
00310    20    CONTINUE
00311 *
00312 *        Handle first block separately
00313 *
00314          JB = MIN( N, NB - IOFF )
00315          IV = IOFF + 1
00316          JV = IOFF + 1
00317 *
00318 *        Compute upper triangular matrix TQ from TAUQ.
00319 *
00320          CALL PZLARFT( 'Forward', 'Columnwise', M, JB, A, IA, JA, DESCA,
00321      $                 TAUQ, WORK( IPTQ ), WORK( IPWK ) )
00322 *
00323 *        Copy Householder vectors into workspace.
00324 *
00325          CALL PZLACPY( 'Lower', M, JB, A, IA, JA, DESCA, WORK( IPV ),
00326      $                 IV, JV, DESCV )
00327          CALL PZLASET( 'Upper', M, JB, ZERO, ONE, WORK( IPV ), IV, JV,
00328      $                 DESCV )
00329 *
00330 *        Zero out the strict lower triangular part of A.
00331 *
00332          CALL PZLASET( 'Lower', M-1, JB, ZERO, ZERO, A, IA+1, JA,
00333      $                 DESCA )
00334 *
00335 *        Compute upper triangular matrix TP from TAUP.
00336 *
00337          CALL PZLARFT( 'Forward', 'Rowwise', N-1, JB, A, IA, JA+1,
00338      $                 DESCA, TAUP, WORK( IPTP ), WORK( IPWK ) )
00339 *
00340 *        Copy Householder vectors into workspace.
00341 *
00342          CALL PZLACPY( 'Upper', JB, N-1, A, IA, JA+1, DESCA,
00343      $                 WORK( IPW ), IV, JV+1, DESCW )
00344          CALL PZLASET( 'Lower', JB, N-1, ZERO, ONE, WORK( IPW ), IV,
00345      $                 JV+1, DESCW )
00346 *
00347 *        Zero out the strict+1 upper triangular part of A.
00348 *
00349          CALL PZLASET( 'Upper', JB, N-2, ZERO, ZERO, A, IA, JA+2,
00350      $                 DESCA )
00351 *
00352 *        Apply block Householder transformation from left.
00353 *
00354          CALL PZLARFB( 'Left', 'No transpose', 'Forward', 'Columnwise',
00355      $                 M, N, JB, WORK( IPV ), IV, JV, DESCV,
00356      $                 WORK( IPTQ ), A, IA, JA, DESCA, WORK( IPWK ) )
00357 *
00358 *        Apply block Householder transformation from right.
00359 *
00360          CALL PZLARFB( 'Right', 'Conjugate transpose', 'Forward',
00361      $                 'Rowwise', M, N-1, JB, WORK( IPW ), IV, JV+1,
00362      $                 DESCW, WORK( IPTP ), A, IA, JA+1, DESCA,
00363      $                 WORK( IPWK ) )
00364 *
00365       ELSE
00366 *
00367          CALL DESCSET( DESCD, IA+MN-1, 1, DESCA( MB_ ), 1,
00368      $                 DESCA( RSRC_ ), MYCOL, DESCA( CTXT_ ),
00369      $                 DESCA( LLD_ ) )
00370          CALL DESCSET( DESCE, 1, JA+MN-2, 1, DESCA( NB_ ), MYROW,
00371      $                 DESCA( CSRC_ ), DESCA( CTXT_ ), 1 )
00372 *
00373          DO 30 J = 0, MN-1
00374             D1 = ZERO
00375             E1 = ZERO
00376             D2 = RZERO
00377             E2 = RZERO
00378             CALL PDELGET( ' ', ' ', D2, D, IA+J, 1, DESCD )
00379             CALL PZELGET( 'Rowwise', ' ', D1, A, IA+J, JA+J, DESCA )
00380             IF( J.LT.(MN-1) ) THEN
00381                CALL PDELGET( ' ', ' ', E2, E, 1, JA+J, DESCE )
00382                CALL PZELGET( 'Columnwise', ' ', E1, A, IA+J+1, JA+J,
00383      $                       DESCA )
00384             END IF
00385 *
00386             IF( ( ABS( D1-DCMPLX( D2 ) ).GT.( ABS( D2 )*ADDBND ) ) .OR.
00387      $          ( ABS( E1-DCMPLX( E2 ) ).GT.( ABS( E2 )*ADDBND ) ) )
00388      $         INFO = INFO + 1
00389    30    CONTINUE
00390 *
00391          DO 40 I = IL, IA+NB-IOFF, -NB
00392             JB = MIN( IA+M-I, NB )
00393             J  = JA + I - IA
00394             K  = J - JA + 1
00395 *
00396 *           Compute upper triangular matrix TQ from TAUQ.
00397 *
00398             CALL PZLARFT( 'Forward', 'Columnwise', M-K, JB, A, I+1, J,
00399      $                    DESCA, TAUQ, WORK( IPTQ ), WORK( IPWK ) )
00400 *
00401 *           Copy Householder vectors into workspace.
00402 *
00403             CALL PZLACPY( 'Lower', M-K, JB, A, I+1, J, DESCA,
00404      $                    WORK( IPV ), IV+1, JV, DESCV )
00405             CALL PZLASET( 'Upper', M-K, JB, ZERO, ONE, WORK( IPV ),
00406      $                    IV+1, JV, DESCV )
00407 *
00408 *           Zero out the strict lower triangular part of A.
00409 *
00410             CALL PZLASET( 'Lower', M-K-1, JB, ZERO, ZERO, A, I+2, J,
00411      $                    DESCA )
00412 *
00413 *           Compute upper triangular matrix TP from TAUP.
00414 *
00415             CALL PZLARFT( 'Forward', 'Rowwise', N-K+1, JB, A, I, J,
00416      $                    DESCA, TAUP, WORK( IPTP ), WORK( IPWK ) )
00417 *
00418 *           Copy Householder vectors into workspace.
00419 *
00420             CALL PZLACPY( 'Upper', JB, N-K+1, A, I, J, DESCA,
00421      $                    WORK( IPW ), IV, JV, DESCW )
00422             CALL PZLASET( 'Lower', JB, N-K+1, ZERO, ONE, WORK( IPW ),
00423      $                    IV, JV, DESCW )
00424 *
00425 *           Zero out the strict+1 upper triangular part of A.
00426 *
00427             CALL PZLASET( 'Upper', JB, N-K, ZERO, ZERO, A, I, J+1,
00428      $                    DESCA )
00429 *
00430 *           Apply block Householder transformation from Left.
00431 *
00432             CALL PZLARFB( 'Left', 'No transpose', 'Forward',
00433      $                    'Columnwise', M-K, N-K+1, JB, WORK( IPV ),
00434      $                    IV+1, JV, DESCV, WORK( IPTQ ), A, I+1, J,
00435      $                    DESCA, WORK( IPWK ) )
00436 *
00437 *           Apply block Householder transformation from Right.
00438 *
00439             CALL PZLARFB( 'Right', 'Conjugate transpose', 'Forward',
00440      $                    'Rowwise', M-K+1, N-K+1, JB, WORK( IPW ), IV,
00441      $                    JV, DESCW, WORK( IPTP ), A, I, J, DESCA,
00442      $                    WORK( IPWK ) )
00443 *
00444             DESCV( M_ ) = DESCV( M_ ) + NB
00445             DESCV( RSRC_ ) = MOD( DESCV( RSRC_ ) + NPROW - 1, NPROW )
00446             DESCV( CSRC_ ) = MOD( DESCV( CSRC_ ) + NPCOL - 1, NPCOL )
00447             DESCW( N_ ) = DESCW( N_ ) + NB
00448             DESCW( RSRC_ ) = DESCV( RSRC_ )
00449             DESCW( CSRC_ ) = DESCV( CSRC_ )
00450 *
00451    40    CONTINUE
00452 *
00453 *        Handle first block separately
00454 *
00455          JB = MIN( M, NB - IOFF )
00456          IV = IOFF + 1
00457          JV = IOFF + 1
00458 *
00459 *        Compute upper triangular matrix TQ from TAUQ.
00460 *
00461          CALL PZLARFT( 'Forward', 'Columnwise', M-1, JB, A, IA+1, JA,
00462      $                 DESCA, TAUQ, WORK( IPTQ ), WORK( IPWK ) )
00463 *
00464 *        Copy Householder vectors into workspace.
00465 *
00466          CALL PZLACPY( 'Lower', M-1, JB, A, IA+1, JA, DESCA,
00467      $                 WORK( IPV ), IV+1, JV, DESCV )
00468          CALL PZLASET( 'Upper', M-1, JB, ZERO, ONE, WORK( IPV ), IV+1,
00469      $                 JV, DESCV )
00470 *
00471 *        Zero out the strict lower triangular part of A.
00472 *
00473          CALL PZLASET( 'Lower', M-2, JB, ZERO, ZERO, A, IA+2, JA,
00474      $                 DESCA )
00475 *
00476 *        Compute upper triangular matrix TP from TAUP.
00477 *
00478          CALL PZLARFT( 'Forward', 'Rowwise', N, JB, A, IA, JA, DESCA,
00479      $                 TAUP, WORK( IPTP ), WORK( IPWK ) )
00480 *
00481 *        Copy Householder vectors into workspace.
00482 *
00483          CALL PZLACPY( 'Upper', JB, N, A, IA, JA, DESCA, WORK( IPW ),
00484      $                 IV, JV, DESCW )
00485          CALL PZLASET( 'Lower', JB, N, ZERO, ONE, WORK( IPW ), IV, JV,
00486      $                 DESCW )
00487 *
00488 *        Zero out the strict+1 upper triangular part of A.
00489 *
00490          CALL PZLASET( 'Upper', JB, N-1, ZERO, ZERO, A, IA, JA+1,
00491      $                 DESCA )
00492 *
00493 *        Apply block Householder transformation from left
00494 *
00495          CALL PZLARFB( 'Left', 'No transpose', 'Forward', 'Columnwise',
00496      $                 M-1, N, JB, WORK( IPV ), IV+1, JV, DESCV,
00497      $                 WORK( IPTQ ), A, IA+1, JA, DESCA, WORK( IPWK ) )
00498 *
00499 *        Apply block Householder transformation from right
00500 *
00501          CALL PZLARFB( 'Right', 'Conjugate transpose', 'Forward',
00502      $                 'Rowwise', M, N, JB, WORK( IPW ), IV, JV, DESCW,
00503      $                 WORK( IPTP ), A, IA, JA, DESCA, WORK( IPWK ) )
00504       END IF
00505 *
00506       CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, INFO, 1, -1, 0 )
00507 *
00508       RETURN
00509 *
00510 *     End of PZGEBDRV
00511 *
00512       END