Routine: PDLAHRD()  File: SRC\pdlahrd.f

 
 
# lines: 287
  # code: 287
  # comment: 0
  # blank:0
# Variables:42
# Callers:1
# Callings:1
# Words:108
# Keywords:56
 

 

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