src/LAPACK/dspgst.f

   1       SUBROUTINE DSPGST( ITYPE, UPLO, N, AP, BP, INFO )
   2 *
   3 *  -- LAPACK routine (version 3.1) --
   4 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
   5 *     November 2006
   6 *
   7 *     .. Scalar Arguments ..
   8       CHARACTER          UPLO
   9       INTEGER            INFO, ITYPE, N
  10 *     ..
  11 *     .. Array Arguments ..
  12       DOUBLE PRECISION   AP( * ), BP( * )
  13 *     ..
  14 *
  15 *  Purpose
  16 *  =======
  17 *
  18 *  DSPGST reduces a real symmetric-definite generalized eigenproblem
  19 *  to standard form, using packed storage.
  20 *
  21 *  If ITYPE = 1, the problem is A*x = lambda*B*x,
  22 *  and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
  23 *
  24 *  If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
  25 *  B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
  26 *
  27 *  B must have been previously factorized as U**T*U or L*L**T by DPPTRF.
  28 *
  29 *  Arguments
  30 *  =========
  31 *
  32 *  ITYPE   (input) INTEGER
  33 *          = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
  34 *          = 2 or 3: compute U*A*U**T or L**T*A*L.
  35 *
  36 *  UPLO    (input) CHARACTER*1
  37 *          = 'U':  Upper triangle of A is stored and B is factored as
  38 *                  U**T*U;
  39 *          = 'L':  Lower triangle of A is stored and B is factored as
  40 *                  L*L**T.
  41 *
  42 *  N       (input) INTEGER
  43 *          The order of the matrices A and B.  N >= 0.
  44 *
  45 *  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
  46 *          On entry, the upper or lower triangle of the symmetric matrix
  47 *          A, packed columnwise in a linear array.  The j-th column of A
  48 *          is stored in the array AP as follows:
  49 *          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  50 *          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  51 *
  52 *          On exit, if INFO = 0, the transformed matrix, stored in the
  53 *          same format as A.
  54 *
  55 *  BP      (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
  56 *          The triangular factor from the Cholesky factorization of B,
  57 *          stored in the same format as A, as returned by DPPTRF.
  58 *
  59 *  INFO    (output) INTEGER
  60 *          = 0:  successful exit
  61 *          < 0:  if INFO = -i, the i-th argument had an illegal value
  62 *
  63 *  =====================================================================
  64 *
  65 *     .. Parameters ..
  66       DOUBLE PRECISION   ONE, HALF
  67       PARAMETER          ( ONE = 1.0D0, HALF = 0.5D0 )
  68 *     ..
  69 *     .. Local Scalars ..
  70       LOGICAL            UPPER
  71       INTEGER            J, J1, J1J1, JJ, K, K1, K1K1, KK
  72       DOUBLE PRECISION   AJJ, AKK, BJJ, BKK, CT
  73 *     ..
  74 *     .. External Subroutines ..
  75       EXTERNAL           DAXPY, DSCAL, DSPMV, DSPR2, DTPMV, DTPSV,
  76      $                   XERBLA
  77 *     ..
  78 *     .. External Functions ..
  79       LOGICAL            LSAME
  80       DOUBLE PRECISION   DDOT
  81       EXTERNAL           LSAME, DDOT
  82 *     ..
  83 *     .. Executable Statements ..
  84 *
  85 *     Test the input parameters.
  86 *
  87       INFO = 0
  88       UPPER = LSAME( UPLO, 'U' )
  89       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
  90          INFO = -1
  91       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  92          INFO = -2
  93       ELSE IF( N.LT.0 ) THEN
  94          INFO = -3
  95       END IF
  96       IF( INFO.NE.0 ) THEN
  97          CALL XERBLA( 'DSPGST', -INFO )
  98          RETURN
  99       END IF
 100 *
 101       IF( ITYPE.EQ.1 ) THEN
 102          IF( UPPER ) THEN
 103 *
 104 *           Compute inv(U')*A*inv(U)
 105 *
 106 *           J1 and JJ are the indices of A(1,j) and A(j,j)
 107 *
 108             JJ = 0
 109             DO 10 J = 1, N
 110                J1 = JJ + 1
 111                JJ = JJ + J
 112 *
 113 *              Compute the j-th column of the upper triangle of A
 114 *
 115                BJJ = BP( JJ )
 116                CALL DTPSV( UPLO, 'Transpose', 'Nonunit', J, BP,
 117      $                     AP( J1 ), 1 )
 118                CALL DSPMV( UPLO, J-1, -ONE, AP, BP( J1 ), 1, ONE,
 119      $                     AP( J1 ), 1 )
 120                CALL DSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
 121                AP( JJ ) = ( AP( JJ )-DDOT( J-1, AP( J1 ), 1, BP( J1 ),
 122      $                    1 ) ) / BJJ
 123    10       CONTINUE
 124          ELSE
 125 *
 126 *           Compute inv(L)*A*inv(L')
 127 *
 128 *           KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
 129 *
 130             KK = 1
 131             DO 20 K = 1, N
 132                K1K1 = KK + N - K + 1
 133 *
 134 *              Update the lower triangle of A(k:n,k:n)
 135 *
 136                AKK = AP( KK )
 137                BKK = BP( KK )
 138                AKK = AKK / BKK**2
 139                AP( KK ) = AKK
 140                IF( K.LT.N ) THEN
 141                   CALL DSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
 142                   CT = -HALF*AKK
 143                   CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
 144                   CALL DSPR2( UPLO, N-K, -ONE, AP( KK+1 ), 1,
 145      $                        BP( KK+1 ), 1, AP( K1K1 ) )
 146                   CALL DAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
 147                   CALL DTPSV( UPLO, 'No transpose', 'Non-unit', N-K,
 148      $                        BP( K1K1 ), AP( KK+1 ), 1 )
 149                END IF
 150                KK = K1K1
 151    20       CONTINUE
 152          END IF
 153       ELSE
 154          IF( UPPER ) THEN
 155 *
 156 *           Compute U*A*U'
 157 *
 158 *           K1 and KK are the indices of A(1,k) and A(k,k)
 159 *
 160             KK = 0
 161             DO 30 K = 1, N
 162                K1 = KK + 1
 163                KK = KK + K
 164 *
 165 *              Update the upper triangle of A(1:k,1:k)
 166 *
 167                AKK = AP( KK )
 168                BKK = BP( KK )
 169                CALL DTPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
 170      $                     AP( K1 ), 1 )
 171                CT = HALF*AKK
 172                CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
 173                CALL DSPR2( UPLO, K-1, ONE, AP( K1 ), 1, BP( K1 ), 1,
 174      $                     AP )
 175                CALL DAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
 176                CALL DSCAL( K-1, BKK, AP( K1 ), 1 )
 177                AP( KK ) = AKK*BKK**2
 178    30       CONTINUE
 179          ELSE
 180 *
 181 *           Compute L'*A*L
 182 *
 183 *           JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
 184 *
 185             JJ = 1
 186             DO 40 J = 1, N
 187                J1J1 = JJ + N - J + 1
 188 *
 189 *              Compute the j-th column of the lower triangle of A
 190 *
 191                AJJ = AP( JJ )
 192                BJJ = BP( JJ )
 193                AP( JJ ) = AJJ*BJJ + DDOT( N-J, AP( JJ+1 ), 1,
 194      $                    BP( JJ+1 ), 1 )
 195                CALL DSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
 196                CALL DSPMV( UPLO, N-J, ONE, AP( J1J1 ), BP( JJ+1 ), 1,
 197      $                     ONE, AP( JJ+1 ), 1 )
 198                CALL DTPMV( UPLO, 'Transpose', 'Non-unit', N-J+1,
 199      $                     BP( JJ ), AP( JJ ), 1 )
 200                JJ = J1J1
 201    40       CONTINUE
 202          END IF
 203       END IF
 204       RETURN
 205 *
 206 *     End of DSPGST
 207 *
 208       END