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 November 2017
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.8.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: * November 2017
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 ILAENV2STAGE
246: DOUBLE PRECISION DLAMCH, ZLANHB
247: EXTERNAL LSAME, DLAMCH, ZLANHB, ILAENV2STAGE
248: * ..
249: * .. External Subroutines ..
250: EXTERNAL DSCAL, DSTERF, XERBLA, ZLASCL, ZSTEQR,
251: $ ZHETRD_2STAGE, ZHETRD_HB2ST
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 = ILAENV2STAGE( 2, 'ZHETRD_HB2ST', JOBZ,
285: $ N, KD, -1, -1 )
286: LHTRD = ILAENV2STAGE( 3, 'ZHETRD_HB2ST', JOBZ,
287: $ N, KD, IB, -1 )
288: LWTRD = ILAENV2STAGE( 4, 'ZHETRD_HB2ST', JOBZ,
289: $ N, KD, IB, -1 )
290: LWMIN = LHTRD + LWTRD
291: WORK( 1 ) = LWMIN
292: ENDIF
293: *
294: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY )
295: $ INFO = -11
296: END IF
297: *
298: IF( INFO.NE.0 ) THEN
299: CALL XERBLA( 'ZHBEV_2STAGE ', -INFO )
300: RETURN
301: ELSE IF( LQUERY ) THEN
302: RETURN
303: END IF
304: *
305: * Quick return if possible
306: *
307: IF( N.EQ.0 )
308: $ RETURN
309: *
310: IF( N.EQ.1 ) THEN
311: IF( LOWER ) THEN
312: W( 1 ) = DBLE( AB( 1, 1 ) )
313: ELSE
314: W( 1 ) = DBLE( AB( KD+1, 1 ) )
315: END IF
316: IF( WANTZ )
317: $ Z( 1, 1 ) = ONE
318: RETURN
319: END IF
320: *
321: * Get machine constants.
322: *
323: SAFMIN = DLAMCH( 'Safe minimum' )
324: EPS = DLAMCH( 'Precision' )
325: SMLNUM = SAFMIN / EPS
326: BIGNUM = ONE / SMLNUM
327: RMIN = SQRT( SMLNUM )
328: RMAX = SQRT( BIGNUM )
329: *
330: * Scale matrix to allowable range, if necessary.
331: *
332: ANRM = ZLANHB( 'M', UPLO, N, KD, AB, LDAB, RWORK )
333: ISCALE = 0
334: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
335: ISCALE = 1
336: SIGMA = RMIN / ANRM
337: ELSE IF( ANRM.GT.RMAX ) THEN
338: ISCALE = 1
339: SIGMA = RMAX / ANRM
340: END IF
341: IF( ISCALE.EQ.1 ) THEN
342: IF( LOWER ) THEN
343: CALL ZLASCL( 'B', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
344: ELSE
345: CALL ZLASCL( 'Q', KD, KD, ONE, SIGMA, N, N, AB, LDAB, INFO )
346: END IF
347: END IF
348: *
349: * Call ZHBTRD_HB2ST to reduce Hermitian band matrix to tridiagonal form.
350: *
351: INDE = 1
352: INDHOUS = 1
353: INDWRK = INDHOUS + LHTRD
354: LLWORK = LWORK - INDWRK + 1
355: *
356: CALL ZHETRD_HB2ST( "N", JOBZ, UPLO, N, KD, AB, LDAB, W,
357: $ RWORK( INDE ), WORK( INDHOUS ), LHTRD,
358: $ WORK( INDWRK ), LLWORK, IINFO )
359: *
360: * For eigenvalues only, call DSTERF. For eigenvectors, call ZSTEQR.
361: *
362: IF( .NOT.WANTZ ) THEN
363: CALL DSTERF( N, W, RWORK( INDE ), INFO )
364: ELSE
365: INDRWK = INDE + N
366: CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), Z, LDZ,
367: $ RWORK( INDRWK ), INFO )
368: END IF
369: *
370: * If matrix was scaled, then rescale eigenvalues appropriately.
371: *
372: IF( ISCALE.EQ.1 ) THEN
373: IF( INFO.EQ.0 ) THEN
374: IMAX = N
375: ELSE
376: IMAX = INFO - 1
377: END IF
378: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
379: END IF
380: *
381: * Set WORK(1) to optimal workspace size.
382: *
383: WORK( 1 ) = LWMIN
384: *
385: RETURN
386: *
387: * End of ZHBEV_2STAGE
388: *
389: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zhbev_2stage.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack