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