![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DSYGS2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DSYGS2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsygs2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsygs2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsygs2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DSYGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER UPLO 00025 * INTEGER INFO, ITYPE, LDA, LDB, N 00026 * .. 00027 * .. Array Arguments .. 00028 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> DSYGS2 reduces a real symmetric-definite generalized eigenproblem 00038 *> to standard form. 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 DPOTRF. 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 *> Specifies whether the upper or lower triangular part of the 00063 *> symmetric matrix A is stored, and how B has been factorized. 00064 *> = 'U': Upper triangular 00065 *> = 'L': Lower triangular 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] A 00075 *> \verbatim 00076 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00077 *> On entry, the symmetric matrix A. If UPLO = 'U', the leading 00078 *> n by n upper triangular part of A contains the upper 00079 *> triangular part of the matrix A, and the strictly lower 00080 *> triangular part of A is not referenced. If UPLO = 'L', the 00081 *> leading n by n lower triangular part of A contains the lower 00082 *> triangular part of the matrix A, and the strictly upper 00083 *> triangular part of A is not referenced. 00084 *> 00085 *> On exit, if INFO = 0, the transformed matrix, stored in the 00086 *> same format as A. 00087 *> \endverbatim 00088 *> 00089 *> \param[in] LDA 00090 *> \verbatim 00091 *> LDA is INTEGER 00092 *> The leading dimension of the array A. LDA >= max(1,N). 00093 *> \endverbatim 00094 *> 00095 *> \param[in] B 00096 *> \verbatim 00097 *> B is DOUBLE PRECISION array, dimension (LDB,N) 00098 *> The triangular factor from the Cholesky factorization of B, 00099 *> as returned by DPOTRF. 00100 *> \endverbatim 00101 *> 00102 *> \param[in] LDB 00103 *> \verbatim 00104 *> LDB is INTEGER 00105 *> The leading dimension of the array B. LDB >= max(1,N). 00106 *> \endverbatim 00107 *> 00108 *> \param[out] INFO 00109 *> \verbatim 00110 *> INFO is INTEGER 00111 *> = 0: successful exit. 00112 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00113 *> \endverbatim 00114 * 00115 * Authors: 00116 * ======== 00117 * 00118 *> \author Univ. of Tennessee 00119 *> \author Univ. of California Berkeley 00120 *> \author Univ. of Colorado Denver 00121 *> \author NAG Ltd. 00122 * 00123 *> \date November 2011 00124 * 00125 *> \ingroup doubleSYcomputational 00126 * 00127 * ===================================================================== 00128 SUBROUTINE DSYGS2( ITYPE, UPLO, N, A, LDA, B, LDB, INFO ) 00129 * 00130 * -- LAPACK computational routine (version 3.4.0) -- 00131 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00132 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00133 * November 2011 00134 * 00135 * .. Scalar Arguments .. 00136 CHARACTER UPLO 00137 INTEGER INFO, ITYPE, LDA, LDB, N 00138 * .. 00139 * .. Array Arguments .. 00140 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 00141 * .. 00142 * 00143 * ===================================================================== 00144 * 00145 * .. Parameters .. 00146 DOUBLE PRECISION ONE, HALF 00147 PARAMETER ( ONE = 1.0D0, HALF = 0.5D0 ) 00148 * .. 00149 * .. Local Scalars .. 00150 LOGICAL UPPER 00151 INTEGER K 00152 DOUBLE PRECISION AKK, BKK, CT 00153 * .. 00154 * .. External Subroutines .. 00155 EXTERNAL DAXPY, DSCAL, DSYR2, DTRMV, DTRSV, XERBLA 00156 * .. 00157 * .. Intrinsic Functions .. 00158 INTRINSIC MAX 00159 * .. 00160 * .. External Functions .. 00161 LOGICAL LSAME 00162 EXTERNAL LSAME 00163 * .. 00164 * .. Executable Statements .. 00165 * 00166 * Test the input parameters. 00167 * 00168 INFO = 0 00169 UPPER = LSAME( UPLO, 'U' ) 00170 IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN 00171 INFO = -1 00172 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00173 INFO = -2 00174 ELSE IF( N.LT.0 ) THEN 00175 INFO = -3 00176 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00177 INFO = -5 00178 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00179 INFO = -7 00180 END IF 00181 IF( INFO.NE.0 ) THEN 00182 CALL XERBLA( 'DSYGS2', -INFO ) 00183 RETURN 00184 END IF 00185 * 00186 IF( ITYPE.EQ.1 ) THEN 00187 IF( UPPER ) THEN 00188 * 00189 * Compute inv(U**T)*A*inv(U) 00190 * 00191 DO 10 K = 1, N 00192 * 00193 * Update the upper triangle of A(k:n,k:n) 00194 * 00195 AKK = A( K, K ) 00196 BKK = B( K, K ) 00197 AKK = AKK / BKK**2 00198 A( K, K ) = AKK 00199 IF( K.LT.N ) THEN 00200 CALL DSCAL( N-K, ONE / BKK, A( K, K+1 ), LDA ) 00201 CT = -HALF*AKK 00202 CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ), 00203 $ LDA ) 00204 CALL DSYR2( UPLO, N-K, -ONE, A( K, K+1 ), LDA, 00205 $ B( K, K+1 ), LDB, A( K+1, K+1 ), LDA ) 00206 CALL DAXPY( N-K, CT, B( K, K+1 ), LDB, A( K, K+1 ), 00207 $ LDA ) 00208 CALL DTRSV( UPLO, 'Transpose', 'Non-unit', N-K, 00209 $ B( K+1, K+1 ), LDB, A( K, K+1 ), LDA ) 00210 END IF 00211 10 CONTINUE 00212 ELSE 00213 * 00214 * Compute inv(L)*A*inv(L**T) 00215 * 00216 DO 20 K = 1, N 00217 * 00218 * Update the lower triangle of A(k:n,k:n) 00219 * 00220 AKK = A( K, K ) 00221 BKK = B( K, K ) 00222 AKK = AKK / BKK**2 00223 A( K, K ) = AKK 00224 IF( K.LT.N ) THEN 00225 CALL DSCAL( N-K, ONE / BKK, A( K+1, K ), 1 ) 00226 CT = -HALF*AKK 00227 CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 ) 00228 CALL DSYR2( UPLO, N-K, -ONE, A( K+1, K ), 1, 00229 $ B( K+1, K ), 1, A( K+1, K+1 ), LDA ) 00230 CALL DAXPY( N-K, CT, B( K+1, K ), 1, A( K+1, K ), 1 ) 00231 CALL DTRSV( UPLO, 'No transpose', 'Non-unit', N-K, 00232 $ B( K+1, K+1 ), LDB, A( K+1, K ), 1 ) 00233 END IF 00234 20 CONTINUE 00235 END IF 00236 ELSE 00237 IF( UPPER ) THEN 00238 * 00239 * Compute U*A*U**T 00240 * 00241 DO 30 K = 1, N 00242 * 00243 * Update the upper triangle of A(1:k,1:k) 00244 * 00245 AKK = A( K, K ) 00246 BKK = B( K, K ) 00247 CALL DTRMV( UPLO, 'No transpose', 'Non-unit', K-1, B, 00248 $ LDB, A( 1, K ), 1 ) 00249 CT = HALF*AKK 00250 CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 ) 00251 CALL DSYR2( UPLO, K-1, ONE, A( 1, K ), 1, B( 1, K ), 1, 00252 $ A, LDA ) 00253 CALL DAXPY( K-1, CT, B( 1, K ), 1, A( 1, K ), 1 ) 00254 CALL DSCAL( K-1, BKK, A( 1, K ), 1 ) 00255 A( K, K ) = AKK*BKK**2 00256 30 CONTINUE 00257 ELSE 00258 * 00259 * Compute L**T *A*L 00260 * 00261 DO 40 K = 1, N 00262 * 00263 * Update the lower triangle of A(1:k,1:k) 00264 * 00265 AKK = A( K, K ) 00266 BKK = B( K, K ) 00267 CALL DTRMV( UPLO, 'Transpose', 'Non-unit', K-1, B, LDB, 00268 $ A( K, 1 ), LDA ) 00269 CT = HALF*AKK 00270 CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA ) 00271 CALL DSYR2( UPLO, K-1, ONE, A( K, 1 ), LDA, B( K, 1 ), 00272 $ LDB, A, LDA ) 00273 CALL DAXPY( K-1, CT, B( K, 1 ), LDB, A( K, 1 ), LDA ) 00274 CALL DSCAL( K-1, BKK, A( K, 1 ), LDA ) 00275 A( K, K ) = AKK*BKK**2 00276 40 CONTINUE 00277 END IF 00278 END IF 00279 RETURN 00280 * 00281 * End of DSYGS2 00282 * 00283 END