1: *> \brief <b> ZHBEV_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 ZHBEV_2STAGE + dependencies
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhbev_2stage.f">
13: *> [TGZ]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhbev_2stage.f">
15: *> [ZIP]</a>
16: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhbev_2stage.f">
17: *> [TXT]</a>
18: *> \endhtmlonly
19: *
20: * Definition:
21: * ===========
22: *
23: * SUBROUTINE ZHBEV_2STAGE( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
24: * WORK, LWORK, RWORK, INFO )
25: *
26: * IMPLICIT NONE
27: *
28: * .. Scalar Arguments ..
29: * CHARACTER JOBZ, UPLO
30: * INTEGER INFO, KD, LDAB, LDZ, N, LWORK
31: * ..
32: * .. Array Arguments ..
33: * DOUBLE PRECISION RWORK( * ), W( * )
34: * COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
35: * ..
36: *
37: *
38: *> \par Purpose:
39: * =============
40: *>
41: *> \verbatim
42: *>
43: *> ZHBEV_2STAGE computes all the eigenvalues and, optionally, eigenvectors of
44: *> a complex Hermitian band matrix A using the 2stage technique for
45: *> the reduction to tridiagonal.
46: *> \endverbatim
47: *
48: * Arguments:
49: * ==========
50: *
51: *> \param[in] JOBZ
52: *> \verbatim
53: *> JOBZ is CHARACTER*1
54: *> = 'N': Compute eigenvalues only;
55: *> = 'V': Compute eigenvalues and eigenvectors.
56: *> Not available in this release.
57: *> \endverbatim
58: *>
59: *> \param[in] UPLO
60: *> \verbatim
61: *> UPLO is CHARACTER*1
62: *> = 'U': Upper triangle of A is stored;
63: *> = 'L': Lower triangle of A is stored.
64: *> \endverbatim
65: *>
66: *> \param[in] N
67: *> \verbatim
68: *> N is INTEGER
69: *> The order of the matrix A. N >= 0.
70: *> \endverbatim
71: *>
72: *> \param[in] KD
73: *> \verbatim
74: *> KD is INTEGER
75: *> The number of superdiagonals of the matrix A if UPLO = 'U',
76: *> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
77: *> \endverbatim
78: *>
79: *> \param[in,out] AB
80: *> \verbatim
81: *> AB is COMPLEX*16 array, dimension (LDAB, N)
82: *> On entry, the upper or lower triangle of the Hermitian band
83: *> matrix A, stored in the first KD+1 rows of the array. The
84: *> j-th column of A is stored in the j-th column of the array AB
85: *> as follows:
86: *> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
87: *> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
88: *>
89: *> On exit, AB is overwritten by values generated during the
90: *> reduction to tridiagonal form. If UPLO = 'U', the first
91: *> superdiagonal and the diagonal of the tridiagonal matrix T
92: *> are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
93: *> the diagonal and first subdiagonal of T are returned in the
94: *> first two rows of AB.
95: *> \endverbatim
96: *>
97: *> \param[in] LDAB
98: *> \verbatim
99: *> LDAB is INTEGER
100: *> The leading dimension of the array AB. LDAB >= KD + 1.
101: *> \endverbatim
102: *>
103: *> \param[out] W
104: *> \verbatim
105: *> W is DOUBLE PRECISION array, dimension (N)
106: *> If INFO = 0, the eigenvalues in ascending order.
107: *> \endverbatim
108: *>
109: *> \param[out] Z
110: *> \verbatim
111: *> Z is COMPLEX*16 array, dimension (LDZ, N)
112: *> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
113: *> eigenvectors of the matrix A, with the i-th column of Z
114: *> holding the eigenvector associated with W(i).
115: *> If JOBZ = 'N', then Z is not referenced.
116: *> \endverbatim
117: *>
118: *> \param[in] LDZ
119: *> \verbatim
120: *> LDZ is INTEGER
121: *> The leading dimension of the array Z. LDZ >= 1, and if
122: *> JOBZ = 'V', LDZ >= max(1,N).
123: *> \endverbatim
124: *>
125: *> \param[out] WORK
126: *> \verbatim
127: *> WORK is COMPLEX*16 array, dimension LWORK
128: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
129: *> \endverbatim
130: *>
131: *> \param[in] LWORK
132: *> \verbatim
133: *> LWORK is INTEGER
134: *> The length of the array WORK. LWORK >= 1, when N <= 1;
135: *> otherwise
136: *> If JOBZ = 'N' and N > 1, LWORK must be queried.
137: *> LWORK = MAX(1, dimension) where
138: *> dimension = (2KD+1)*N + KD*NTHREADS
139: *> where KD is the size of the band.
140: *> NTHREADS is the number of threads used when
141: *> openMP compilation is enabled, otherwise =1.
142: *> If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available.
143: *>
144: *> If LWORK = -1, then a workspace query is assumed; the routine
145: *> only calculates the optimal sizes of the WORK, RWORK and
146: *> IWORK arrays, returns these values as the first entries of
147: *> the WORK, RWORK and IWORK arrays, and no error message
148: *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
149: *> \endverbatim
150: *>
151: *> \param[out] RWORK
152: *> \verbatim
153: *> RWORK is DOUBLE PRECISION array, dimension (max(1,3*N-2))
154: *> \endverbatim
155: *>
156: *> \param[out] INFO
157: *> \verbatim
158: *> INFO is INTEGER
159: *> = 0: successful exit.
160: *> < 0: if INFO = -i, the i-th argument had an illegal value.
161: *> > 0: if INFO = i, the algorithm failed to converge; i
162: *> off-diagonal elements of an intermediate tridiagonal
163: *> form did not converge to zero.
164: *> \endverbatim
165: *
166: * Authors:
167: * ========
168: *
169: *> \author Univ. of Tennessee
170: *> \author Univ. of California Berkeley
171: *> \author Univ. of Colorado Denver
172: *> \author NAG Ltd.
173: *
174: *> \date December 2016
175: *
176: *> \ingroup complex16OTHEReigen
177: *
178: *> \par Further Details:
179: * =====================
180: *>
181: *> \verbatim
182: *>
183: *> All details about the 2stage techniques are available in:
184: *>
185: *> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
186: *> Parallel reduction to condensed forms for symmetric eigenvalue problems
187: *> using aggregated fine-grained and memory-aware kernels. In Proceedings
188: *> of 2011 International Conference for High Performance Computing,
189: *> Networking, Storage and Analysis (SC '11), New York, NY, USA,
190: *> Article 8 , 11 pages.
191: *> http://doi.acm.org/10.1145/2063384.2063394
192: *>
193: *> A. Haidar, J. Kurzak, P. Luszczek, 2013.
194: *> An improved parallel singular value algorithm and its implementation
195: *> for multicore hardware, In Proceedings of 2013 International Conference
196: *> for High Performance Computing, Networking, Storage and Analysis (SC '13).
197: *> Denver, Colorado, USA, 2013.
198: *> Article 90, 12 pages.
199: *> http://doi.acm.org/10.1145/2503210.2503292
200: *>
201: *> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
202: *> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
203: *> calculations based on fine-grained memory aware tasks.
204: *> International Journal of High Performance Computing Applications.
205: *> Volume 28 Issue 2, Pages 196-209, May 2014.
206: *> http://hpc.sagepub.com/content/28/2/196
207: *>
208: *> \endverbatim
209: *
210: * =====================================================================
211: SUBROUTINE ZHBEV_2STAGE( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
212: $ WORK, LWORK, RWORK, INFO )
213: *
214: IMPLICIT NONE
215: *
216: * -- LAPACK driver routine (version 3.7.0) --
217: * -- LAPACK is a software package provided by Univ. of Tennessee, --
218: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
219: * December 2016
220: *
221: * .. Scalar Arguments ..
222: CHARACTER JOBZ, UPLO
223: INTEGER INFO, KD, LDAB, LDZ, N, LWORK
224: * ..
225: * .. Array Arguments ..
226: DOUBLE PRECISION RWORK( * ), W( * )
227: COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
228: * ..
229: *
230: * =====================================================================
231: *
232: * .. Parameters ..
233: DOUBLE PRECISION ZERO, ONE
234: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
235: * ..
236: * .. Local Scalars ..
237: LOGICAL LOWER, WANTZ, LQUERY
238: INTEGER IINFO, IMAX, INDE, INDWRK, INDRWK, ISCALE,
239: $ LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS
240: DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
241: $ SMLNUM
242: * ..
243: * .. External Functions ..
244: LOGICAL LSAME
245: INTEGER ILAENV
246: DOUBLE PRECISION DLAMCH, ZLANHB
247: EXTERNAL LSAME, DLAMCH, ZLANHB, ILAENV
248: * ..
249: * .. External Subroutines ..
250: EXTERNAL DSCAL, DSTERF, XERBLA, ZLASCL, ZSTEQR
251: $ ZHETRD_2STAGE
252: * ..
253: * .. Intrinsic Functions ..
254: INTRINSIC DBLE, SQRT
255: * ..
256: * .. Executable Statements ..
257: *
258: * Test the input parameters.
259: *
260: WANTZ = LSAME( JOBZ, 'V' )
261: LOWER = LSAME( UPLO, 'L' )
262: LQUERY = ( LWORK.EQ.-1 )
263: *
264: INFO = 0
265: IF( .NOT.( LSAME( JOBZ, 'N' ) ) ) THEN
266: INFO = -1
267: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
268: INFO = -2
269: ELSE IF( N.LT.0 ) THEN
270: INFO = -3
271: ELSE IF( KD.LT.0 ) THEN
272: INFO = -4
273: ELSE IF( LDAB.LT.KD+1 ) THEN
274: INFO = -6
275: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
276: INFO = -9
277: END IF
278: *
279: IF( INFO.EQ.0 ) THEN
280: IF( N.LE.1 ) THEN
281: LWMIN = 1
282: WORK( 1 ) = LWMIN
283: ELSE
284: IB = ILAENV( 18, 'ZHETRD_HB2ST', JOBZ, N, KD, -1, -1 )
285: LHTRD = ILAENV( 19, 'ZHETRD_HB2ST', JOBZ, N, KD, IB, -1 )
286: LWTRD = ILAENV( 20, 'ZHETRD_HB2ST', JOBZ, N, KD, IB, -1 )
287: LWMIN = LHTRD + LWTRD
288: WORK( 1 ) = LWMIN
289: ENDIF
290: *
291: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY )
292: $ INFO = -11
293: END IF
294: *
295: IF( INFO.NE.0 ) THEN
296: CALL XERBLA( 'ZHBEV_2STAGE ', -INFO )
297: RETURN
298: ELSE IF( LQUERY ) THEN
299: RETURN
300: END IF
301: *
302: * Quick return if possible
303: *
304: IF( N.EQ.0 )
305: $ RETURN
306: *
307: IF( N.EQ.1 ) THEN
308: IF( LOWER ) THEN
309: W( 1 ) = DBLE( AB( 1, 1 ) )
310: ELSE
311: W( 1 ) = DBLE( AB( KD+1, 1 ) )
312: END IF
313: IF( WANTZ )
314: $ Z( 1, 1 ) = ONE
315: RETURN
316: END IF
317: *
318: * Get machine constants.
319: *
320: SAFMIN = DLAMCH( 'Safe minimum' )
321: EPS = DLAMCH( 'Precision' )
322: SMLNUM = SAFMIN / EPS
323: BIGNUM = ONE / SMLNUM
324: RMIN = SQRT( SMLNUM )
325: RMAX = SQRT( BIGNUM )
326: *
327: * Scale matrix to allowable range, if necessary.
328: *
329: ANRM = ZLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK )
330: ISCALE = 0
331: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
332: ISCALE = 1
333: SIGMA = RMIN / ANRM
334: ELSE IF( ANRM.GT.RMAX ) THEN
335: ISCALE = 1
336: SIGMA = RMAX / ANRM
337: END IF
338: IF( ISCALE.EQ.1 ) THEN
339: IF( LOWER ) THEN
340: CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
341: ELSE
342: CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
343: END IF
344: END IF
345: *
346: * Call ZHBTRD_HB2ST to reduce Hermitian band matrix to tridiagonal form.
347: *
348: INDE = 1
349: INDHOUS = 1
350: INDWRK = INDHOUS + LHTRD
351: LLWORK = LWORK - INDWRK + 1
352: *
353: CALL ZHETRD_HB2ST( "N", JOBZ, UPLO, N, KD, AB, LDAB, W,
354: $ RWORK( INDE ), WORK( INDHOUS ), LHTRD,
355: $ WORK( INDWRK ), LLWORK, IINFO )
356: *
357: * For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEQR.
358: *
359: IF( .NOT.WANTZ ) THEN
360: CALL DSTERF( N, W, RWORK( INDE ), INFO )
361: ELSE
362: INDRWK = INDE + N
363: CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), Z, LDZ,
364: $ RWORK( INDRWK ), INFO )
365: END IF
366: *
367: * If matrix was scaled, then rescale eigenvalues appropriately.
368: *
369: IF( ISCALE.EQ.1 ) THEN
370: IF( INFO.EQ.0 ) THEN
371: IMAX = N
372: ELSE
373: IMAX = INFO - 1
374: END IF
375: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
376: END IF
377: *
378: * Set WORK(1) to optimal workspace size.
379: *
380: WORK( 1 ) = LWMIN
381: *
382: RETURN
383: *
384: * End of ZHBEV_2STAGE
385: *
386: END
CVSweb interface <joel.bertrand@systella.fr>