1*> \brief \b DSBGVD
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DSBGVD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbgvd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbgvd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbgvd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DSBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
22*                          Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          JOBZ, UPLO
26*       INTEGER            INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LWORK, N
27*       ..
28*       .. Array Arguments ..
29*       INTEGER            IWORK( * )
30*       DOUBLE PRECISION   AB( LDAB, * ), BB( LDBB, * ), W( * ),
31*      $                   WORK( * ), Z( LDZ, * )
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*> \verbatim
39*>
40*> DSBGVD computes all the eigenvalues, and optionally, the eigenvectors
41*> of a real generalized symmetric-definite banded eigenproblem, of the
42*> form A*x=(lambda)*B*x.  Here A and B are assumed to be symmetric and
43*> banded, and B is also positive definite.  If eigenvectors are
44*> desired, it uses a divide and conquer algorithm.
45*>
46*> The divide and conquer algorithm makes very mild assumptions about
47*> floating point arithmetic. It will work on machines with a guard
48*> digit in add/subtract, or on those binary machines without guard
49*> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
50*> Cray-2. It could conceivably fail on hexadecimal or decimal machines
51*> without guard digits, but we know of none.
52*> \endverbatim
53*
54*  Arguments:
55*  ==========
56*
57*> \param[in] JOBZ
58*> \verbatim
59*>          JOBZ is CHARACTER*1
60*>          = 'N':  Compute eigenvalues only;
61*>          = 'V':  Compute eigenvalues and eigenvectors.
62*> \endverbatim
63*>
64*> \param[in] UPLO
65*> \verbatim
66*>          UPLO is CHARACTER*1
67*>          = 'U':  Upper triangles of A and B are stored;
68*>          = 'L':  Lower triangles of A and B are stored.
69*> \endverbatim
70*>
71*> \param[in] N
72*> \verbatim
73*>          N is INTEGER
74*>          The order of the matrices A and B.  N >= 0.
75*> \endverbatim
76*>
77*> \param[in] KA
78*> \verbatim
79*>          KA is INTEGER
80*>          The number of superdiagonals of the matrix A if UPLO = 'U',
81*>          or the number of subdiagonals if UPLO = 'L'.  KA >= 0.
82*> \endverbatim
83*>
84*> \param[in] KB
85*> \verbatim
86*>          KB is INTEGER
87*>          The number of superdiagonals of the matrix B if UPLO = 'U',
88*>          or the number of subdiagonals if UPLO = 'L'.  KB >= 0.
89*> \endverbatim
90*>
91*> \param[in,out] AB
92*> \verbatim
93*>          AB is DOUBLE PRECISION array, dimension (LDAB, N)
94*>          On entry, the upper or lower triangle of the symmetric band
95*>          matrix A, stored in the first ka+1 rows of the array.  The
96*>          j-th column of A is stored in the j-th column of the array AB
97*>          as follows:
98*>          if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
99*>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+ka).
100*>
101*>          On exit, the contents of AB are destroyed.
102*> \endverbatim
103*>
104*> \param[in] LDAB
105*> \verbatim
106*>          LDAB is INTEGER
107*>          The leading dimension of the array AB.  LDAB >= KA+1.
108*> \endverbatim
109*>
110*> \param[in,out] BB
111*> \verbatim
112*>          BB is DOUBLE PRECISION array, dimension (LDBB, N)
113*>          On entry, the upper or lower triangle of the symmetric band
114*>          matrix B, stored in the first kb+1 rows of the array.  The
115*>          j-th column of B is stored in the j-th column of the array BB
116*>          as follows:
117*>          if UPLO = 'U', BB(ka+1+i-j,j) = B(i,j) for max(1,j-kb)<=i<=j;
118*>          if UPLO = 'L', BB(1+i-j,j)    = B(i,j) for j<=i<=min(n,j+kb).
119*>
120*>          On exit, the factor S from the split Cholesky factorization
121*>          B = S**T*S, as returned by DPBSTF.
122*> \endverbatim
123*>
124*> \param[in] LDBB
125*> \verbatim
126*>          LDBB is INTEGER
127*>          The leading dimension of the array BB.  LDBB >= KB+1.
128*> \endverbatim
129*>
130*> \param[out] W
131*> \verbatim
132*>          W is DOUBLE PRECISION array, dimension (N)
133*>          If INFO = 0, the eigenvalues in ascending order.
134*> \endverbatim
135*>
136*> \param[out] Z
137*> \verbatim
138*>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
139*>          If JOBZ = 'V', then if INFO = 0, Z contains the matrix Z of
140*>          eigenvectors, with the i-th column of Z holding the
141*>          eigenvector associated with W(i).  The eigenvectors are
142*>          normalized so Z**T*B*Z = I.
143*>          If JOBZ = 'N', then Z is not referenced.
144*> \endverbatim
145*>
146*> \param[in] LDZ
147*> \verbatim
148*>          LDZ is INTEGER
149*>          The leading dimension of the array Z.  LDZ >= 1, and if
150*>          JOBZ = 'V', LDZ >= max(1,N).
151*> \endverbatim
152*>
153*> \param[out] WORK
154*> \verbatim
155*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
156*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
157*> \endverbatim
158*>
159*> \param[in] LWORK
160*> \verbatim
161*>          LWORK is INTEGER
162*>          The dimension of the array WORK.
163*>          If N <= 1,               LWORK >= 1.
164*>          If JOBZ = 'N' and N > 1, LWORK >= 3*N.
165*>          If JOBZ = 'V' and N > 1, LWORK >= 1 + 5*N + 2*N**2.
166*>
167*>          If LWORK = -1, then a workspace query is assumed; the routine
168*>          only calculates the optimal sizes of the WORK and IWORK
169*>          arrays, returns these values as the first entries of the WORK
170*>          and IWORK arrays, and no error message related to LWORK or
171*>          LIWORK is issued by XERBLA.
172*> \endverbatim
173*>
174*> \param[out] IWORK
175*> \verbatim
176*>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
177*>          On exit, if LIWORK > 0, IWORK(1) returns the optimal LIWORK.
178*> \endverbatim
179*>
180*> \param[in] LIWORK
181*> \verbatim
182*>          LIWORK is INTEGER
183*>          The dimension of the array IWORK.
184*>          If JOBZ  = 'N' or N <= 1, LIWORK >= 1.
185*>          If JOBZ  = 'V' and N > 1, LIWORK >= 3 + 5*N.
186*>
187*>          If LIWORK = -1, then a workspace query is assumed; the
188*>          routine only calculates the optimal sizes of the WORK and
189*>          IWORK arrays, returns these values as the first entries of
190*>          the WORK and IWORK arrays, and no error message related to
191*>          LWORK or LIWORK is issued by XERBLA.
192*> \endverbatim
193*>
194*> \param[out] INFO
195*> \verbatim
196*>          INFO is INTEGER
197*>          = 0:  successful exit
198*>          < 0:  if INFO = -i, the i-th argument had an illegal value
199*>          > 0:  if INFO = i, and i is:
200*>             <= N:  the algorithm failed to converge:
201*>                    i off-diagonal elements of an intermediate
202*>                    tridiagonal form did not converge to zero;
203*>             > N:   if INFO = N + i, for 1 <= i <= N, then DPBSTF
204*>                    returned INFO = i: B is not positive definite.
205*>                    The factorization of B could not be completed and
206*>                    no eigenvalues or eigenvectors were computed.
207*> \endverbatim
208*
209*  Authors:
210*  ========
211*
212*> \author Univ. of Tennessee
213*> \author Univ. of California Berkeley
214*> \author Univ. of Colorado Denver
215*> \author NAG Ltd.
216*
217*> \date November 2015
218*
219*> \ingroup doubleOTHEReigen
220*
221*> \par Contributors:
222*  ==================
223*>
224*>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
225*
226*  =====================================================================
227      SUBROUTINE DSBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
228     $                   Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO )
229*
230*  -- LAPACK driver routine (version 3.6.0) --
231*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
232*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
233*     November 2015
234*
235*     .. Scalar Arguments ..
236      CHARACTER          JOBZ, UPLO
237      INTEGER            INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LWORK, N
238*     ..
239*     .. Array Arguments ..
240      INTEGER            IWORK( * )
241      DOUBLE PRECISION   AB( LDAB, * ), BB( LDBB, * ), W( * ),
242     $                   WORK( * ), Z( LDZ, * )
243*     ..
244*
245*  =====================================================================
246*
247*     .. Parameters ..
248      DOUBLE PRECISION   ONE, ZERO
249      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
250*     ..
251*     .. Local Scalars ..
252      LOGICAL            LQUERY, UPPER, WANTZ
253      CHARACTER          VECT
254      INTEGER            IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLWRK2,
255     $                   LWMIN
256*     ..
257*     .. External Functions ..
258      LOGICAL            LSAME
259      EXTERNAL           LSAME
260*     ..
261*     .. External Subroutines ..
262      EXTERNAL           DGEMM, DLACPY, DPBSTF, DSBGST, DSBTRD, DSTEDC,
263     $                   DSTERF, XERBLA
264*     ..
265*     .. Executable Statements ..
266*
267*     Test the input parameters.
268*
269      WANTZ = LSAME( JOBZ, 'V' )
270      UPPER = LSAME( UPLO, 'U' )
271      LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
272*
273      INFO = 0
274      IF( N.LE.1 ) THEN
275         LIWMIN = 1
276         LWMIN = 1
277      ELSE IF( WANTZ ) THEN
278         LIWMIN = 3 + 5*N
279         LWMIN = 1 + 5*N + 2*N**2
280      ELSE
281         LIWMIN = 1
282         LWMIN = 2*N
283      END IF
284*
285      IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
286         INFO = -1
287      ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
288         INFO = -2
289      ELSE IF( N.LT.0 ) THEN
290         INFO = -3
291      ELSE IF( KA.LT.0 ) THEN
292         INFO = -4
293      ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
294         INFO = -5
295      ELSE IF( LDAB.LT.KA+1 ) THEN
296         INFO = -7
297      ELSE IF( LDBB.LT.KB+1 ) THEN
298         INFO = -9
299      ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
300         INFO = -12
301      END IF
302*
303      IF( INFO.EQ.0 ) THEN
304         WORK( 1 ) = LWMIN
305         IWORK( 1 ) = LIWMIN
306*
307         IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
308            INFO = -14
309         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
310            INFO = -16
311         END IF
312      END IF
313*
314      IF( INFO.NE.0 ) THEN
315         CALL XERBLA( 'DSBGVD', -INFO )
316         RETURN
317      ELSE IF( LQUERY ) THEN
318         RETURN
319      END IF
320*
321*     Quick return if possible
322*
323      IF( N.EQ.0 )
324     $   RETURN
325*
326*     Form a split Cholesky factorization of B.
327*
328      CALL DPBSTF( UPLO, N, KB, BB, LDBB, INFO )
329      IF( INFO.NE.0 ) THEN
330         INFO = N + INFO
331         RETURN
332      END IF
333*
334*     Transform problem to standard eigenvalue problem.
335*
336      INDE = 1
337      INDWRK = INDE + N
338      INDWK2 = INDWRK + N*N
339      LLWRK2 = LWORK - INDWK2 + 1
340      CALL DSBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ,
341     $             WORK( INDWRK ), IINFO )
342*
343*     Reduce to tridiagonal form.
344*
345      IF( WANTZ ) THEN
346         VECT = 'U'
347      ELSE
348         VECT = 'N'
349      END IF
350      CALL DSBTRD( VECT, UPLO, N, KA, AB, LDAB, W, WORK( INDE ), Z, LDZ,
351     $             WORK( INDWRK ), IINFO )
352*
353*     For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
354*
355      IF( .NOT.WANTZ ) THEN
356         CALL DSTERF( N, W, WORK( INDE ), INFO )
357      ELSE
358         CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
359     $                WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
360         CALL DGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N,
361     $               ZERO, WORK( INDWK2 ), N )
362         CALL DLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
363      END IF
364*
365      WORK( 1 ) = LWMIN
366      IWORK( 1 ) = LIWMIN
367*
368      RETURN
369*
370*     End of DSBGVD
371*
372      END
373