Routine: PSSYGST()  File: SRC\pssygst.f

 
 
# lines: 438
  # code: 438
  # comment: 0
  # blank:0
# Variables:44
# Callers:1
# Callings:1
# Words:231
# Keywords:142
 

 

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