Routine: PDSYNGST()  File: SRC\pdsyngst.f

 
 
# lines: 423
  # code: 423
  # comment: 0
  # blank:0
# Variables:61
# Callers:1
# Callings:2
# Words:222
# Keywords:134
 

 

 
..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PDSYNGST reduces a complex Hermitian-definite generalized
  eigenproblem to standard form.
  PDSYNGST performs the same function as PDHEGST, but is based on
  rank 2K updates, which are faster and more scalable than
  triangular solves (the basis of PDSYNGST).
  PDSYNGST calls PDHEGST when UPLO='U', hence PDHENGST provides
  improved performance only when UPLO='L', IBTYPE=1.
  PDSYNGST also calls PDHEGST when insufficient workspace is
  provided,  hence PDSYNGST 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
  PDPOTRF.
  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) DOUBLE PRECISION 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) DOUBLE PRECISION 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
          PDPOTRF.
  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) DOUBLE PRECISION 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', PDSYNGST 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 PDSYNGST( 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        DOUBLE PRECISION ONEHALF , ONE , MONE
014        PARAMETER( ONEHALF = 0.5D0 , ONE = 1.0D0 , MONE = - 1.0D0 )
015        INTEGER DLEN_ , CTXT_ , MB_ , NB_ , RSRC_ , CSRC_ , LLD_
016        PARAMETER( DLEN_ = 9 , CTXT_ = 2 , MB_ = 5 , NB_ = 6 ,
017       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
018  *     ..
019  *     .. Local Scalars ..
020        LOGICAL LQUERY , UPPER
021        INTEGER I , IACOL , IAROW , IBCOL , IBROW , ICOFFA , ICOFFB ,
022       $ICTXT , INDAA , INDG , INDR , INDRT , IROFFA ,
023       $IROFFB , J , K , KB , LWMIN , LWOPT , MYCOL , MYROW ,
024       $NB , NP0 , NPCOL , NPK , NPROW , NQ0 , POSTK
025  *     ..
026  *     .. Local Arrays ..
027        INTEGER DESCAA( DLEN_ ) , DESCG( DLEN_ ) ,
028       $DESCR( DLEN_ ) , DESCRT( DLEN_ ) , IDUM1( 2 ) ,
029       $IDUM2( 2 )
030  *     ..
031  *     .. External Functions ..
032        LOGICAL LSAME
033        INTEGER INDXG2P , NUMROC
034        EXTERNAL LSAME , INDXG2P , NUMROC
035  *     ..
036  *     .. External Subroutines ..
037        EXTERNAL BLACS_GRIDINFO , CHK1MAT , DESCSET , PCHK2MAT ,
038       $PDGEMM , PDLACPY , PDSYGST , PDSYMM , PDSYR2K ,
039       $PDTRSM , PXERBLA
040  *     ..
041  *     .. Intrinsic Functions ..
042        INTRINSIC DBLE , ICHAR , MAX , MIN , MOD
043  *     ..
044  *     .. Executable Statements ..
045        ICTXT = DESCA( CTXT_ )
046        CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
047        SCALE = 1.0D0
048  
049        NB = DESCA( MB_ )
050  
051  *     Test the input parameters
052  
053        INFO = 0
054        IF( NPROW.EQ. - 1 ) THEN
055            INFO = - ( 700 + CTXT_ )
056        ELSE
057            UPPER = LSAME( UPLO , 'U' )
058            CALL CHK1MAT( N , 3 , N , 3 , IA , JA , DESCA , 7 , INFO )
059            CALL CHK1MAT( N , 3 , N , 3 , IB , JB , DESCB , 11 , INFO )
060            IF( INFO.EQ.0 ) THEN
061                IAROW = INDXG2P( IA , DESCA( MB_ ) , MYROW , DESCA( RSRC_ ) ,
062       $        NPROW )
063                IBROW = INDXG2P( IB , DESCB( MB_ ) , MYROW , DESCB( RSRC_ ) ,
064       $        NPROW )
065                IACOL = INDXG2P( JA , DESCA( NB_ ) , MYCOL , DESCA( CSRC_ ) ,
066       $        NPCOL )
067                IBCOL = INDXG2P( JB , DESCB( NB_ ) , MYCOL , DESCB( CSRC_ ) ,
068       $        NPCOL )
069                IROFFA = MOD( IA - 1 , DESCA( MB_ ) )
070                ICOFFA = MOD( JA - 1 , DESCA( NB_ ) )
071                IROFFB = MOD( IB - 1 , DESCB( MB_ ) )
072                ICOFFB = MOD( JB - 1 , DESCB( NB_ ) )
073                NP0 = NUMROC( N , NB , 0 , 0 , NPROW )
074                NQ0 = NUMROC( N , NB , 0 , 0 , NPCOL )
075                LWMIN = MAX( NB*( NP0 + 1 ) , 3*NB )
076                IF( IBTYPE.EQ.1 .AND. .NOT.UPPER ) THEN
077                    LWOPT = 2*NP0*NB + NQ0*NB + NB*NB
078                ELSE
079                    LWOPT = LWMIN
080                END IF
081                WORK( 1 ) = DBLE( LWOPT )
082                LQUERY =( LWORK.EQ. - 1 )
083                IF( IBTYPE.LT.1 .OR. IBTYPE.GT.3 ) THEN
084                    INFO = - 1
085                ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO , 'L' ) ) THEN
086                    INFO = - 2
087                ELSE IF( N.LT.0 ) THEN
088                    INFO = - 3
089                ELSE IF( IROFFA.NE.0 ) THEN
090                    INFO = - 5
091                ELSE IF( ICOFFA.NE.0 ) THEN
092                    INFO = - 6
093                ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
094                    INFO = - ( 700 + NB_ )
095                ELSE IF( IROFFB.NE.0 .OR. IBROW.NE.IAROW ) THEN
096                    INFO = - 9
097                ELSE IF( ICOFFB.NE.0 .OR. IBCOL.NE.IACOL ) THEN
098                    INFO = - 10
099                ELSE IF( DESCB( MB_ ).NE.DESCA( MB_ ) ) THEN
100                    INFO = - ( 1100 + MB_ )
101                ELSE IF( DESCB( NB_ ).NE.DESCA( NB_ ) ) THEN
102                    INFO = - ( 1100 + NB_ )
103                ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN
104                    INFO = - ( 1100 + CTXT_ )
105                ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
106                    INFO = - 13
107                END IF
108            END IF
109            IDUM1( 1 ) = IBTYPE
110            IDUM2( 1 ) = 1
111            IF( UPPER ) THEN
112                IDUM1( 2 ) = ICHAR( 'U' )
113            ELSE
114                IDUM1( 2 ) = ICHAR( 'L' )
115            END IF
116            IDUM2( 2 ) = 2
117            CALL PCHK2MAT( N , 3 , N , 3 , IA , JA , DESCA , 7 , N , 3 , N , 3 , IB ,
118       $    JB , DESCB , 11 , 2 , IDUM1 , IDUM2 , INFO )
119        END IF
120  
121        IF( INFO.NE.0 ) THEN
122            CALL PXERBLA( ICTXT , 'PDSYNGST' , - INFO )
123            RETURN
124        ELSE IF( LQUERY ) THEN
125            RETURN
126        END IF
127  
128  *     Quick return if possible
129  
130        IF( N.EQ.0 )
131       $    RETURN
132  
133            IF( IBTYPE.NE.1 .OR. UPPER .OR. LWORK.LT.LWOPT ) THEN
134                CALL PDSYGST ( IBTYPE , UPLO , N , A , IA , JA , DESCA , B , IB , JB ,
135       $        DESCB , SCALE , INFO )
136                RETURN
137            END IF
138  
139            CALL DESCSET( DESCG , N , NB , NB , NB , IAROW , IACOL , ICTXT , NP0 )
140            CALL DESCSET( DESCR , N , NB , NB , NB , IAROW , IACOL , ICTXT , NP0 )
141            CALL DESCSET( DESCRT , NB , N , NB , NB , IAROW , IACOL , ICTXT , NB )
142            CALL DESCSET( DESCAA , NB , NB , NB , NB , IAROW , IACOL , ICTXT , NB )
143  
144            INDG = 1
145            INDR = INDG + DESCG( LLD_ )*NB
146            INDAA = INDR + DESCR( LLD_ )*NB
147            INDRT = INDAA + DESCAA( LLD_ )*NB
148  
149            DO 30 K = 1 , N , NB
150  
151                KB = MIN( N - K + 1 , NB )
152                POSTK = K + KB
153                NPK = N - POSTK + 1
154  
155                CALL PDLACPY ( 'A' , N - POSTK + 1 , KB , B , POSTK + IB - 1 , K + JB - 1 , DESCB ,
156       $        WORK( INDG ) , POSTK , 1 , DESCG )
157                CALL PDLACPY ( 'A' , N - POSTK + 1 , KB , A , POSTK + IA - 1 , K + JA - 1 , DESCA ,
158       $        WORK( INDR ) , POSTK , 1 , DESCR )
159                CALL PDLACPY ( 'A' , KB , K - 1 , A , K + IA - 1 , JA , DESCA ,
160       $        WORK( INDRT ) , 1 , 1 , DESCRT )
161  
162                CALL PDLACPY ( 'L' , KB , KB , A , K + IA - 1 , K + JA - 1 , DESCA ,
163       $        WORK( INDR ) , K , 1 , DESCR )
164                CALL PDTRSM( 'Right' , 'L' , 'N' , 'N' , NPK , KB , MONE , B , K + IB - 1 ,
165       $        K + JB - 1 , DESCB , WORK( INDG ) , POSTK , 1 , DESCG )
166  
167                CALL PDSYMM( 'Right' , 'L' , NPK , KB , ONEHALF , A , K + IA - 1 , K + JA - 1 ,
168       $        DESCA , WORK( INDG ) , POSTK , 1 , DESCG , ONE ,
169       $        WORK( INDR ) , POSTK , 1 , DESCR )
170  
171                CALL PDSYR2K( 'Lower' , 'No T' , NPK , KB , ONE , WORK( INDG ) ,
172       $        POSTK , 1 , DESCG , WORK( INDR ) , POSTK , 1 , DESCR ,
173       $        ONE , A , POSTK + IA - 1 , POSTK + JA - 1 , DESCA )
174  
175                CALL PDGEMM( 'No T' , 'No Conj' , NPK , K - 1 , KB , ONE ,
176       $        WORK( INDG ) , POSTK , 1 , DESCG , WORK( INDRT ) , 1 ,
177       $        1 , DESCRT , ONE , A , POSTK + IA - 1 , JA , DESCA )
178  
179                CALL PDSYMM( 'Right' , 'L' , NPK , KB , ONE , WORK( INDR ) , K , 1 ,
180       $        DESCR , WORK( INDG ) , POSTK , 1 , DESCG , ONE , A ,
181       $        POSTK + IA - 1 , K + JA - 1 , DESCA )
182  
183                CALL PDTRSM( 'Left' , 'Lower' , 'No Conj' , 'Non - unit' , KB , K - 1 ,
184       $        ONE , B , K + IB - 1 , K + JB - 1 , DESCB , A , K + IA - 1 , JA ,
185       $        DESCA )
186  
187                CALL PDLACPY ( 'L' , KB , KB , A , K + IA - 1 , K + JA - 1 , DESCA ,
188       $        WORK( INDAA ) , 1 , 1 , DESCAA )
189  
190                IF( MYROW.EQ.DESCAA( RSRC_ ) .AND. MYCOL.EQ.DESCAA( CSRC_ ) )
191       $            THEN
192                    DO 20 I = 1 , KB
193                        DO 10 J = 1 , I
194                            WORK( INDAA + J - 1 + ( I - 1 )*DESCAA( LLD_ ) )
195       $                    = WORK( INDAA + I - 1 + ( J - 1 )*DESCAA( LLD_ ) )
196     10                 CONTINUE
197     20             CONTINUE
198                END IF
199  
200                CALL PDTRSM( 'Left' , 'Lower' , 'No Conj' , 'Non - unit' , KB , KB ,
201       $        ONE , B , K + IB - 1 , K + JB - 1 , DESCB , WORK( INDAA ) , 1 ,
202       $        1 , DESCAA )
203  
204                CALL PDTRSM( 'Right' , 'Lower' , 'Conj' , 'Non - unit' , KB , KB , ONE ,
205       $        B , K + IB - 1 , K + JB - 1 , DESCB , WORK( INDAA ) , 1 , 1 ,
206       $        DESCAA )
207  
208                CALL PDLACPY ( 'L' , KB , KB , WORK( INDAA ) , 1 , 1 , DESCAA , A ,
209       $        K + IA - 1 , K + JA - 1 , DESCA )
210  
211                CALL PDTRSM( 'Right' , 'Lower' , 'Conj' , 'Non - unit' , NPK , KB ,
212       $        ONE , B , K + IB - 1 , K + JB - 1 , DESCB , A , POSTK + IA - 1 ,
213       $        K + JA - 1 , DESCA )
214  
215                DESCR( CSRC_ ) = MOD( DESCR( CSRC_ ) + 1 , NPCOL )
216                DESCG( CSRC_ ) = MOD( DESCG( CSRC_ ) + 1 , NPCOL )
217                DESCRT( RSRC_ ) = MOD( DESCRT( RSRC_ ) + 1 , NPROW )
218                DESCAA( RSRC_ ) = MOD( DESCAA( RSRC_ ) + 1 , NPROW )
219                DESCAA( CSRC_ ) = MOD( DESCAA( CSRC_ ) + 1 , NPCOL )
220     30     CONTINUE
221  
222            WORK( 1 ) = DBLE( LWOPT )
223  
224            RETURN
225        END