1 SUBROUTINE DPTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
2 *
3 * -- LAPACK 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 COMPZ
10 INTEGER INFO, LDZ, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DPTEQR computes all eigenvalues and, optionally, eigenvectors of a
20 * symmetric positive definite tridiagonal matrix by first factoring the
21 * matrix using DPTTRF, and then calling DBDSQR to compute the singular
22 * values of the bidiagonal factor.
23 *
24 * This routine computes the eigenvalues of the positive definite
25 * tridiagonal matrix to high relative accuracy. This means that if the
26 * eigenvalues range over many orders of magnitude in size, then the
27 * small eigenvalues and corresponding eigenvectors will be computed
28 * more accurately than, for example, with the standard QR method.
29 *
30 * The eigenvectors of a full or band symmetric positive definite matrix
31 * can also be found if DSYTRD, DSPTRD, or DSBTRD has been used to
32 * reduce this matrix to tridiagonal form. (The reduction to tridiagonal
33 * form, however, may preclude the possibility of obtaining high
34 * relative accuracy in the small eigenvalues of the original matrix, if
35 * these eigenvalues range over many orders of magnitude.)
36 *
37 * Arguments
38 * =========
39 *
40 * COMPZ (input) CHARACTER*1
41 * = 'N': Compute eigenvalues only.
42 * = 'V': Compute eigenvectors of original symmetric
43 * matrix also. Array Z contains the orthogonal
44 * matrix used to reduce the original matrix to
45 * tridiagonal form.
46 * = 'I': Compute eigenvectors of tridiagonal matrix also.
47 *
48 * N (input) INTEGER
49 * The order of the matrix. N >= 0.
50 *
51 * D (input/output) DOUBLE PRECISION array, dimension (N)
52 * On entry, the n diagonal elements of the tridiagonal
53 * matrix.
54 * On normal exit, D contains the eigenvalues, in descending
55 * order.
56 *
57 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
58 * On entry, the (n-1) subdiagonal elements of the tridiagonal
59 * matrix.
60 * On exit, E has been destroyed.
61 *
62 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
63 * On entry, if COMPZ = 'V', the orthogonal matrix used in the
64 * reduction to tridiagonal form.
65 * On exit, if COMPZ = 'V', the orthonormal eigenvectors of the
66 * original symmetric matrix;
67 * if COMPZ = 'I', the orthonormal eigenvectors of the
68 * tridiagonal matrix.
69 * If INFO > 0 on exit, Z contains the eigenvectors associated
70 * with only the stored eigenvalues.
71 * If COMPZ = 'N', then Z is not referenced.
72 *
73 * LDZ (input) INTEGER
74 * The leading dimension of the array Z. LDZ >= 1, and if
75 * COMPZ = 'V' or 'I', LDZ >= max(1,N).
76 *
77 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
78 *
79 * INFO (output) INTEGER
80 * = 0: successful exit.
81 * < 0: if INFO = -i, the i-th argument had an illegal value.
82 * > 0: if INFO = i, and i is:
83 * <= N the Cholesky factorization of the matrix could
84 * not be performed because the i-th principal minor
85 * was not positive definite.
86 * > N the SVD algorithm failed to converge;
87 * if INFO = N+i, i off-diagonal elements of the
88 * bidiagonal factor did not converge to zero.
89 *
90 * =====================================================================
91 *
92 * .. Parameters ..
93 DOUBLE PRECISION ZERO, ONE
94 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
95 * ..
96 * .. External Functions ..
97 LOGICAL LSAME
98 EXTERNAL LSAME
99 * ..
100 * .. External Subroutines ..
101 EXTERNAL DBDSQR, DLASET, DPTTRF, XERBLA
102 * ..
103 * .. Local Arrays ..
104 DOUBLE PRECISION C( 1, 1 ), VT( 1, 1 )
105 * ..
106 * .. Local Scalars ..
107 INTEGER I, ICOMPZ, NRU
108 * ..
109 * .. Intrinsic Functions ..
110 INTRINSIC MAX, SQRT
111 * ..
112 * .. Executable Statements ..
113 *
114 * Test the input parameters.
115 *
116 INFO = 0
117 *
118 IF( LSAME( COMPZ, 'N' ) ) THEN
119 ICOMPZ = 0
120 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
121 ICOMPZ = 1
122 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
123 ICOMPZ = 2
124 ELSE
125 ICOMPZ = -1
126 END IF
127 IF( ICOMPZ.LT.0 ) THEN
128 INFO = -1
129 ELSE IF( N.LT.0 ) THEN
130 INFO = -2
131 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
132 $ N ) ) ) THEN
133 INFO = -6
134 END IF
135 IF( INFO.NE.0 ) THEN
136 CALL XERBLA( 'DPTEQR', -INFO )
137 RETURN
138 END IF
139 *
140 * Quick return if possible
141 *
142 IF( N.EQ.0 )
143 $ RETURN
144 *
145 IF( N.EQ.1 ) THEN
146 IF( ICOMPZ.GT.0 )
147 $ Z( 1, 1 ) = ONE
148 RETURN
149 END IF
150 IF( ICOMPZ.EQ.2 )
151 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
152 *
153 * Call DPTTRF to factor the matrix.
154 *
155 CALL DPTTRF( N, D, E, INFO )
156 IF( INFO.NE.0 )
157 $ RETURN
158 DO 10 I = 1, N
159 D( I ) = SQRT( D( I ) )
160 10 CONTINUE
161 DO 20 I = 1, N - 1
162 E( I ) = E( I )*D( I )
163 20 CONTINUE
164 *
165 * Call DBDSQR to compute the singular values/vectors of the
166 * bidiagonal factor.
167 *
168 IF( ICOMPZ.GT.0 ) THEN
169 NRU = N
170 ELSE
171 NRU = 0
172 END IF
173 CALL DBDSQR( 'Lower', N, 0, NRU, 0, D, E, VT, 1, Z, LDZ, C, 1,
174 $ WORK, INFO )
175 *
176 * Square the singular values.
177 *
178 IF( INFO.EQ.0 ) THEN
179 DO 30 I = 1, N
180 D( I ) = D( I )*D( I )
181 30 CONTINUE
182 ELSE
183 INFO = N + INFO
184 END IF
185 *
186 RETURN
187 *
188 * End of DPTEQR
189 *
190 END
2 *
3 * -- LAPACK 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 COMPZ
10 INTEGER INFO, LDZ, N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION D( * ), E( * ), WORK( * ), Z( LDZ, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DPTEQR computes all eigenvalues and, optionally, eigenvectors of a
20 * symmetric positive definite tridiagonal matrix by first factoring the
21 * matrix using DPTTRF, and then calling DBDSQR to compute the singular
22 * values of the bidiagonal factor.
23 *
24 * This routine computes the eigenvalues of the positive definite
25 * tridiagonal matrix to high relative accuracy. This means that if the
26 * eigenvalues range over many orders of magnitude in size, then the
27 * small eigenvalues and corresponding eigenvectors will be computed
28 * more accurately than, for example, with the standard QR method.
29 *
30 * The eigenvectors of a full or band symmetric positive definite matrix
31 * can also be found if DSYTRD, DSPTRD, or DSBTRD has been used to
32 * reduce this matrix to tridiagonal form. (The reduction to tridiagonal
33 * form, however, may preclude the possibility of obtaining high
34 * relative accuracy in the small eigenvalues of the original matrix, if
35 * these eigenvalues range over many orders of magnitude.)
36 *
37 * Arguments
38 * =========
39 *
40 * COMPZ (input) CHARACTER*1
41 * = 'N': Compute eigenvalues only.
42 * = 'V': Compute eigenvectors of original symmetric
43 * matrix also. Array Z contains the orthogonal
44 * matrix used to reduce the original matrix to
45 * tridiagonal form.
46 * = 'I': Compute eigenvectors of tridiagonal matrix also.
47 *
48 * N (input) INTEGER
49 * The order of the matrix. N >= 0.
50 *
51 * D (input/output) DOUBLE PRECISION array, dimension (N)
52 * On entry, the n diagonal elements of the tridiagonal
53 * matrix.
54 * On normal exit, D contains the eigenvalues, in descending
55 * order.
56 *
57 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
58 * On entry, the (n-1) subdiagonal elements of the tridiagonal
59 * matrix.
60 * On exit, E has been destroyed.
61 *
62 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
63 * On entry, if COMPZ = 'V', the orthogonal matrix used in the
64 * reduction to tridiagonal form.
65 * On exit, if COMPZ = 'V', the orthonormal eigenvectors of the
66 * original symmetric matrix;
67 * if COMPZ = 'I', the orthonormal eigenvectors of the
68 * tridiagonal matrix.
69 * If INFO > 0 on exit, Z contains the eigenvectors associated
70 * with only the stored eigenvalues.
71 * If COMPZ = 'N', then Z is not referenced.
72 *
73 * LDZ (input) INTEGER
74 * The leading dimension of the array Z. LDZ >= 1, and if
75 * COMPZ = 'V' or 'I', LDZ >= max(1,N).
76 *
77 * WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
78 *
79 * INFO (output) INTEGER
80 * = 0: successful exit.
81 * < 0: if INFO = -i, the i-th argument had an illegal value.
82 * > 0: if INFO = i, and i is:
83 * <= N the Cholesky factorization of the matrix could
84 * not be performed because the i-th principal minor
85 * was not positive definite.
86 * > N the SVD algorithm failed to converge;
87 * if INFO = N+i, i off-diagonal elements of the
88 * bidiagonal factor did not converge to zero.
89 *
90 * =====================================================================
91 *
92 * .. Parameters ..
93 DOUBLE PRECISION ZERO, ONE
94 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
95 * ..
96 * .. External Functions ..
97 LOGICAL LSAME
98 EXTERNAL LSAME
99 * ..
100 * .. External Subroutines ..
101 EXTERNAL DBDSQR, DLASET, DPTTRF, XERBLA
102 * ..
103 * .. Local Arrays ..
104 DOUBLE PRECISION C( 1, 1 ), VT( 1, 1 )
105 * ..
106 * .. Local Scalars ..
107 INTEGER I, ICOMPZ, NRU
108 * ..
109 * .. Intrinsic Functions ..
110 INTRINSIC MAX, SQRT
111 * ..
112 * .. Executable Statements ..
113 *
114 * Test the input parameters.
115 *
116 INFO = 0
117 *
118 IF( LSAME( COMPZ, 'N' ) ) THEN
119 ICOMPZ = 0
120 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
121 ICOMPZ = 1
122 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
123 ICOMPZ = 2
124 ELSE
125 ICOMPZ = -1
126 END IF
127 IF( ICOMPZ.LT.0 ) THEN
128 INFO = -1
129 ELSE IF( N.LT.0 ) THEN
130 INFO = -2
131 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
132 $ N ) ) ) THEN
133 INFO = -6
134 END IF
135 IF( INFO.NE.0 ) THEN
136 CALL XERBLA( 'DPTEQR', -INFO )
137 RETURN
138 END IF
139 *
140 * Quick return if possible
141 *
142 IF( N.EQ.0 )
143 $ RETURN
144 *
145 IF( N.EQ.1 ) THEN
146 IF( ICOMPZ.GT.0 )
147 $ Z( 1, 1 ) = ONE
148 RETURN
149 END IF
150 IF( ICOMPZ.EQ.2 )
151 $ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
152 *
153 * Call DPTTRF to factor the matrix.
154 *
155 CALL DPTTRF( N, D, E, INFO )
156 IF( INFO.NE.0 )
157 $ RETURN
158 DO 10 I = 1, N
159 D( I ) = SQRT( D( I ) )
160 10 CONTINUE
161 DO 20 I = 1, N - 1
162 E( I ) = E( I )*D( I )
163 20 CONTINUE
164 *
165 * Call DBDSQR to compute the singular values/vectors of the
166 * bidiagonal factor.
167 *
168 IF( ICOMPZ.GT.0 ) THEN
169 NRU = N
170 ELSE
171 NRU = 0
172 END IF
173 CALL DBDSQR( 'Lower', N, 0, NRU, 0, D, E, VT, 1, Z, LDZ, C, 1,
174 $ WORK, INFO )
175 *
176 * Square the singular values.
177 *
178 IF( INFO.EQ.0 ) THEN
179 DO 30 I = 1, N
180 D( I ) = D( I )*D( I )
181 30 CONTINUE
182 ELSE
183 INFO = N + INFO
184 END IF
185 *
186 RETURN
187 *
188 * End of DPTEQR
189 *
190 END