001:       SUBROUTINE CPTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          COMPZ
010:       INTEGER            INFO, LDZ, N
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               D( * ), E( * ), WORK( * )
014:       COMPLEX            Z( LDZ, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  CPTEQR computes all eigenvalues and, optionally, eigenvectors of a
021: *  symmetric positive definite tridiagonal matrix by first factoring the
022: *  matrix using SPTTRF and then calling CBDSQR to compute the singular
023: *  values of the bidiagonal factor.
024: *
025: *  This routine computes the eigenvalues of the positive definite
026: *  tridiagonal matrix to high relative accuracy.  This means that if the
027: *  eigenvalues range over many orders of magnitude in size, then the
028: *  small eigenvalues and corresponding eigenvectors will be computed
029: *  more accurately than, for example, with the standard QR method.
030: *
031: *  The eigenvectors of a full or band positive definite Hermitian matrix
032: *  can also be found if CHETRD, CHPTRD, or CHBTRD has been used to
033: *  reduce this matrix to tridiagonal form.  (The reduction to
034: *  tridiagonal form, however, may preclude the possibility of obtaining
035: *  high relative accuracy in the small eigenvalues of the original
036: *  matrix, if these eigenvalues range over many orders of magnitude.)
037: *
038: *  Arguments
039: *  =========
040: *
041: *  COMPZ   (input) CHARACTER*1
042: *          = 'N':  Compute eigenvalues only.
043: *          = 'V':  Compute eigenvectors of original Hermitian
044: *                  matrix also.  Array Z contains the unitary matrix
045: *                  used to reduce the original matrix to tridiagonal
046: *                  form.
047: *          = 'I':  Compute eigenvectors of tridiagonal matrix also.
048: *
049: *  N       (input) INTEGER
050: *          The order of the matrix.  N >= 0.
051: *
052: *  D       (input/output) REAL array, dimension (N)
053: *          On entry, the n diagonal elements of the tridiagonal matrix.
054: *          On normal exit, D contains the eigenvalues, in descending
055: *          order.
056: *
057: *  E       (input/output) REAL array, dimension (N-1)
058: *          On entry, the (n-1) subdiagonal elements of the tridiagonal
059: *          matrix.
060: *          On exit, E has been destroyed.
061: *
062: *  Z       (input/output) COMPLEX array, dimension (LDZ, N)
063: *          On entry, if COMPZ = 'V', the unitary matrix used in the
064: *          reduction to tridiagonal form.
065: *          On exit, if COMPZ = 'V', the orthonormal eigenvectors of the
066: *          original Hermitian matrix;
067: *          if COMPZ = 'I', the orthonormal eigenvectors of the
068: *          tridiagonal matrix.
069: *          If INFO > 0 on exit, Z contains the eigenvectors associated
070: *          with only the stored eigenvalues.
071: *          If  COMPZ = 'N', then Z is not referenced.
072: *
073: *  LDZ     (input) INTEGER
074: *          The leading dimension of the array Z.  LDZ >= 1, and if
075: *          COMPZ = 'V' or 'I', LDZ >= max(1,N).
076: *
077: *  WORK    (workspace) REAL array, dimension (4*N)
078: *
079: *  INFO    (output) INTEGER
080: *          = 0:  successful exit.
081: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
082: *          > 0:  if INFO = i, and i is:
083: *                <= N  the Cholesky factorization of the matrix could
084: *                      not be performed because the i-th principal minor
085: *                      was not positive definite.
086: *                > N   the SVD algorithm failed to converge;
087: *                      if INFO = N+i, i off-diagonal elements of the
088: *                      bidiagonal factor did not converge to zero.
089: *
090: *  ====================================================================
091: *
092: *     .. Parameters ..
093:       COMPLEX            CZERO, CONE
094:       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
095:      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
096: *     ..
097: *     .. External Functions ..
098:       LOGICAL            LSAME
099:       EXTERNAL           LSAME
100: *     ..
101: *     .. External Subroutines ..
102:       EXTERNAL           CBDSQR, CLASET, SPTTRF, XERBLA
103: *     ..
104: *     .. Local Arrays ..
105:       COMPLEX            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( 'CPTEQR', -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 CLASET( 'Full', N, N, CZERO, CONE, Z, LDZ )
153: *
154: *     Call SPTTRF to factor the matrix.
155: *
156:       CALL SPTTRF( 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 CBDSQR 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 CBDSQR( '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 CPTEQR
190: *
191:       END
192: