Lots of trailing whitespaces in the files of Syd. Cleaning this. No big deal.
[platform/upstream/lapack.git] / SRC / cheevx.f
1 *> \brief <b> CHEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices</b>
2 *
3 *  =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 *            http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CHEEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cheevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cheevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cheevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE CHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
22 *                          ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK,
23 *                          IWORK, IFAIL, INFO )
24 *
25 *       .. Scalar Arguments ..
26 *       CHARACTER          JOBZ, RANGE, UPLO
27 *       INTEGER            IL, INFO, IU, LDA, LDZ, LWORK, M, N
28 *       REAL               ABSTOL, VL, VU
29 *       ..
30 *       .. Array Arguments ..
31 *       INTEGER            IFAIL( * ), IWORK( * )
32 *       REAL               RWORK( * ), W( * )
33 *       COMPLEX            A( LDA, * ), WORK( * ), Z( LDZ, * )
34 *       ..
35 *
36 *
37 *> \par Purpose:
38 *  =============
39 *>
40 *> \verbatim
41 *>
42 *> CHEEVX computes selected eigenvalues and, optionally, eigenvectors
43 *> of a complex Hermitian matrix A.  Eigenvalues and eigenvectors can
44 *> be selected by specifying either a range of values or a range of
45 *> indices for the desired eigenvalues.
46 *> \endverbatim
47 *
48 *  Arguments:
49 *  ==========
50 *
51 *> \param[in] JOBZ
52 *> \verbatim
53 *>          JOBZ is CHARACTER*1
54 *>          = 'N':  Compute eigenvalues only;
55 *>          = 'V':  Compute eigenvalues and eigenvectors.
56 *> \endverbatim
57 *>
58 *> \param[in] RANGE
59 *> \verbatim
60 *>          RANGE is CHARACTER*1
61 *>          = 'A': all eigenvalues will be found.
62 *>          = 'V': all eigenvalues in the half-open interval (VL,VU]
63 *>                 will be found.
64 *>          = 'I': the IL-th through IU-th eigenvalues will be found.
65 *> \endverbatim
66 *>
67 *> \param[in] UPLO
68 *> \verbatim
69 *>          UPLO is CHARACTER*1
70 *>          = 'U':  Upper triangle of A is stored;
71 *>          = 'L':  Lower triangle of A is stored.
72 *> \endverbatim
73 *>
74 *> \param[in] N
75 *> \verbatim
76 *>          N is INTEGER
77 *>          The order of the matrix A.  N >= 0.
78 *> \endverbatim
79 *>
80 *> \param[in,out] A
81 *> \verbatim
82 *>          A is COMPLEX array, dimension (LDA, N)
83 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the
84 *>          leading N-by-N upper triangular part of A contains the
85 *>          upper triangular part of the matrix A.  If UPLO = 'L',
86 *>          the leading N-by-N lower triangular part of A contains
87 *>          the lower triangular part of the matrix A.
88 *>          On exit, the lower triangle (if UPLO='L') or the upper
89 *>          triangle (if UPLO='U') of A, including the diagonal, is
90 *>          destroyed.
91 *> \endverbatim
92 *>
93 *> \param[in] LDA
94 *> \verbatim
95 *>          LDA is INTEGER
96 *>          The leading dimension of the array A.  LDA >= max(1,N).
97 *> \endverbatim
98 *>
99 *> \param[in] VL
100 *> \verbatim
101 *>          VL is REAL
102 *>          If RANGE='V', the lower bound of the interval to
103 *>          be searched for eigenvalues. VL < VU.
104 *>          Not referenced if RANGE = 'A' or 'I'.
105 *> \endverbatim
106 *>
107 *> \param[in] VU
108 *> \verbatim
109 *>          VU is REAL
110 *>          If RANGE='V', the upper bound of the interval to
111 *>          be searched for eigenvalues. VL < VU.
112 *>          Not referenced if RANGE = 'A' or 'I'.
113 *> \endverbatim
114 *>
115 *> \param[in] IL
116 *> \verbatim
117 *>          IL is INTEGER
118 *>          If RANGE='I', the index of the
119 *>          smallest eigenvalue to be returned.
120 *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
121 *>          Not referenced if RANGE = 'A' or 'V'.
122 *> \endverbatim
123 *>
124 *> \param[in] IU
125 *> \verbatim
126 *>          IU is INTEGER
127 *>          If RANGE='I', the index of the
128 *>          largest eigenvalue to be returned.
129 *>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
130 *>          Not referenced if RANGE = 'A' or 'V'.
131 *> \endverbatim
132 *>
133 *> \param[in] ABSTOL
134 *> \verbatim
135 *>          ABSTOL is REAL
136 *>          The absolute error tolerance for the eigenvalues.
137 *>          An approximate eigenvalue is accepted as converged
138 *>          when it is determined to lie in an interval [a,b]
139 *>          of width less than or equal to
140 *>
141 *>                  ABSTOL + EPS *   max( |a|,|b| ) ,
142 *>
143 *>          where EPS is the machine precision.  If ABSTOL is less than
144 *>          or equal to zero, then  EPS*|T|  will be used in its place,
145 *>          where |T| is the 1-norm of the tridiagonal matrix obtained
146 *>          by reducing A to tridiagonal form.
147 *>
148 *>          Eigenvalues will be computed most accurately when ABSTOL is
149 *>          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
150 *>          If this routine returns with INFO>0, indicating that some
151 *>          eigenvectors did not converge, try setting ABSTOL to
152 *>          2*SLAMCH('S').
153 *>
154 *>          See "Computing Small Singular Values of Bidiagonal Matrices
155 *>          with Guaranteed High Relative Accuracy," by Demmel and
156 *>          Kahan, LAPACK Working Note #3.
157 *> \endverbatim
158 *>
159 *> \param[out] M
160 *> \verbatim
161 *>          M is INTEGER
162 *>          The total number of eigenvalues found.  0 <= M <= N.
163 *>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
164 *> \endverbatim
165 *>
166 *> \param[out] W
167 *> \verbatim
168 *>          W is REAL array, dimension (N)
169 *>          On normal exit, the first M elements contain the selected
170 *>          eigenvalues in ascending order.
171 *> \endverbatim
172 *>
173 *> \param[out] Z
174 *> \verbatim
175 *>          Z is COMPLEX array, dimension (LDZ, max(1,M))
176 *>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
177 *>          contain the orthonormal eigenvectors of the matrix A
178 *>          corresponding to the selected eigenvalues, with the i-th
179 *>          column of Z holding the eigenvector associated with W(i).
180 *>          If an eigenvector fails to converge, then that column of Z
181 *>          contains the latest approximation to the eigenvector, and the
182 *>          index of the eigenvector is returned in IFAIL.
183 *>          If JOBZ = 'N', then Z is not referenced.
184 *>          Note: the user must ensure that at least max(1,M) columns are
185 *>          supplied in the array Z; if RANGE = 'V', the exact value of M
186 *>          is not known in advance and an upper bound must be used.
187 *> \endverbatim
188 *>
189 *> \param[in] LDZ
190 *> \verbatim
191 *>          LDZ is INTEGER
192 *>          The leading dimension of the array Z.  LDZ >= 1, and if
193 *>          JOBZ = 'V', LDZ >= max(1,N).
194 *> \endverbatim
195 *>
196 *> \param[out] WORK
197 *> \verbatim
198 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
199 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
200 *> \endverbatim
201 *>
202 *> \param[in] LWORK
203 *> \verbatim
204 *>          LWORK is INTEGER
205 *>          The length of the array WORK.  LWORK >= 1, when N <= 1;
206 *>          otherwise 2*N.
207 *>          For optimal efficiency, LWORK >= (NB+1)*N,
208 *>          where NB is the max of the blocksize for CHETRD and for
209 *>          CUNMTR as returned by ILAENV.
210 *>
211 *>          If LWORK = -1, then a workspace query is assumed; the routine
212 *>          only calculates the optimal size of the WORK array, returns
213 *>          this value as the first entry of the WORK array, and no error
214 *>          message related to LWORK is issued by XERBLA.
215 *> \endverbatim
216 *>
217 *> \param[out] RWORK
218 *> \verbatim
219 *>          RWORK is REAL array, dimension (7*N)
220 *> \endverbatim
221 *>
222 *> \param[out] IWORK
223 *> \verbatim
224 *>          IWORK is INTEGER array, dimension (5*N)
225 *> \endverbatim
226 *>
227 *> \param[out] IFAIL
228 *> \verbatim
229 *>          IFAIL is INTEGER array, dimension (N)
230 *>          If JOBZ = 'V', then if INFO = 0, the first M elements of
231 *>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
232 *>          indices of the eigenvectors that failed to converge.
233 *>          If JOBZ = 'N', then IFAIL is not referenced.
234 *> \endverbatim
235 *>
236 *> \param[out] INFO
237 *> \verbatim
238 *>          INFO is INTEGER
239 *>          = 0:  successful exit
240 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
241 *>          > 0:  if INFO = i, then i eigenvectors failed to converge.
242 *>                Their indices are stored in array IFAIL.
243 *> \endverbatim
244 *
245 *  Authors:
246 *  ========
247 *
248 *> \author Univ. of Tennessee
249 *> \author Univ. of California Berkeley
250 *> \author Univ. of Colorado Denver
251 *> \author NAG Ltd.
252 *
253 *> \date June 2016
254 *
255 *> \ingroup complexHEeigen
256 *
257 *  =====================================================================
258       SUBROUTINE CHEEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU,
259      $                   ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK,
260      $                   IWORK, IFAIL, INFO )
261 *
262 *  -- LAPACK driver routine (version 3.6.1) --
263 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
264 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
265 *     June 2016
266 *
267 *     .. Scalar Arguments ..
268       CHARACTER          JOBZ, RANGE, UPLO
269       INTEGER            IL, INFO, IU, LDA, LDZ, LWORK, M, N
270       REAL               ABSTOL, VL, VU
271 *     ..
272 *     .. Array Arguments ..
273       INTEGER            IFAIL( * ), IWORK( * )
274       REAL               RWORK( * ), W( * )
275       COMPLEX            A( LDA, * ), WORK( * ), Z( LDZ, * )
276 *     ..
277 *
278 *  =====================================================================
279 *
280 *     .. Parameters ..
281       REAL               ZERO, ONE
282       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
283       COMPLEX            CONE
284       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
285 *     ..
286 *     .. Local Scalars ..
287       LOGICAL            ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
288      $                   WANTZ
289       CHARACTER          ORDER
290       INTEGER            I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
291      $                   INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
292      $                   ITMP1, J, JJ, LLWORK, LWKMIN, LWKOPT, NB,
293      $                   NSPLIT
294       REAL               ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
295      $                   SIGMA, SMLNUM, TMP1, VLL, VUU
296 *     ..
297 *     .. External Functions ..
298       LOGICAL            LSAME
299       INTEGER            ILAENV
300       REAL               SLAMCH, CLANHE
301       EXTERNAL           LSAME, ILAENV, SLAMCH, CLANHE
302 *     ..
303 *     .. External Subroutines ..
304       EXTERNAL           SCOPY, SSCAL, SSTEBZ, SSTERF, XERBLA, CSSCAL,
305      $                   CHETRD, CLACPY, CSTEIN, CSTEQR, CSWAP, CUNGTR,
306      $                   CUNMTR
307 *     ..
308 *     .. Intrinsic Functions ..
309       INTRINSIC          REAL, MAX, MIN, SQRT
310 *     ..
311 *     .. Executable Statements ..
312 *
313 *     Test the input parameters.
314 *
315       LOWER = LSAME( UPLO, 'L' )
316       WANTZ = LSAME( JOBZ, 'V' )
317       ALLEIG = LSAME( RANGE, 'A' )
318       VALEIG = LSAME( RANGE, 'V' )
319       INDEIG = LSAME( RANGE, 'I' )
320       LQUERY = ( LWORK.EQ.-1 )
321 *
322       INFO = 0
323       IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
324          INFO = -1
325       ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
326          INFO = -2
327       ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
328          INFO = -3
329       ELSE IF( N.LT.0 ) THEN
330          INFO = -4
331       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
332          INFO = -6
333       ELSE
334          IF( VALEIG ) THEN
335             IF( N.GT.0 .AND. VU.LE.VL )
336      $         INFO = -8
337          ELSE IF( INDEIG ) THEN
338             IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
339                INFO = -9
340             ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
341                INFO = -10
342             END IF
343          END IF
344       END IF
345       IF( INFO.EQ.0 ) THEN
346          IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
347             INFO = -15
348          END IF
349       END IF
350 *
351       IF( INFO.EQ.0 ) THEN
352          IF( N.LE.1 ) THEN
353             LWKMIN = 1
354             WORK( 1 ) = LWKMIN
355          ELSE
356             LWKMIN = 2*N
357             NB = ILAENV( 1, 'CHETRD', UPLO, N, -1, -1, -1 )
358             NB = MAX( NB, ILAENV( 1, 'CUNMTR', UPLO, N, -1, -1, -1 ) )
359             LWKOPT = MAX( 1, ( NB + 1 )*N )
360             WORK( 1 ) = LWKOPT
361          END IF
362 *
363          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
364      $      INFO = -17
365       END IF
366 *
367       IF( INFO.NE.0 ) THEN
368          CALL XERBLA( 'CHEEVX', -INFO )
369          RETURN
370       ELSE IF( LQUERY ) THEN
371          RETURN
372       END IF
373 *
374 *     Quick return if possible
375 *
376       M = 0
377       IF( N.EQ.0 ) THEN
378          RETURN
379       END IF
380 *
381       IF( N.EQ.1 ) THEN
382          IF( ALLEIG .OR. INDEIG ) THEN
383             M = 1
384             W( 1 ) = A( 1, 1 )
385          ELSE IF( VALEIG ) THEN
386             IF( VL.LT.REAL( A( 1, 1 ) ) .AND. VU.GE.REAL( A( 1, 1 ) ) )
387      $           THEN
388                M = 1
389                W( 1 ) = A( 1, 1 )
390             END IF
391          END IF
392          IF( WANTZ )
393      $      Z( 1, 1 ) = CONE
394          RETURN
395       END IF
396 *
397 *     Get machine constants.
398 *
399       SAFMIN = SLAMCH( 'Safe minimum' )
400       EPS = SLAMCH( 'Precision' )
401       SMLNUM = SAFMIN / EPS
402       BIGNUM = ONE / SMLNUM
403       RMIN = SQRT( SMLNUM )
404       RMAX = MIN( SQRT( BIGNUM ), ONE / SQRT( SQRT( SAFMIN ) ) )
405 *
406 *     Scale matrix to allowable range, if necessary.
407 *
408       ISCALE = 0
409       ABSTLL = ABSTOL
410       IF( VALEIG ) THEN
411          VLL = VL
412          VUU = VU
413       END IF
414       ANRM = CLANHE( 'M', UPLO, N, A, LDA, RWORK )
415       IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
416          ISCALE = 1
417          SIGMA = RMIN / ANRM
418       ELSE IF( ANRM.GT.RMAX ) THEN
419          ISCALE = 1
420          SIGMA = RMAX / ANRM
421       END IF
422       IF( ISCALE.EQ.1 ) THEN
423          IF( LOWER ) THEN
424             DO 10 J = 1, N
425                CALL CSSCAL( N-J+1, SIGMA, A( J, J ), 1 )
426    10       CONTINUE
427          ELSE
428             DO 20 J = 1, N
429                CALL CSSCAL( J, SIGMA, A( 1, J ), 1 )
430    20       CONTINUE
431          END IF
432          IF( ABSTOL.GT.0 )
433      $      ABSTLL = ABSTOL*SIGMA
434          IF( VALEIG ) THEN
435             VLL = VL*SIGMA
436             VUU = VU*SIGMA
437          END IF
438       END IF
439 *
440 *     Call CHETRD to reduce Hermitian matrix to tridiagonal form.
441 *
442       INDD = 1
443       INDE = INDD + N
444       INDRWK = INDE + N
445       INDTAU = 1
446       INDWRK = INDTAU + N
447       LLWORK = LWORK - INDWRK + 1
448       CALL CHETRD( UPLO, N, A, LDA, RWORK( INDD ), RWORK( INDE ),
449      $             WORK( INDTAU ), WORK( INDWRK ), LLWORK, IINFO )
450 *
451 *     If all eigenvalues are desired and ABSTOL is less than or equal to
452 *     zero, then call SSTERF or CUNGTR and CSTEQR.  If this fails for
453 *     some eigenvalue, then try SSTEBZ.
454 *
455       TEST = .FALSE.
456       IF( INDEIG ) THEN
457          IF( IL.EQ.1 .AND. IU.EQ.N ) THEN
458             TEST = .TRUE.
459          END IF
460       END IF
461       IF( ( ALLEIG .OR. TEST ) .AND. ( ABSTOL.LE.ZERO ) ) THEN
462          CALL SCOPY( N, RWORK( INDD ), 1, W, 1 )
463          INDEE = INDRWK + 2*N
464          IF( .NOT.WANTZ ) THEN
465             CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 )
466             CALL SSTERF( N, W, RWORK( INDEE ), INFO )
467          ELSE
468             CALL CLACPY( 'A', N, N, A, LDA, Z, LDZ )
469             CALL CUNGTR( UPLO, N, Z, LDZ, WORK( INDTAU ),
470      $                   WORK( INDWRK ), LLWORK, IINFO )
471             CALL SCOPY( N-1, RWORK( INDE ), 1, RWORK( INDEE ), 1 )
472             CALL CSTEQR( JOBZ, N, W, RWORK( INDEE ), Z, LDZ,
473      $                   RWORK( INDRWK ), INFO )
474             IF( INFO.EQ.0 ) THEN
475                DO 30 I = 1, N
476                   IFAIL( I ) = 0
477    30          CONTINUE
478             END IF
479          END IF
480          IF( INFO.EQ.0 ) THEN
481             M = N
482             GO TO 40
483          END IF
484          INFO = 0
485       END IF
486 *
487 *     Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
488 *
489       IF( WANTZ ) THEN
490          ORDER = 'B'
491       ELSE
492          ORDER = 'E'
493       END IF
494       INDIBL = 1
495       INDISP = INDIBL + N
496       INDIWK = INDISP + N
497       CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL,
498      $             RWORK( INDD ), RWORK( INDE ), M, NSPLIT, W,
499      $             IWORK( INDIBL ), IWORK( INDISP ), RWORK( INDRWK ),
500      $             IWORK( INDIWK ), INFO )
501 *
502       IF( WANTZ ) THEN
503          CALL CSTEIN( N, RWORK( INDD ), RWORK( INDE ), M, W,
504      $                IWORK( INDIBL ), IWORK( INDISP ), Z, LDZ,
505      $                RWORK( INDRWK ), IWORK( INDIWK ), IFAIL, INFO )
506 *
507 *        Apply unitary matrix used in reduction to tridiagonal
508 *        form to eigenvectors returned by CSTEIN.
509 *
510          CALL CUNMTR( 'L', UPLO, 'N', N, M, A, LDA, WORK( INDTAU ), Z,
511      $                LDZ, WORK( INDWRK ), LLWORK, IINFO )
512       END IF
513 *
514 *     If matrix was scaled, then rescale eigenvalues appropriately.
515 *
516    40 CONTINUE
517       IF( ISCALE.EQ.1 ) THEN
518          IF( INFO.EQ.0 ) THEN
519             IMAX = M
520          ELSE
521             IMAX = INFO - 1
522          END IF
523          CALL SSCAL( IMAX, ONE / SIGMA, W, 1 )
524       END IF
525 *
526 *     If eigenvalues are not in order, then sort them, along with
527 *     eigenvectors.
528 *
529       IF( WANTZ ) THEN
530          DO 60 J = 1, M - 1
531             I = 0
532             TMP1 = W( J )
533             DO 50 JJ = J + 1, M
534                IF( W( JJ ).LT.TMP1 ) THEN
535                   I = JJ
536                   TMP1 = W( JJ )
537                END IF
538    50       CONTINUE
539 *
540             IF( I.NE.0 ) THEN
541                ITMP1 = IWORK( INDIBL+I-1 )
542                W( I ) = W( J )
543                IWORK( INDIBL+I-1 ) = IWORK( INDIBL+J-1 )
544                W( J ) = TMP1
545                IWORK( INDIBL+J-1 ) = ITMP1
546                CALL CSWAP( N, Z( 1, I ), 1, Z( 1, J ), 1 )
547                IF( INFO.NE.0 ) THEN
548                   ITMP1 = IFAIL( I )
549                   IFAIL( I ) = IFAIL( J )
550                   IFAIL( J ) = ITMP1
551                END IF
552             END IF
553    60    CONTINUE
554       END IF
555 *
556 *     Set WORK(1) to optimal complex workspace size.
557 *
558       WORK( 1 ) = LWKOPT
559 *
560       RETURN
561 *
562 *     End of CHEEVX
563 *
564       END