LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
chbevd_2stage.f
Go to the documentation of this file.
1 *> \brief <b> CHBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices</b>
2 *
3 * @generated from zhbevd_2stage.f, fortran z -> c, Sat Nov 5 23:18:17 2016
4 *
5 * =========== DOCUMENTATION ===========
6 *
7 * Online html documentation available at
8 * http://www.netlib.org/lapack/explore-html/
9 *
10 *> \htmlonly
11 *> Download CHBEVD_2STAGE + dependencies
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chbevd_2stage.f">
13 *> [TGZ]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chbevd_2stage.f">
15 *> [ZIP]</a>
16 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chbevd_2stage.f">
17 *> [TXT]</a>
18 *> \endhtmlonly
19 *
20 * Definition:
21 * ===========
22 *
23 * SUBROUTINE CHBEVD_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 * REAL RWORK( * ), W( * )
36 * COMPLEX AB( LDAB, * ), WORK( * ), Z( LDZ, * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> CHBEVD_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 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 REAL 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 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 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 REAL 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 *> \date November 2017
223 *
224 *> \ingroup complexOTHEReigen
225 *
226 *> \par Further Details:
227 * =====================
228 *>
229 *> \verbatim
230 *>
231 *> All details about the 2stage techniques are available in:
232 *>
233 *> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
234 *> Parallel reduction to condensed forms for symmetric eigenvalue problems
235 *> using aggregated fine-grained and memory-aware kernels. In Proceedings
236 *> of 2011 International Conference for High Performance Computing,
237 *> Networking, Storage and Analysis (SC '11), New York, NY, USA,
238 *> Article 8 , 11 pages.
239 *> http://doi.acm.org/10.1145/2063384.2063394
240 *>
241 *> A. Haidar, J. Kurzak, P. Luszczek, 2013.
242 *> An improved parallel singular value algorithm and its implementation
243 *> for multicore hardware, In Proceedings of 2013 International Conference
244 *> for High Performance Computing, Networking, Storage and Analysis (SC '13).
245 *> Denver, Colorado, USA, 2013.
246 *> Article 90, 12 pages.
247 *> http://doi.acm.org/10.1145/2503210.2503292
248 *>
249 *> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
250 *> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
251 *> calculations based on fine-grained memory aware tasks.
252 *> International Journal of High Performance Computing Applications.
253 *> Volume 28 Issue 2, Pages 196-209, May 2014.
254 *> http://hpc.sagepub.com/content/28/2/196
255 *>
256 *> \endverbatim
257 *
258 * =====================================================================
259  SUBROUTINE chbevd_2stage( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
260  $ WORK, LWORK, RWORK, LRWORK, IWORK,
261  $ LIWORK, INFO )
262 *
263  IMPLICIT NONE
264 *
265 * -- LAPACK driver routine (version 3.8.0) --
266 * -- LAPACK is a software package provided by Univ. of Tennessee, --
267 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
268 * November 2017
269 *
270 * .. Scalar Arguments ..
271  CHARACTER JOBZ, UPLO
272  INTEGER INFO, KD, LDAB, LDZ, LIWORK, LRWORK, LWORK, N
273 * ..
274 * .. Array Arguments ..
275  INTEGER IWORK( * )
276  REAL RWORK( * ), W( * )
277  COMPLEX AB( ldab, * ), WORK( * ), Z( ldz, * )
278 * ..
279 *
280 * =====================================================================
281 *
282 * .. Parameters ..
283  REAL ZERO, ONE
284  parameter( zero = 0.0e0, one = 1.0e0 )
285  COMPLEX CZERO, CONE
286  parameter( czero = ( 0.0e0, 0.0e0 ),
287  $ cone = ( 1.0e0, 0.0e0 ) )
288 * ..
289 * .. Local Scalars ..
290  LOGICAL LOWER, LQUERY, WANTZ
291  INTEGER IINFO, IMAX, INDE, INDWK2, INDRWK, ISCALE,
292  $ llwork, indwk, lhtrd, lwtrd, ib, indhous,
293  $ liwmin, llrwk, llwk2, lrwmin, lwmin
294  REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
295  $ smlnum
296 * ..
297 * .. External Functions ..
298  LOGICAL LSAME
299  INTEGER ILAENV2STAGE
300  REAL SLAMCH, CLANHB
301  EXTERNAL lsame, slamch, clanhb, ilaenv2stage
302 * ..
303 * .. External Subroutines ..
304  EXTERNAL sscal, ssterf, xerbla, cgemm, clacpy,
305  $ clascl, cstedc, chetrd_hb2st
306 * ..
307 * .. Intrinsic Functions ..
308  INTRINSIC REAL, SQRT
309 * ..
310 * .. Executable Statements ..
311 *
312 * Test the input parameters.
313 *
314  wantz = lsame( jobz, 'V' )
315  lower = lsame( uplo, 'L' )
316  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 .OR. lrwork.EQ.-1 )
317 *
318  info = 0
319  IF( n.LE.1 ) THEN
320  lwmin = 1
321  lrwmin = 1
322  liwmin = 1
323  ELSE
324  ib = ilaenv2stage( 2, 'CHETRD_HB2ST', jobz, n, kd, -1, -1 )
325  lhtrd = ilaenv2stage( 3, 'CHETRD_HB2ST', jobz, n, kd, ib, -1 )
326  lwtrd = ilaenv2stage( 4, 'CHETRD_HB2ST', jobz, n, kd, ib, -1 )
327  IF( wantz ) THEN
328  lwmin = 2*n**2
329  lrwmin = 1 + 5*n + 2*n**2
330  liwmin = 3 + 5*n
331  ELSE
332  lwmin = max( n, lhtrd + lwtrd )
333  lrwmin = n
334  liwmin = 1
335  END IF
336  END IF
337  IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
338  info = -1
339  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
340  info = -2
341  ELSE IF( n.LT.0 ) THEN
342  info = -3
343  ELSE IF( kd.LT.0 ) THEN
344  info = -4
345  ELSE IF( ldab.LT.kd+1 ) THEN
346  info = -6
347  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
348  info = -9
349  END IF
350 *
351  IF( info.EQ.0 ) THEN
352  work( 1 ) = lwmin
353  rwork( 1 ) = lrwmin
354  iwork( 1 ) = liwmin
355 *
356  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
357  info = -11
358  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
359  info = -13
360  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
361  info = -15
362  END IF
363  END IF
364 *
365  IF( info.NE.0 ) THEN
366  CALL xerbla( 'CHBEVD_2STAGE', -info )
367  RETURN
368  ELSE IF( lquery ) THEN
369  RETURN
370  END IF
371 *
372 * Quick return if possible
373 *
374  IF( n.EQ.0 )
375  $ RETURN
376 *
377  IF( n.EQ.1 ) THEN
378  w( 1 ) = REAL( AB( 1, 1 ) )
379  IF( wantz )
380  $ z( 1, 1 ) = cone
381  RETURN
382  END IF
383 *
384 * Get machine constants.
385 *
386  safmin = slamch( 'Safe minimum' )
387  eps = slamch( 'Precision' )
388  smlnum = safmin / eps
389  bignum = one / smlnum
390  rmin = sqrt( smlnum )
391  rmax = sqrt( bignum )
392 *
393 * Scale matrix to allowable range, if necessary.
394 *
395  anrm = clanhb( 'M', uplo, n, kd, ab, ldab, rwork )
396  iscale = 0
397  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
398  iscale = 1
399  sigma = rmin / anrm
400  ELSE IF( anrm.GT.rmax ) THEN
401  iscale = 1
402  sigma = rmax / anrm
403  END IF
404  IF( iscale.EQ.1 ) THEN
405  IF( lower ) THEN
406  CALL clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
407  ELSE
408  CALL clascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
409  END IF
410  END IF
411 *
412 * Call CHBTRD_HB2ST to reduce Hermitian band matrix to tridiagonal form.
413 *
414  inde = 1
415  indrwk = inde + n
416  llrwk = lrwork - indrwk + 1
417  indhous = 1
418  indwk = indhous + lhtrd
419  llwork = lwork - indwk + 1
420  indwk2 = indwk + n*n
421  llwk2 = lwork - indwk2 + 1
422 *
423  CALL chetrd_hb2st( "N", jobz, uplo, n, kd, ab, ldab, w,
424  $ rwork( inde ), work( indhous ), lhtrd,
425  $ work( indwk ), llwork, iinfo )
426 *
427 * For eigenvalues only, call SSTERF. For eigenvectors, call CSTEDC.
428 *
429  IF( .NOT.wantz ) THEN
430  CALL ssterf( n, w, rwork( inde ), info )
431  ELSE
432  CALL cstedc( 'I', n, w, rwork( inde ), work, n, work( indwk2 ),
433  $ llwk2, rwork( indrwk ), llrwk, iwork, liwork,
434  $ info )
435  CALL cgemm( 'N', 'N', n, n, n, cone, z, ldz, work, n, czero,
436  $ work( indwk2 ), n )
437  CALL clacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
438  END IF
439 *
440 * If matrix was scaled, then rescale eigenvalues appropriately.
441 *
442  IF( iscale.EQ.1 ) THEN
443  IF( info.EQ.0 ) THEN
444  imax = n
445  ELSE
446  imax = info - 1
447  END IF
448  CALL sscal( imax, one / sigma, w, 1 )
449  END IF
450 *
451  work( 1 ) = lwmin
452  rwork( 1 ) = lrwmin
453  iwork( 1 ) = liwmin
454  RETURN
455 *
456 * End of CHBEVD_2STAGE
457 *
458  END
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:145
subroutine chbevd_2stage(JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CHBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER...
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine cstedc(COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, INFO)
CSTEDC
Definition: cstedc.f:214
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189