LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zhbevd_2stage.f
Go to the documentation of this file.
1*> \brief <b> ZHBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
2*
3* @precisions fortran z -> s d c
4*
5* =========== DOCUMENTATION ===========
6*
7* Online html documentation available at
8* http://www.netlib.org/lapack/explore-html/
9*
10*> \htmlonly
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbevd_2stage.f">
13*> [TGZ]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbevd_2stage.f">
15*> [ZIP]</a>
16*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbevd_2stage.f">
17*> [TXT]</a>
18*> \endhtmlonly
19*
20* Definition:
21* ===========
22*
23* SUBROUTINE ZHBEVD_2STAGE( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
24* WORK, LWORK, RWORK, LRWORK, IWORK,
25* LIWORK, INFO )
26*
27* IMPLICIT NONE
28*
29* .. Scalar Arguments ..
30* CHARACTER JOBZ, UPLO
31* INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
32* ..
33* .. Array Arguments ..
34* INTEGER IWORK( * )
35* DOUBLE PRECISION RWORK( * ), W( * )
36* COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
37* ..
38*
39*
40*> \par Purpose:
41* =============
42*>
43*> \verbatim
44*>
45*> ZHBEVD_2STAGE computes all the eigenvalues and, optionally, eigenvectors of
46*> a complex Hermitian band matrix A using the 2stage technique for
47*> the reduction to tridiagonal. If eigenvectors are desired, it
48*> uses a divide and conquer algorithm.
49*>
50*> \endverbatim
51*
52* Arguments:
53* ==========
54*
55*> \param[in] JOBZ
56*> \verbatim
57*> JOBZ is CHARACTER*1
58*> = 'N': Compute eigenvalues only;
59*> = 'V': Compute eigenvalues and eigenvectors.
60*> Not available in this release.
61*> \endverbatim
62*>
63*> \param[in] UPLO
64*> \verbatim
65*> UPLO is CHARACTER*1
66*> = 'U': Upper triangle of A is stored;
67*> = 'L': Lower triangle of A is stored.
68*> \endverbatim
69*>
70*> \param[in] N
71*> \verbatim
72*> N is INTEGER
73*> The order of the matrix A. N >= 0.
74*> \endverbatim
75*>
76*> \param[in] KD
77*> \verbatim
78*> KD is INTEGER
79*> The number of superdiagonals of the matrix A if UPLO = 'U',
80*> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
81*> \endverbatim
82*>
83*> \param[in,out] AB
84*> \verbatim
85*> AB is COMPLEX*16 array, dimension (LDAB, N)
86*> On entry, the upper or lower triangle of the Hermitian band
87*> matrix A, stored in the first KD+1 rows of the array. The
88*> j-th column of A is stored in the j-th column of the array AB
89*> as follows:
90*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
91*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
92*>
93*> On exit, AB is overwritten by values generated during the
94*> reduction to tridiagonal form. If UPLO = 'U', the first
95*> superdiagonal and the diagonal of the tridiagonal matrix T
96*> are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
97*> the diagonal and first subdiagonal of T are returned in the
98*> first two rows of AB.
99*> \endverbatim
100*>
101*> \param[in] LDAB
102*> \verbatim
103*> LDAB is INTEGER
104*> The leading dimension of the array AB. LDAB >= KD + 1.
105*> \endverbatim
106*>
107*> \param[out] W
108*> \verbatim
109*> W is DOUBLE PRECISION array, dimension (N)
110*> If INFO = 0, the eigenvalues in ascending order.
111*> \endverbatim
112*>
113*> \param[out] Z
114*> \verbatim
115*> Z is COMPLEX*16 array, dimension (LDZ, N)
116*> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
117*> eigenvectors of the matrix A, with the i-th column of Z
118*> holding the eigenvector associated with W(i).
119*> If JOBZ = 'N', then Z is not referenced.
120*> \endverbatim
121*>
122*> \param[in] LDZ
123*> \verbatim
124*> LDZ is INTEGER
125*> The leading dimension of the array Z. LDZ >= 1, and if
126*> JOBZ = 'V', LDZ >= max(1,N).
127*> \endverbatim
128*>
129*> \param[out] WORK
130*> \verbatim
131*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
132*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
133*> \endverbatim
134*>
135*> \param[in] LWORK
136*> \verbatim
137*> LWORK is INTEGER
138*> The length of the array WORK. LWORK >= 1, when N <= 1;
139*> otherwise
140*> If JOBZ = 'N' and N > 1, LWORK must be queried.
141*> LWORK = MAX(1, dimension) where
142*> dimension = (2KD+1)*N + KD*NTHREADS
143*> where KD is the size of the band.
145*> openMP compilation is enabled, otherwise =1.
146*> If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available.
147*>
148*> If LWORK = -1, then a workspace query is assumed; the routine
149*> only calculates the optimal sizes of the WORK, RWORK and
150*> IWORK arrays, returns these values as the first entries of
151*> the WORK, RWORK and IWORK arrays, and no error message
152*> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
153*> \endverbatim
154*>
155*> \param[out] RWORK
156*> \verbatim
157*> RWORK is DOUBLE PRECISION array,
158*> dimension (LRWORK)
159*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
160*> \endverbatim
161*>
162*> \param[in] LRWORK
163*> \verbatim
164*> LRWORK is INTEGER
165*> The dimension of array RWORK.
166*> If N <= 1, LRWORK must be at least 1.
167*> If JOBZ = 'N' and N > 1, LRWORK must be at least N.
168*> If JOBZ = 'V' and N > 1, LRWORK must be at least
169*> 1 + 5*N + 2*N**2.
170*>
171*> If LRWORK = -1, then a workspace query is assumed; the
172*> routine only calculates the optimal sizes of the WORK, RWORK
173*> and IWORK arrays, returns these values as the first entries
174*> of the WORK, RWORK and IWORK arrays, and no error message
175*> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
176*> \endverbatim
177*>
178*> \param[out] IWORK
179*> \verbatim
180*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
181*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
182*> \endverbatim
183*>
184*> \param[in] LIWORK
185*> \verbatim
186*> LIWORK is INTEGER
187*> The dimension of array IWORK.
188*> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
189*> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N .
190*>
191*> If LIWORK = -1, then a workspace query is assumed; the
192*> routine only calculates the optimal sizes of the WORK, RWORK
193*> and IWORK arrays, returns these values as the first entries
194*> of the WORK, RWORK and IWORK arrays, and no error message
195*> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
196*> \endverbatim
197*>
198*> \param[out] INFO
199*> \verbatim
200*> INFO is INTEGER
201*> = 0: successful exit.
202*> < 0: if INFO = -i, the i-th argument had an illegal value.
203*> > 0: if INFO = i, the algorithm failed to converge; i
204*> off-diagonal elements of an intermediate tridiagonal
205*> form did not converge to zero.
206*> \endverbatim
207*
208* Authors:
209* ========
210*
211*> \author Univ. of Tennessee
212*> \author Univ. of California Berkeley
213*> \author Univ. of Colorado Denver
214*> \author NAG Ltd.
215*
216*> \ingroup hbevd_2stage
217*
218*> \par Further Details:
219* =====================
220*>
221*> \verbatim
222*>
223*> All details about the 2stage techniques are available in:
224*>
225*> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
226*> Parallel reduction to condensed forms for symmetric eigenvalue problems
227*> using aggregated fine-grained and memory-aware kernels. In Proceedings
228*> of 2011 International Conference for High Performance Computing,
229*> Networking, Storage and Analysis (SC '11), New York, NY, USA,
230*> Article 8 , 11 pages.
231*> http://doi.acm.org/10.1145/2063384.2063394
232*>
233*> A. Haidar, J. Kurzak, P. Luszczek, 2013.
234*> An improved parallel singular value algorithm and its implementation
235*> for multicore hardware, In Proceedings of 2013 International Conference
236*> for High Performance Computing, Networking, Storage and Analysis (SC '13).
238*> Article 90, 12 pages.
239*> http://doi.acm.org/10.1145/2503210.2503292
240*>
241*> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
242*> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
243*> calculations based on fine-grained memory aware tasks.
244*> International Journal of High Performance Computing Applications.
245*> Volume 28 Issue 2, Pages 196-209, May 2014.
246*> http://hpc.sagepub.com/content/28/2/196
247*>
248*> \endverbatim
249*
250* =====================================================================
251 SUBROUTINE zhbevd_2stage( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
252 \$ WORK, LWORK, RWORK, LRWORK, IWORK,
253 \$ LIWORK, INFO )
254*
255 IMPLICIT NONE
256*
257* -- LAPACK driver routine --
258* -- LAPACK is a software package provided by Univ. of Tennessee, --
259* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
260*
261* .. Scalar Arguments ..
262 CHARACTER JOBZ, UPLO
263 INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
264* ..
265* .. Array Arguments ..
266 INTEGER IWORK( * )
267 DOUBLE PRECISION RWORK( * ), W( * )
268 COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
269* ..
270*
271* =====================================================================
272*
273* .. Parameters ..
274 DOUBLE PRECISION ZERO, ONE
275 PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
276 COMPLEX*16 CZERO, CONE
277 parameter( czero = ( 0.0d0, 0.0d0 ),
278 \$ cone = ( 1.0d0, 0.0d0 ) )
279* ..
280* .. Local Scalars ..
281 LOGICAL LOWER, LQUERY, WANTZ
282 INTEGER IINFO, IMAX, INDE, INDWK2, INDRWK, ISCALE,
283 \$ llwork, indwk, lhtrd, lwtrd, ib, indhous,
284 \$ liwmin, llrwk, llwk2, lrwmin, lwmin
285 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
286 \$ SMLNUM
287* ..
288* .. External Functions ..
289 LOGICAL LSAME
290 INTEGER ILAENV2STAGE
291 DOUBLE PRECISION DLAMCH, ZLANHB
292 EXTERNAL lsame, dlamch, zlanhb, ilaenv2stage
293* ..
294* .. External Subroutines ..
295 EXTERNAL dscal, dsterf, xerbla, zgemm, zlacpy,
297* ..
298* .. Intrinsic Functions ..
299 INTRINSIC dble, sqrt
300* ..
301* .. Executable Statements ..
302*
303* Test the input parameters.
304*
305 wantz = lsame( jobz, 'V' )
306 lower = lsame( uplo, 'L' )
307 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 .OR. lrwork.EQ.-1 )
308*
309 info = 0
310 IF( n.LE.1 ) THEN
311 lwmin = 1
312 lrwmin = 1
313 liwmin = 1
314 ELSE
315 ib = ilaenv2stage( 2, 'ZHETRD_HB2ST', jobz, n, kd, -1, -1 )
316 lhtrd = ilaenv2stage( 3, 'ZHETRD_HB2ST', jobz, n, kd, ib, -1 )
317 lwtrd = ilaenv2stage( 4, 'ZHETRD_HB2ST', jobz, n, kd, ib, -1 )
318 IF( wantz ) THEN
319 lwmin = 2*n**2
320 lrwmin = 1 + 5*n + 2*n**2
321 liwmin = 3 + 5*n
322 ELSE
323 lwmin = max( n, lhtrd + lwtrd )
324 lrwmin = n
325 liwmin = 1
326 END IF
327 END IF
328 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
329 info = -1
330 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
331 info = -2
332 ELSE IF( n.LT.0 ) THEN
333 info = -3
334 ELSE IF( kd.LT.0 ) THEN
335 info = -4
336 ELSE IF( ldab.LT.kd+1 ) THEN
337 info = -6
338 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
339 info = -9
340 END IF
341*
342 IF( info.EQ.0 ) THEN
343 work( 1 ) = lwmin
344 rwork( 1 ) = lrwmin
345 iwork( 1 ) = liwmin
346*
347 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
348 info = -11
349 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
350 info = -13
351 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
352 info = -15
353 END IF
354 END IF
355*
356 IF( info.NE.0 ) THEN
357 CALL xerbla( 'ZHBEVD_2STAGE', -info )
358 RETURN
359 ELSE IF( lquery ) THEN
360 RETURN
361 END IF
362*
363* Quick return if possible
364*
365 IF( n.EQ.0 )
366 \$ RETURN
367*
368 IF( n.EQ.1 ) THEN
369 w( 1 ) = dble( ab( 1, 1 ) )
370 IF( wantz )
371 \$ z( 1, 1 ) = cone
372 RETURN
373 END IF
374*
375* Get machine constants.
376*
377 safmin = dlamch( 'Safe minimum' )
378 eps = dlamch( 'Precision' )
379 smlnum = safmin / eps
380 bignum = one / smlnum
381 rmin = sqrt( smlnum )
382 rmax = sqrt( bignum )
383*
384* Scale matrix to allowable range, if necessary.
385*
386 anrm = zlanhb( 'M', uplo, n, kd, ab, ldab, rwork )
387 iscale = 0
388 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
389 iscale = 1
390 sigma = rmin / anrm
391 ELSE IF( anrm.GT.rmax ) THEN
392 iscale = 1
393 sigma = rmax / anrm
394 END IF
395 IF( iscale.EQ.1 ) THEN
396 IF( lower ) THEN
397 CALL zlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
398 ELSE
399 CALL zlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
400 END IF
401 END IF
402*
403* Call ZHBTRD_HB2ST to reduce Hermitian band matrix to tridiagonal form.
404*
405 inde = 1
406 indrwk = inde + n
407 llrwk = lrwork - indrwk + 1
408 indhous = 1
409 indwk = indhous + lhtrd
410 llwork = lwork - indwk + 1
411 indwk2 = indwk + n*n
412 llwk2 = lwork - indwk2 + 1
413*
414 CALL zhetrd_hb2st( "N", jobz, uplo, n, kd, ab, ldab, w,
415 \$ rwork( inde ), work( indhous ), lhtrd,
416 \$ work( indwk ), llwork, iinfo )
417*
418* For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEDC.
419*
420 IF( .NOT.wantz ) THEN
421 CALL dsterf( n, w, rwork( inde ), info )
422 ELSE
423 CALL zstedc( 'I', n, w, rwork( inde ), work, n, work( indwk2 ),
424 \$ llwk2, rwork( indrwk ), llrwk, iwork, liwork,
425 \$ info )
426 CALL zgemm( 'N', 'N', n, n, n, cone, z, ldz, work, n, czero,
427 \$ work( indwk2 ), n )
428 CALL zlacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
429 END IF
430*
431* If matrix was scaled, then rescale eigenvalues appropriately.
432*
433 IF( iscale.EQ.1 ) THEN
434 IF( info.EQ.0 ) THEN
435 imax = n
436 ELSE
437 imax = info - 1
438 END IF
439 CALL dscal( imax, one / sigma, w, 1 )
440 END IF
441*
442 work( 1 ) = lwmin
443 rwork( 1 ) = lrwmin
444 iwork( 1 ) = liwmin
445 RETURN
446*
447* End of ZHBEVD_2STAGE
448*
449 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
subroutine zhbevd_2stage(jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
ZHBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER ...
subroutine zhetrd_hb2st(stage1, vect, uplo, n, kd, ab, ldab, d, e, hous, lhous, work, lwork, info)
ZHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric tridiagonal form T
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
subroutine zlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition zlascl.f:143
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine zstedc(compz, n, d, e, z, ldz, work, lwork, rwork, lrwork, iwork, liwork, info)
ZSTEDC
Definition zstedc.f:206
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86