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: *> \ingroup complex16OTHEReigen
175: *
176: *> \par Further Details:
177: * =====================
178: *>
179: *> \verbatim
180: *>
181: *> All details about the 2stage techniques are available in:
182: *>
183: *> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
184: *> Parallel reduction to condensed forms for symmetric eigenvalue problems
185: *> using aggregated fine-grained and memory-aware kernels. In Proceedings
186: *> of 2011 International Conference for High Performance Computing,
187: *> Networking, Storage and Analysis (SC '11), New York, NY, USA,
188: *> Article 8 , 11 pages.
189: *> http://doi.acm.org/10.1145/2063384.2063394
190: *>
191: *> A. Haidar, J. Kurzak, P. Luszczek, 2013.
192: *> An improved parallel singular value algorithm and its implementation
193: *> for multicore hardware, In Proceedings of 2013 International Conference
194: *> for High Performance Computing, Networking, Storage and Analysis (SC '13).
195: *> Denver, Colorado, USA, 2013.
196: *> Article 90, 12 pages.
197: *> http://doi.acm.org/10.1145/2503210.2503292
198: *>
199: *> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
200: *> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
201: *> calculations based on fine-grained memory aware tasks.
202: *> International Journal of High Performance Computing Applications.
203: *> Volume 28 Issue 2, Pages 196-209, May 2014.
204: *> http://hpc.sagepub.com/content/28/2/196
205: *>
206: *> \endverbatim
207: *
208: * =====================================================================
209: SUBROUTINE ZHBEV_2STAGE( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
210: $ WORK, LWORK, RWORK, INFO )
211: *
212: IMPLICIT NONE
213: *
214: * -- LAPACK driver routine --
215: * -- LAPACK is a software package provided by Univ. of Tennessee, --
216: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217: *
218: * .. Scalar Arguments ..
219: CHARACTER JOBZ, UPLO
220: INTEGER INFO, KD, LDAB, LDZ, N, LWORK
221: * ..
222: * .. Array Arguments ..
223: DOUBLE PRECISION RWORK( * ), W( * )
224: COMPLEX*16 AB( LDAB, * ), WORK( * ), Z( LDZ, * )
225: * ..
226: *
227: * =====================================================================
228: *
229: * .. Parameters ..
230: DOUBLE PRECISION ZERO, ONE
231: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
232: * ..
233: * .. Local Scalars ..
234: LOGICAL LOWER, WANTZ, LQUERY
235: INTEGER IINFO, IMAX, INDE, INDWRK, INDRWK, ISCALE,
236: $ LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS
237: DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
238: $ SMLNUM
239: * ..
240: * .. External Functions ..
241: LOGICAL LSAME
242: INTEGER ILAENV2STAGE
243: DOUBLE PRECISION DLAMCH, ZLANHB
244: EXTERNAL LSAME, DLAMCH, ZLANHB, ILAENV2STAGE
245: * ..
246: * .. External Subroutines ..
247: EXTERNAL DSCAL, DSTERF, XERBLA, ZLASCL, ZSTEQR,
248: $ ZHETRD_2STAGE, ZHETRD_HB2ST
249: * ..
250: * .. Intrinsic Functions ..
251: INTRINSIC DBLE, SQRT
252: * ..
253: * .. Executable Statements ..
254: *
255: * Test the input parameters.
256: *
257: WANTZ = LSAME( JOBZ, 'V' )
258: LOWER = LSAME( UPLO, 'L' )
259: LQUERY = ( LWORK.EQ.-1 )
260: *
261: INFO = 0
262: IF( .NOT.( LSAME( JOBZ, 'N' ) ) ) THEN
263: INFO = -1
264: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
265: INFO = -2
266: ELSE IF( N.LT.0 ) THEN
267: INFO = -3
268: ELSE IF( KD.LT.0 ) THEN
269: INFO = -4
270: ELSE IF( LDAB.LT.KD+1 ) THEN
271: INFO = -6
272: ELSE IF( LDZ.LT.1 .OR. ( WANTZ .AND. LDZ.LT.N ) ) THEN
273: INFO = -9
274: END IF
275: *
276: IF( INFO.EQ.0 ) THEN
277: IF( N.LE.1 ) THEN
278: LWMIN = 1
279: WORK( 1 ) = LWMIN
280: ELSE
281: IB = ILAENV2STAGE( 2, 'ZHETRD_HB2ST', JOBZ,
282: $ N, KD, -1, -1 )
283: LHTRD = ILAENV2STAGE( 3, 'ZHETRD_HB2ST', JOBZ,
284: $ N, KD, IB, -1 )
285: LWTRD = ILAENV2STAGE( 4, 'ZHETRD_HB2ST', JOBZ,
286: $ 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>
![[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