ScaLAPACK  2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
pdsygst.f
Go to the documentation of this file.
00001 *
00002 *
00003       SUBROUTINE PDSYGST( IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB,
00004      $                    DESCB, SCALE, INFO )
00005 *
00006 *  -- ScaLAPACK routine (version 1.7) --
00007 *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
00008 *     and University of California, Berkeley.
00009 *     May 1, 1997
00010 *
00011 *     .. Scalar Arguments ..
00012       CHARACTER          UPLO
00013       INTEGER            IA, IB, IBTYPE, INFO, JA, JB, N
00014       DOUBLE PRECISION   SCALE
00015 *     ..
00016 *     .. Array Arguments ..
00017       INTEGER            DESCA( * ), DESCB( * )
00018       DOUBLE PRECISION   A( * ), B( * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  PDSYGST reduces a real symmetric-definite generalized eigenproblem
00025 *  to standard form.
00026 *
00027 *  In the following sub( A ) denotes A( IA:IA+N-1, JA:JA+N-1 ) and
00028 *  sub( B ) denotes B( IB:IB+N-1, JB:JB+N-1 ).
00029 *
00030 *  If IBTYPE = 1, the problem is sub( A )*x = lambda*sub( B )*x,
00031 *  and sub( A ) is overwritten by inv(U**T)*sub( A )*inv(U) or
00032 *  inv(L)*sub( A )*inv(L**T)
00033 *
00034 *  If IBTYPE = 2 or 3, the problem is sub( A )*sub( B )*x = lambda*x or
00035 *  sub( B )*sub( A )*x = lambda*x, and sub( A ) is overwritten by
00036 *  U*sub( A )*U**T or L**T*sub( A )*L.
00037 *
00038 *  sub( B ) must have been previously factorized as U**T*U or L*L**T by
00039 *  PDPOTRF.
00040 *
00041 *  Notes
00042 *  =====
00043 *
00044 *  Each global data object is described by an associated description
00045 *  vector.  This vector stores the information required to establish
00046 *  the mapping between an object element and its corresponding process
00047 *  and memory location.
00048 *
00049 *  Let A be a generic term for any 2D block cyclicly distributed array.
00050 *  Such a global array has an associated description vector DESCA.
00051 *  In the following comments, the character _ should be read as
00052 *  "of the global array".
00053 *
00054 *  NOTATION        STORED IN      EXPLANATION
00055 *  --------------- -------------- --------------------------------------
00056 *  DTYPE_A(global) DESCA( DTYPE_ )The descriptor type.  In this case,
00057 *                                 DTYPE_A = 1.
00058 *  CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
00059 *                                 the BLACS process grid A is distribu-
00060 *                                 ted over. The context itself is glo-
00061 *                                 bal, but the handle (the integer
00062 *                                 value) may vary.
00063 *  M_A    (global) DESCA( M_ )    The number of rows in the global
00064 *                                 array A.
00065 *  N_A    (global) DESCA( N_ )    The number of columns in the global
00066 *                                 array A.
00067 *  MB_A   (global) DESCA( MB_ )   The blocking factor used to distribute
00068 *                                 the rows of the array.
00069 *  NB_A   (global) DESCA( NB_ )   The blocking factor used to distribute
00070 *                                 the columns of the array.
00071 *  RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
00072 *                                 row of the array A is distributed.
00073 *  CSRC_A (global) DESCA( CSRC_ ) The process column over which the
00074 *                                 first column of the array A is
00075 *                                 distributed.
00076 *  LLD_A  (local)  DESCA( LLD_ )  The leading dimension of the local
00077 *                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
00078 *
00079 *  Let K be the number of rows or columns of a distributed matrix,
00080 *  and assume that its process grid has dimension p x q.
00081 *  LOCr( K ) denotes the number of elements of K that a process
00082 *  would receive if K were distributed over the p processes of its
00083 *  process column.
00084 *  Similarly, LOCc( K ) denotes the number of elements of K that a
00085 *  process would receive if K were distributed over the q processes of
00086 *  its process row.
00087 *  The values of LOCr() and LOCc() may be determined via a call to the
00088 *  ScaLAPACK tool function, NUMROC:
00089 *          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
00090 *          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
00091 *  An upper bound for these quantities may be computed by:
00092 *          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
00093 *          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
00094 *
00095 *  Arguments
00096 *  =========
00097 *
00098 *  IBTYPE   (global input) INTEGER
00099 *          = 1: compute inv(U**T)*sub( A )*inv(U) or
00100 *               inv(L)*sub( A )*inv(L**T);
00101 *          = 2 or 3: compute U*sub( A )*U**T or L**T*sub( A )*L.
00102 *
00103 *  UPLO    (global input) CHARACTER
00104 *          = 'U':  Upper triangle of sub( A ) is stored and sub( B ) is
00105 *                  factored as U**T*U;
00106 *          = 'L':  Lower triangle of sub( A ) is stored and sub( B ) is
00107 *                  factored as L*L**T.
00108 *
00109 *  N       (global input) INTEGER
00110 *          The order of the matrices sub( A ) and sub( B ).  N >= 0.
00111 *
00112 *  A       (local input/local output) DOUBLE PRECISION pointer into the
00113 *          local memory to an array of dimension (LLD_A, LOCc(JA+N-1)).
00114 *          On entry, this array contains the local pieces of the
00115 *          N-by-N symmetric distributed matrix sub( A ). If UPLO = 'U',
00116 *          the leading N-by-N upper triangular part of sub( A ) contains
00117 *          the upper triangular part of the matrix, and its strictly
00118 *          lower triangular part is not referenced.  If UPLO = 'L', the
00119 *          leading N-by-N lower triangular part of sub( A ) contains
00120 *          the lower triangular part of the matrix, and its strictly
00121 *          upper triangular part is not referenced.
00122 *
00123 *          On exit, if INFO = 0, the transformed matrix, stored in the
00124 *          same format as sub( A ).
00125 *
00126 *  IA      (global input) INTEGER
00127 *          A's global row index, which points to the beginning of the
00128 *          submatrix which is to be operated on.
00129 *
00130 *  JA      (global input) INTEGER
00131 *          A's global column index, which points to the beginning of
00132 *          the submatrix which is to be operated on.
00133 *
00134 *  DESCA   (global and local input) INTEGER array of dimension DLEN_.
00135 *          The array descriptor for the distributed matrix A.
00136 *
00137 *  B       (local input) DOUBLE PRECISION pointer into the local memory
00138 *          to an array of dimension (LLD_B, LOCc(JB+N-1)). On entry,
00139 *          this array contains the local pieces of the triangular factor
00140 *          from the Cholesky factorization of sub( B ), as returned by
00141 *          PDPOTRF.
00142 *
00143 *  IB      (global input) INTEGER
00144 *          B's global row index, which points to the beginning of the
00145 *          submatrix which is to be operated on.
00146 *
00147 *  JB      (global input) INTEGER
00148 *          B's global column index, which points to the beginning of
00149 *          the submatrix which is to be operated on.
00150 *
00151 *  DESCB   (global and local input) INTEGER array of dimension DLEN_.
00152 *          The array descriptor for the distributed matrix B.
00153 *
00154 *  SCALE   (global output) DOUBLE PRECISION
00155 *          Amount by which the eigenvalues should be scaled to
00156 *          compensate for the scaling performed in this routine.
00157 *          At present, SCALE is always returned as 1.0, it is
00158 *          returned here to allow for future enhancement.
00159 *
00160 *  INFO    (global output) INTEGER
00161 *          = 0:  successful exit
00162 *          < 0:  If the i-th argument is an array and the j-entry had
00163 *                an illegal value, then INFO = -(i*100+j), if the i-th
00164 *                argument is a scalar and had an illegal value, then
00165 *                INFO = -i.
00166 *
00167 *  =====================================================================
00168 *
00169 *     .. Parameters ..
00170       INTEGER            BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
00171      $                   MB_, NB_, RSRC_, CSRC_, LLD_
00172       PARAMETER          ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1,
00173      $                   CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6,
00174      $                   RSRC_ = 7, CSRC_ = 8, LLD_ = 9 )
00175       DOUBLE PRECISION   ONE, HALF
00176       PARAMETER          ( ONE = 1.0D+0, HALF = 0.5D+0 )
00177 *     ..
00178 *     .. Local Scalars ..
00179       LOGICAL            UPPER
00180       INTEGER            IACOL, IAROW, IBCOL, IBROW, ICOFFA, ICOFFB,
00181      $                   ICTXT, IROFFA, IROFFB, K, KB, MYCOL, MYROW, NB,
00182      $                   NPCOL, NPROW
00183 *     ..
00184 *     .. Local Arrays ..
00185       INTEGER            IDUM1( 2 ), IDUM2( 2 )
00186 *     ..
00187 *     .. External Subroutines ..
00188       EXTERNAL           BLACS_GRIDINFO, CHK1MAT, PCHK2MAT, PDSYGS2,
00189      $                   PDSYMM, PDSYR2K, PDTRMM, PDTRSM, PXERBLA
00190 *     ..
00191 *     .. Intrinsic Functions ..
00192       INTRINSIC          ICHAR, MIN, MOD
00193 *     ..
00194 *     .. External Functions ..
00195       LOGICAL            LSAME
00196       INTEGER            ICEIL, INDXG2P
00197       EXTERNAL           LSAME, ICEIL, INDXG2P
00198 *     ..
00199 *     .. Executable Statements ..
00200 *       This is just to keep ftnchek happy
00201       IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_*
00202      $    RSRC_.LT.0 )RETURN
00203 *
00204 *     Get grid parameters
00205 *
00206       SCALE = ONE
00207 *
00208       ICTXT = DESCA( CTXT_ )
00209       CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
00210 *
00211 *     Test the input parameters
00212 *
00213       INFO = 0
00214       IF( NPROW.EQ.-1 ) THEN
00215          INFO = -( 700+CTXT_ )
00216       ELSE
00217          UPPER = LSAME( UPLO, 'U' )
00218          CALL CHK1MAT( N, 3, N, 3, IA, JA, DESCA, 7, INFO )
00219          CALL CHK1MAT( N, 3, N, 3, IB, JB, DESCB, 11, INFO )
00220          IF( INFO.EQ.0 ) THEN
00221             IAROW = INDXG2P( IA, DESCA( MB_ ), MYROW, DESCA( RSRC_ ),
00222      $              NPROW )
00223             IBROW = INDXG2P( IB, DESCB( MB_ ), MYROW, DESCB( RSRC_ ),
00224      $              NPROW )
00225             IACOL = INDXG2P( JA, DESCA( NB_ ), MYCOL, DESCA( CSRC_ ),
00226      $              NPCOL )
00227             IBCOL = INDXG2P( JB, DESCB( NB_ ), MYCOL, DESCB( CSRC_ ),
00228      $              NPCOL )
00229             IROFFA = MOD( IA-1, DESCA( MB_ ) )
00230             ICOFFA = MOD( JA-1, DESCA( NB_ ) )
00231             IROFFB = MOD( IB-1, DESCB( MB_ ) )
00232             ICOFFB = MOD( JB-1, DESCB( NB_ ) )
00233             IF( IBTYPE.LT.1 .OR. IBTYPE.GT.3 ) THEN
00234                INFO = -1
00235             ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00236                INFO = -2
00237             ELSE IF( N.LT.0 ) THEN
00238                INFO = -3
00239             ELSE IF( IROFFA.NE.0 ) THEN
00240                INFO = -5
00241             ELSE IF( ICOFFA.NE.0 ) THEN
00242                INFO = -6
00243             ELSE IF( DESCA( MB_ ).NE.DESCA( NB_ ) ) THEN
00244                INFO = -( 700+NB_ )
00245             ELSE IF( IROFFB.NE.0 .OR. IBROW.NE.IAROW ) THEN
00246                INFO = -9
00247             ELSE IF( ICOFFB.NE.0 .OR. IBCOL.NE.IACOL ) THEN
00248                INFO = -10
00249             ELSE IF( DESCB( MB_ ).NE.DESCA( MB_ ) ) THEN
00250                INFO = -( 1100+MB_ )
00251             ELSE IF( DESCB( NB_ ).NE.DESCA( NB_ ) ) THEN
00252                INFO = -( 1100+NB_ )
00253             ELSE IF( ICTXT.NE.DESCB( CTXT_ ) ) THEN
00254                INFO = -( 1100+CTXT_ )
00255             END IF
00256          END IF
00257          IDUM1( 1 ) = IBTYPE
00258          IDUM2( 1 ) = 1
00259          IF( UPPER ) THEN
00260             IDUM1( 2 ) = ICHAR( 'U' )
00261          ELSE
00262             IDUM1( 2 ) = ICHAR( 'L' )
00263          END IF
00264          IDUM2( 2 ) = 2
00265          CALL PCHK2MAT( N, 3, N, 3, IA, JA, DESCA, 7, N, 3, N, 3, IB,
00266      $                  JB, DESCB, 11, 2, IDUM1, IDUM2, INFO )
00267       END IF
00268 *
00269       IF( INFO.NE.0 ) THEN
00270          CALL PXERBLA( ICTXT, 'PDSYGST', -INFO )
00271          RETURN
00272       END IF
00273 *
00274 *     Quick return if possible
00275 *
00276       IF( N.EQ.0 )
00277      $   RETURN
00278 *
00279       IF( IBTYPE.EQ.1 ) THEN
00280          IF( UPPER ) THEN
00281 *
00282 *           Compute inv(U')*sub( A )*inv(U)
00283 *
00284             K = 1
00285             NB = DESCA( NB_ )
00286             KB = MIN( ICEIL( JA, NB )*NB, JA+N-1 ) - JA + 1
00287 *
00288    10       CONTINUE
00289 *
00290 *           Update the upper triangle of A(ia+k-1:ia+n-1,ja+k-1:ja+n-1)
00291 *
00292             CALL PDSYGS2( IBTYPE, UPLO, KB, A, IA+K-1, JA+K-1, DESCA, B,
00293      $                    IB+K-1, IB+K-1, DESCB, INFO )
00294             IF( K+KB.LE.N ) THEN
00295                CALL PDTRSM( 'Left', UPLO, 'Transpose', 'Non-unit', KB,
00296      $                      N-K-KB+1, ONE, B, IB+K-1, JB+K-1, DESCB, A,
00297      $                      IA+K-1, JA+K+KB-1, DESCA )
00298                CALL PDSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF, A,
00299      $                      IA+K-1, JA+K-1, DESCA, B, IB+K-1, JB+K+KB-1,
00300      $                      DESCB, ONE, A, IA+K-1, JA+K+KB-1, DESCA )
00301                CALL PDSYR2K( UPLO, 'Transpose', N-K-KB+1, KB, -ONE, A,
00302      $                       IA+K-1, JA+K+KB-1, DESCA, B, IB+K-1,
00303      $                       JB+K+KB-1, DESCB, ONE, A, IA+K+KB-1,
00304      $                       JA+K+KB-1, DESCA )
00305                CALL PDSYMM( 'Left', UPLO, KB, N-K-KB+1, -HALF, A,
00306      $                      IA+K-1, JA+K-1, DESCA, B, IB+K-1, JB+K+KB-1,
00307      $                      DESCB, ONE, A, IA+K-1, JA+K+KB-1, DESCA )
00308                CALL PDTRSM( 'Right', UPLO, 'No transpose', 'Non-unit',
00309      $                      KB, N-K-KB+1, ONE, B, IB+K+KB-1, JB+K+KB-1,
00310      $                      DESCB, A, IA+K-1, JA+K+KB-1, DESCA )
00311             END IF
00312             K = K + KB
00313             KB = MIN( N-K+1, NB )
00314 *
00315             IF( K.LE.N )
00316      $         GO TO 10
00317 *
00318          ELSE
00319 *
00320 *           Compute inv(L)*sub( A )*inv(L')
00321 *
00322             K = 1
00323             NB = DESCA( MB_ )
00324             KB = MIN( ICEIL( IA, NB )*NB, IA+N-1 ) - IA + 1
00325 *
00326    20       CONTINUE
00327 *
00328 *           Update the lower triangle of A(ia+k-1:ia+n-1,ja+k-1:ja+n-1)
00329 *
00330             CALL PDSYGS2( IBTYPE, UPLO, KB, A, IA+K-1, JA+K-1, DESCA, B,
00331      $                    IB+K-1, JB+K-1, DESCB, INFO )
00332             IF( K+KB.LE.N ) THEN
00333                CALL PDTRSM( 'Right', UPLO, 'Transpose', 'Non-unit',
00334      $                      N-K-KB+1, KB, ONE, B, IB+K-1, JB+K-1, DESCB,
00335      $                      A, IA+K+KB-1, JA+K-1, DESCA )
00336                CALL PDSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF, A,
00337      $                      IA+K-1, JA+K-1, DESCA, B, IB+K+KB-1, JB+K-1,
00338      $                      DESCB, ONE, A, IA+K+KB-1, JA+K-1, DESCA )
00339                CALL PDSYR2K( UPLO, 'No transpose', N-K-KB+1, KB, -ONE,
00340      $                       A, IA+K+KB-1, JA+K-1, DESCA, B, IB+K+KB-1,
00341      $                       JB+K-1, DESCB, ONE, A, IA+K+KB-1,
00342      $                       JA+K+KB-1, DESCA )
00343                CALL PDSYMM( 'Right', UPLO, N-K-KB+1, KB, -HALF, A,
00344      $                      IA+K-1, JA+K-1, DESCA, B, IB+K+KB-1, JB+K-1,
00345      $                      DESCB, ONE, A, IA+K+KB-1, JA+K-1, DESCA )
00346                CALL PDTRSM( 'Left', UPLO, 'No transpose', 'Non-unit',
00347      $                      N-K-KB+1, KB, ONE, B, IB+K+KB-1, JB+K+KB-1,
00348      $                      DESCB, A, IA+K+KB-1, JA+K-1, DESCA )
00349             END IF
00350             K = K + KB
00351             KB = MIN( N-K+1, NB )
00352 *
00353             IF( K.LE.N )
00354      $         GO TO 20
00355 *
00356          END IF
00357 *
00358       ELSE
00359 *
00360          IF( UPPER ) THEN
00361 *
00362 *           Compute U*sub( A )*U'
00363 *
00364             K = 1
00365             NB = DESCA( NB_ )
00366             KB = MIN( ICEIL( JA, NB )*NB, JA+N-1 ) - JA + 1
00367 *
00368    30       CONTINUE
00369 *
00370 *           Update the upper triangle of A(ia:ia+k+kb-2,ja:ja+k+kb-2)
00371 *
00372             CALL PDTRMM( 'Left', UPLO, 'No transpose', 'Non-unit', K-1,
00373      $                   KB, ONE, B, IB, JB, DESCB, A, IA, JA+K-1,
00374      $                   DESCA )
00375             CALL PDSYMM( 'Right', UPLO, K-1, KB, HALF, A, IA+K-1,
00376      $                   JA+K-1, DESCA, B, IB, JB+K-1, DESCB, ONE, A,
00377      $                   IA, JA+K-1, DESCA )
00378             CALL PDSYR2K( UPLO, 'No transpose', K-1, KB, ONE, A, IA,
00379      $                    JA+K-1, DESCA, B, IB, JB+K-1, DESCB, ONE, A,
00380      $                    IA, JA, DESCA )
00381             CALL PDSYMM( 'Right', UPLO, K-1, KB, HALF, A, IA+K-1,
00382      $                   JA+K-1, DESCA, B, IB, JB+K-1, DESCB, ONE, A,
00383      $                   IA, JA+K-1, DESCA )
00384             CALL PDTRMM( 'Right', UPLO, 'Transpose', 'Non-unit', K-1,
00385      $                   KB, ONE, B, IB+K-1, JB+K-1, DESCB, A, IA,
00386      $                   JA+K-1, DESCA )
00387             CALL PDSYGS2( IBTYPE, UPLO, KB, A, IA+K-1, JA+K-1, DESCA, B,
00388      $                    IB+K-1, JB+K-1, DESCB, INFO )
00389 *
00390             K = K + KB
00391             KB = MIN( N-K+1, NB )
00392 *
00393             IF( K.LE.N )
00394      $         GO TO 30
00395 *
00396          ELSE
00397 *
00398 *           Compute L'*sub( A )*L
00399 *
00400             K = 1
00401             NB = DESCA( MB_ )
00402             KB = MIN( ICEIL( IA, NB )*NB, IA+N-1 ) - IA + 1
00403 *
00404    40       CONTINUE
00405 *
00406 *           Update the lower triangle of A(ia:ia+k+kb-2,ja:ja+k+kb-2)
00407 *
00408             CALL PDTRMM( 'Right', UPLO, 'No transpose', 'Non-unit', KB,
00409      $                   K-1, ONE, B, IB, JB, DESCB, A, IA+K-1, JA,
00410      $                   DESCA )
00411             CALL PDSYMM( 'Left', UPLO, KB, K-1, HALF, A, IA+K-1, JA+K-1,
00412      $                   DESCA, B, IB+K-1, JB, DESCB, ONE, A, IA+K-1,
00413      $                   JA, DESCA )
00414             CALL PDSYR2K( UPLO, 'Transpose', K-1, KB, ONE, A, IA+K-1,
00415      $                    JA, DESCA, B, IB+K-1, JB, DESCB, ONE, A, IA,
00416      $                    JA, DESCA )
00417             CALL PDSYMM( 'Left', UPLO, KB, K-1, HALF, A, IA+K-1, JA+K-1,
00418      $                   DESCA, B, IB+K-1, JB, DESCB, ONE, A, IA+K-1,
00419      $                   JA, DESCA )
00420             CALL PDTRMM( 'Left', UPLO, 'Transpose', 'Non-unit', KB, K-1,
00421      $                   ONE, B, IB+K-1, JB+K-1, DESCB, A, IA+K-1, JA,
00422      $                   DESCA )
00423             CALL PDSYGS2( IBTYPE, UPLO, KB, A, IA+K-1, JA+K-1, DESCA, B,
00424      $                    IB+K-1, JB+K-1, DESCB, INFO )
00425 *
00426             K = K + KB
00427             KB = MIN( N-K+1, NB )
00428 *
00429             IF( K.LE.N )
00430      $         GO TO 40
00431 *
00432          END IF
00433 *
00434       END IF
00435 *
00436       RETURN
00437 *
00438 *     End of PDSYGST
00439 *
00440       END