LAPACK: SRC/dspgst.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DSPGST( ITYPE, UPLO, N, AP, BP, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          UPLO
00010       INTEGER            INFO, ITYPE, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   AP( * ), BP( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DSPGST reduces a real symmetric-definite generalized eigenproblem
00020 *  to standard form, using packed storage.
00021 *
00022 *  If ITYPE = 1, the problem is A*x = lambda*B*x,
00023 *  and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
00024 *
00025 *  If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
00026 *  B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
00027 *
00028 *  B must have been previously factorized as U**T*U or L*L**T by DPPTRF.
00029 *
00030 *  Arguments
00031 *  =========
00032 *
00033 *  ITYPE   (input) INTEGER
00034 *          = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
00035 *          = 2 or 3: compute U*A*U**T or L**T*A*L.
00036 *
00037 *  UPLO    (input) CHARACTER*1
00038 *          = 'U':  Upper triangle of A is stored and B is factored as
00039 *                  U**T*U;
00040 *          = 'L':  Lower triangle of A is stored and B is factored as
00041 *                  L*L**T.
00042 *
00043 *  N       (input) INTEGER
00044 *          The order of the matrices A and B.  N >= 0.
00045 *
00046 *  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00047 *          On entry, the upper or lower triangle of the symmetric matrix
00048 *          A, packed columnwise in a linear array.  The j-th column of A
00049 *          is stored in the array AP as follows:
00050 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
00051 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
00052 *
00053 *          On exit, if INFO = 0, the transformed matrix, stored in the
00054 *          same format as A.
00055 *
00056 *  BP      (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
00057 *          The triangular factor from the Cholesky factorization of B,
00058 *          stored in the same format as A, as returned by DPPTRF.
00059 *
00060 *  INFO    (output) INTEGER
00061 *          = 0:  successful exit
00062 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00063 *
00064 *  =====================================================================
00065 *
00066 *     .. Parameters ..
00067       DOUBLE PRECISION   ONE, HALF
00068       PARAMETER          ( ONE = 1.0D0, HALF = 0.5D0 )
00069 *     ..
00070 *     .. Local Scalars ..
00071       LOGICAL            UPPER
00072       INTEGER            J, J1, J1J1, JJ, K, K1, K1K1, KK
00073       DOUBLE PRECISION   AJJ, AKK, BJJ, BKK, CT
00074 *     ..
00075 *     .. External Subroutines ..
00076       EXTERNAL           DAXPY, DSCAL, DSPMV, DSPR2, DTPMV, DTPSV,
00077      $                   XERBLA
00078 *     ..
00079 *     .. External Functions ..
00080       LOGICAL            LSAME
00081       DOUBLE PRECISION   DDOT
00082       EXTERNAL           LSAME, DDOT
00083 *     ..
00084 *     .. Executable Statements ..
00085 *
00086 *     Test the input parameters.
00087 *
00088       INFO = 0
00089       UPPER = LSAME( UPLO, 'U' )
00090       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
00091          INFO = -1
00092       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00093          INFO = -2
00094       ELSE IF( N.LT.0 ) THEN
00095          INFO = -3
00096       END IF
00097       IF( INFO.NE.0 ) THEN
00098          CALL XERBLA( 'DSPGST', -INFO )
00099          RETURN
00100       END IF
00101 *
00102       IF( ITYPE.EQ.1 ) THEN
00103          IF( UPPER ) THEN
00104 *
00105 *           Compute inv(U')*A*inv(U)
00106 *
00107 *           J1 and JJ are the indices of A(1,j) and A(j,j)
00108 *
00109             JJ = 0
00110             DO 10 J = 1, N
00111                J1 = JJ + 1
00112                JJ = JJ + J
00113 *
00114 *              Compute the j-th column of the upper triangle of A
00115 *
00116                BJJ = BP( JJ )
00117                CALL DTPSV( UPLO, 'Transpose', 'Nonunit', J, BP,
00118      $                     AP( J1 ), 1 )
00119                CALL DSPMV( UPLO, J-1, -ONE, AP, BP( J1 ), 1, ONE,
00120      $                     AP( J1 ), 1 )
00121                CALL DSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
00122                AP( JJ ) = ( AP( JJ )-DDOT( J-1, AP( J1 ), 1, BP( J1 ),
00123      $                    1 ) ) / BJJ
00124    10       CONTINUE
00125          ELSE
00126 *
00127 *           Compute inv(L)*A*inv(L')
00128 *
00129 *           KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
00130 *
00131             KK = 1
00132             DO 20 K = 1, N
00133                K1K1 = KK + N - K + 1
00134 *
00135 *              Update the lower triangle of A(k:n,k:n)
00136 *
00137                AKK = AP( KK )
00138                BKK = BP( KK )
00139                AKK = AKK / BKK**2
00140                AP( KK ) = AKK
00141                IF( K.LT.N ) THEN
00142                   CALL DSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
00143                   CT = -HALF*AKK
00144                   CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
00145                   CALL DSPR2( UPLO, N-K, -ONE, AP( KK+1 ), 1,
00146      $                        BP( KK+1 ), 1, AP( K1K1 ) )
00147                   CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
00148                   CALL DTPSV( UPLO, 'No transpose', 'Non-unit', N-K,
00149      $                        BP( K1K1 ), AP( KK+1 ), 1 )
00150                END IF
00151                KK = K1K1
00152    20       CONTINUE
00153          END IF
00154       ELSE
00155          IF( UPPER ) THEN
00156 *
00157 *           Compute U*A*U'
00158 *
00159 *           K1 and KK are the indices of A(1,k) and A(k,k)
00160 *
00161             KK = 0
00162             DO 30 K = 1, N
00163                K1 = KK + 1
00164                KK = KK + K
00165 *
00166 *              Update the upper triangle of A(1:k,1:k)
00167 *
00168                AKK = AP( KK )
00169                BKK = BP( KK )
00170                CALL DTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
00171      $                     AP( K1 ), 1 )
00172                CT = HALF*AKK
00173                CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
00174                CALL DSPR2( UPLO, K-1, ONE, AP( K1 ), 1, BP( K1 ), 1,
00175      $                     AP )
00176                CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
00177                CALL DSCAL( K-1, BKK, AP( K1 ), 1 )
00178                AP( KK ) = AKK*BKK**2
00179    30       CONTINUE
00180          ELSE
00181 *
00182 *           Compute L'*A*L
00183 *
00184 *           JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
00185 *
00186             JJ = 1
00187             DO 40 J = 1, N
00188                J1J1 = JJ + N - J + 1
00189 *
00190 *              Compute the j-th column of the lower triangle of A
00191 *
00192                AJJ = AP( JJ )
00193                BJJ = BP( JJ )
00194                AP( JJ ) = AJJ*BJJ + DDOT( N-J, AP( JJ+1 ), 1,
00195      $                    BP( JJ+1 ), 1 )
00196                CALL DSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
00197                CALL DSPMV( UPLO, N-J, ONE, AP( J1J1 ), BP( JJ+1 ), 1,
00198      $                     ONE, AP( JJ+1 ), 1 )
00199                CALL DTPMV( UPLO, 'Transpose', 'Non-unit', N-J+1,
00200      $                     BP( JJ ), AP( JJ ), 1 )
00201                JJ = J1J1
00202    40       CONTINUE
00203          END IF
00204       END IF
00205       RETURN
00206 *
00207 *     End of DSPGST
00208 *
00209       END