LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
zhegv_2stage.f
Go to the documentation of this file.
1*> \brief \b ZHEGV_2STAGE
2*
3* @precisions fortran z -> c
4*
5* =========== DOCUMENTATION ===========
6*
7* Online html documentation available at
8* http://www.netlib.org/lapack/explore-html/
9*
10*> Download ZHEGV_2STAGE + dependencies
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhegv_2stage.f">
12*> [TGZ]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhegv_2stage.f">
14*> [ZIP]</a>
15*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhegv_2stage.f">
16*> [TXT]</a>
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZHEGV_2STAGE( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W,
22* WORK, LWORK, RWORK, INFO )
23*
24* IMPLICIT NONE
25*
26* .. Scalar Arguments ..
27* CHARACTER JOBZ, UPLO
28* INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
29* ..
30* .. Array Arguments ..
31* DOUBLE PRECISION RWORK( * ), W( * )
32* COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> ZHEGV_2STAGE computes all the eigenvalues, and optionally, the eigenvectors
42*> of a complex generalized Hermitian-definite eigenproblem, of the form
43*> A*x=(lambda)*B*x, A*Bx=(lambda)*x, or B*A*x=(lambda)*x.
44*> Here A and B are assumed to be Hermitian and B is also
45*> positive definite.
46*> This routine use the 2stage technique for the reduction to tridiagonal
47*> which showed higher performance on recent architecture and for large
48*> sizes N>2000.
49*> \endverbatim
50*
51* Arguments:
52* ==========
53*
54*> \param[in] ITYPE
55*> \verbatim
56*> ITYPE is INTEGER
57*> Specifies the problem type to be solved:
58*> = 1: A*x = (lambda)*B*x
59*> = 2: A*B*x = (lambda)*x
60*> = 3: B*A*x = (lambda)*x
61*> \endverbatim
62*>
63*> \param[in] JOBZ
64*> \verbatim
65*> JOBZ is CHARACTER*1
66*> = 'N': Compute eigenvalues only;
67*> = 'V': Compute eigenvalues and eigenvectors.
68*> Not available in this release.
69*> \endverbatim
70*>
71*> \param[in] UPLO
72*> \verbatim
73*> UPLO is CHARACTER*1
74*> = 'U': Upper triangles of A and B are stored;
75*> = 'L': Lower triangles of A and B are stored.
76*> \endverbatim
77*>
78*> \param[in] N
79*> \verbatim
80*> N is INTEGER
81*> The order of the matrices A and B. N >= 0.
82*> \endverbatim
83*>
84*> \param[in,out] A
85*> \verbatim
86*> A is COMPLEX*16 array, dimension (LDA, N)
87*> On entry, the Hermitian matrix A. If UPLO = 'U', the
88*> leading N-by-N upper triangular part of A contains the
89*> upper triangular part of the matrix A. If UPLO = 'L',
90*> the leading N-by-N lower triangular part of A contains
91*> the lower triangular part of the matrix A.
92*>
93*> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
94*> matrix Z of eigenvectors. The eigenvectors are normalized
95*> as follows:
96*> if ITYPE = 1 or 2, Z**H*B*Z = I;
97*> if ITYPE = 3, Z**H*inv(B)*Z = I.
98*> If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
99*> or the lower triangle (if UPLO='L') of A, including the
100*> diagonal, is destroyed.
101*> \endverbatim
102*>
103*> \param[in] LDA
104*> \verbatim
105*> LDA is INTEGER
106*> The leading dimension of the array A. LDA >= max(1,N).
107*> \endverbatim
108*>
109*> \param[in,out] B
110*> \verbatim
111*> B is COMPLEX*16 array, dimension (LDB, N)
112*> On entry, the Hermitian positive definite matrix B.
113*> If UPLO = 'U', the leading N-by-N upper triangular part of B
114*> contains the upper triangular part of the matrix B.
115*> If UPLO = 'L', the leading N-by-N lower triangular part of B
116*> contains the lower triangular part of the matrix B.
117*>
118*> On exit, if INFO <= N, the part of B containing the matrix is
119*> overwritten by the triangular factor U or L from the Cholesky
120*> factorization B = U**H*U or B = L*L**H.
121*> \endverbatim
122*>
123*> \param[in] LDB
124*> \verbatim
125*> LDB is INTEGER
126*> The leading dimension of the array B. LDB >= max(1,N).
127*> \endverbatim
128*>
129*> \param[out] W
130*> \verbatim
131*> W is DOUBLE PRECISION array, dimension (N)
132*> If INFO = 0, the eigenvalues in ascending order.
133*> \endverbatim
134*>
135*> \param[out] WORK
136*> \verbatim
137*> WORK is COMPLEX*16 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 = max(stage1,stage2) + (KD+1)*N + N
149*> = N*KD + N*max(KD+1,FACTOPTNB)
150*> + max(2*KD*KD, KD*NTHREADS)
151*> + (KD+1)*N + N
152*> where KD is the blocking size of the reduction,
153*> FACTOPTNB is the blocking used by the QR or LQ
154*> algorithm, usually FACTOPTNB=128 is a good choice
155*> NTHREADS is the number of threads used when
156*> openMP compilation is enabled, otherwise =1.
157*> If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available
158*>
159*> If LWORK = -1, then a workspace query is assumed; the routine
160*> only calculates the optimal size of the WORK array, returns
161*> this value as the first entry of the WORK array, and no error
162*> message related to LWORK is issued by XERBLA.
163*> \endverbatim
164*>
165*> \param[out] RWORK
166*> \verbatim
167*> RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2))
168*> \endverbatim
169*>
170*> \param[out] INFO
171*> \verbatim
172*> INFO is INTEGER
173*> = 0: successful exit
174*> < 0: if INFO = -i, the i-th argument had an illegal value
175*> > 0: ZPOTRF or ZHEEV returned an error code:
176*> <= N: if INFO = i, ZHEEV failed to converge;
177*> i off-diagonal elements of an intermediate
178*> tridiagonal form did not converge to zero;
179*> > N: if INFO = N + i, for 1 <= i <= N, then the leading
180*> principal minor of order i of B is not positive.
181*> The factorization of B could not be completed and
182*> no eigenvalues or eigenvectors were computed.
183*> \endverbatim
184*
185* Authors:
186* ========
187*
188*> \author Univ. of Tennessee
189*> \author Univ. of California Berkeley
190*> \author Univ. of Colorado Denver
191*> \author NAG Ltd.
192*
193*> \ingroup hegv_2stage
194*
195*> \par Further Details:
196* =====================
197*>
198*> \verbatim
199*>
200*> All details about the 2stage techniques are available in:
201*>
202*> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
203*> Parallel reduction to condensed forms for symmetric eigenvalue problems
204*> using aggregated fine-grained and memory-aware kernels. In Proceedings
205*> of 2011 International Conference for High Performance Computing,
206*> Networking, Storage and Analysis (SC '11), New York, NY, USA,
207*> Article 8 , 11 pages.
208*> http://doi.acm.org/10.1145/2063384.2063394
209*>
210*> A. Haidar, J. Kurzak, P. Luszczek, 2013.
211*> An improved parallel singular value algorithm and its implementation
212*> for multicore hardware, In Proceedings of 2013 International Conference
213*> for High Performance Computing, Networking, Storage and Analysis (SC '13).
214*> Denver, Colorado, USA, 2013.
215*> Article 90, 12 pages.
216*> http://doi.acm.org/10.1145/2503210.2503292
217*>
218*> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
219*> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
220*> calculations based on fine-grained memory aware tasks.
221*> International Journal of High Performance Computing Applications.
222*> Volume 28 Issue 2, Pages 196-209, May 2014.
223*> http://hpc.sagepub.com/content/28/2/196
224*>
225*> \endverbatim
226*
227* =====================================================================
228 SUBROUTINE zhegv_2stage( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB,
229 $ W,
230 $ WORK, LWORK, RWORK, INFO )
231*
232 IMPLICIT NONE
233*
234* -- LAPACK driver routine --
235* -- LAPACK is a software package provided by Univ. of Tennessee, --
236* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
237*
238* .. Scalar Arguments ..
239 CHARACTER JOBZ, UPLO
240 INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
241* ..
242* .. Array Arguments ..
243 DOUBLE PRECISION RWORK( * ), W( * )
244 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
245* ..
246*
247* =====================================================================
248*
249* .. Parameters ..
250 COMPLEX*16 ONE
251 PARAMETER ( ONE = ( 1.0d+0, 0.0d+0 ) )
252* ..
253* .. Local Scalars ..
254 LOGICAL LQUERY, UPPER, WANTZ
255 CHARACTER TRANS
256 INTEGER NEIG, LWMIN, LHTRD, LWTRD, KD, IB
257* ..
258* .. External Functions ..
259 LOGICAL LSAME
260 INTEGER ILAENV2STAGE
261 EXTERNAL lsame, ilaenv2stage
262* ..
263* .. External Subroutines ..
264 EXTERNAL xerbla, zhegst, zpotrf, ztrmm,
265 $ ztrsm,
267* ..
268* .. Intrinsic Functions ..
269 INTRINSIC max
270* ..
271* .. Executable Statements ..
272*
273* Test the input parameters.
274*
275 wantz = lsame( jobz, 'V' )
276 upper = lsame( uplo, 'U' )
277 lquery = ( lwork.EQ.-1 )
278*
279 info = 0
280 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
281 info = -1
282 ELSE IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
283 info = -2
284 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
285 info = -3
286 ELSE IF( n.LT.0 ) THEN
287 info = -4
288 ELSE IF( lda.LT.max( 1, n ) ) THEN
289 info = -6
290 ELSE IF( ldb.LT.max( 1, n ) ) THEN
291 info = -8
292 END IF
293*
294 IF( info.EQ.0 ) THEN
295 kd = ilaenv2stage( 1, 'ZHETRD_2STAGE', jobz, n, -1, -1,
296 $ -1 )
297 ib = ilaenv2stage( 2, 'ZHETRD_2STAGE', jobz, n, kd, -1,
298 $ -1 )
299 lhtrd = ilaenv2stage( 3, 'ZHETRD_2STAGE', jobz, n, kd, ib,
300 $ -1 )
301 lwtrd = ilaenv2stage( 4, 'ZHETRD_2STAGE', jobz, n, kd, ib,
302 $ -1 )
303 lwmin = n + lhtrd + lwtrd
304 work( 1 ) = lwmin
305*
306 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
307 info = -11
308 END IF
309 END IF
310*
311 IF( info.NE.0 ) THEN
312 CALL xerbla( 'ZHEGV_2STAGE ', -info )
313 RETURN
314 ELSE IF( lquery ) THEN
315 RETURN
316 END IF
317*
318* Quick return if possible
319*
320 IF( n.EQ.0 )
321 $ RETURN
322*
323* Form a Cholesky factorization of B.
324*
325 CALL zpotrf( uplo, n, b, ldb, info )
326 IF( info.NE.0 ) THEN
327 info = n + info
328 RETURN
329 END IF
330*
331* Transform problem to standard eigenvalue problem and solve.
332*
333 CALL zhegst( itype, uplo, n, a, lda, b, ldb, info )
334 CALL zheev_2stage( jobz, uplo, n, a, lda, w,
335 $ work, lwork, rwork, info )
336*
337 IF( wantz ) THEN
338*
339* Backtransform eigenvectors to the original problem.
340*
341 neig = n
342 IF( info.GT.0 )
343 $ neig = info - 1
344 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
345*
346* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
347* backtransform eigenvectors: x = inv(L)**H *y or inv(U)*y
348*
349 IF( upper ) THEN
350 trans = 'N'
351 ELSE
352 trans = 'C'
353 END IF
354*
355 CALL ztrsm( 'Left', uplo, trans, 'Non-unit', n, neig,
356 $ one,
357 $ b, ldb, a, lda )
358*
359 ELSE IF( itype.EQ.3 ) THEN
360*
361* For B*A*x=(lambda)*x;
362* backtransform eigenvectors: x = L*y or U**H *y
363*
364 IF( upper ) THEN
365 trans = 'C'
366 ELSE
367 trans = 'N'
368 END IF
369*
370 CALL ztrmm( 'Left', uplo, trans, 'Non-unit', n, neig,
371 $ one,
372 $ b, ldb, a, lda )
373 END IF
374 END IF
375*
376 work( 1 ) = lwmin
377*
378 RETURN
379*
380* End of ZHEGV_2STAGE
381*
382 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 zhegst(itype, uplo, n, a, lda, b, ldb, info)
ZHEGST
Definition zhegst.f:126
subroutine zhegv_2stage(itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, rwork, info)
ZHEGV_2STAGE
subroutine zpotrf(uplo, n, a, lda, info)
ZPOTRF
Definition zpotrf.f:105
subroutine ztrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRMM
Definition ztrmm.f:177
subroutine ztrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
ZTRSM
Definition ztrsm.f:180