SRC/sspgst.f

   1 *> \brief \b SSPGST
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SSPGST + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sspgst.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sspgst.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sspgst.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SSPGST( ITYPE, UPLO, N, AP, BP, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, ITYPE, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       REAL               AP( * ), BP( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> SSPGST reduces a real symmetric-definite generalized eigenproblem
  38 *> to standard form, using packed storage.
  39 *>
  40 *> If ITYPE = 1, the problem is A*x = lambda*B*x,
  41 *> and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T)
  42 *>
  43 *> If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or
  44 *> B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L.
  45 *>
  46 *> B must have been previously factorized as U**T*U or L*L**T by SPPTRF.
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] ITYPE
  53 *> \verbatim
  54 *>          ITYPE is INTEGER
  55 *>          = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T);
  56 *>          = 2 or 3: compute U*A*U**T or L**T*A*L.
  57 *> \endverbatim
  58 *>
  59 *> \param[in] UPLO
  60 *> \verbatim
  61 *>          UPLO is CHARACTER*1
  62 *>          = 'U':  Upper triangle of A is stored and B is factored as
  63 *>                  U**T*U;
  64 *>          = 'L':  Lower triangle of A is stored and B is factored as
  65 *>                  L*L**T.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] N
  69 *> \verbatim
  70 *>          N is INTEGER
  71 *>          The order of the matrices A and B.  N >= 0.
  72 *> \endverbatim
  73 *>
  74 *> \param[in,out] AP
  75 *> \verbatim
  76 *>          AP is REAL array, dimension (N*(N+1)/2)
  77 *>          On entry, the upper or lower triangle of the symmetric matrix
  78 *>          A, packed columnwise in a linear array.  The j-th column of A
  79 *>          is stored in the array AP as follows:
  80 *>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
  81 *>          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
  82 *>
  83 *>          On exit, if INFO = 0, the transformed matrix, stored in the
  84 *>          same format as A.
  85 *> \endverbatim
  86 *>
  87 *> \param[in] BP
  88 *> \verbatim
  89 *>          BP is REAL array, dimension (N*(N+1)/2)
  90 *>          The triangular factor from the Cholesky factorization of B,
  91 *>          stored in the same format as A, as returned by SPPTRF.
  92 *> \endverbatim
  93 *>
  94 *> \param[out] INFO
  95 *> \verbatim
  96 *>          INFO is INTEGER
  97 *>          = 0:  successful exit
  98 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  99 *> \endverbatim
 100 *
 101 *  Authors:
 102 *  ========
 103 *
 104 *> \author Univ. of Tennessee
 105 *> \author Univ. of California Berkeley
 106 *> \author Univ. of Colorado Denver
 107 *> \author NAG Ltd.
 108 *
 109 *> \date November 2011
 110 *
 111 *> \ingroup realOTHERcomputational
 112 *
 113 *  =====================================================================
 114       SUBROUTINE SSPGST( ITYPE, UPLO, N, AP, BP, INFO )
 115 *
 116 *  -- LAPACK computational routine (version 3.4.0) --
 117 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 118 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 119 *     November 2011
 120 *
 121 *     .. Scalar Arguments ..
 122       CHARACTER          UPLO
 123       INTEGER            INFO, ITYPE, N
 124 *     ..
 125 *     .. Array Arguments ..
 126       REAL               AP( * ), BP( * )
 127 *     ..
 128 *
 129 *  =====================================================================
 130 *
 131 *     .. Parameters ..
 132       REAL               ONE, HALF
 133       PARAMETER          ( ONE = 1.0, HALF = 0.5 )
 134 *     ..
 135 *     .. Local Scalars ..
 136       LOGICAL            UPPER
 137       INTEGER            J, J1, J1J1, JJ, K, K1, K1K1, KK
 138       REAL               AJJ, AKK, BJJ, BKK, CT
 139 *     ..
 140 *     .. External Subroutines ..
 141       EXTERNAL           SAXPY, SSCAL, SSPMV, SSPR2, STPMV, STPSV,
 142      $                   XERBLA
 143 *     ..
 144 *     .. External Functions ..
 145       LOGICAL            LSAME
 146       REAL               SDOT
 147       EXTERNAL           LSAME, SDOT
 148 *     ..
 149 *     .. Executable Statements ..
 150 *
 151 *     Test the input parameters.
 152 *
 153       INFO = 0
 154       UPPER = LSAME( UPLO, 'U' )
 155       IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
 156          INFO = -1
 157       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 158          INFO = -2
 159       ELSE IF( N.LT.0 ) THEN
 160          INFO = -3
 161       END IF
 162       IF( INFO.NE.0 ) THEN
 163          CALL XERBLA( 'SSPGST', -INFO )
 164          RETURN
 165       END IF
 166 *
 167       IF( ITYPE.EQ.1 ) THEN
 168          IF( UPPER ) THEN
 169 *
 170 *           Compute inv(U**T)*A*inv(U)
 171 *
 172 *           J1 and JJ are the indices of A(1,j) and A(j,j)
 173 *
 174             JJ = 0
 175             DO 10 J = 1, N
 176                J1 = JJ + 1
 177                JJ = JJ + J
 178 *
 179 *              Compute the j-th column of the upper triangle of A
 180 *
 181                BJJ = BP( JJ )
 182                CALL STPSV( UPLO, 'Transpose', 'Nonunit', J, BP,
 183      $                     AP( J1 ), 1 )
 184                CALL SSPMV( UPLO, J-1, -ONE, AP, BP( J1 ), 1, ONE,
 185      $                     AP( J1 ), 1 )
 186                CALL SSCAL( J-1, ONE / BJJ, AP( J1 ), 1 )
 187                AP( JJ ) = ( AP( JJ )-SDOT( J-1, AP( J1 ), 1, BP( J1 ),
 188      $                    1 ) ) / BJJ
 189    10       CONTINUE
 190          ELSE
 191 *
 192 *           Compute inv(L)*A*inv(L**T)
 193 *
 194 *           KK and K1K1 are the indices of A(k,k) and A(k+1,k+1)
 195 *
 196             KK = 1
 197             DO 20 K = 1, N
 198                K1K1 = KK + N - K + 1
 199 *
 200 *              Update the lower triangle of A(k:n,k:n)
 201 *
 202                AKK = AP( KK )
 203                BKK = BP( KK )
 204                AKK = AKK / BKK**2
 205                AP( KK ) = AKK
 206                IF( K.LT.N ) THEN
 207                   CALL SSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 )
 208                   CT = -HALF*AKK
 209                   CALL SAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
 210                   CALL SSPR2( UPLO, N-K, -ONE, AP( KK+1 ), 1,
 211      $                        BP( KK+1 ), 1, AP( K1K1 ) )
 212                   CALL SAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 )
 213                   CALL STPSV( UPLO, 'No transpose', 'Non-unit', N-K,
 214      $                        BP( K1K1 ), AP( KK+1 ), 1 )
 215                END IF
 216                KK = K1K1
 217    20       CONTINUE
 218          END IF
 219       ELSE
 220          IF( UPPER ) THEN
 221 *
 222 *           Compute U*A*U**T
 223 *
 224 *           K1 and KK are the indices of A(1,k) and A(k,k)
 225 *
 226             KK = 0
 227             DO 30 K = 1, N
 228                K1 = KK + 1
 229                KK = KK + K
 230 *
 231 *              Update the upper triangle of A(1:k,1:k)
 232 *
 233                AKK = AP( KK )
 234                BKK = BP( KK )
 235                CALL STPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP,
 236      $                     AP( K1 ), 1 )
 237                CT = HALF*AKK
 238                CALL SAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
 239                CALL SSPR2( UPLO, K-1, ONE, AP( K1 ), 1, BP( K1 ), 1,
 240      $                     AP )
 241                CALL SAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 )
 242                CALL SSCAL( K-1, BKK, AP( K1 ), 1 )
 243                AP( KK ) = AKK*BKK**2
 244    30       CONTINUE
 245          ELSE
 246 *
 247 *           Compute L**T *A*L
 248 *
 249 *           JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1)
 250 *
 251             JJ = 1
 252             DO 40 J = 1, N
 253                J1J1 = JJ + N - J + 1
 254 *
 255 *              Compute the j-th column of the lower triangle of A
 256 *
 257                AJJ = AP( JJ )
 258                BJJ = BP( JJ )
 259                AP( JJ ) = AJJ*BJJ + SDOT( N-J, AP( JJ+1 ), 1,
 260      $                    BP( JJ+1 ), 1 )
 261                CALL SSCAL( N-J, BJJ, AP( JJ+1 ), 1 )
 262                CALL SSPMV( UPLO, N-J, ONE, AP( J1J1 ), BP( JJ+1 ), 1,
 263      $                     ONE, AP( JJ+1 ), 1 )
 264                CALL STPMV( UPLO, 'Transpose', 'Non-unit', N-J+1,
 265      $                     BP( JJ ), AP( JJ ), 1 )
 266                JJ = J1J1
 267    40       CONTINUE
 268          END IF
 269       END IF
 270       RETURN
 271 *
 272 *     End of SSPGST
 273 *
 274       END