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