Routine: PZLAHRD()  File: SRC\pzlahrd.f

 
 
# lines: 290
  # code: 290
  # comment: 0
  # blank:0
# Variables:42
# Callers:1
# Callings:2
# Words:111
# Keywords:58
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PZLAHRD reduces the first NB columns of a complex general
  N-by-(N-K+1) distributed matrix A(IA:IA+N-1,JA:JA+N-K) so that
  elements below the k-th subdiagonal are zero. The reduction is
  performed by an unitary similarity transformation Q' * A * Q. The
  routine returns the matrices V and T which determine Q as a block
  reflector I - V*T*V', and also the matrix Y = A * V * T.
  This is an auxiliary routine called by PZGEHRD. In the following
  comments sub( A ) denotes A(IA:IA+N-1,JA:JA+N-1).
  Arguments
  =========
  N       (global input) INTEGER
          The number of rows and columns to be operated on, i.e. the
          order of the distributed submatrix sub( A ).
          N >= 0.
  K       (global input) INTEGER
          The offset for the reduction. Elements below the k-th
          subdiagonal in the first NB columns are reduced to zero.
  NB      (global input) INTEGER
          The number of columns to be reduced.
  A       (local input/local output) COMPLEX*16 pointer into
          the local memory to an array of dimension (LLD_A,
          LOCc(JA+N-K)). On entry, this array contains the the local
          pieces of the N-by-(N-K+1) general distributed matrix
          A(IA:IA+N-1,JA:JA+N-K). On exit, the elements on and above
          the k-th subdiagonal in the first NB columns are overwritten
          with the corresponding elements of the reduced distributed
          matrix; the elements below the k-th subdiagonal, with the
          array TAU, represent the matrix Q as a product of elementary
          reflectors. The other columns of A(IA:IA+N-1,JA:JA+N-K) are
          unchanged. See Further Details.
  IA      (global input) INTEGER
          The row index in the global array A indicating the first
          row of sub( A ).
  JA      (global input) INTEGER
          The column index in the global array A indicating the
          first column of sub( A ).
  DESCA   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix A.
  TAU     (local output) COMPLEX*16 array, dimension LOCc(JA+N-2)
          The scalar factors of the elementary reflectors (see Further
          Details). TAU is tied to the distributed matrix A.
  T       (local output) COMPLEX*16 array, dimension (NB_A,NB_A)
          The upper triangular matrix T.
  Y       (local output) COMPLEX*16 pointer into the local memory
          to an array of dimension (LLD_Y,NB_A). On exit, this array
          contains the local pieces of the N-by-NB distributed
          matrix Y. LLD_Y >= LOCr(IA+N-1).
  IY      (global input) INTEGER
          The row index in the global array Y indicating the first
          row of sub( Y ).
  JY      (global input) INTEGER
          The column index in the global array Y indicating the
          first column of sub( Y ).
  DESCY   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix Y.
  WORK    (local workspace) COMPLEX*16 array, dimension (NB)
  Further Details
  ===============
  The matrix Q is represented as a product of nb elementary reflectors
     Q = H(1) H(2) . . . H(nb).
  Each H(i) has the form
     H(i) = I - tau * v * v'
  where tau is a complex scalar, and v is a complex vector with
  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
  A(ia+i+k:ia+n-1,ja+i-1), and tau in TAU(ja+i-1).
  The elements of the vectors v together form the (n-k+1)-by-nb matrix
  V which is needed, with T and Y, to apply the transformation to the
  unreduced part of the matrix, using an update of the form:
  A(ia:ia+n-1,ja:ja+n-k) := (I-V*T*V')*(A(ia:ia+n-1,ja:ja+n-k)-Y*V').
  The contents of A(ia:ia+n-1,ja:ja+n-k) on exit are illustrated by the
  following example with n = 7, k = 3 and nb = 2:
     ( a   h   a   a   a )
     ( a   h   a   a   a )
     ( a   h   a   a   a )
     ( h   h   a   a   a )
     ( v1  h   a   a   a )
     ( v1  v2  a   a   a )
     ( v1  v2  a   a   a )
  where a denotes an element of the original matrix
  A(ia:ia+n-1,ja:ja+n-k), h denotes a modified element of the upper
  Hessenberg matrix H, and vi denotes an element of the vector
  defining H(i).
  =====================================================================
     .. Parameters ..

 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

 
001        SUBROUTINE PZLAHRD( N , K , NB , A , IA , JA , DESCA , TAU , T , Y , IY , JY ,
002       $DESCY , 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        INTEGER IA , IY , JA , JY , K , N , NB
011        INTEGER BLOCK_CYCLIC_2D , CSRC_ , CTXT_ , DLEN_ , DTYPE_ ,
012       $LLD_ , MB_ , M_ , NB_ , N_ , RSRC_
013        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
014       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
015       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
016        COMPLEX*16 ONE , ZERO
017        PARAMETER( ONE =( 1.0D + 0 , 0.0D + 0 ) ,
018       $ZERO =( 0.0D + 0 , 0.0D + 0 ) )
019  *     ..
020  *     .. Local Scalars ..
021        LOGICAL IPROC
022        INTEGER I , IACOL , IAROW , ICTXT , IOFF , II , J , JJ , JL ,
023       $JT , JW , L , MYROW , MYCOL , NPCOL , NPROW , NQ
024        COMPLEX*16 EI
025  *     ..
026  *     .. Local Arrays ..
027        INTEGER DESCW( DLEN_ )
028  *     ..
029  *     .. External Functions ..
030        INTEGER NUMROC
031        EXTERNAL NUMROC
032  *     ..
033  *     .. External Subroutines ..
034        EXTERNAL BLACS_GRIDINFO , DESCSET , INFOG2L , PZELSET ,
035       $PZGEMV , PZLACGV , PZLARFG , PZSCAL ,
036       $ZAXPY , ZCOPY , ZSCAL , ZTRMV
037  *     ..
038  *     .. Intrinsic Functions ..
039        INTRINSIC MIN , MOD
040  *     ..
041  *     .. Executable Statements ..
042  
043  *     Quick return if possible
044  
045        IF( N.LE.1 )
046       $    RETURN
047  
048            ICTXT = DESCA( CTXT_ )
049            CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
050  
051            IOFF = MOD( JA - 1 , DESCA( NB_ ) )
052            CALL INFOG2L( IA + K , JA , DESCA , NPROW , NPCOL , MYROW , MYCOL , II ,
053       $    JJ , IAROW , IACOL )
054  
055            IPROC =( MYROW.EQ.IAROW .AND. MYCOL.EQ.IACOL )
056            NQ = NUMROC( N + JA - 1 , DESCA( NB_ ) , MYCOL , IACOL , NPCOL )
057            IF( MYCOL.EQ.IACOL )
058       $        NQ = NQ - IOFF
059  
060                EI = ZERO
061  
062                JW = IOFF + 1
063                CALL DESCSET( DESCW , 1 , DESCA( MB_ ) , 1 , DESCA( MB_ ) , IAROW ,
064       $        IACOL , ICTXT , 1 )
065  
066                DO 10 L = 1 , NB
067                    I = IA + K + L - 2
068                    J = JA + L - 1
069  
070                    IF( L.GT.1 ) THEN
071  
072  *                     Update A(ia : ia + n - 1 , j)
073  
074  *                     Compute i - th column of A - Y * V'
075  
076                        CALL PZLACGV ( L - 1 , A , I , JA , DESCA , DESCA( M_ ) )
077                        CALL PZGEMV( 'No transpose' , N , L - 1 , - ONE , Y , IY , JY , DESCY ,
078       $                A , I , JA , DESCA , DESCA( M_ ) , ONE , A , IA , J ,
079       $                DESCA , 1 )
080                        CALL PZLACGV ( L - 1 , A , I , JA , DESCA , DESCA( M_ ) )
081  
082  *                     Apply I - V * T' * V' to this column(call it b) from the
083  *                     left , using the last column of T as workspace
084  
085  *                     Let V =( V1 ) and b =( b1 )(first I - 1 rows)
086  *                     ( V2 )( b2 )
087  
088  *                     where V1 is unit lower triangular
089  
090  *                     w := V1' * b1
091  
092                        IF( IPROC ) THEN
093                            CALL ZCOPY( L - 1 , A((JJ + L - 2)*DESCA( LLD_ ) + II ) , 1 ,
094       $                    WORK( JW ) , 1 )
095                            CALL ZTRMV( 'Lower' , 'Conjugate transpose' , 'Unit' , L - 1 ,
096       $                    A((JJ - 1)*DESCA( LLD_ ) + II ) , DESCA( LLD_ ) ,
097       $                    WORK( JW ) , 1 )
098                        END IF
099  
100  *                     w := w + V2'*b2
101  
102                        CALL PZGEMV( 'Conjugate transpose' , N - K - L + 1 , L - 1 , ONE , A ,
103       $                I + 1 , JA , DESCA , A , I + 1 , J , DESCA , 1 , ONE , WORK ,
104       $                1 , JW , DESCW , DESCW( M_ ) )
105  
106  *                     w := T'*w
107  
108                        IF( IPROC )
109       $                    CALL ZTRMV( 'Upper' , 'Conjugate transpose' , 'Non - unit' ,
110       $                    L - 1 , T , DESCA( NB_ ) , WORK( JW ) , 1 )
111  
112  *                         b2 := b2 - V2*w
113  
114                            CALL PZGEMV( 'No transpose' , N - K - L + 1 , L - 1 , - ONE , A , I + 1 , JA ,
115       $                    DESCA , WORK , 1 , JW , DESCW , DESCW( M_ ) , ONE ,
116       $                    A , I + 1 , J , DESCA , 1 )
117  
118  *                         b1 := b1 - V1*w
119  
120                            IF( IPROC ) THEN
121                                CALL ZTRMV( 'Lower' , 'No transpose' , 'Unit' , L - 1 ,
122       $                        A((JJ - 1)*DESCA( LLD_ ) + II ) , DESCA( LLD_ ) ,
123       $                        WORK( JW ) , 1 )
124                                CALL ZAXPY( L - 1 , - ONE , WORK( JW ) , 1 ,
125       $                        A(( JJ + L - 2 )*DESCA( LLD_ ) + II ) , 1 )
126                            END IF
127                            CALL PZELSET( A , I , J - 1 , DESCA , EI )
128                        END IF
129  
130  *                     Generate the elementary reflector H(i) to annihilate
131  *                     A(ia + k + i : ia + n - 1 , j)
132  
133                        CALL PZLARFG ( N - K - L + 1 , EI , I + 1 , J , A , MIN( I + 2 , N + IA - 1 ) , J ,
134       $                DESCA , 1 , TAU )
135                        CALL PZELSET( A , I + 1 , J , DESCA , ONE )
136  
137  *                     Compute Y(iy : y + n - 1 , jy + l - 1)
138  
139                        CALL PZGEMV( 'No transpose' , N , N - K - L + 1 , ONE , A , IA , J + 1 ,
140       $                DESCA , A , I + 1 , J , DESCA , 1 , ZERO , Y , IY , JY + L - 1 ,
141       $                DESCY , 1 )
142                        CALL PZGEMV( 'Conjugate transpose' , N - K - L + 1 , L - 1 , ONE , A , I + 1 ,
143       $                JA , DESCA , A , I + 1 , J , DESCA , 1 , ZERO , WORK , 1 , JW ,
144       $                DESCW , DESCW( M_ ) )
145                        CALL PZGEMV( 'No transpose' , N , L - 1 , - ONE , Y , IY , JY , DESCY ,
146       $                WORK , 1 , JW , DESCW , DESCW( M_ ) , ONE , Y , IY ,
147       $                JY + L - 1 , DESCY , 1 )
148                        JL = MIN( JJ + L - 1 , JA + NQ - 1 )
149                        CALL PZSCAL( N , TAU( JL ) , Y , IY , JY + L - 1 , DESCY , 1 )
150  
151  *                     Compute T(1 : i , i)
152  
153                        IF( IPROC ) THEN
154                            JT =( L - 1 ) * DESCA( NB_ )
155                            CALL ZSCAL( L - 1 , - TAU( JL ) , WORK( JW ) , 1 )
156                            CALL ZCOPY( L - 1 , WORK( JW ) , 1 , T( JT + 1 ) , 1 )
157                            CALL ZTRMV( 'Upper' , 'No transpose' , 'Non - unit' , L - 1 , T ,
158       $                    DESCA( NB_ ) , T( JT + 1 ) , 1 )
159                            T( JT + L ) = TAU( JL )
160                        END IF
161     10         CONTINUE
162  
163                CALL PZELSET( A , K + NB + IA - 1 , J , DESCA , EI )
164  
165                RETURN
166  
167  *             End of PZLAHRD
168  
169            END