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