```001:       SUBROUTINE DPTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
005: *     November 2006
006: *
007: *     .. Scalar Arguments ..
008:       CHARACTER          COMPZ
009:       INTEGER            INFO, LDZ, N
010: *     ..
011: *     .. Array Arguments ..
012:       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
013: *     ..
014: *
015: *  Purpose
016: *  =======
017: *
018: *  DPTEQR computes all eigenvalues and, optionally, eigenvectors of a
019: *  symmetric positive definite tridiagonal matrix by first factoring the
020: *  matrix using DPTTRF, and then calling DBDSQR to compute the singular
021: *  values of the bidiagonal factor.
022: *
023: *  This routine computes the eigenvalues of the positive definite
024: *  tridiagonal matrix to high relative accuracy.  This means that if the
025: *  eigenvalues range over many orders of magnitude in size, then the
026: *  small eigenvalues and corresponding eigenvectors will be computed
027: *  more accurately than, for example, with the standard QR method.
028: *
029: *  The eigenvectors of a full or band symmetric positive definite matrix
030: *  can also be found if DSYTRD, DSPTRD, or DSBTRD has been used to
031: *  reduce this matrix to tridiagonal form. (The reduction to tridiagonal
032: *  form, however, may preclude the possibility of obtaining high
033: *  relative accuracy in the small eigenvalues of the original matrix, if
034: *  these eigenvalues range over many orders of magnitude.)
035: *
036: *  Arguments
037: *  =========
038: *
039: *  COMPZ   (input) CHARACTER*1
040: *          = 'N':  Compute eigenvalues only.
041: *          = 'V':  Compute eigenvectors of original symmetric
042: *                  matrix also.  Array Z contains the orthogonal
043: *                  matrix used to reduce the original matrix to
044: *                  tridiagonal form.
045: *          = 'I':  Compute eigenvectors of tridiagonal matrix also.
046: *
047: *  N       (input) INTEGER
048: *          The order of the matrix.  N >= 0.
049: *
050: *  D       (input/output) DOUBLE PRECISION array, dimension (N)
051: *          On entry, the n diagonal elements of the tridiagonal
052: *          matrix.
053: *          On normal exit, D contains the eigenvalues, in descending
054: *          order.
055: *
056: *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
057: *          On entry, the (n-1) subdiagonal elements of the tridiagonal
058: *          matrix.
059: *          On exit, E has been destroyed.
060: *
061: *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
062: *          On entry, if COMPZ = 'V', the orthogonal matrix used in the
063: *          reduction to tridiagonal form.
064: *          On exit, if COMPZ = 'V', the orthonormal eigenvectors of the
065: *          original symmetric matrix;
066: *          if COMPZ = 'I', the orthonormal eigenvectors of the
067: *          tridiagonal matrix.
068: *          If INFO > 0 on exit, Z contains the eigenvectors associated
069: *          with only the stored eigenvalues.
070: *          If  COMPZ = 'N', then Z is not referenced.
071: *
072: *  LDZ     (input) INTEGER
073: *          The leading dimension of the array Z.  LDZ >= 1, and if
074: *          COMPZ = 'V' or 'I', LDZ >= max(1,N).
075: *
076: *  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)
077: *
078: *  INFO    (output) INTEGER
079: *          = 0:  successful exit.
080: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
081: *          > 0:  if INFO = i, and i is:
082: *                <= N  the Cholesky factorization of the matrix could
083: *                      not be performed because the i-th principal minor
084: *                      was not positive definite.
085: *                > N   the SVD algorithm failed to converge;
086: *                      if INFO = N+i, i off-diagonal elements of the
087: *                      bidiagonal factor did not converge to zero.
088: *
089: *  =====================================================================
090: *
091: *     .. Parameters ..
092:       DOUBLE PRECISION   ZERO, ONE
093:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
094: *     ..
095: *     .. External Functions ..
096:       LOGICAL            LSAME
097:       EXTERNAL           LSAME
098: *     ..
099: *     .. External Subroutines ..
100:       EXTERNAL           DBDSQR, DLASET, DPTTRF, XERBLA
101: *     ..
102: *     .. Local Arrays ..
103:       DOUBLE PRECISION   C( 1, 1 ), VT( 1, 1 )
104: *     ..
105: *     .. Local Scalars ..
106:       INTEGER            I, ICOMPZ, NRU
107: *     ..
108: *     .. Intrinsic Functions ..
109:       INTRINSIC          MAX, SQRT
110: *     ..
111: *     .. Executable Statements ..
112: *
113: *     Test the input parameters.
114: *
115:       INFO = 0
116: *
117:       IF( LSAME( COMPZ, 'N' ) ) THEN
118:          ICOMPZ = 0
119:       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
120:          ICOMPZ = 1
121:       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
122:          ICOMPZ = 2
123:       ELSE
124:          ICOMPZ = -1
125:       END IF
126:       IF( ICOMPZ.LT.0 ) THEN
127:          INFO = -1
128:       ELSE IF( N.LT.0 ) THEN
129:          INFO = -2
130:       ELSE IF( ( LDZ.LT.1 ) .OR. ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1,
131:      \$         N ) ) ) THEN
132:          INFO = -6
133:       END IF
134:       IF( INFO.NE.0 ) THEN
135:          CALL XERBLA( 'DPTEQR', -INFO )
136:          RETURN
137:       END IF
138: *
139: *     Quick return if possible
140: *
141:       IF( N.EQ.0 )
142:      \$   RETURN
143: *
144:       IF( N.EQ.1 ) THEN
145:          IF( ICOMPZ.GT.0 )
146:      \$      Z( 1, 1 ) = ONE
147:          RETURN
148:       END IF
149:       IF( ICOMPZ.EQ.2 )
150:      \$   CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
151: *
152: *     Call DPTTRF to factor the matrix.
153: *
154:       CALL DPTTRF( N, D, E, INFO )
155:       IF( INFO.NE.0 )
156:      \$   RETURN
157:       DO 10 I = 1, N
158:          D( I ) = SQRT( D( I ) )
159:    10 CONTINUE
160:       DO 20 I = 1, N - 1
161:          E( I ) = E( I )*D( I )
162:    20 CONTINUE
163: *
164: *     Call DBDSQR to compute the singular values/vectors of the
165: *     bidiagonal factor.
166: *
167:       IF( ICOMPZ.GT.0 ) THEN
168:          NRU = N
169:       ELSE
170:          NRU = 0
171:       END IF
172:       CALL DBDSQR( 'Lower', N, 0, NRU, 0, D, E, VT, 1, Z, LDZ, C, 1,
173:      \$             WORK, INFO )
174: *
175: *     Square the singular values.
176: *
177:       IF( INFO.EQ.0 ) THEN
178:          DO 30 I = 1, N
179:             D( I ) = D( I )*D( I )
180:    30    CONTINUE
181:       ELSE
182:          INFO = N + INFO
183:       END IF
184: *
185:       RETURN
186: *
187: *     End of DPTEQR
188: *
189:       END
190: ```