Routine: PCHEGST()  File: SRC\pchegst.f

 
 
# lines: 441
  # code: 441
  # comment: 0
  # blank:0
# Variables:45
# Callers:1
# Callings:1
# Words:233
# Keywords:143
 

 

..
     .. Array Arguments ..
     ..
  Purpose
  =======
  PCHEGST reduces a complex Hermitian-definite generalized eigenproblem
  to standard form.
  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.
  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 PCHEGST( IBTYPE , UPLO , N , A , IA , JA , DESCA , B , IB , JB ,
002       $DESCB , SCALE , INFO )
003  
004  *     -- ScaLAPACK 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 , IB , IBTYPE , INFO , JA , JB , N
012        REAL SCALE
013        INTEGER BLOCK_CYCLIC_2D , DLEN_ , DTYPE_ , CTXT_ , M_ , N_ ,
014       $MB_ , NB_ , RSRC_ , CSRC_ , LLD_
015        PARAMETER( BLOCK_CYCLIC_2D = 1 , DLEN_ = 9 , DTYPE_ = 1 ,
016       $CTXT_ = 2 , M_ = 3 , N_ = 4 , MB_ = 5 , NB_ = 6 ,
017       $RSRC_ = 7 , CSRC_ = 8 , LLD_ = 9 )
018        REAL ONE
019        PARAMETER( ONE = 1.0E + 0 )
020        COMPLEX CONE , HALF
021        PARAMETER( CONE =( 1.0E + 0 , 0.0E + 0 ) ,
022       $HALF =( 0.5E + 0 , 0.0E + 0 ) )
023  *     ..
024  *     .. Local Scalars ..
025        LOGICAL UPPER
026        INTEGER IACOL , IAROW , IBCOL , IBROW , ICOFFA , ICOFFB ,
027       $ICTXT , IROFFA , IROFFB , K , KB , MYCOL , MYROW , NB ,
028       $NPCOL , NPROW
029  *     ..
030  *     .. Local Arrays ..
031        INTEGER IDUM1( 2 ) , IDUM2( 2 )
032  *     ..
033  *     .. External Subroutines ..
034        EXTERNAL BLACS_GRIDINFO , CHK1MAT , PCHEGS2 , PCHEMM ,
035       $PCHER2K , PCHK2MAT , PCTRMM , PCTRSM , PXERBLA
036  *     ..
037  *     .. Intrinsic Functions ..
038        INTRINSIC ICHAR , MIN , MOD
039  *     ..
040  *     .. External Functions ..
041        LOGICAL LSAME
042        INTEGER ICEIL , INDXG2P
043        EXTERNAL LSAME , ICEIL , INDXG2P
044  *     ..
045  *     .. Executable Statements ..
046  *     This is just to keep ftnchek happy
047        IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_*
048       $    RSRC_.LT.0 )RETURN
049  
050  *         Get grid parameters
051  
052            SCALE = ONE
053  
054            ICTXT = DESCA( CTXT_ )
055            CALL BLACS_GRIDINFO( ICTXT , NPROW , NPCOL , MYROW , MYCOL )
056  
057  *         Test the input parameters
058  
059            INFO = 0
060            IF( NPROW.EQ. - 1 ) THEN
061                INFO = - ( 700 + CTXT_ )
062            ELSE
063                UPPER = LSAME( UPLO , 'U' )
064                CALL CHK1MAT( N , 3 , N , 3 , IA , JA , DESCA , 7 , INFO )
065                CALL CHK1MAT( N , 3 , N , 3 , IB , JB , DESCB , 11 , INFO )
066                IF( INFO.EQ.0 ) THEN
067                    IAROW = INDXG2P( IA , DESCA( MB_ ) , MYROW , DESCA( RSRC_ ) ,
068       $            NPROW )
069                    IBROW = INDXG2P( IB , DESCB( MB_ ) , MYROW , DESCB( RSRC_ ) ,
070       $            NPROW )
071                    IACOL = INDXG2P( JA , DESCA( NB_ ) , MYCOL , DESCA( CSRC_ ) ,
072       $            NPCOL )
073                    IBCOL = INDXG2P( JB , DESCB( NB_ ) , MYCOL , DESCB( CSRC_ ) ,
074       $            NPCOL )
075                    IROFFA = MOD( IA - 1 , DESCA( MB_ ) )
076                    ICOFFA = MOD( JA - 1 , DESCA( NB_ ) )
077                    IROFFB = MOD( IB - 1 , DESCB( MB_ ) )
078                    ICOFFB = MOD( JB - 1 , DESCB( NB_ ) )
079                    IF( IBTYPE.LT.1 .OR. IBTYPE.GT.3 ) THEN
080                        INFO = - 1
081                    ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO , 'L' ) ) THEN
082                        INFO = - 2
083                    ELSE IF( N.LT.0 ) THEN
084                        INFO = - 3
085                    ELSE IF( IROFFA.NE.0 ) THEN
086                        INFO = - 5
087                    ELSE IF( ICOFFA.NE.0 ) THEN
088                        INFO = - 6
089                    ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
090                        INFO = - ( 700 + NB_ )
091                    ELSE IF( IROFFB.NE.0 .OR. IBROW.NE.IAROW ) THEN
092                        INFO = - 9
093                    ELSE IF( ICOFFB.NE.0 .OR. IBCOL.NE.IACOL ) THEN
094                        INFO = - 10
095                    ELSE IF( DESCB( MB_ ).NE.DESCA( MB_ ) ) THEN
096                        INFO = - ( 1100 + MB_ )
097                    ELSE IF( DESCB( NB_ ).NE.DESCA( NB_ ) ) THEN
098                        INFO = - ( 1100 + NB_ )
099                    ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN
100                        INFO = - ( 1100 + CTXT_ )
101                    END IF
102                END IF
103                IDUM1( 1 ) = IBTYPE
104                IDUM2( 1 ) = 1
105                IF( UPPER ) THEN
106                    IDUM1( 2 ) = ICHAR( 'U' )
107                ELSE
108                    IDUM1( 2 ) = ICHAR( 'L' )
109                END IF
110                IDUM2( 2 ) = 2
111                CALL PCHK2MAT( N , 3 , N , 3 , IA , JA , DESCA , 7 , N , 3 , N , 3 , IB ,
112       $        JB , DESCB , 11 , 2 , IDUM1 , IDUM2 , INFO )
113            END IF
114  
115            IF( INFO.NE.0 ) THEN
116                CALL PXERBLA( ICTXT , 'PCHEGST' , - INFO )
117                RETURN
118            END IF
119  
120  *         Quick return if possible
121  
122            IF( N.EQ.0 )
123       $        RETURN
124  
125                IF( IBTYPE.EQ.1 ) THEN
126                    IF( UPPER ) THEN
127  
128  *                     Compute inv(U')*sub( A )*inv(U)
129  
130                        K = 1
131                        NB = DESCA( NB_ )
132                        KB = MIN( ICEIL( JA , NB )*NB , JA + N - 1 ) - JA + 1
133  
134     10 CONTINUE
135  
136  *     Update the upper triangle of A(ia + k - 1 : ia + n - 1 , ja + k - 1 : ja + n - 1)
137  
138        CALL PCHEGS2 ( IBTYPE , UPLO , KB , A , IA + K - 1 , JA + K - 1 , DESCA , B ,
139       $IB + K - 1 , IB + K - 1 , DESCB , INFO )
140        IF( K + KB.LE.N ) THEN
141            CALL PCTRSM( 'Left' , UPLO , 'Conjugate Transpose' ,
142       $    'Non - unit' , KB , N - K - KB + 1 , CONE , B , IB + K - 1 ,
143       $    JB + K - 1 , DESCB , A , IA + K - 1 , JA + K + KB - 1 , DESCA )
144            CALL PCHEMM( 'Left' , UPLO , KB , N - K - KB + 1 , - HALF , A ,
145       $    IA + K - 1 , JA + K - 1 , DESCA , B , IB + K - 1 , JB + K + KB - 1 ,
146       $    DESCB , CONE , A , IA + K - 1 , JA + K + KB - 1 , DESCA )
147            CALL PCHER2K( UPLO , 'Conjugate Transpose' , N - K - KB + 1 , KB ,
148       $    - CONE , A , IA + K - 1 , JA + K + KB - 1 , DESCA , B ,
149       $    IB + K - 1 , JB + K + KB - 1 , DESCB , ONE , A ,
150       $    IA + K + KB - 1 , JA + K + KB - 1 , DESCA )
151            CALL PCHEMM( 'Left' , UPLO , KB , N - K - KB + 1 , - HALF , A ,
152       $    IA + K - 1 , JA + K - 1 , DESCA , B , IB + K - 1 , JB + K + KB - 1 ,
153       $    DESCB , CONE , A , IA + K - 1 , JA + K + KB - 1 , DESCA )
154            CALL PCTRSM( 'Right' , UPLO , 'No transpose' , 'Non - unit' ,
155       $    KB , N - K - KB + 1 , CONE , B , IB + K + KB - 1 , JB + K + KB - 1 ,
156       $    DESCB , A , IA + K - 1 , JA + K + KB - 1 , DESCA )
157        END IF
158        K = K + KB
159        KB = MIN( N - K + 1 , NB )
160  
161        IF( K.LE.N )
162       $    GO TO 10
163  
164        ELSE
165  
166  *         Compute inv(L)*sub( A )*inv(L')
167  
168            K = 1
169            NB = DESCA( MB_ )
170            KB = MIN( ICEIL( IA , NB )*NB , IA + N - 1 ) - IA + 1
171  
172     20 CONTINUE
173  
174  *     Update the lower triangle of A(ia + k - 1 : ia + n - 1 , ja + k - 1 : ja + n - 1)
175  
176        CALL PCHEGS2 ( IBTYPE , UPLO , KB , A , IA + K - 1 , JA + K - 1 , DESCA , B ,
177       $IB + K - 1 , JB + K - 1 , DESCB , INFO )
178        IF( K + KB.LE.N ) THEN
179            CALL PCTRSM( 'Right' , UPLO , 'Conjugate transpose' ,
180       $    'Non - unit' , N - K - KB + 1 , KB , CONE , B , IB + K - 1 ,
181       $    JB + K - 1 , DESCB , A , IA + K + KB - 1 , JA + K - 1 , DESCA )
182            CALL PCHEMM( 'Right' , UPLO , N - K - KB + 1 , KB , - HALF , A ,
183       $    IA + K - 1 , JA + K - 1 , DESCA , B , IB + K + KB - 1 , JB + K - 1 ,
184       $    DESCB , CONE , A , IA + K + KB - 1 , JA + K - 1 , DESCA )
185            CALL PCHER2K( UPLO , 'No transpose' , N - K - KB + 1 , KB , - CONE ,
186       $    A , IA + K + KB - 1 , JA + K - 1 , DESCA , B , IB + K + KB - 1 ,
187       $    JB + K - 1 , DESCB , ONE , A , IA + K + KB - 1 ,
188       $    JA + K + KB - 1 , DESCA )
189            CALL PCHEMM( 'Right' , UPLO , N - K - KB + 1 , KB , - HALF , A ,
190       $    IA + K - 1 , JA + K - 1 , DESCA , B , IB + K + KB - 1 , JB + K - 1 ,
191       $    DESCB , CONE , A , IA + K + KB - 1 , JA + K - 1 , DESCA )
192            CALL PCTRSM( 'Left' , UPLO , 'No transpose' , 'Non - unit' ,
193       $    N - K - KB + 1 , KB , CONE , B , IB + K + KB - 1 , JB + K + KB - 1 ,
194       $    DESCB , A , IA + K + KB - 1 , JA + K - 1 , DESCA )
195        END IF
196        K = K + KB
197        KB = MIN( N - K + 1 , NB )
198  
199        IF( K.LE.N )
200       $    GO TO 20
201  
202        END IF
203  
204        ELSE
205  
206            IF( UPPER ) THEN
207  
208  *             Compute U*sub( A )*U'
209  
210                K = 1
211                NB = DESCA( NB_ )
212                KB = MIN( ICEIL( JA , NB )*NB , JA + N - 1 ) - JA + 1
213  
214     30 CONTINUE
215  
216  *     Update the upper triangle of A(ia : ia + k + kb - 2 , ja : ja + k + kb - 2)
217  
218        CALL PCTRMM( 'Left' , UPLO , 'No transpose' , 'Non - unit' , K - 1 ,
219       $KB , CONE , B , IB , JB , DESCB , A , IA , JA + K - 1 ,
220       $DESCA )
221        CALL PCHEMM( 'Right' , UPLO , K - 1 , KB , HALF , A , IA + K - 1 ,
222       $JA + K - 1 , DESCA , B , IB , JB + K - 1 , DESCB , CONE , A ,
223       $IA , JA + K - 1 , DESCA )
224        CALL PCHER2K( UPLO , 'No transpose' , K - 1 , KB , CONE , A , IA ,
225       $JA + K - 1 , DESCA , B , IB , JB + K - 1 , DESCB , ONE , A ,
226       $IA , JA , DESCA )
227        CALL PCHEMM( 'Right' , UPLO , K - 1 , KB , HALF , A , IA + K - 1 ,
228       $JA + K - 1 , DESCA , B , IB , JB + K - 1 , DESCB , CONE , A ,
229       $IA , JA + K - 1 , DESCA )
230        CALL PCTRMM( 'Right' , UPLO , 'Conjugate transpose' ,
231       $'Non - unit' , K - 1 , KB , CONE , B , IB + K - 1 , JB + K - 1 ,
232       $DESCB , A , IA , JA + K - 1 , DESCA )
233        CALL PCHEGS2 ( IBTYPE , UPLO , KB , A , IA + K - 1 , JA + K - 1 , DESCA , B ,
234       $IB + K - 1 , JB + K - 1 , DESCB , INFO )
235  
236        K = K + KB
237        KB = MIN( N - K + 1 , NB )
238  
239        IF( K.LE.N )
240       $    GO TO 30
241  
242        ELSE
243  
244  *         Compute L'*sub( A )*L
245  
246            K = 1
247            NB = DESCA( MB_ )
248            KB = MIN( ICEIL( IA , NB )*NB , IA + N - 1 ) - IA + 1
249  
250     40 CONTINUE
251  
252  *     Update the lower triangle of A(ia : ia + k + kb - 2 , ja : ja + k + kb - 2)
253  
254        CALL PCTRMM( 'Right' , UPLO , 'No transpose' , 'Non - unit' , KB ,
255       $K - 1 , CONE , B , IB , JB , DESCB , A , IA + K - 1 , JA ,
256       $DESCA )
257        CALL PCHEMM( 'Left' , UPLO , KB , K - 1 , HALF , A , IA + K - 1 , JA + K - 1 ,
258       $DESCA , B , IB + K - 1 , JB , DESCB , CONE , A , IA + K - 1 ,
259       $JA , DESCA )
260        CALL PCHER2K( UPLO , 'Conjugate transpose' , K - 1 , KB , CONE , A ,
261       $IA + K - 1 , JA , DESCA , B , IB + K - 1 , JB , DESCB , ONE ,
262       $A , IA , JA , DESCA )
263        CALL PCHEMM( 'Left' , UPLO , KB , K - 1 , HALF , A , IA + K - 1 , JA + K - 1 ,
264       $DESCA , B , IB + K - 1 , JB , DESCB , CONE , A , IA + K - 1 ,
265       $JA , DESCA )
266        CALL PCTRMM( 'Left' , UPLO , 'Conjugate transpose' ,
267       $'Non - unit' , KB , K - 1 , CONE , B , IB + K - 1 , JB + K - 1 ,
268       $DESCB , A , IA + K - 1 , JA , DESCA )
269        CALL PCHEGS2 ( IBTYPE , UPLO , KB , A , IA + K - 1 , JA + K - 1 , DESCA , B ,
270       $IB + K - 1 , JB + K - 1 , DESCB , INFO )
271  
272        K = K + KB
273        KB = MIN( N - K + 1 , NB )
274  
275        IF( K.LE.N )
276       $    GO TO 40
277  
278        END IF
279  
280        END IF
281  
282        RETURN
283  
284  *     End of PCHEGST
285  
286        END