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