 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DSPGST( ITYPE, UPLO, N, AP, BP, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.3.1) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * -- April 2011 -- 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**T)*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**T) 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**T 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**T *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
1.7.3 
