1 SUBROUTINE ZPTEQR( 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( * )
14 COMPLEX*16 Z( LDZ, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZPTEQR computes all eigenvalues and, optionally, eigenvectors of a
21 * symmetric positive definite tridiagonal matrix by first factoring the
22 * matrix using DPTTRF and then calling ZBDSQR to compute the singular
23 * values of the bidiagonal factor.
24 *
25 * This routine computes the eigenvalues of the positive definite
26 * tridiagonal matrix to high relative accuracy. This means that if the
27 * eigenvalues range over many orders of magnitude in size, then the
28 * small eigenvalues and corresponding eigenvectors will be computed
29 * more accurately than, for example, with the standard QR method.
30 *
31 * The eigenvectors of a full or band positive definite Hermitian matrix
32 * can also be found if ZHETRD, ZHPTRD, or ZHBTRD has been used to
33 * reduce this matrix to tridiagonal form. (The reduction to
34 * tridiagonal form, however, may preclude the possibility of obtaining
35 * high relative accuracy in the small eigenvalues of the original
36 * matrix, if these eigenvalues range over many orders of magnitude.)
37 *
38 * Arguments
39 * =========
40 *
41 * COMPZ (input) CHARACTER*1
42 * = 'N': Compute eigenvalues only.
43 * = 'V': Compute eigenvectors of original Hermitian
44 * matrix also. Array Z contains the unitary matrix
45 * used to reduce the original matrix to tridiagonal
46 * form.
47 * = 'I': Compute eigenvectors of tridiagonal matrix also.
48 *
49 * N (input) INTEGER
50 * The order of the matrix. N >= 0.
51 *
52 * D (input/output) DOUBLE PRECISION array, dimension (N)
53 * On entry, the n diagonal elements of the tridiagonal 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) COMPLEX*16 array, dimension (LDZ, N)
63 * On entry, if COMPZ = 'V', the unitary matrix used in the
64 * reduction to tridiagonal form.
65 * On exit, if COMPZ = 'V', the orthonormal eigenvectors of the
66 * original Hermitian 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 COMPLEX*16 CZERO, CONE
94 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
95 $ CONE = ( 1.0D+0, 0.0D+0 ) )
96 * ..
97 * .. External Functions ..
98 LOGICAL LSAME
99 EXTERNAL LSAME
100 * ..
101 * .. External Subroutines ..
102 EXTERNAL DPTTRF, XERBLA, ZBDSQR, ZLASET
103 * ..
104 * .. Local Arrays ..
105 COMPLEX*16 C( 1, 1 ), VT( 1, 1 )
106 * ..
107 * .. Local Scalars ..
108 INTEGER I, ICOMPZ, NRU
109 * ..
110 * .. Intrinsic Functions ..
111 INTRINSIC MAX, SQRT
112 * ..
113 * .. Executable Statements ..
114 *
115 * Test the input parameters.
116 *
117 INFO = 0
118 *
119 IF( LSAME( COMPZ, 'N' ) ) THEN
120 ICOMPZ = 0
121 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
122 ICOMPZ = 1
123 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
124 ICOMPZ = 2
125 ELSE
126 ICOMPZ = -1
127 END IF
128 IF( ICOMPZ.LT.0 ) THEN
129 INFO = -1
130 ELSE IF( N.LT.0 ) THEN
131 INFO = -2
132 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
133 $ N ) ) ) THEN
134 INFO = -6
135 END IF
136 IF( INFO.NE.0 ) THEN
137 CALL XERBLA( 'ZPTEQR', -INFO )
138 RETURN
139 END IF
140 *
141 * Quick return if possible
142 *
143 IF( N.EQ.0 )
144 $ RETURN
145 *
146 IF( N.EQ.1 ) THEN
147 IF( ICOMPZ.GT.0 )
148 $ Z( 1, 1 ) = CONE
149 RETURN
150 END IF
151 IF( ICOMPZ.EQ.2 )
152 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
153 *
154 * Call DPTTRF to factor the matrix.
155 *
156 CALL DPTTRF( N, D, E, INFO )
157 IF( INFO.NE.0 )
158 $ RETURN
159 DO 10 I = 1, N
160 D( I ) = SQRT( D( I ) )
161 10 CONTINUE
162 DO 20 I = 1, N - 1
163 E( I ) = E( I )*D( I )
164 20 CONTINUE
165 *
166 * Call ZBDSQR to compute the singular values/vectors of the
167 * bidiagonal factor.
168 *
169 IF( ICOMPZ.GT.0 ) THEN
170 NRU = N
171 ELSE
172 NRU = 0
173 END IF
174 CALL ZBDSQR( 'Lower', N, 0, NRU, 0, D, E, VT, 1, Z, LDZ, C, 1,
175 $ WORK, INFO )
176 *
177 * Square the singular values.
178 *
179 IF( INFO.EQ.0 ) THEN
180 DO 30 I = 1, N
181 D( I ) = D( I )*D( I )
182 30 CONTINUE
183 ELSE
184 INFO = N + INFO
185 END IF
186 *
187 RETURN
188 *
189 * End of ZPTEQR
190 *
191 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( * )
14 COMPLEX*16 Z( LDZ, * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZPTEQR computes all eigenvalues and, optionally, eigenvectors of a
21 * symmetric positive definite tridiagonal matrix by first factoring the
22 * matrix using DPTTRF and then calling ZBDSQR to compute the singular
23 * values of the bidiagonal factor.
24 *
25 * This routine computes the eigenvalues of the positive definite
26 * tridiagonal matrix to high relative accuracy. This means that if the
27 * eigenvalues range over many orders of magnitude in size, then the
28 * small eigenvalues and corresponding eigenvectors will be computed
29 * more accurately than, for example, with the standard QR method.
30 *
31 * The eigenvectors of a full or band positive definite Hermitian matrix
32 * can also be found if ZHETRD, ZHPTRD, or ZHBTRD has been used to
33 * reduce this matrix to tridiagonal form. (The reduction to
34 * tridiagonal form, however, may preclude the possibility of obtaining
35 * high relative accuracy in the small eigenvalues of the original
36 * matrix, if these eigenvalues range over many orders of magnitude.)
37 *
38 * Arguments
39 * =========
40 *
41 * COMPZ (input) CHARACTER*1
42 * = 'N': Compute eigenvalues only.
43 * = 'V': Compute eigenvectors of original Hermitian
44 * matrix also. Array Z contains the unitary matrix
45 * used to reduce the original matrix to tridiagonal
46 * form.
47 * = 'I': Compute eigenvectors of tridiagonal matrix also.
48 *
49 * N (input) INTEGER
50 * The order of the matrix. N >= 0.
51 *
52 * D (input/output) DOUBLE PRECISION array, dimension (N)
53 * On entry, the n diagonal elements of the tridiagonal 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) COMPLEX*16 array, dimension (LDZ, N)
63 * On entry, if COMPZ = 'V', the unitary matrix used in the
64 * reduction to tridiagonal form.
65 * On exit, if COMPZ = 'V', the orthonormal eigenvectors of the
66 * original Hermitian 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 COMPLEX*16 CZERO, CONE
94 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
95 $ CONE = ( 1.0D+0, 0.0D+0 ) )
96 * ..
97 * .. External Functions ..
98 LOGICAL LSAME
99 EXTERNAL LSAME
100 * ..
101 * .. External Subroutines ..
102 EXTERNAL DPTTRF, XERBLA, ZBDSQR, ZLASET
103 * ..
104 * .. Local Arrays ..
105 COMPLEX*16 C( 1, 1 ), VT( 1, 1 )
106 * ..
107 * .. Local Scalars ..
108 INTEGER I, ICOMPZ, NRU
109 * ..
110 * .. Intrinsic Functions ..
111 INTRINSIC MAX, SQRT
112 * ..
113 * .. Executable Statements ..
114 *
115 * Test the input parameters.
116 *
117 INFO = 0
118 *
119 IF( LSAME( COMPZ, 'N' ) ) THEN
120 ICOMPZ = 0
121 ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
122 ICOMPZ = 1
123 ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
124 ICOMPZ = 2
125 ELSE
126 ICOMPZ = -1
127 END IF
128 IF( ICOMPZ.LT.0 ) THEN
129 INFO = -1
130 ELSE IF( N.LT.0 ) THEN
131 INFO = -2
132 ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
133 $ N ) ) ) THEN
134 INFO = -6
135 END IF
136 IF( INFO.NE.0 ) THEN
137 CALL XERBLA( 'ZPTEQR', -INFO )
138 RETURN
139 END IF
140 *
141 * Quick return if possible
142 *
143 IF( N.EQ.0 )
144 $ RETURN
145 *
146 IF( N.EQ.1 ) THEN
147 IF( ICOMPZ.GT.0 )
148 $ Z( 1, 1 ) = CONE
149 RETURN
150 END IF
151 IF( ICOMPZ.EQ.2 )
152 $ CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
153 *
154 * Call DPTTRF to factor the matrix.
155 *
156 CALL DPTTRF( N, D, E, INFO )
157 IF( INFO.NE.0 )
158 $ RETURN
159 DO 10 I = 1, N
160 D( I ) = SQRT( D( I ) )
161 10 CONTINUE
162 DO 20 I = 1, N - 1
163 E( I ) = E( I )*D( I )
164 20 CONTINUE
165 *
166 * Call ZBDSQR to compute the singular values/vectors of the
167 * bidiagonal factor.
168 *
169 IF( ICOMPZ.GT.0 ) THEN
170 NRU = N
171 ELSE
172 NRU = 0
173 END IF
174 CALL ZBDSQR( 'Lower', N, 0, NRU, 0, D, E, VT, 1, Z, LDZ, C, 1,
175 $ WORK, INFO )
176 *
177 * Square the singular values.
178 *
179 IF( INFO.EQ.0 ) THEN
180 DO 30 I = 1, N
181 D( I ) = D( I )*D( I )
182 30 CONTINUE
183 ELSE
184 INFO = N + INFO
185 END IF
186 *
187 RETURN
188 *
189 * End of ZPTEQR
190 *
191 END