Routine: PCHENGST()  File: SRC\pchengst.f

 
 
# lines: 426
  # code: 426
  # comment: 0
  # blank:0
# Variables:62
# Callers:1
# Callings:2
# Words:226
# Keywords:136
 

 

 
..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PCHENGST reduces a complex Hermitian-definite generalized
  eigenproblem to standard form.
  PCHENGST performs the same function as PCHEGST, but is based on
  rank 2K updates, which are faster and more scalable than
  triangular solves (the basis of PCHENGST).
  PCHENGST calls PCHEGST when UPLO='U', hence PCHENGST provides
  improved performance only when UPLO='L', IBTYPE=1.
  PCHENGST also calls PCHEGST when insufficient workspace is
  provided,  hence PCHENGST 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
  PCPOTRF.
  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 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 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
          PCPOTRF.
  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) REAL
          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 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', PCHENGST 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 PCHENGST( 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        REAL SCALE
013        COMPLEX ONEHALF , ONE , MONE
014        REAL RONE
015        PARAMETER( ONEHALF =( 0.5E0 , 0.0E0 ) ,
016       $ONE =( 1.0E0 , 0.0E0 ) ,
017       $MONE =( - 1.0E0 , 0.0E0 ) , RONE = 1.0E0 )
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 , PCGEMM ,
041       $PCHEGST , PCHEMM , PCHER2K , PCHK2MAT , PCLACPY ,
042       $PCTRSM , PXERBLA
043  *     ..
044  *     .. Intrinsic Functions ..
045        INTRINSIC CMPLX , CONJG , ICHAR , MAX , MIN , MOD , REAL
046  *     ..
047  *     .. Executable Statements ..
048        ICTXT = DESCA( CTXT_ )
049        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
050        SCALE = 1.0E0
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 ) = CMPLX( REAL( 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 , 'PCHENGST' , - 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 PCHEGST ( 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 PCLACPY ( 'A' , N - POSTK + 1 , KB , B , POSTK + IB - 1 , K + JB - 1 , DESCB ,
159       $        WORK( INDG ) , POSTK , 1 , DESCG )
160                CALL PCLACPY ( 'A' , N - POSTK + 1 , KB , A , POSTK + IA - 1 , K + JA - 1 , DESCA ,
161       $        WORK( INDR ) , POSTK , 1 , DESCR )
162                CALL PCLACPY ( 'A' , KB , K - 1 , A , K + IA - 1 , JA , DESCA ,
163       $        WORK( INDRT ) , 1 , 1 , DESCRT )
164  
165                CALL PCLACPY ( 'L' , KB , KB , A , K + IA - 1 , K + JA - 1 , DESCA ,
166       $        WORK( INDR ) , K , 1 , DESCR )
167                CALL PCTRSM( 'Right' , 'L' , 'N' , 'N' , NPK , KB , MONE , B , K + IB - 1 ,
168       $        K + JB - 1 , DESCB , WORK( INDG ) , POSTK , 1 , DESCG )
169  
170                CALL PCHEMM( '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 PCHER2K( '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 PCGEMM( '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 PCHEMM( '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 PCTRSM( '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 PCLACPY ( '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       $                    = CONJG( WORK( INDAA + I - 1 + ( J - 1 )*DESCAA( LLD_ ) ) )
199     10                 CONTINUE
200     20             CONTINUE
201                END IF
202  
203                CALL PCTRSM( 'Left' , 'Lower' , 'No Conj' , 'Non - unit' , KB , KB ,
204       $        ONE , B , K + IB - 1 , K + JB - 1 , DESCB , WORK( INDAA ) , 1 ,
205       $        1 , DESCAA )
206  
207                CALL PCTRSM( 'Right' , 'Lower' , 'Conj' , 'Non - unit' , KB , KB , ONE ,
208       $        B , K + IB - 1 , K + JB - 1 , DESCB , WORK( INDAA ) , 1 , 1 ,
209       $        DESCAA )
210  
211                CALL PCLACPY ( 'L' , KB , KB , WORK( INDAA ) , 1 , 1 , DESCAA , A ,
212       $        K + IA - 1 , K + JA - 1 , DESCA )
213  
214                CALL PCTRSM( 'Right' , 'Lower' , 'Conj' , 'Non - unit' , NPK , KB ,
215       $        ONE , B , K + IB - 1 , K + JB - 1 , DESCB , A , POSTK + IA - 1 ,
216       $        K + JA - 1 , DESCA )
217  
218                DESCR( CSRC_ ) = MOD( DESCR( CSRC_ ) + 1 , NPCOL )
219                DESCG( CSRC_ ) = MOD( DESCG( CSRC_ ) + 1 , NPCOL )
220                DESCRT( RSRC_ ) = MOD( DESCRT( RSRC_ ) + 1 , NPROW )
221                DESCAA( RSRC_ ) = MOD( DESCAA( RSRC_ ) + 1 , NPROW )
222                DESCAA( CSRC_ ) = MOD( DESCAA( CSRC_ ) + 1 , NPCOL )
223     30     CONTINUE
224  
225            WORK( 1 ) = CMPLX( REAL( LWOPT ) )
226  
227            RETURN
228        END