![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SSPGST 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SSPGST + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sspgst.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sspgst.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sspgst.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SSPGST( ITYPE, UPLO, N, AP, BP, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, ITYPE, N 00026 * .. 00027 * .. Array Arguments .. 00028 * REAL AP( * ), BP( * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> SSPGST reduces a real symmetric-definite generalized eigenproblem 00038 *> to standard form, using packed storage. 00039 *> 00040 *> If ITYPE = 1, the problem is A*x = lambda*B*x, 00041 *> and A is overwritten by inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T) 00042 *> 00043 *> If ITYPE = 2 or 3, the problem is A*B*x = lambda*x or 00044 *> B*A*x = lambda*x, and A is overwritten by U*A*U**T or L**T*A*L. 00045 *> 00046 *> B must have been previously factorized as U**T*U or L*L**T by SPPTRF. 00047 *> \endverbatim 00048 * 00049 * Arguments: 00050 * ========== 00051 * 00052 *> \param[in] ITYPE 00053 *> \verbatim 00054 *> ITYPE is INTEGER 00055 *> = 1: compute inv(U**T)*A*inv(U) or inv(L)*A*inv(L**T); 00056 *> = 2 or 3: compute U*A*U**T or L**T*A*L. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] UPLO 00060 *> \verbatim 00061 *> UPLO is CHARACTER*1 00062 *> = 'U': Upper triangle of A is stored and B is factored as 00063 *> U**T*U; 00064 *> = 'L': Lower triangle of A is stored and B is factored as 00065 *> L*L**T. 00066 *> \endverbatim 00067 *> 00068 *> \param[in] N 00069 *> \verbatim 00070 *> N is INTEGER 00071 *> The order of the matrices A and B. N >= 0. 00072 *> \endverbatim 00073 *> 00074 *> \param[in,out] AP 00075 *> \verbatim 00076 *> AP is REAL array, dimension (N*(N+1)/2) 00077 *> On entry, the upper or lower triangle of the symmetric matrix 00078 *> A, packed columnwise in a linear array. The j-th column of A 00079 *> is stored in the array AP as follows: 00080 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j; 00081 *> if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n. 00082 *> 00083 *> On exit, if INFO = 0, the transformed matrix, stored in the 00084 *> same format as A. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] BP 00088 *> \verbatim 00089 *> BP is REAL array, dimension (N*(N+1)/2) 00090 *> The triangular factor from the Cholesky factorization of B, 00091 *> stored in the same format as A, as returned by SPPTRF. 00092 *> \endverbatim 00093 *> 00094 *> \param[out] INFO 00095 *> \verbatim 00096 *> INFO is INTEGER 00097 *> = 0: successful exit 00098 *> < 0: if INFO = -i, the i-th argument had an illegal value 00099 *> \endverbatim 00100 * 00101 * Authors: 00102 * ======== 00103 * 00104 *> \author Univ. of Tennessee 00105 *> \author Univ. of California Berkeley 00106 *> \author Univ. of Colorado Denver 00107 *> \author NAG Ltd. 00108 * 00109 *> \date November 2011 00110 * 00111 *> \ingroup realOTHERcomputational 00112 * 00113 * ===================================================================== 00114 SUBROUTINE SSPGST( ITYPE, UPLO, N, AP, BP, INFO ) 00115 * 00116 * -- LAPACK computational routine (version 3.4.0) -- 00117 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00118 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00119 * November 2011 00120 * 00121 * .. Scalar Arguments .. 00122 CHARACTER UPLO 00123 INTEGER INFO, ITYPE, N 00124 * .. 00125 * .. Array Arguments .. 00126 REAL AP( * ), BP( * ) 00127 * .. 00128 * 00129 * ===================================================================== 00130 * 00131 * .. Parameters .. 00132 REAL ONE, HALF 00133 PARAMETER ( ONE = 1.0, HALF = 0.5 ) 00134 * .. 00135 * .. Local Scalars .. 00136 LOGICAL UPPER 00137 INTEGER J, J1, J1J1, JJ, K, K1, K1K1, KK 00138 REAL AJJ, AKK, BJJ, BKK, CT 00139 * .. 00140 * .. External Subroutines .. 00141 EXTERNAL SAXPY, SSCAL, SSPMV, SSPR2, STPMV, STPSV, 00142 $ XERBLA 00143 * .. 00144 * .. External Functions .. 00145 LOGICAL LSAME 00146 REAL SDOT 00147 EXTERNAL LSAME, SDOT 00148 * .. 00149 * .. Executable Statements .. 00150 * 00151 * Test the input parameters. 00152 * 00153 INFO = 0 00154 UPPER = LSAME( UPLO, 'U' ) 00155 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00156 INFO = -1 00157 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00158 INFO = -2 00159 ELSE IF( N.LT.0 ) THEN 00160 INFO = -3 00161 END IF 00162 IF( INFO.NE.0 ) THEN 00163 CALL XERBLA( 'SSPGST', -INFO ) 00164 RETURN 00165 END IF 00166 * 00167 IF( ITYPE.EQ.1 ) THEN 00168 IF( UPPER ) THEN 00169 * 00170 * Compute inv(U**T)*A*inv(U) 00171 * 00172 * J1 and JJ are the indices of A(1,j) and A(j,j) 00173 * 00174 JJ = 0 00175 DO 10 J = 1, N 00176 J1 = JJ + 1 00177 JJ = JJ + J 00178 * 00179 * Compute the j-th column of the upper triangle of A 00180 * 00181 BJJ = BP( JJ ) 00182 CALL STPSV( UPLO, 'Transpose', 'Nonunit', J, BP, 00183 $ AP( J1 ), 1 ) 00184 CALL SSPMV( UPLO, J-1, -ONE, AP, BP( J1 ), 1, ONE, 00185 $ AP( J1 ), 1 ) 00186 CALL SSCAL( J-1, ONE / BJJ, AP( J1 ), 1 ) 00187 AP( JJ ) = ( AP( JJ )-SDOT( J-1, AP( J1 ), 1, BP( J1 ), 00188 $ 1 ) ) / BJJ 00189 10 CONTINUE 00190 ELSE 00191 * 00192 * Compute inv(L)*A*inv(L**T) 00193 * 00194 * KK and K1K1 are the indices of A(k,k) and A(k+1,k+1) 00195 * 00196 KK = 1 00197 DO 20 K = 1, N 00198 K1K1 = KK + N - K + 1 00199 * 00200 * Update the lower triangle of A(k:n,k:n) 00201 * 00202 AKK = AP( KK ) 00203 BKK = BP( KK ) 00204 AKK = AKK / BKK**2 00205 AP( KK ) = AKK 00206 IF( K.LT.N ) THEN 00207 CALL SSCAL( N-K, ONE / BKK, AP( KK+1 ), 1 ) 00208 CT = -HALF*AKK 00209 CALL SAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) 00210 CALL SSPR2( UPLO, N-K, -ONE, AP( KK+1 ), 1, 00211 $ BP( KK+1 ), 1, AP( K1K1 ) ) 00212 CALL SAXPY( N-K, CT, BP( KK+1 ), 1, AP( KK+1 ), 1 ) 00213 CALL STPSV( UPLO, 'No transpose', 'Non-unit', N-K, 00214 $ BP( K1K1 ), AP( KK+1 ), 1 ) 00215 END IF 00216 KK = K1K1 00217 20 CONTINUE 00218 END IF 00219 ELSE 00220 IF( UPPER ) THEN 00221 * 00222 * Compute U*A*U**T 00223 * 00224 * K1 and KK are the indices of A(1,k) and A(k,k) 00225 * 00226 KK = 0 00227 DO 30 K = 1, N 00228 K1 = KK + 1 00229 KK = KK + K 00230 * 00231 * Update the upper triangle of A(1:k,1:k) 00232 * 00233 AKK = AP( KK ) 00234 BKK = BP( KK ) 00235 CALL STPMV( UPLO, 'No transpose', 'Non-unit', K-1, BP, 00236 $ AP( K1 ), 1 ) 00237 CT = HALF*AKK 00238 CALL SAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) 00239 CALL SSPR2( UPLO, K-1, ONE, AP( K1 ), 1, BP( K1 ), 1, 00240 $ AP ) 00241 CALL SAXPY( K-1, CT, BP( K1 ), 1, AP( K1 ), 1 ) 00242 CALL SSCAL( K-1, BKK, AP( K1 ), 1 ) 00243 AP( KK ) = AKK*BKK**2 00244 30 CONTINUE 00245 ELSE 00246 * 00247 * Compute L**T *A*L 00248 * 00249 * JJ and J1J1 are the indices of A(j,j) and A(j+1,j+1) 00250 * 00251 JJ = 1 00252 DO 40 J = 1, N 00253 J1J1 = JJ + N - J + 1 00254 * 00255 * Compute the j-th column of the lower triangle of A 00256 * 00257 AJJ = AP( JJ ) 00258 BJJ = BP( JJ ) 00259 AP( JJ ) = AJJ*BJJ + SDOT( N-J, AP( JJ+1 ), 1, 00260 $ BP( JJ+1 ), 1 ) 00261 CALL SSCAL( N-J, BJJ, AP( JJ+1 ), 1 ) 00262 CALL SSPMV( UPLO, N-J, ONE, AP( J1J1 ), BP( JJ+1 ), 1, 00263 $ ONE, AP( JJ+1 ), 1 ) 00264 CALL STPMV( UPLO, 'Transpose', 'Non-unit', N-J+1, 00265 $ BP( JJ ), AP( JJ ), 1 ) 00266 JJ = J1J1 00267 40 CONTINUE 00268 END IF 00269 END IF 00270 RETURN 00271 * 00272 * End of SSPGST 00273 * 00274 END