Routine: PCLATRD()  File: SRC\pclatrd.f

 
 
# lines: 435
  # code: 435
  # comment: 0
  # blank:0
# Variables:45
# Callers:2
# Callings:2
# Words:203
# Keywords:95
 

 

 
..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PCLATRD reduces NB rows and columns of a complex Hermitian
  distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1) to complex
  tridiagonal form by an unitary similarity transformation
  Q' * sub( A ) * Q, and returns the matrices V and W which are
  needed to apply the transformation to the unreduced part of sub( A ).
  If UPLO = 'U', PCLATRD reduces the last NB rows and columns of a
  matrix, of which the upper triangle is supplied;
  if UPLO = 'L', PCLATRD reduces the first NB rows and columns of a
  matrix, of which the lower triangle is supplied.
  This is an auxiliary routine called by PCHETRD.
  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
  =========
  UPLO    (global input) CHARACTER
          Specifies whether the upper or lower triangular part of the
          Hermitian matrix sub( A ) is stored:
          = 'U': Upper triangular
          = 'L': Lower triangular
  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.
  NB      (global input) INTEGER
          The number of rows and columns to be reduced.
  A       (local input/local output) COMPLEX pointer into the
          local memory to an array of dimension (LLD_A,LOCc(JA+N-1)).
          On entry, this array contains the local pieces of the
          Hermitian distributed matrix sub( A ).  If UPLO = 'U', the
          leading N-by-N upper triangular part of sub( A ) contains
          the upper triangular part of the matrix, and its strictly
          lower triangular part is not referenced. If UPLO = 'L', the
          leading N-by-N lower triangular part of sub( A ) contains the
          lower triangular part of the matrix, and its strictly upper
          triangular part is not referenced.
          On exit, if UPLO = 'U', the last NB columns have been reduced
          to tridiagonal form, with the diagonal elements overwriting
          the diagonal elements of sub( A ); the elements above the
          diagonal with the array TAU, represent the unitary matrix Q
          as a product of elementary reflectors. If UPLO = 'L', the
          first NB columns have been reduced to tridiagonal form, with
          the diagonal elements overwriting the diagonal elements of
          sub( A ); the elements below the diagonal with the array TAU,
          represent the unitary matrix Q as a product of elementary
          reflectors; 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.
  D       (local output) REAL array, dimension LOCc(JA+N-1)
          The diagonal elements of the tridiagonal matrix T:
          D(i) = A(i,i). D is tied to the distributed matrix A.
  E       (local output) REAL array, dimension LOCc(JA+N-1)
          if UPLO = 'U', LOCc(JA+N-2) otherwise. The off-diagonal
          elements of the tridiagonal matrix T: E(i) = A(i,i+1) if
          UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. E is tied to the
          distributed matrix A.
  TAU     (local output) COMPLEX, array, dimension
          LOCc(JA+N-1). This array contains the scalar factors TAU of
          the elementary reflectors. TAU is tied to the distributed
          matrix A.
  W       (local output) COMPLEX pointer into the local memory
          to an array of dimension (LLD_W,NB_W), This array contains
          the local pieces of the N-by-NB_W matrix W required to
          update the unreduced part of sub( A ).
  IW      (global input) INTEGER
          The row index in the global array W indicating the first
          row of sub( W ).
  JW      (global input) INTEGER
          The column index in the global array W indicating the
          first column of sub( W ).
  DESCW   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix W.
  WORK    (local workspace) COMPLEX array, dimension (NB_A)
  Further Details
  ===============
  If UPLO = 'U', the matrix Q is represented as a product of elementary
  reflectors
     Q = H(n) H(n-1) . . . H(n-nb+1).
  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(i:n) = 0 and v(i-1) = 1; v(1:i-1) is stored on exit in
  A(ia:ia+i-2,ja+i), and tau in TAU(ja+i-1).
  If UPLO = 'L', the matrix Q is represented as a product of 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) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in
  A(ia+i+1:ia+n-1,ja+i-1), and tau in TAU(ja+i-1).
  The elements of the vectors v together form the N-by-NB matrix V
  which is needed, with W, to apply the transformation to the unreduced
  part of the matrix, using a Hermitian rank-2k update of the form:
  sub( A ) := sub( A ) - V*W' - W*V'.
  The contents of A on exit are illustrated by the following examples
  with n = 5 and nb = 2:
  if UPLO = 'U':                       if UPLO = 'L':
    (  a   a   a   v4  v5 )              (  d                  )
    (      a   a   v4  v5 )              (  1   d              )
    (          a   1   v5 )              (  v1  1   a          )
    (              d   1  )              (  v1  v2  a   a      )
    (                  d  )              (  v1  v2  a   a   a  )
  where d denotes a diagonal element of the reduced matrix, a denotes
  an element of the original matrix that is unchanged, 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 PCLATRD( UPLO , N , NB , A , IA , JA , DESCA , D , E , TAU , W ,
002       $IW , JW , DESCW , 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 UPLO
011        INTEGER IA , IW , JA , JW , N , NB
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        COMPLEX HALF , ONE , ZERO
018        PARAMETER( HALF =( 0.5E + 0 , 0.0E + 0 ) ,
019       $ONE =( 1.0E + 0 , 0.0E + 0 ) ,
020       $ZERO =( 0.0E + 0 , 0.0E + 0 ) )
021  *     ..
022  *     .. Local Scalars ..
023        INTEGER I , IACOL , IAROW , ICTXT , II , J , JJ , JP , JWK , K ,
024       $KW , MYCOL , MYROW , NPCOL , NPROW , NQ
025        COMPLEX AII , ALPHA , BETA
026  *     ..
027  *     .. Local Arrays ..
028        INTEGER DESCD( DLEN_ ) , DESCE( DLEN_ ) , DESCWK( DLEN_ )
029  *     ..
030  *     .. External Subroutines ..
031        EXTERNAL BLACS_GRIDINFO , DESCSET , INFOG2L , PCAXPY ,
032       $PCDOTC , PCELGET , PCELSET , PCGEMV ,
033       $PCHEMV , PCLACGV , PCLARFG , PCSCAL ,
034       $PSELSET , SGEBR2D , SGEBS2D
035  *     ..
036  *     .. External Functions ..
037        LOGICAL LSAME
038        INTEGER NUMROC
039        EXTERNAL LSAME , NUMROC
040  *     ..
041  *     .. Intrinsic Functions ..
042        INTRINSIC CMPLX , MIN , REAL
043  *     ..
044  *     .. Executable Statements ..
045  
046  *     Quick return if possible
047  
048        IF( N.LE.0 )
049       $    RETURN
050  
051            ICTXT = DESCA( CTXT_ )
052            CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
053            NQ = MAX( 1 , NUMROC( JA + N - 1 , DESCA( NB_ ) , MYCOL , DESCA( CSRC_ ) ,
054       $    NPCOL ) )
055            CALL DESCSET( DESCD , 1 , JA + N - 1 , 1 , DESCA( NB_ ) , MYROW ,
056       $    DESCA( CSRC_ ) , DESCA( CTXT_ ) , 1 )
057            AII = ZERO
058            BETA = ZERO
059  
060            IF( LSAME( UPLO , 'U' ) ) THEN
061  
062                CALL INFOG2L( N + IA - NB , N + JA - NB , DESCA , NPROW , NPCOL , MYROW ,
063       $        MYCOL , II , JJ , IAROW , IACOL )
064                CALL DESCSET( DESCWK , 1 , DESCW( NB_ ) , 1 , DESCW( NB_ ) , IAROW ,
065       $        IACOL , ICTXT , 1 )
066                CALL DESCSET( DESCE , 1 , JA + N - 1 , 1 , DESCA( NB_ ) , MYROW ,
067       $        DESCA( CSRC_ ) , DESCA( CTXT_ ) , 1 )
068  
069  *             Reduce last NB columns of upper triangle
070  
071                DO 10 J = JA + N - 1 , JA + N - NB , - 1
072                    I = IA + J - JA
073                    K = J - JA + 1
074                    KW = MOD( K - 1 , DESCA( MB_ ) ) + 1
075  
076  *                 Update A(IA : I , I)
077  
078                    CALL PCELGET( 'E' , ' ' , AII , A , I , J , DESCA )
079                    CALL PCELSET( A , I , J , DESCA , CMPLX( REAL( AII ) ) )
080                    CALL PCLACGV ( N - K , W , IW + K - 1 , JW + KW , DESCW , DESCW( M_ ) )
081                    CALL PCGEMV( 'No transpose' , K , N - K , - ONE , A , IA , J + 1 ,
082       $            DESCA , W , IW + K - 1 , JW + KW , DESCW , DESCW( M_ ) ,
083       $            ONE , A , IA , J , DESCA , 1 )
084                    CALL PCLACGV ( N - K , W , IW + K - 1 , JW + KW , DESCW , DESCW( M_ ) )
085                    CALL PCLACGV ( N - K , A , I , J + 1 , DESCA , DESCA( M_ ) )
086                    CALL PCGEMV( 'No transpose' , K , N - K , - ONE , W , IW , JW + KW ,
087       $            DESCW , A , I , J + 1 , DESCA , DESCA( M_ ) , ONE , A ,
088       $            IA , J , DESCA , 1 )
089                    CALL PCLACGV ( N - K , A , I , J + 1 , DESCA , DESCA( M_ ) )
090                    CALL PCELGET( 'E' , ' ' , AII , A , I , J , DESCA )
091                    CALL PCELSET( A , I , J , DESCA , CMPLX( REAL( AII ) ) )
092                    IF( N - K.GT.0 )
093       $                CALL PCELSET( A , I , J + 1 , DESCA , CMPLX( E( JP ) ) )
094  
095  *                     Generate elementary reflector H(i) to annihilate
096  *                     A(IA : I - 2 , I)
097  
098                        JP = MIN( JJ + KW - 1 , NQ )
099                        CALL PCLARFG ( K - 1 , BETA , I - 1 , J , A , IA , J , DESCA , 1 ,
100       $                TAU )
101                        CALL PSELSET( E , 1 , J , DESCE , REAL( BETA ) )
102                        CALL PCELSET( A , I - 1 , J , DESCA , ONE )
103  
104  *                     Compute W(IW : IW + K - 2 , JW + KW - 1)
105  
106                        CALL PCHEMV( 'Upper' , K - 1 , ONE , A , IA , JA , DESCA , A , IA , J ,
107       $                DESCA , 1 , ZERO , W , IW , JW + KW - 1 , DESCW , 1 )
108  
109                        JWK = MOD( K - 1 , DESCWK( NB_ ) ) + 2
110                        CALL PCGEMV( 'Conjugate transpose' , K - 1 , N - K , ONE , W , IW ,
111       $                JW + KW , DESCW , A , IA , J , DESCA , 1 , ZERO , WORK ,
112       $                1 , JWK , DESCWK , DESCWK( M_ ) )
113                        CALL PCGEMV( 'No transpose' , K - 1 , N - K , - ONE , A , IA , J + 1 ,
114       $                DESCA , WORK , 1 , JWK , DESCWK , DESCWK( M_ ) , ONE ,
115       $                W , IW , JW + KW - 1 , DESCW , 1 )
116                        CALL PCGEMV( 'Conjugate transpose' , K - 1 , N - K , ONE , A , IA ,
117       $                J + 1 , DESCA , A , IA , J , DESCA , 1 , ZERO , WORK , 1 ,
118       $                JWK , DESCWK , DESCWK( M_ ) )
119                        CALL PCGEMV( 'No transpose' , K - 1 , N - K , - ONE , W , IW , JW + KW ,
120       $                DESCW , WORK , 1 , JWK , DESCWK , DESCWK( M_ ) , ONE ,
121       $                W , IW , JW + KW - 1 , DESCW , 1 )
122                        CALL PCSCAL( K - 1 , TAU( JP ) , W , IW , JW + KW - 1 , DESCW , 1 )
123  
124                        CALL PCDOTC( K - 1 , ALPHA , W , IW , JW + KW - 1 , DESCW , 1 , A , IA , J ,
125       $                DESCA , 1 )
126                        IF( MYCOL.EQ.IACOL )
127       $                    ALPHA = - HALF*TAU( JP )*ALPHA
128                            CALL PCAXPY( K - 1 , ALPHA , A , IA , J , DESCA , 1 , W , IW , JW + KW - 1 ,
129       $                    DESCW , 1 )
130                            CALL PCELGET( 'E' , ' ' , BETA , A , I , J , DESCA )
131                            CALL PSELSET( D , 1 , J , DESCD , REAL( BETA ) )
132  
133     10         CONTINUE
134  
135            ELSE
136  
137                CALL INFOG2L( IA , JA , DESCA , NPROW , NPCOL , MYROW , MYCOL , II ,
138       $        JJ , IAROW , IACOL )
139                CALL DESCSET( DESCWK , 1 , DESCW( NB_ ) , 1 , DESCW( NB_ ) , IAROW ,
140       $        IACOL , ICTXT , 1 )
141                CALL DESCSET( DESCE , 1 , JA + N - 2 , 1 , DESCA( NB_ ) , MYROW ,
142       $        DESCA( CSRC_ ) , DESCA( CTXT_ ) , 1 )
143  
144  *             Reduce first NB columns of lower triangle
145  
146                DO 20 J = JA , JA + NB - 1
147                    I = IA + J - JA
148                    K = J - JA + 1
149  
150  *                 Update A(J : JA + N - 1 , J)
151  
152                    CALL PCELGET( 'E' , ' ' , AII , A , I , J , DESCA )
153                    CALL PCELSET( A , I , J , DESCA , CMPLX( REAL( AII ) ) )
154                    CALL PCLACGV ( K - 1 , W , IW + K - 1 , JW , DESCW , DESCW( M_ ) )
155                    CALL PCGEMV( 'No transpose' , N - K + 1 , K - 1 , - ONE , A , I , JA ,
156       $            DESCA , W , IW + K - 1 , JW , DESCW , DESCW( M_ ) , ONE ,
157       $            A , I , J , DESCA , 1 )
158                    CALL PCLACGV ( K - 1 , W , IW + K - 1 , JW , DESCW , DESCW( M_ ) )
159                    CALL PCLACGV ( K - 1 , A , I , JA , DESCA , DESCA( M_ ) )
160                    CALL PCGEMV( 'No transpose' , N - K + 1 , K - 1 , - ONE , W , IW + K - 1 ,
161       $            JW , DESCW , A , I , JA , DESCA , DESCA( M_ ) , ONE ,
162       $            A , I , J , DESCA , 1 )
163                    CALL PCLACGV ( K - 1 , A , I , JA , DESCA , DESCA( M_ ) )
164                    CALL PCELGET( 'E' , ' ' , AII , A , I , J , DESCA )
165                    CALL PCELSET( A , I , J , DESCA , CMPLX( REAL( AII ) ) )
166                    IF( K.GT.1 )
167       $                CALL PCELSET( A , I , J - 1 , DESCA , CMPLX( E( JP ) ) )
168  
169  *                     Generate elementary reflector H(i) to annihilate
170  *                     A(I + 2 : IA + N - 1 , I)
171  
172                        JP = MIN( JJ + K - 1 , NQ )
173                        CALL PCLARFG ( N - K , BETA , I + 1 , J , A , I + 2 , J , DESCA , 1 ,
174       $                TAU )
175                        CALL PSELSET( E , 1 , J , DESCE , REAL( BETA ) )
176                        CALL PCELSET( A , I + 1 , J , DESCA , ONE )
177  
178  *                     Compute W(IW + K : IW + N - 1 , JW + K - 1)
179  
180                        CALL PCHEMV( 'Lower' , N - K , ONE , A , I + 1 , J + 1 , DESCA , A , I + 1 ,
181       $                J , DESCA , 1 , ZERO , W , IW + K , JW + K - 1 , DESCW , 1 )
182  
183                        CALL PCGEMV( 'Conjugate Transpose' , N - K , K - 1 , ONE , W , IW + K ,
184       $                JW , DESCW , A , I + 1 , J , DESCA , 1 , ZERO , WORK , 1 ,
185       $                1 , DESCWK , DESCWK( M_ ) )
186                        CALL PCGEMV( 'No transpose' , N - K , K - 1 , - ONE , A , I + 1 , JA ,
187       $                DESCA , WORK , 1 , 1 , DESCWK , DESCWK( M_ ) , ONE , W ,
188       $                IW + K , JW + K - 1 , DESCW , 1 )
189                        CALL PCGEMV( 'Conjugate transpose' , N - K , K - 1 , ONE , A , I + 1 ,
190       $                JA , DESCA , A , I + 1 , J , DESCA , 1 , ZERO , WORK , 1 ,
191       $                1 , DESCWK , DESCWK( M_ ) )
192                        CALL PCGEMV( 'No transpose' , N - K , K - 1 , - ONE , W , IW + K , JW ,
193       $                DESCW , WORK , 1 , 1 , DESCWK , DESCWK( M_ ) , ONE , W ,
194       $                IW + K , JW + K - 1 , DESCW , 1 )
195                        CALL PCSCAL( N - K , TAU( JP ) , W , IW + K , JW + K - 1 , DESCW , 1 )
196                        CALL PCDOTC( N - K , ALPHA , W , IW + K , JW + K - 1 , DESCW , 1 , A , I + 1 ,
197       $                J , DESCA , 1 )
198                        IF( MYCOL.EQ.IACOL )
199       $                    ALPHA = - HALF*TAU( JP )*ALPHA
200                            CALL PCAXPY( N - K , ALPHA , A , I + 1 , J , DESCA , 1 , W , IW + K ,
201       $                    JW + K - 1 , DESCW , 1 )
202                            CALL PCELGET( 'E' , ' ' , BETA , A , I , J , DESCA )
203                            CALL PSELSET( D , 1 , J , DESCD , REAL( BETA ) )
204  
205     20         CONTINUE
206  
207            END IF
208  
209  *         Broadcast columnwise the diagonal elements into D.
210  
211            IF( MYCOL.EQ.IACOL ) THEN
212                IF( MYROW.EQ.IAROW ) THEN
213                    CALL SGEBS2D( ICTXT , 'Columnwise' , ' ' , 1 , NB , D( JJ ) , 1 )
214                ELSE
215                    CALL SGEBR2D( ICTXT , 'Columnwise' , ' ' , 1 , NB , D( JJ ) , 1 ,
216       $            IAROW , MYCOL )
217                END IF
218            END IF
219  
220            RETURN
221  
222  *         End of PCLATRD
223  
224        END