1 SUBROUTINE DSPEV( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, INFO )
2 *
3 * -- LAPACK driver routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER JOBZ, UPLO
10 INTEGER INFO, LDZ, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSPEV computes all the eigenvalues and, optionally, eigenvectors of a
20 * real symmetric matrix A in packed storage.
21 *
22 * Arguments
23 * =========
24 *
25 * JOBZ (input) CHARACTER*1
26 * = 'N': Compute eigenvalues only;
27 * = 'V': Compute eigenvalues and eigenvectors.
28 *
29 * UPLO (input) CHARACTER*1
30 * = 'U': Upper triangle of A is stored;
31 * = 'L': Lower triangle of A is stored.
32 *
33 * N (input) INTEGER
34 * The order of the matrix A. N >= 0.
35 *
36 * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
37 * On entry, the upper or lower triangle of the symmetric matrix
38 * A, packed columnwise in a linear array. The j-th column of A
39 * is stored in the array AP as follows:
40 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
41 * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
42 *
43 * On exit, AP is overwritten by values generated during the
44 * reduction to tridiagonal form. If UPLO = 'U', the diagonal
45 * and first superdiagonal of the tridiagonal matrix T overwrite
46 * the corresponding elements of A, and if UPLO = 'L', the
47 * diagonal and first subdiagonal of T overwrite the
48 * corresponding elements of A.
49 *
50 * W (output) DOUBLE PRECISION array, dimension (N)
51 * If INFO = 0, the eigenvalues in ascending order.
52 *
53 * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
54 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
55 * eigenvectors of the matrix A, with the i-th column of Z
56 * holding the eigenvector associated with W(i).
57 * If JOBZ = 'N', then Z is not referenced.
58 *
59 * LDZ (input) INTEGER
60 * The leading dimension of the array Z. LDZ >= 1, and if
61 * JOBZ = 'V', LDZ >= max(1,N).
62 *
63 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
64 *
65 * INFO (output) INTEGER
66 * = 0: successful exit.
67 * < 0: if INFO = -i, the i-th argument had an illegal value.
68 * > 0: if INFO = i, the algorithm failed to converge; i
69 * off-diagonal elements of an intermediate tridiagonal
70 * form did not converge to zero.
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 DOUBLE PRECISION ZERO, ONE
76 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
77 * ..
78 * .. Local Scalars ..
79 LOGICAL WANTZ
80 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE
81 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
82 $ SMLNUM
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 DOUBLE PRECISION DLAMCH, DLANSP
87 EXTERNAL LSAME, DLAMCH, DLANSP
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL DOPGTR, DSCAL, DSPTRD, DSTEQR, DSTERF, XERBLA
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC SQRT
94 * ..
95 * .. Executable Statements ..
96 *
97 * Test the input parameters.
98 *
99 WANTZ = LSAME( JOBZ, 'V' )
100 *
101 INFO = 0
102 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
103 INFO = -1
104 ELSE IF( .NOT.( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) )
105 $ THEN
106 INFO = -2
107 ELSE IF( N.LT.0 ) THEN
108 INFO = -3
109 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
110 INFO = -7
111 END IF
112 *
113 IF( INFO.NE.0 ) THEN
114 CALL XERBLA( 'DSPEV ', -INFO )
115 RETURN
116 END IF
117 *
118 * Quick return if possible
119 *
120 IF( N.EQ.0 )
121 $ RETURN
122 *
123 IF( N.EQ.1 ) THEN
124 W( 1 ) = AP( 1 )
125 IF( WANTZ )
126 $ Z( 1, 1 ) = ONE
127 RETURN
128 END IF
129 *
130 * Get machine constants.
131 *
132 SAFMIN = DLAMCH( 'Safe minimum' )
133 EPS = DLAMCH( 'Precision' )
134 SMLNUM = SAFMIN / EPS
135 BIGNUM = ONE / SMLNUM
136 RMIN = SQRT( SMLNUM )
137 RMAX = SQRT( BIGNUM )
138 *
139 * Scale matrix to allowable range, if necessary.
140 *
141 ANRM = DLANSP( 'M', UPLO, N, AP, WORK )
142 ISCALE = 0
143 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
144 ISCALE = 1
145 SIGMA = RMIN / ANRM
146 ELSE IF( ANRM.GT.RMAX ) THEN
147 ISCALE = 1
148 SIGMA = RMAX / ANRM
149 END IF
150 IF( ISCALE.EQ.1 ) THEN
151 CALL DSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
152 END IF
153 *
154 * Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
155 *
156 INDE = 1
157 INDTAU = INDE + N
158 CALL DSPTRD( UPLO, N, AP, W, WORK( INDE ), WORK( INDTAU ), IINFO )
159 *
160 * For eigenvalues only, call DSTERF. For eigenvectors, first call
161 * DOPGTR to generate the orthogonal matrix, then call DSTEQR.
162 *
163 IF( .NOT.WANTZ ) THEN
164 CALL DSTERF( N, W, WORK( INDE ), INFO )
165 ELSE
166 INDWRK = INDTAU + N
167 CALL DOPGTR( UPLO, N, AP, WORK( INDTAU ), Z, LDZ,
168 $ WORK( INDWRK ), IINFO )
169 CALL DSTEQR( JOBZ, N, W, WORK( INDE ), Z, LDZ, WORK( INDTAU ),
170 $ INFO )
171 END IF
172 *
173 * If matrix was scaled, then rescale eigenvalues appropriately.
174 *
175 IF( ISCALE.EQ.1 ) THEN
176 IF( INFO.EQ.0 ) THEN
177 IMAX = N
178 ELSE
179 IMAX = INFO - 1
180 END IF
181 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
182 END IF
183 *
184 RETURN
185 *
186 * End of DSPEV
187 *
188 END
2 *
3 * -- LAPACK driver routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER JOBZ, UPLO
10 INTEGER INFO, LDZ, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AP( * ), W( * ), WORK( * ), Z( LDZ, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DSPEV computes all the eigenvalues and, optionally, eigenvectors of a
20 * real symmetric matrix A in packed storage.
21 *
22 * Arguments
23 * =========
24 *
25 * JOBZ (input) CHARACTER*1
26 * = 'N': Compute eigenvalues only;
27 * = 'V': Compute eigenvalues and eigenvectors.
28 *
29 * UPLO (input) CHARACTER*1
30 * = 'U': Upper triangle of A is stored;
31 * = 'L': Lower triangle of A is stored.
32 *
33 * N (input) INTEGER
34 * The order of the matrix A. N >= 0.
35 *
36 * AP (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
37 * On entry, the upper or lower triangle of the symmetric matrix
38 * A, packed columnwise in a linear array. The j-th column of A
39 * is stored in the array AP as follows:
40 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
41 * if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
42 *
43 * On exit, AP is overwritten by values generated during the
44 * reduction to tridiagonal form. If UPLO = 'U', the diagonal
45 * and first superdiagonal of the tridiagonal matrix T overwrite
46 * the corresponding elements of A, and if UPLO = 'L', the
47 * diagonal and first subdiagonal of T overwrite the
48 * corresponding elements of A.
49 *
50 * W (output) DOUBLE PRECISION array, dimension (N)
51 * If INFO = 0, the eigenvalues in ascending order.
52 *
53 * Z (output) DOUBLE PRECISION array, dimension (LDZ, N)
54 * If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
55 * eigenvectors of the matrix A, with the i-th column of Z
56 * holding the eigenvector associated with W(i).
57 * If JOBZ = 'N', then Z is not referenced.
58 *
59 * LDZ (input) INTEGER
60 * The leading dimension of the array Z. LDZ >= 1, and if
61 * JOBZ = 'V', LDZ >= max(1,N).
62 *
63 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
64 *
65 * INFO (output) INTEGER
66 * = 0: successful exit.
67 * < 0: if INFO = -i, the i-th argument had an illegal value.
68 * > 0: if INFO = i, the algorithm failed to converge; i
69 * off-diagonal elements of an intermediate tridiagonal
70 * form did not converge to zero.
71 *
72 * =====================================================================
73 *
74 * .. Parameters ..
75 DOUBLE PRECISION ZERO, ONE
76 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
77 * ..
78 * .. Local Scalars ..
79 LOGICAL WANTZ
80 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE
81 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
82 $ SMLNUM
83 * ..
84 * .. External Functions ..
85 LOGICAL LSAME
86 DOUBLE PRECISION DLAMCH, DLANSP
87 EXTERNAL LSAME, DLAMCH, DLANSP
88 * ..
89 * .. External Subroutines ..
90 EXTERNAL DOPGTR, DSCAL, DSPTRD, DSTEQR, DSTERF, XERBLA
91 * ..
92 * .. Intrinsic Functions ..
93 INTRINSIC SQRT
94 * ..
95 * .. Executable Statements ..
96 *
97 * Test the input parameters.
98 *
99 WANTZ = LSAME( JOBZ, 'V' )
100 *
101 INFO = 0
102 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
103 INFO = -1
104 ELSE IF( .NOT.( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) )
105 $ THEN
106 INFO = -2
107 ELSE IF( N.LT.0 ) THEN
108 INFO = -3
109 ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
110 INFO = -7
111 END IF
112 *
113 IF( INFO.NE.0 ) THEN
114 CALL XERBLA( 'DSPEV ', -INFO )
115 RETURN
116 END IF
117 *
118 * Quick return if possible
119 *
120 IF( N.EQ.0 )
121 $ RETURN
122 *
123 IF( N.EQ.1 ) THEN
124 W( 1 ) = AP( 1 )
125 IF( WANTZ )
126 $ Z( 1, 1 ) = ONE
127 RETURN
128 END IF
129 *
130 * Get machine constants.
131 *
132 SAFMIN = DLAMCH( 'Safe minimum' )
133 EPS = DLAMCH( 'Precision' )
134 SMLNUM = SAFMIN / EPS
135 BIGNUM = ONE / SMLNUM
136 RMIN = SQRT( SMLNUM )
137 RMAX = SQRT( BIGNUM )
138 *
139 * Scale matrix to allowable range, if necessary.
140 *
141 ANRM = DLANSP( 'M', UPLO, N, AP, WORK )
142 ISCALE = 0
143 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
144 ISCALE = 1
145 SIGMA = RMIN / ANRM
146 ELSE IF( ANRM.GT.RMAX ) THEN
147 ISCALE = 1
148 SIGMA = RMAX / ANRM
149 END IF
150 IF( ISCALE.EQ.1 ) THEN
151 CALL DSCAL( ( N*( N+1 ) ) / 2, SIGMA, AP, 1 )
152 END IF
153 *
154 * Call DSPTRD to reduce symmetric packed matrix to tridiagonal form.
155 *
156 INDE = 1
157 INDTAU = INDE + N
158 CALL DSPTRD( UPLO, N, AP, W, WORK( INDE ), WORK( INDTAU ), IINFO )
159 *
160 * For eigenvalues only, call DSTERF. For eigenvectors, first call
161 * DOPGTR to generate the orthogonal matrix, then call DSTEQR.
162 *
163 IF( .NOT.WANTZ ) THEN
164 CALL DSTERF( N, W, WORK( INDE ), INFO )
165 ELSE
166 INDWRK = INDTAU + N
167 CALL DOPGTR( UPLO, N, AP, WORK( INDTAU ), Z, LDZ,
168 $ WORK( INDWRK ), IINFO )
169 CALL DSTEQR( JOBZ, N, W, WORK( INDE ), Z, LDZ, WORK( INDTAU ),
170 $ INFO )
171 END IF
172 *
173 * If matrix was scaled, then rescale eigenvalues appropriately.
174 *
175 IF( ISCALE.EQ.1 ) THEN
176 IF( INFO.EQ.0 ) THEN
177 IMAX = N
178 ELSE
179 IMAX = INFO - 1
180 END IF
181 CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
182 END IF
183 *
184 RETURN
185 *
186 * End of DSPEV
187 *
188 END