Routine: PZHENGST()  File: SRC\pzhengst.f

 
 
# lines: 427
  # code: 427
  # comment: 0
  # blank:0
# Variables:62
# Callers:1
# Callings:2
# Words:228
# Keywords:135
 

 

 
..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PZHENGST reduces a complex Hermitian-definite generalized
  eigenproblem to standard form.
  PZHENGST performs the same function as PZHEGST, but is based on
  rank 2K updates, which are faster and more scalable than
  triangular solves (the basis of PZHENGST).
  PZHENGST calls PZHEGST when UPLO='U', hence PZHENGST provides
  improved performance only when UPLO='L', IBTYPE=1.
  PZHENGST also calls PZHEGST when insufficient workspace is
  provided,  hence PZHENGST provides improved
  performance only when LWORK >= 2 * NP0 * NB + NQ0 * NB + NB * NB
  In the following sub( A ) denotes A( IA:IA+N-1, JA:JA+N-1 ) and
  sub( B ) denotes B( IB:IB+N-1, JB:JB+N-1 ).
  If IBTYPE = 1, the problem is sub( A )*x = lambda*sub( B )*x,
  and sub( A ) is overwritten by inv(U**H)*sub( A )*inv(U) or
  inv(L)*sub( A )*inv(L**H)
  If IBTYPE = 2 or 3, the problem is sub( A )*sub( B )*x = lambda*x or
  sub( B )*sub( A )*x = lambda*x, and sub( A ) is overwritten by
  U*sub( A )*U**H or L**H*sub( A )*L.
  sub( B ) must have been previously factorized as U**H*U or L*L**H by
  PZPOTRF.
  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
  =========
  IBTYPE   (global input) INTEGER
          = 1: compute inv(U**H)*sub( A )*inv(U) or
               inv(L)*sub( A )*inv(L**H);
          = 2 or 3: compute U*sub( A )*U**H or L**H*sub( A )*L.
  UPLO    (global input) CHARACTER
          = 'U':  Upper triangle of sub( A ) is stored and sub( B ) is
                  factored as U**H*U;
          = 'L':  Lower triangle of sub( A ) is stored and sub( B ) is
                  factored as L*L**H.
  N       (global input) INTEGER
          The order of the matrices sub( A ) and sub( B ).  N >= 0.
  A       (local input/local output) COMPLEX*16 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
          N-by-N 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 INFO = 0, the transformed matrix, stored in the
          same format as sub( A ).
  IA      (global input) INTEGER
          A's global row index, which points to the beginning of the
          submatrix which is to be operated on.
  JA      (global input) INTEGER
          A's global column index, which points to the beginning of
          the submatrix which is to be operated on.
  DESCA   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix A.
  B       (local input) COMPLEX*16 pointer into the local memory
          to an array of dimension (LLD_B, LOCc(JB+N-1)). On entry,
          this array contains the local pieces of the triangular factor
          from the Cholesky factorization of sub( B ), as returned by
          PZPOTRF.
  IB      (global input) INTEGER
          B's global row index, which points to the beginning of the
          submatrix which is to be operated on.
  JB      (global input) INTEGER
          B's global column index, which points to the beginning of
          the submatrix which is to be operated on.
  DESCB   (global and local input) INTEGER array of dimension DLEN_.
          The array descriptor for the distributed matrix B.
  SCALE   (global output) DOUBLE PRECISION
          Amount by which the eigenvalues should be scaled to
          compensate for the scaling performed in this routine.
          At present, SCALE is always returned as 1.0, it is
          returned here to allow for future enhancement.
  WORK    (local workspace/local output) COMPLEX*16 array,
                                                  dimension (LWORK)
          On exit, WORK( 1 ) returns the minimal and optimal LWORK.
  LWORK   (local or global input) INTEGER
          The dimension of the array WORK.
          LWORK is local input and must be at least
          LWORK >= MAX( NB * ( NP0 +1 ), 3 * NB )
          When IBTYPE = 1 and UPLO = 'L', PZHENGST provides improved
          performance when LWORK >= 2 * NP0 * NB + NQ0 * NB + NB * NB
          where NB = MB_A = NB_A,
          NP0 = NUMROC( N, NB, 0, 0, NPROW ),
          NQ0 = NUMROC( N, NB, 0, 0, NPROW ),
          NUMROC ia a ScaLAPACK tool functions
          MYROW, MYCOL, NPROW and NPCOL can be determined by calling
          the subroutine BLACS_GRIDINFO.
          If LWORK = -1, then LWORK is global input and a workspace
          query is assumed; the routine only calculates the
          optimal size for all work arrays. Each of these
          values is returned in the first entry of the corresponding
          work array, and no error message is issued by PXERBLA.
  INFO    (global output) INTEGER
          = 0:  successful exit
          < 0:  If the i-th argument is an array and the j-entry had
                an illegal value, then INFO = -(i*100+j), if the i-th
                argument is a scalar and had an illegal value, then
                INFO = -i.
  =====================================================================
     .. Parameters ..
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE PZHENGST( IBTYPE , UPLO , N , A , IA , JA , DESCA , B , IB , JB ,
002       $DESCB , SCALE , WORK , LWORK , INFO )
003  
004  *     -- ScaLAPACK routine(version 1.7) --
005  *     University of Tennessee , Knoxville , Oak Ridge National Laboratory ,
006  *     and University of California , Berkeley.
007  *     October 15 , 1999
008  
009  *     .. Scalar Arguments ..
010        CHARACTER UPLO
011        INTEGER IA , IB , IBTYPE , INFO , JA , JB , LWORK , N
012        DOUBLE PRECISION SCALE
013        COMPLEX*16 ONEHALF , ONE , MONE
014        DOUBLE PRECISION RONE
015        PARAMETER( ONEHALF =( 0.5D0 , 0.0D0 ) ,
016       $ONE =( 1.0D0 , 0.0D0 ) ,
017       $MONE =( - 1.0D0 , 0.0D0 ) , RONE = 1.0D0 )
018        INTEGER DLEN_ , CTXT_ , MB_ , NB_ , RSRC_ , CSRC_ , LLD_
019        PARAMETER( DLEN_ = 9 , CTXT_ = 2 , MB_ = 5 , NB_ = 6 ,
020       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
021  *     ..
022  *     .. Local Scalars ..
023        LOGICAL LQUERY , UPPER
024        INTEGER I , IACOL , IAROW , IBCOL , IBROW , ICOFFA , ICOFFB ,
025       $ICTXT , INDAA , INDG , INDR , INDRT , IROFFA ,
026       $IROFFB , J , K , KB , LWMIN , LWOPT , MYCOL , MYROW ,
027       $NB , NP0 , NPCOL , NPK , NPROW , NQ0 , POSTK
028  *     ..
029  *     .. Local Arrays ..
030        INTEGER DESCAA( DLEN_ ) , DESCG( DLEN_ ) ,
031       $DESCR( DLEN_ ) , DESCRT( DLEN_ ) , IDUM1( 2 ) ,
032       $IDUM2( 2 )
033  *     ..
034  *     .. External Functions ..
035        LOGICAL LSAME
036        INTEGER INDXG2P , NUMROC
037        EXTERNAL LSAME , INDXG2P , NUMROC
038  *     ..
039  *     .. External Subroutines ..
040        EXTERNAL BLACS_GRIDINFO , CHK1MAT , DESCSET , PCHK2MAT ,
041       $PXERBLA , PZGEMM , PZHEGST , PZHEMM , PZHER2K ,
042       $PZLACPY , PZTRSM
043  *     ..
044  *     .. Intrinsic Functions ..
045        INTRINSIC DBLE , DCMPLX , DCONJG , ICHAR , MAX , MIN , MOD
046  *     ..
047  *     .. Executable Statements ..
048        ICTXT = DESCA( CTXT_ )
049        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
050        SCALE = 1.0D0
051  
052        NB = DESCA( MB_ )
053  
054  *     Test the input parameters
055  
056        INFO = 0
057        IF( NPROW.EQ. - 1 ) THEN
058            INFO = - ( 700 + CTXT_ )
059        ELSE
060            UPPER = LSAME( UPLO , 'U' )
061            CALL CHK1MAT( N , 3 , N , 3 , IA , JA , DESCA , 7 , INFO )
062            CALL CHK1MAT( N , 3 , N , 3 , IB , JB , DESCB , 11 , INFO )
063            IF( INFO.EQ.0 ) THEN
064                IAROW = INDXG2P( IA , DESCA( MB_ ) , MYROW , DESCA( RSRC_ ) ,
065       $        NPROW )
066                IBROW = INDXG2P( IB , DESCB( MB_ ) , MYROW , DESCB( RSRC_ ) ,
067       $        NPROW )
068                IACOL = INDXG2P( JA , DESCA( NB_ ) , MYCOL , DESCA( CSRC_ ) ,
069       $        NPCOL )
070                IBCOL = INDXG2P( JB , DESCB( NB_ ) , MYCOL , DESCB( CSRC_ ) ,
071       $        NPCOL )
072                IROFFA = MOD( IA - 1 , DESCA( MB_ ) )
073                ICOFFA = MOD( JA - 1 , DESCA( NB_ ) )
074                IROFFB = MOD( IB - 1 , DESCB( MB_ ) )
075                ICOFFB = MOD( JB - 1 , DESCB( NB_ ) )
076                NP0 = NUMROC( N , NB , 0 , 0 , NPROW )
077                NQ0 = NUMROC( N , NB , 0 , 0 , NPCOL )
078                LWMIN = MAX( NB*( NP0 + 1 ) , 3*NB )
079                IF( IBTYPE.EQ.1 .AND. .NOT.UPPER ) THEN
080                    LWOPT = 2*NP0*NB + NQ0*NB + NB*NB
081                ELSE
082                    LWOPT = LWMIN
083                END IF
084                WORK( 1 ) = DCMPLX( DBLE( LWOPT ) )
085                LQUERY =( LWORK.EQ. - 1 )
086                IF( IBTYPE.LT.1 .OR. IBTYPE.GT.3 ) THEN
087                    INFO = - 1
088                ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO , 'L' ) ) THEN
089                    INFO = - 2
090                ELSE IF( N.LT.0 ) THEN
091                    INFO = - 3
092                ELSE IF( IROFFA.NE.0 ) THEN
093                    INFO = - 5
094                ELSE IF( ICOFFA.NE.0 ) THEN
095                    INFO = - 6
096                ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
097                    INFO = - ( 700 + NB_ )
098                ELSE IF( IROFFB.NE.0 .OR. IBROW.NE.IAROW ) THEN
099                    INFO = - 9
100                ELSE IF( ICOFFB.NE.0 .OR. IBCOL.NE.IACOL ) THEN
101                    INFO = - 10
102                ELSE IF( DESCB( MB_ ).NE.DESCA( MB_ ) ) THEN
103                    INFO = - ( 1100 + MB_ )
104                ELSE IF( DESCB( NB_ ).NE.DESCA( NB_ ) ) THEN
105                    INFO = - ( 1100 + NB_ )
106                ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN
107                    INFO = - ( 1100 + CTXT_ )
108                ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
109                    INFO = - 13
110                END IF
111            END IF
112            IDUM1( 1 ) = IBTYPE
113            IDUM2( 1 ) = 1
114            IF( UPPER ) THEN
115                IDUM1( 2 ) = ICHAR( 'U' )
116            ELSE
117                IDUM1( 2 ) = ICHAR( 'L' )
118            END IF
119            IDUM2( 2 ) = 2
120            CALL PCHK2MAT( N , 3 , N , 3 , IA , JA , DESCA , 7 , N , 3 , N , 3 , IB ,
121       $    JB , DESCB , 11 , 2 , IDUM1 , IDUM2 , INFO )
122        END IF
123  
124        IF( INFO.NE.0 ) THEN
125            CALL PXERBLA( ICTXT , 'PZHENGST' , - INFO )
126            RETURN
127        ELSE IF( LQUERY ) THEN
128            RETURN
129        END IF
130  
131  *     Quick return if possible
132  
133        IF( N.EQ.0 )
134       $    RETURN
135  
136            IF( IBTYPE.NE.1 .OR. UPPER .OR. LWORK.LT.LWOPT ) THEN
137                CALL PZHEGST ( IBTYPE , UPLO , N , A , IA , JA , DESCA , B , IB , JB ,
138       $        DESCB , SCALE , INFO )
139                RETURN
140            END IF
141  
142            CALL DESCSET( DESCG , N , NB , NB , NB , IAROW , IACOL , ICTXT , NP0 )
143            CALL DESCSET( DESCR , N , NB , NB , NB , IAROW , IACOL , ICTXT , NP0 )
144            CALL DESCSET( DESCRT , NB , N , NB , NB , IAROW , IACOL , ICTXT , NB )
145            CALL DESCSET( DESCAA , NB , NB , NB , NB , IAROW , IACOL , ICTXT , NB )
146  
147            INDG = 1
148            INDR = INDG + DESCG( LLD_ )*NB
149            INDAA = INDR + DESCR( LLD_ )*NB
150            INDRT = INDAA + DESCAA( LLD_ )*NB
151  
152            DO 30 K = 1 , N , NB
153  
154                KB = MIN( N - K + 1 , NB )
155                POSTK = K + KB
156                NPK = N - POSTK + 1
157  
158                CALL PZLACPY ( 'A' , N - POSTK + 1 , KB , B , POSTK + IB - 1 , K + JB - 1 , DESCB ,
159       $        WORK( INDG ) , POSTK , 1 , DESCG )
160                CALL PZLACPY ( 'A' , N - POSTK + 1 , KB , A , POSTK + IA - 1 , K + JA - 1 , DESCA ,
161       $        WORK( INDR ) , POSTK , 1 , DESCR )
162                CALL PZLACPY ( 'A' , KB , K - 1 , A , K + IA - 1 , JA , DESCA ,
163       $        WORK( INDRT ) , 1 , 1 , DESCRT )
164  
165                CALL PZLACPY ( 'L' , KB , KB , A , K + IA - 1 , K + JA - 1 , DESCA ,
166       $        WORK( INDR ) , K , 1 , DESCR )
167                CALL PZTRSM( 'Right' , 'L' , 'N' , 'N' , NPK , KB , MONE , B , K + IB - 1 ,
168       $        K + JB - 1 , DESCB , WORK( INDG ) , POSTK , 1 , DESCG )
169  
170                CALL PZHEMM( 'Right' , 'L' , NPK , KB , ONEHALF , A , K + IA - 1 , K + JA - 1 ,
171       $        DESCA , WORK( INDG ) , POSTK , 1 , DESCG , ONE ,
172       $        WORK( INDR ) , POSTK , 1 , DESCR )
173  
174                CALL PZHER2K( 'Lower' , 'No T' , NPK , KB , ONE , WORK( INDG ) ,
175       $        POSTK , 1 , DESCG , WORK( INDR ) , POSTK , 1 , DESCR ,
176       $        RONE , A , POSTK + IA - 1 , POSTK + JA - 1 , DESCA )
177  
178                CALL PZGEMM( 'No T' , 'No Conj' , NPK , K - 1 , KB , ONE ,
179       $        WORK( INDG ) , POSTK , 1 , DESCG , WORK( INDRT ) , 1 ,
180       $        1 , DESCRT , ONE , A , POSTK + IA - 1 , JA , DESCA )
181  
182                CALL PZHEMM( 'Right' , 'L' , NPK , KB , ONE , WORK( INDR ) , K , 1 ,
183       $        DESCR , WORK( INDG ) , POSTK , 1 , DESCG , ONE , A ,
184       $        POSTK + IA - 1 , K + JA - 1 , DESCA )
185  
186                CALL PZTRSM( 'Left' , 'Lower' , 'No Conj' , 'Non - unit' , KB , K - 1 ,
187       $        ONE , B , K + IB - 1 , K + JB - 1 , DESCB , A , K + IA - 1 , JA ,
188       $        DESCA )
189  
190                CALL PZLACPY ( 'L' , KB , KB , A , K + IA - 1 , K + JA - 1 , DESCA ,
191       $        WORK( INDAA ) , 1 , 1 , DESCAA )
192  
193                IF( MYROW.EQ.DESCAA( RSRC_ ) .AND. MYCOL.EQ.DESCAA( CSRC_ ) )
194       $            THEN
195                    DO 20 I = 1 , KB
196                        DO 10 J = 1 , I
197                            WORK( INDAA + J - 1 + ( I - 1 )*DESCAA( LLD_ ) )
198       $                    = DCONJG( WORK( INDAA + I - 1 + ( J - 1 )*
199       $                    DESCAA( LLD_ ) ) )
200     10                 CONTINUE
201     20             CONTINUE
202                END IF
203  
204                CALL PZTRSM( 'Left' , 'Lower' , 'No Conj' , 'Non - unit' , KB , KB ,
205       $        ONE , B , K + IB - 1 , K + JB - 1 , DESCB , WORK( INDAA ) , 1 ,
206       $        1 , DESCAA )
207  
208                CALL PZTRSM( 'Right' , 'Lower' , 'Conj' , 'Non - unit' , KB , KB , ONE ,
209       $        B , K + IB - 1 , K + JB - 1 , DESCB , WORK( INDAA ) , 1 , 1 ,
210       $        DESCAA )
211  
212                CALL PZLACPY ( 'L' , KB , KB , WORK( INDAA ) , 1 , 1 , DESCAA , A ,
213       $        K + IA - 1 , K + JA - 1 , DESCA )
214  
215                CALL PZTRSM( 'Right' , 'Lower' , 'Conj' , 'Non - unit' , NPK , KB ,
216       $        ONE , B , K + IB - 1 , K + JB - 1 , DESCB , A , POSTK + IA - 1 ,
217       $        K + JA - 1 , DESCA )
218  
219                DESCR( CSRC_ ) = MOD( DESCR( CSRC_ ) + 1 , NPCOL )
220                DESCG( CSRC_ ) = MOD( DESCG( CSRC_ ) + 1 , NPCOL )
221                DESCRT( RSRC_ ) = MOD( DESCRT( RSRC_ ) + 1 , NPROW )
222                DESCAA( RSRC_ ) = MOD( DESCAA( RSRC_ ) + 1 , NPROW )
223                DESCAA( CSRC_ ) = MOD( DESCAA( CSRC_ ) + 1 , NPCOL )
224     30     CONTINUE
225  
226            WORK( 1 ) = DCMPLX( DBLE( LWOPT ) )
227  
228            RETURN
229        END