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 >= 2*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*> \ingroup doubleOTHEReigen
218*
219*> \par Contributors:
220*  ==================
221*>
222*>     Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA
223*
224*  =====================================================================
225      SUBROUTINE DSBGVD( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, W,
226     $                   Z, LDZ, WORK, LWORK, IWORK, LIWORK, INFO )
227*
228*  -- LAPACK driver routine --
229*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
230*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231*
232*     .. Scalar Arguments ..
233      CHARACTER          JOBZ, UPLO
234      INTEGER            INFO, KA, KB, LDAB, LDBB, LDZ, LIWORK, LWORK, N
235*     ..
236*     .. Array Arguments ..
237      INTEGER            IWORK( * )
238      DOUBLE PRECISION   AB( LDAB, * ), BB( LDBB, * ), W( * ),
239     $                   WORK( * ), Z( LDZ, * )
240*     ..
241*
242*  =====================================================================
243*
244*     .. Parameters ..
245      DOUBLE PRECISION   ONE, ZERO
246      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
247*     ..
248*     .. Local Scalars ..
249      LOGICAL            LQUERY, UPPER, WANTZ
250      CHARACTER          VECT
251      INTEGER            IINFO, INDE, INDWK2, INDWRK, LIWMIN, LLWRK2,
252     $                   LWMIN
253*     ..
254*     .. External Functions ..
255      LOGICAL            LSAME
256      EXTERNAL           LSAME
257*     ..
258*     .. External Subroutines ..
259      EXTERNAL           DGEMM, DLACPY, DPBSTF, DSBGST, DSBTRD, DSTEDC,
260     $                   DSTERF, XERBLA
261*     ..
262*     .. Executable Statements ..
263*
264*     Test the input parameters.
265*
266      WANTZ = LSAME( JOBZ, 'V' )
267      UPPER = LSAME( UPLO, 'U' )
268      LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
269*
270      INFO = 0
271      IF( N.LE.1 ) THEN
272         LIWMIN = 1
273         LWMIN = 1
274      ELSE IF( WANTZ ) THEN
275         LIWMIN = 3 + 5*N
276         LWMIN = 1 + 5*N + 2*N**2
277      ELSE
278         LIWMIN = 1
279         LWMIN = 2*N
280      END IF
281*
282      IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
283         INFO = -1
284      ELSE IF( .NOT.( UPPER .OR. LSAME( UPLO, 'L' ) ) ) THEN
285         INFO = -2
286      ELSE IF( N.LT.0 ) THEN
287         INFO = -3
288      ELSE IF( KA.LT.0 ) THEN
289         INFO = -4
290      ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
291         INFO = -5
292      ELSE IF( LDAB.LT.KA+1 ) THEN
293         INFO = -7
294      ELSE IF( LDBB.LT.KB+1 ) THEN
295         INFO = -9
296      ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
297         INFO = -12
298      END IF
299*
300      IF( INFO.EQ.0 ) THEN
301         WORK( 1 ) = LWMIN
302         IWORK( 1 ) = LIWMIN
303*
304         IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
305            INFO = -14
306         ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
307            INFO = -16
308         END IF
309      END IF
310*
311      IF( INFO.NE.0 ) THEN
312         CALL XERBLA( 'DSBGVD', -INFO )
313         RETURN
314      ELSE IF( LQUERY ) THEN
315         RETURN
316      END IF
317*
318*     Quick return if possible
319*
320      IF( N.EQ.0 )
321     $   RETURN
322*
323*     Form a split Cholesky factorization of B.
324*
325      CALL DPBSTF( UPLO, N, KB, BB, LDBB, INFO )
326      IF( INFO.NE.0 ) THEN
327         INFO = N + INFO
328         RETURN
329      END IF
330*
331*     Transform problem to standard eigenvalue problem.
332*
333      INDE = 1
334      INDWRK = INDE + N
335      INDWK2 = INDWRK + N*N
336      LLWRK2 = LWORK - INDWK2 + 1
337      CALL DSBGST( JOBZ, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, Z, LDZ,
338     $             WORK, IINFO )
339*
340*     Reduce to tridiagonal form.
341*
342      IF( WANTZ ) THEN
343         VECT = 'U'
344      ELSE
345         VECT = 'N'
346      END IF
347      CALL DSBTRD( VECT, UPLO, N, KA, AB, LDAB, W, WORK( INDE ), Z, LDZ,
348     $             WORK( INDWRK ), IINFO )
349*
350*     For eigenvalues only, call DSTERF. For eigenvectors, call SSTEDC.
351*
352      IF( .NOT.WANTZ ) THEN
353         CALL DSTERF( N, W, WORK( INDE ), INFO )
354      ELSE
355         CALL DSTEDC( 'I', N, W, WORK( INDE ), WORK( INDWRK ), N,
356     $                WORK( INDWK2 ), LLWRK2, IWORK, LIWORK, INFO )
357         CALL DGEMM( 'N', 'N', N, N, N, ONE, Z, LDZ, WORK( INDWRK ), N,
358     $               ZERO, WORK( INDWK2 ), N )
359         CALL DLACPY( 'A', N, N, WORK( INDWK2 ), N, Z, LDZ )
360      END IF
361*
362      WORK( 1 ) = LWMIN
363      IWORK( 1 ) = LIWMIN
364*
365      RETURN
366*
367*     End of DSBGVD
368*
369      END
370