LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zheev_2stage.f
Go to the documentation of this file.
1*> \brief <b> ZHEEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE 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 ZHEEV_2STAGE + dependencies
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheev_2stage.f">
13*> [TGZ]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheev_2stage.f">
15*> [ZIP]</a>
16*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheev_2stage.f">
17*> [TXT]</a>
18*> \endhtmlonly
19*
20* Definition:
21* ===========
22*
23* SUBROUTINE ZHEEV_2STAGE( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK,
24* RWORK, INFO )
25*
26* IMPLICIT NONE
27*
28* .. Scalar Arguments ..
29* CHARACTER JOBZ, UPLO
30* INTEGER INFO, LDA, LWORK, N
31* ..
32* .. Array Arguments ..
33* DOUBLE PRECISION RWORK( * ), W( * )
34* COMPLEX*16 A( LDA, * ), WORK( * )
35* ..
36*
37*
38*> \par Purpose:
39* =============
40*>
41*> \verbatim
42*>
43*> ZHEEV_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
44*> complex Hermitian 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,out] A
73*> \verbatim
74*> A is COMPLEX*16 array, dimension (LDA, N)
75*> On entry, the Hermitian matrix A. If UPLO = 'U', the
76*> leading N-by-N upper triangular part of A contains the
77*> upper triangular part of the matrix A. If UPLO = 'L',
78*> the leading N-by-N lower triangular part of A contains
79*> the lower triangular part of the matrix A.
80*> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
81*> orthonormal eigenvectors of the matrix A.
82*> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
83*> or the upper triangle (if UPLO='U') of A, including the
84*> diagonal, is destroyed.
85*> \endverbatim
86*>
87*> \param[in] LDA
88*> \verbatim
89*> LDA is INTEGER
90*> The leading dimension of the array A. LDA >= max(1,N).
91*> \endverbatim
92*>
93*> \param[out] W
94*> \verbatim
95*> W is DOUBLE PRECISION array, dimension (N)
96*> If INFO = 0, the eigenvalues in ascending order.
97*> \endverbatim
98*>
99*> \param[out] WORK
100*> \verbatim
101*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
102*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
103*> \endverbatim
104*>
105*> \param[in] LWORK
106*> \verbatim
107*> LWORK is INTEGER
108*> The length of the array WORK. LWORK >= 1, when N <= 1;
109*> otherwise
110*> If JOBZ = 'N' and N > 1, LWORK must be queried.
111*> LWORK = MAX(1, dimension) where
112*> dimension = max(stage1,stage2) + (KD+1)*N + N
113*> = N*KD + N*max(KD+1,FACTOPTNB)
114*> + max(2*KD*KD, KD*NTHREADS)
115*> + (KD+1)*N + N
116*> where KD is the blocking size of the reduction,
117*> FACTOPTNB is the blocking used by the QR or LQ
118*> algorithm, usually FACTOPTNB=128 is a good choice
119*> NTHREADS is the number of threads used when
120*> openMP compilation is enabled, otherwise =1.
121*> If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available
122*>
123*> If LWORK = -1, then a workspace query is assumed; the routine
124*> only calculates the optimal size of the WORK array, returns
125*> this value as the first entry of the WORK array, and no error
126*> message related to LWORK is issued by XERBLA.
127*> \endverbatim
128*>
129*> \param[out] RWORK
130*> \verbatim
131*> RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2))
132*> \endverbatim
133*>
134*> \param[out] INFO
135*> \verbatim
136*> INFO is INTEGER
137*> = 0: successful exit
138*> < 0: if INFO = -i, the i-th argument had an illegal value
139*> > 0: if INFO = i, the algorithm failed to converge; i
140*> off-diagonal elements of an intermediate tridiagonal
141*> form did not converge to zero.
142*> \endverbatim
143*
144* Authors:
145* ========
146*
147*> \author Univ. of Tennessee
148*> \author Univ. of California Berkeley
149*> \author Univ. of Colorado Denver
150*> \author NAG Ltd.
151*
152*> \ingroup heev_2stage
153*
154*> \par Further Details:
155* =====================
156*>
157*> \verbatim
158*>
159*> All details about the 2stage techniques are available in:
160*>
161*> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
162*> Parallel reduction to condensed forms for symmetric eigenvalue problems
163*> using aggregated fine-grained and memory-aware kernels. In Proceedings
164*> of 2011 International Conference for High Performance Computing,
165*> Networking, Storage and Analysis (SC '11), New York, NY, USA,
166*> Article 8 , 11 pages.
167*> http://doi.acm.org/10.1145/2063384.2063394
168*>
169*> A. Haidar, J. Kurzak, P. Luszczek, 2013.
170*> An improved parallel singular value algorithm and its implementation
171*> for multicore hardware, In Proceedings of 2013 International Conference
172*> for High Performance Computing, Networking, Storage and Analysis (SC '13).
173*> Denver, Colorado, USA, 2013.
174*> Article 90, 12 pages.
175*> http://doi.acm.org/10.1145/2503210.2503292
176*>
177*> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
178*> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
179*> calculations based on fine-grained memory aware tasks.
180*> International Journal of High Performance Computing Applications.
181*> Volume 28 Issue 2, Pages 196-209, May 2014.
182*> http://hpc.sagepub.com/content/28/2/196
183*>
184*> \endverbatim
185*
186* =====================================================================
187 SUBROUTINE zheev_2stage( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK,
188 $ RWORK, INFO )
189*
190 IMPLICIT NONE
191*
192* -- LAPACK driver routine --
193* -- LAPACK is a software package provided by Univ. of Tennessee, --
194* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
195*
196* .. Scalar Arguments ..
197 CHARACTER JOBZ, UPLO
198 INTEGER INFO, LDA, LWORK, N
199* ..
200* .. Array Arguments ..
201 DOUBLE PRECISION RWORK( * ), W( * )
202 COMPLEX*16 A( LDA, * ), WORK( * )
203* ..
204*
205* =====================================================================
206*
207* .. Parameters ..
208 DOUBLE PRECISION ZERO, ONE
209 parameter( zero = 0.0d0, one = 1.0d0 )
210 COMPLEX*16 CONE
211 parameter( cone = ( 1.0d0, 0.0d0 ) )
212* ..
213* .. Local Scalars ..
214 LOGICAL LOWER, LQUERY, WANTZ
215 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
216 $ llwork, lwmin, lhtrd, lwtrd, kd, ib, indhous
217 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
218 $ smlnum
219* ..
220* .. External Functions ..
221 LOGICAL LSAME
222 INTEGER ILAENV2STAGE
223 DOUBLE PRECISION DLAMCH, ZLANHE
224 EXTERNAL lsame, dlamch, zlanhe, ilaenv2stage
225* ..
226* .. External Subroutines ..
227 EXTERNAL dscal, dsterf, xerbla, zlascl, zsteqr,
229* ..
230* .. Intrinsic Functions ..
231 INTRINSIC dble, max, sqrt
232* ..
233* .. Executable Statements ..
234*
235* Test the input parameters.
236*
237 wantz = lsame( jobz, 'V' )
238 lower = lsame( uplo, 'L' )
239 lquery = ( lwork.EQ.-1 )
240*
241 info = 0
242 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
243 info = -1
244 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
245 info = -2
246 ELSE IF( n.LT.0 ) THEN
247 info = -3
248 ELSE IF( lda.LT.max( 1, n ) ) THEN
249 info = -5
250 END IF
251*
252 IF( info.EQ.0 ) THEN
253 kd = ilaenv2stage( 1, 'ZHETRD_2STAGE', jobz, n, -1, -1, -1 )
254 ib = ilaenv2stage( 2, 'ZHETRD_2STAGE', jobz, n, kd, -1, -1 )
255 lhtrd = ilaenv2stage( 3, 'ZHETRD_2STAGE', jobz, n, kd, ib, -1 )
256 lwtrd = ilaenv2stage( 4, 'ZHETRD_2STAGE', jobz, n, kd, ib, -1 )
257 lwmin = n + lhtrd + lwtrd
258 work( 1 ) = lwmin
259*
260 IF( lwork.LT.lwmin .AND. .NOT.lquery )
261 $ info = -8
262 END IF
263*
264 IF( info.NE.0 ) THEN
265 CALL xerbla( 'ZHEEV_2STAGE ', -info )
266 RETURN
267 ELSE IF( lquery ) THEN
268 RETURN
269 END IF
270*
271* Quick return if possible
272*
273 IF( n.EQ.0 ) THEN
274 RETURN
275 END IF
276*
277 IF( n.EQ.1 ) THEN
278 w( 1 ) = dble( a( 1, 1 ) )
279 work( 1 ) = 1
280 IF( wantz )
281 $ a( 1, 1 ) = cone
282 RETURN
283 END IF
284*
285* Get machine constants.
286*
287 safmin = dlamch( 'Safe minimum' )
288 eps = dlamch( 'Precision' )
289 smlnum = safmin / eps
290 bignum = one / smlnum
291 rmin = sqrt( smlnum )
292 rmax = sqrt( bignum )
293*
294* Scale matrix to allowable range, if necessary.
295*
296 anrm = zlanhe( 'M', uplo, n, a, lda, rwork )
297 iscale = 0
298 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
299 iscale = 1
300 sigma = rmin / anrm
301 ELSE IF( anrm.GT.rmax ) THEN
302 iscale = 1
303 sigma = rmax / anrm
304 END IF
305 IF( iscale.EQ.1 )
306 $ CALL zlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
307*
308* Call ZHETRD_2STAGE to reduce Hermitian matrix to tridiagonal form.
309*
310 inde = 1
311 indtau = 1
312 indhous = indtau + n
313 indwrk = indhous + lhtrd
314 llwork = lwork - indwrk + 1
315*
316 CALL zhetrd_2stage( jobz, uplo, n, a, lda, w, rwork( inde ),
317 $ work( indtau ), work( indhous ), lhtrd,
318 $ work( indwrk ), llwork, iinfo )
319*
320* For eigenvalues only, call DSTERF. For eigenvectors, first call
321* ZUNGTR to generate the unitary matrix, then call ZSTEQR.
322*
323 IF( .NOT.wantz ) THEN
324 CALL dsterf( n, w, rwork( inde ), info )
325 ELSE
326 CALL zungtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
327 $ llwork, iinfo )
328 indwrk = inde + n
329 CALL zsteqr( jobz, n, w, rwork( inde ), a, lda,
330 $ rwork( indwrk ), info )
331 END IF
332*
333* If matrix was scaled, then rescale eigenvalues appropriately.
334*
335 IF( iscale.EQ.1 ) THEN
336 IF( info.EQ.0 ) THEN
337 imax = n
338 ELSE
339 imax = info - 1
340 END IF
341 CALL dscal( imax, one / sigma, w, 1 )
342 END IF
343*
344* Set WORK(1) to optimal complex workspace size.
345*
346 work( 1 ) = lwmin
347*
348 RETURN
349*
350* End of ZHEEV_2STAGE
351*
352 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zheev_2stage(jobz, uplo, n, a, lda, w, work, lwork, rwork, info)
ZHEEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matr...
subroutine zhetrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
ZHETRD_2STAGE
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 dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine zsteqr(compz, n, d, e, z, ldz, work, info)
ZSTEQR
Definition zsteqr.f:132
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
subroutine zungtr(uplo, n, a, lda, tau, work, lwork, info)
ZUNGTR
Definition zungtr.f:123