Mercurial > octave-nkf

comparison libcruft/lapack/ssyev.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 SUBROUTINE SSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
2 *
3 * -- LAPACK driver routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 CHARACTER JOBZ, UPLO
9 INTEGER INFO, LDA, LWORK, N
10 * ..
11 * .. Array Arguments ..
12 REAL A( LDA, * ), W( * ), WORK( * )
13 * ..
14 *
15 * Purpose
16 * =======
17 *
18 * SSYEV computes all eigenvalues and, optionally, eigenvectors of a
19 * real symmetric matrix A.
20 *
21 * Arguments
22 * =========
23 *
24 * JOBZ (input) CHARACTER*1
25 * = 'N': Compute eigenvalues only;
26 * = 'V': Compute eigenvalues and eigenvectors.
27 *
28 * UPLO (input) CHARACTER*1
29 * = 'U': Upper triangle of A is stored;
30 * = 'L': Lower triangle of A is stored.
31 *
32 * N (input) INTEGER
33 * The order of the matrix A. N >= 0.
34 *
35 * A (input/output) REAL array, dimension (LDA, N)
36 * On entry, the symmetric matrix A. If UPLO = 'U', the
37 * leading N-by-N upper triangular part of A contains the
38 * upper triangular part of the matrix A. If UPLO = 'L',
39 * the leading N-by-N lower triangular part of A contains
40 * the lower triangular part of the matrix A.
41 * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
42 * orthonormal eigenvectors of the matrix A.
43 * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
44 * or the upper triangle (if UPLO='U') of A, including the
45 * diagonal, is destroyed.
46 *
47 * LDA (input) INTEGER
48 * The leading dimension of the array A. LDA >= max(1,N).
49 *
50 * W (output) REAL array, dimension (N)
51 * If INFO = 0, the eigenvalues in ascending order.
52 *
53 * WORK (workspace/output) REAL array, dimension (MAX(1,LWORK))
54 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
55 *
56 * LWORK (input) INTEGER
57 * The length of the array WORK. LWORK >= max(1,3*N-1).
58 * For optimal efficiency, LWORK >= (NB+2)*N,
59 * where NB is the blocksize for SSYTRD returned by ILAENV.
60 *
61 * If LWORK = -1, then a workspace query is assumed; the routine
62 * only calculates the optimal size of the WORK array, returns
63 * this value as the first entry of the WORK array, and no error
64 * message related to LWORK is issued by XERBLA.
65 *
66 * INFO (output) INTEGER
67 * = 0: successful exit
68 * < 0: if INFO = -i, the i-th argument had an illegal value
69 * > 0: if INFO = i, the algorithm failed to converge; i
70 * off-diagonal elements of an intermediate tridiagonal
71 * form did not converge to zero.
72 *
73 * =====================================================================
74 *
75 * .. Parameters ..
76 REAL ZERO, ONE
77 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
78 * ..
79 * .. Local Scalars ..
80 LOGICAL LOWER, LQUERY, WANTZ
81 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
82 $ LLWORK, LWKOPT, NB
83 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
84 $ SMLNUM
85 * ..
86 * .. External Functions ..
87 LOGICAL LSAME
88 INTEGER ILAENV
89 REAL SLAMCH, SLANSY
90 EXTERNAL ILAENV, LSAME, SLAMCH, SLANSY
91 * ..
92 * .. External Subroutines ..
93 EXTERNAL SLASCL, SORGTR, SSCAL, SSTEQR, SSTERF, SSYTRD,
94 $ XERBLA
95 * ..
96 * .. Intrinsic Functions ..
97 INTRINSIC MAX, SQRT
98 * ..
99 * .. Executable Statements ..
100 *
101 * Test the input parameters.
102 *
103 WANTZ = LSAME( JOBZ, 'V' )
104 LOWER = LSAME( UPLO, 'L' )
105 LQUERY = ( LWORK.EQ.-1 )
106 *
107 INFO = 0
108 IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
109 INFO = -1
110 ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
111 INFO = -2
112 ELSE IF( N.LT.0 ) THEN
113 INFO = -3
114 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
115 INFO = -5
116 END IF
117 *
118 IF( INFO.EQ.0 ) THEN
119 NB = ILAENV( 1, 'SSYTRD', UPLO, N, -1, -1, -1 )
120 LWKOPT = MAX( 1, ( NB+2 )*N )
121 WORK( 1 ) = LWKOPT
122 *
123 IF( LWORK.LT.MAX( 1, 3*N-1 ) .AND. .NOT.LQUERY )
124 $ INFO = -8
125 END IF
126 *
127 IF( INFO.NE.0 ) THEN
128 CALL XERBLA( 'SSYEV ', -INFO )
129 RETURN
130 ELSE IF( LQUERY ) THEN
131 RETURN
132 END IF
133 *
134 * Quick return if possible
135 *
136 IF( N.EQ.0 ) THEN
137 RETURN
138 END IF
139 *
140 IF( N.EQ.1 ) THEN
141 W( 1 ) = A( 1, 1 )
142 WORK( 1 ) = 2
143 IF( WANTZ )
144 $ A( 1, 1 ) = ONE
145 RETURN
146 END IF
147 *
148 * Get machine constants.
149 *
150 SAFMIN = SLAMCH( 'Safe minimum' )
151 EPS = SLAMCH( 'Precision' )
152 SMLNUM = SAFMIN / EPS
153 BIGNUM = ONE / SMLNUM
154 RMIN = SQRT( SMLNUM )
155 RMAX = SQRT( BIGNUM )
156 *
157 * Scale matrix to allowable range, if necessary.
158 *
159 ANRM = SLANSY( 'M', UPLO, N, A, LDA, WORK )
160 ISCALE = 0
161 IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
162 ISCALE = 1
163 SIGMA = RMIN / ANRM
164 ELSE IF( ANRM.GT.RMAX ) THEN
165 ISCALE = 1
166 SIGMA = RMAX / ANRM
167 END IF
168 IF( ISCALE.EQ.1 )
169 $ CALL SLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
170 *
171 * Call SSYTRD to reduce symmetric matrix to tridiagonal form.
172 *
173 INDE = 1
174 INDTAU = INDE + N
175 INDWRK = INDTAU + N
176 LLWORK = LWORK - INDWRK + 1
177 CALL SSYTRD( UPLO, N, A, LDA, W, WORK( INDE ), WORK( INDTAU ),
178 $ WORK( INDWRK ), LLWORK, IINFO )
179 *
180 * For eigenvalues only, call SSTERF. For eigenvectors, first call
181 * SORGTR to generate the orthogonal matrix, then call SSTEQR.
182 *
183 IF( .NOT.WANTZ ) THEN
184 CALL SSTERF( N, W, WORK( INDE ), INFO )
185 ELSE
186 CALL SORGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
187 $ LLWORK, IINFO )
188 CALL SSTEQR( JOBZ, N, W, WORK( INDE ), A, LDA, WORK( INDTAU ),
189 $ INFO )
190 END IF
191 *
192 * If matrix was scaled, then rescale eigenvalues appropriately.
193 *
194 IF( ISCALE.EQ.1 ) THEN
195 IF( INFO.EQ.0 ) THEN
196 IMAX = N
197 ELSE
198 IMAX = INFO - 1
199 END IF
200 CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
201 END IF
202 *
203 * Set WORK(1) to optimal workspace size.
204 *
205 WORK( 1 ) = LWKOPT
206 *
207 RETURN
208 *
209 * End of SSYEV
210 *
211 END