1*> \brief \b DSYGVX
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSYGVX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsygvx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsygvx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsygvx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB,
22*                          VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK,
23*                          LWORK, IWORK, IFAIL, INFO )
24*
25*       .. Scalar Arguments ..
26*       CHARACTER          JOBZ, RANGE, UPLO
27*       INTEGER            IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
28*       DOUBLE PRECISION   ABSTOL, VL, VU
29*       ..
30*       .. Array Arguments ..
31*       INTEGER            IFAIL( * ), IWORK( * )
32*       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), W( * ), WORK( * ),
33*      $                   Z( LDZ, * )
34*       ..
35*
36*
37*> \par Purpose:
38*  =============
39*>
40*> \verbatim
41*>
42*> DSYGVX computes selected eigenvalues, and optionally, eigenvectors
43*> of a real generalized symmetric-definite eigenproblem, of the form
44*> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A
45*> and B are assumed to be symmetric and B is also positive definite.
46*> Eigenvalues and eigenvectors can be selected by specifying either a
47*> range of values or a range of indices for the desired eigenvalues.
48*> \endverbatim
49*
50*  Arguments:
51*  ==========
52*
53*> \param[in] ITYPE
54*> \verbatim
55*>          ITYPE is INTEGER
56*>          Specifies the problem type to be solved:
57*>          = 1:  A*x = (lambda)*B*x
58*>          = 2:  A*B*x = (lambda)*x
59*>          = 3:  B*A*x = (lambda)*x
60*> \endverbatim
61*>
62*> \param[in] JOBZ
63*> \verbatim
64*>          JOBZ is CHARACTER*1
65*>          = 'N':  Compute eigenvalues only;
66*>          = 'V':  Compute eigenvalues and eigenvectors.
67*> \endverbatim
68*>
69*> \param[in] RANGE
70*> \verbatim
71*>          RANGE is CHARACTER*1
72*>          = 'A': all eigenvalues will be found.
73*>          = 'V': all eigenvalues in the half-open interval (VL,VU]
74*>                 will be found.
75*>          = 'I': the IL-th through IU-th eigenvalues will be found.
76*> \endverbatim
77*>
78*> \param[in] UPLO
79*> \verbatim
80*>          UPLO is CHARACTER*1
81*>          = 'U':  Upper triangle of A and B are stored;
82*>          = 'L':  Lower triangle of A and B are stored.
83*> \endverbatim
84*>
85*> \param[in] N
86*> \verbatim
87*>          N is INTEGER
88*>          The order of the matrix pencil (A,B).  N >= 0.
89*> \endverbatim
90*>
91*> \param[in,out] A
92*> \verbatim
93*>          A is DOUBLE PRECISION array, dimension (LDA, N)
94*>          On entry, the symmetric matrix A.  If UPLO = 'U', the
95*>          leading N-by-N upper triangular part of A contains the
96*>          upper triangular part of the matrix A.  If UPLO = 'L',
97*>          the leading N-by-N lower triangular part of A contains
98*>          the lower triangular part of the matrix A.
99*>
100*>          On exit, the lower triangle (if UPLO='L') or the upper
101*>          triangle (if UPLO='U') of A, including the diagonal, is
102*>          destroyed.
103*> \endverbatim
104*>
105*> \param[in] LDA
106*> \verbatim
107*>          LDA is INTEGER
108*>          The leading dimension of the array A.  LDA >= max(1,N).
109*> \endverbatim
110*>
111*> \param[in,out] B
112*> \verbatim
113*>          B is DOUBLE PRECISION array, dimension (LDB, N)
114*>          On entry, the symmetric matrix B.  If UPLO = 'U', the
115*>          leading N-by-N upper triangular part of B contains the
116*>          upper triangular part of the matrix B.  If UPLO = 'L',
117*>          the leading N-by-N lower triangular part of B contains
118*>          the lower triangular part of the matrix B.
119*>
120*>          On exit, if INFO <= N, the part of B containing the matrix is
121*>          overwritten by the triangular factor U or L from the Cholesky
122*>          factorization B = U**T*U or B = L*L**T.
123*> \endverbatim
124*>
125*> \param[in] LDB
126*> \verbatim
127*>          LDB is INTEGER
128*>          The leading dimension of the array B.  LDB >= max(1,N).
129*> \endverbatim
130*>
131*> \param[in] VL
132*> \verbatim
133*>          VL is DOUBLE PRECISION
134*> \endverbatim
135*>
136*> \param[in] VU
137*> \verbatim
138*>          VU is DOUBLE PRECISION
139*>          If RANGE='V', the lower and upper bounds of the interval to
140*>          be searched for eigenvalues. VL < VU.
141*>          Not referenced if RANGE = 'A' or 'I'.
142*> \endverbatim
143*>
144*> \param[in] IL
145*> \verbatim
146*>          IL is INTEGER
147*> \endverbatim
148*>
149*> \param[in] IU
150*> \verbatim
151*>          IU is INTEGER
152*>          If RANGE='I', the indices (in ascending order) of the
153*>          smallest and largest eigenvalues to be returned.
154*>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
155*>          Not referenced if RANGE = 'A' or 'V'.
156*> \endverbatim
157*>
158*> \param[in] ABSTOL
159*> \verbatim
160*>          ABSTOL is DOUBLE PRECISION
161*>          The absolute error tolerance for the eigenvalues.
162*>          An approximate eigenvalue is accepted as converged
163*>          when it is determined to lie in an interval [a,b]
164*>          of width less than or equal to
165*>
166*>                  ABSTOL + EPS *   max( |a|,|b| ) ,
167*>
168*>          where EPS is the machine precision.  If ABSTOL is less than
169*>          or equal to zero, then  EPS*|T|  will be used in its place,
170*>          where |T| is the 1-norm of the tridiagonal matrix obtained
171*>          by reducing C to tridiagonal form, where C is the symmetric
172*>          matrix of the standard symmetric problem to which the
173*>          generalized problem is transformed.
174*>
175*>          Eigenvalues will be computed most accurately when ABSTOL is
176*>          set to twice the underflow threshold 2*DLAMCH('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*DLAMCH('S').
180*> \endverbatim
181*>
182*> \param[out] M
183*> \verbatim
184*>          M is INTEGER
185*>          The total number of eigenvalues found.  0 <= M <= N.
186*>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
187*> \endverbatim
188*>
189*> \param[out] W
190*> \verbatim
191*>          W is DOUBLE PRECISION array, dimension (N)
192*>          On normal exit, the first M elements contain the selected
193*>          eigenvalues in ascending order.
194*> \endverbatim
195*>
196*> \param[out] Z
197*> \verbatim
198*>          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M))
199*>          If JOBZ = 'N', then Z is not referenced.
200*>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
201*>          contain the orthonormal eigenvectors of the matrix A
202*>          corresponding to the selected eigenvalues, with the i-th
203*>          column of Z holding the eigenvector associated with W(i).
204*>          The eigenvectors are normalized as follows:
205*>          if ITYPE = 1 or 2, Z**T*B*Z = I;
206*>          if ITYPE = 3, Z**T*inv(B)*Z = I.
207*>
208*>          If an eigenvector fails to converge, then that column of Z
209*>          contains the latest approximation to the eigenvector, and the
210*>          index of the eigenvector is returned in IFAIL.
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 DOUBLE PRECISION array, dimension (MAX(1,LWORK))
226*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
227*> \endverbatim
228*>
229*> \param[in] LWORK
230*> \verbatim
231*>          LWORK is INTEGER
232*>          The length of the array WORK.  LWORK >= max(1,8*N).
233*>          For optimal efficiency, LWORK >= (NB+3)*N,
234*>          where NB is the blocksize for DSYTRD returned by ILAENV.
235*>
236*>          If LWORK = -1, then a workspace query is assumed; the routine
237*>          only calculates the optimal size of the WORK array, returns
238*>          this value as the first entry of the WORK array, and no error
239*>          message related to LWORK is issued by XERBLA.
240*> \endverbatim
241*>
242*> \param[out] IWORK
243*> \verbatim
244*>          IWORK is INTEGER array, dimension (5*N)
245*> \endverbatim
246*>
247*> \param[out] IFAIL
248*> \verbatim
249*>          IFAIL is INTEGER array, dimension (N)
250*>          If JOBZ = 'V', then if INFO = 0, the first M elements of
251*>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
252*>          indices of the eigenvectors that failed to converge.
253*>          If JOBZ = 'N', then IFAIL is not referenced.
254*> \endverbatim
255*>
256*> \param[out] INFO
257*> \verbatim
258*>          INFO is INTEGER
259*>          = 0:  successful exit
260*>          < 0:  if INFO = -i, the i-th argument had an illegal value
261*>          > 0:  DPOTRF or DSYEVX returned an error code:
262*>             <= N:  if INFO = i, DSYEVX failed to converge;
263*>                    i eigenvectors failed to converge.  Their indices
264*>                    are stored in array IFAIL.
265*>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
266*>                    minor of order i of B is not positive definite.
267*>                    The factorization of B could not be completed and
268*>                    no eigenvalues or eigenvectors were computed.
269*> \endverbatim
270*
271*  Authors:
272*  ========
273*
274*> \author Univ. of Tennessee
275*> \author Univ. of California Berkeley
276*> \author Univ. of Colorado Denver
277*> \author NAG Ltd.
278*
279*> \date November 2015
280*
281*> \ingroup doubleSYeigen
282*
283*> \par Contributors:
284*  ==================
285*>
286*>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
287*
288*  =====================================================================
289      SUBROUTINE DSYGVX( ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB,
290     $                   VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK,
291     $                   LWORK, IWORK, IFAIL, INFO )
292*
293*  -- LAPACK driver routine (version 3.6.0) --
294*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
295*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
296*     November 2015
297*
298*     .. Scalar Arguments ..
299      CHARACTER          JOBZ, RANGE, UPLO
300      INTEGER            IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
301      DOUBLE PRECISION   ABSTOL, VL, VU
302*     ..
303*     .. Array Arguments ..
304      INTEGER            IFAIL( * ), IWORK( * )
305      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), W( * ), WORK( * ),
306     $                   Z( LDZ, * )
307*     ..
308*
309* =====================================================================
310*
311*     .. Parameters ..
312      DOUBLE PRECISION   ONE
313      PARAMETER          ( ONE = 1.0D+0 )
314*     ..
315*     .. Local Scalars ..
316      LOGICAL            ALLEIG, INDEIG, LQUERY, UPPER, VALEIG, WANTZ
317      CHARACTER          TRANS
318      INTEGER            LWKMIN, LWKOPT, NB
319*     ..
320*     .. External Functions ..
321      LOGICAL            LSAME
322      INTEGER            ILAENV
323      EXTERNAL           LSAME, ILAENV
324*     ..
325*     .. External Subroutines ..
326      EXTERNAL           DPOTRF, DSYEVX, DSYGST, DTRMM, DTRSM, XERBLA
327*     ..
328*     .. Intrinsic Functions ..
329      INTRINSIC          MAX, MIN
330*     ..
331*     .. Executable Statements ..
332*
333*     Test the input parameters.
334*
335      UPPER = LSAME( UPLO, 'U' )
336      WANTZ = LSAME( JOBZ, 'V' )
337      ALLEIG = LSAME( RANGE, 'A' )
338      VALEIG = LSAME( RANGE, 'V' )
339      INDEIG = LSAME( RANGE, 'I' )
340      LQUERY = ( LWORK.EQ.-1 )
341*
342      INFO = 0
343      IF( ITYPE.LT.1 .OR. ITYPE.GT.3 ) THEN
344         INFO = -1
345      ELSE IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
346         INFO = -2
347      ELSE IF( .NOT.( ALLEIG .OR. VALEIG .OR. INDEIG ) ) THEN
348         INFO = -3
349      ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
350         INFO = -4
351      ELSE IF( N.LT.0 ) THEN
352         INFO = -5
353      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
354         INFO = -7
355      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
356         INFO = -9
357      ELSE
358         IF( VALEIG ) THEN
359            IF( N.GT.0 .AND. VU.LE.VL )
360     $         INFO = -11
361         ELSE IF( INDEIG ) THEN
362            IF( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) THEN
363               INFO = -12
364            ELSE IF( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) THEN
365               INFO = -13
366            END IF
367         END IF
368      END IF
369      IF (INFO.EQ.0) THEN
370         IF (LDZ.LT.1 .OR. (WANTZ .AND. LDZ.LT.N)) THEN
371            INFO = -18
372         END IF
373      END IF
374*
375      IF( INFO.EQ.0 ) THEN
376         LWKMIN = MAX( 1, 8*N )
377         NB = ILAENV( 1, 'DSYTRD', UPLO, N, -1, -1, -1 )
378         LWKOPT = MAX( LWKMIN, ( NB + 3 )*N )
379         WORK( 1 ) = LWKOPT
380*
381         IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
382            INFO = -20
383         END IF
384      END IF
385*
386      IF( INFO.NE.0 ) THEN
387         CALL XERBLA( 'DSYGVX', -INFO )
388         RETURN
389      ELSE IF( LQUERY ) THEN
390         RETURN
391      END IF
392*
393*     Quick return if possible
394*
395      M = 0
396      IF( N.EQ.0 ) THEN
397         RETURN
398      END IF
399*
400*     Form a Cholesky factorization of B.
401*
402      CALL DPOTRF( UPLO, N, B, LDB, INFO )
403      IF( INFO.NE.0 ) THEN
404         INFO = N + INFO
405         RETURN
406      END IF
407*
408*     Transform problem to standard eigenvalue problem and solve.
409*
410      CALL DSYGST( ITYPE, UPLO, N, A, LDA, B, LDB, INFO )
411      CALL DSYEVX( JOBZ, RANGE, UPLO, N, A, LDA, VL, VU, IL, IU, ABSTOL,
412     $             M, W, Z, LDZ, WORK, LWORK, IWORK, IFAIL, INFO )
413*
414      IF( WANTZ ) THEN
415*
416*        Backtransform eigenvectors to the original problem.
417*
418         IF( INFO.GT.0 )
419     $      M = INFO - 1
420         IF( ITYPE.EQ.1 .OR. ITYPE.EQ.2 ) THEN
421*
422*           For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
423*           backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
424*
425            IF( UPPER ) THEN
426               TRANS = 'N'
427            ELSE
428               TRANS = 'T'
429            END IF
430*
431            CALL DTRSM( 'Left', UPLO, TRANS, 'Non-unit', N, M, ONE, B,
432     $                  LDB, Z, LDZ )
433*
434         ELSE IF( ITYPE.EQ.3 ) THEN
435*
436*           For B*A*x=(lambda)*x;
437*           backtransform eigenvectors: x = L*y or U**T*y
438*
439            IF( UPPER ) THEN
440               TRANS = 'T'
441            ELSE
442               TRANS = 'N'
443            END IF
444*
445            CALL DTRMM( 'Left', UPLO, TRANS, 'Non-unit', N, M, ONE, B,
446     $                  LDB, Z, LDZ )
447         END IF
448      END IF
449*
450*     Set WORK(1) to optimal workspace size.
451*
452      WORK( 1 ) = LWKOPT
453*
454      RETURN
455*
456*     End of DSYGVX
457*
458      END
459