File:
[local] /
rpl /
lapack /
lapack /
dpteqr.f
Revision
1.7:
download - view:
text,
annotated -
select for diffs -
revision graph
Tue Dec 21 13:53:37 2010 UTC (13 years, 9 months ago) by
bertrand
Branches:
MAIN
CVS tags:
rpl-4_1_3,
rpl-4_1_2,
rpl-4_1_1,
rpl-4_1_0,
rpl-4_0_24,
rpl-4_0_22,
rpl-4_0_21,
rpl-4_0_20,
rpl-4_0,
HEAD
Mise à jour de lapack vers la version 3.3.0.
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
CVSweb interface <joel.bertrand@systella.fr>