LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
cheev_2stage.f
Go to the documentation of this file.
1 *> \brief <b> CHEEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices</b>
2 *
3 * @generated from zheev_2stage.f, fortran z -> c, Sat Nov 5 23:18:06 2016
4 *
5 * =========== DOCUMENTATION ===========
6 *
7 * Online html documentation available at
8 * http://www.netlib.org/lapack/explore-html/
9 *
10 *> \htmlonly
11 *> Download CHEEV_2STAGE + dependencies
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cheev_2stage.f">
13 *> [TGZ]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cheev_2stage.f">
15 *> [ZIP]</a>
16 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cheev_2stage.f">
17 *> [TXT]</a>
18 *> \endhtmlonly
19 *
20 * Definition:
21 * ===========
22 *
23 * SUBROUTINE CHEEV_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 * REAL RWORK( * ), W( * )
34 * COMPLEX A( LDA, * ), WORK( * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> CHEEV_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 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 REAL 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 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 REAL 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 *> \date November 2017
153 *
154 *> \ingroup complexHEeigen
155 *
156 *> \par Further Details:
157 * =====================
158 *>
159 *> \verbatim
160 *>
161 *> All details about the 2stage techniques are available in:
162 *>
163 *> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
164 *> Parallel reduction to condensed forms for symmetric eigenvalue problems
165 *> using aggregated fine-grained and memory-aware kernels. In Proceedings
166 *> of 2011 International Conference for High Performance Computing,
167 *> Networking, Storage and Analysis (SC '11), New York, NY, USA,
168 *> Article 8 , 11 pages.
169 *> http://doi.acm.org/10.1145/2063384.2063394
170 *>
171 *> A. Haidar, J. Kurzak, P. Luszczek, 2013.
172 *> An improved parallel singular value algorithm and its implementation
173 *> for multicore hardware, In Proceedings of 2013 International Conference
174 *> for High Performance Computing, Networking, Storage and Analysis (SC '13).
175 *> Denver, Colorado, USA, 2013.
176 *> Article 90, 12 pages.
177 *> http://doi.acm.org/10.1145/2503210.2503292
178 *>
179 *> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
180 *> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
181 *> calculations based on fine-grained memory aware tasks.
182 *> International Journal of High Performance Computing Applications.
183 *> Volume 28 Issue 2, Pages 196-209, May 2014.
184 *> http://hpc.sagepub.com/content/28/2/196
185 *>
186 *> \endverbatim
187 *
188 * =====================================================================
189  SUBROUTINE cheev_2stage( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK,
190  $ RWORK, INFO )
191 *
192  IMPLICIT NONE
193 *
194 * -- LAPACK driver routine (version 3.8.0) --
195 * -- LAPACK is a software package provided by Univ. of Tennessee, --
196 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
197 * November 2017
198 *
199 * .. Scalar Arguments ..
200  CHARACTER JOBZ, UPLO
201  INTEGER INFO, LDA, LWORK, N
202 * ..
203 * .. Array Arguments ..
204  REAL RWORK( * ), W( * )
205  COMPLEX A( lda, * ), WORK( * )
206 * ..
207 *
208 * =====================================================================
209 *
210 * .. Parameters ..
211  REAL ZERO, ONE
212  parameter( zero = 0.0e0, one = 1.0e0 )
213  COMPLEX CONE
214  parameter( cone = ( 1.0e0, 0.0e0 ) )
215 * ..
216 * .. Local Scalars ..
217  LOGICAL LOWER, LQUERY, WANTZ
218  INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
219  $ llwork, lwmin, lhtrd, lwtrd, kd, ib, indhous
220  REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
221  $ smlnum
222 * ..
223 * .. External Functions ..
224  LOGICAL LSAME
225  INTEGER ILAENV2STAGE
226  REAL SLAMCH, CLANHE
227  EXTERNAL lsame, slamch, clanhe, ilaenv2stage
228 * ..
229 * .. External Subroutines ..
230  EXTERNAL sscal, ssterf, xerbla, clascl, csteqr,
232 * ..
233 * .. Intrinsic Functions ..
234  INTRINSIC REAL, MAX, SQRT
235 * ..
236 * .. Executable Statements ..
237 *
238 * Test the input parameters.
239 *
240  wantz = lsame( jobz, 'V' )
241  lower = lsame( uplo, 'L' )
242  lquery = ( lwork.EQ.-1 )
243 *
244  info = 0
245  IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
246  info = -1
247  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
248  info = -2
249  ELSE IF( n.LT.0 ) THEN
250  info = -3
251  ELSE IF( lda.LT.max( 1, n ) ) THEN
252  info = -5
253  END IF
254 *
255  IF( info.EQ.0 ) THEN
256  kd = ilaenv2stage( 1, 'CHETRD_2STAGE', jobz, n, -1, -1, -1 )
257  ib = ilaenv2stage( 2, 'CHETRD_2STAGE', jobz, n, kd, -1, -1 )
258  lhtrd = ilaenv2stage( 3, 'CHETRD_2STAGE', jobz, n, kd, ib, -1 )
259  lwtrd = ilaenv2stage( 4, 'CHETRD_2STAGE', jobz, n, kd, ib, -1 )
260  lwmin = n + lhtrd + lwtrd
261  work( 1 ) = lwmin
262 *
263  IF( lwork.LT.lwmin .AND. .NOT.lquery )
264  $ info = -8
265  END IF
266 *
267  IF( info.NE.0 ) THEN
268  CALL xerbla( 'CHEEV_2STAGE ', -info )
269  RETURN
270  ELSE IF( lquery ) THEN
271  RETURN
272  END IF
273 *
274 * Quick return if possible
275 *
276  IF( n.EQ.0 ) THEN
277  RETURN
278  END IF
279 *
280  IF( n.EQ.1 ) THEN
281  w( 1 ) = REAL( A( 1, 1 ) )
282  work( 1 ) = 1
283  IF( wantz )
284  $ a( 1, 1 ) = cone
285  RETURN
286  END IF
287 *
288 * Get machine constants.
289 *
290  safmin = slamch( 'Safe minimum' )
291  eps = slamch( 'Precision' )
292  smlnum = safmin / eps
293  bignum = one / smlnum
294  rmin = sqrt( smlnum )
295  rmax = sqrt( bignum )
296 *
297 * Scale matrix to allowable range, if necessary.
298 *
299  anrm = clanhe( 'M', uplo, n, a, lda, rwork )
300  iscale = 0
301  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
302  iscale = 1
303  sigma = rmin / anrm
304  ELSE IF( anrm.GT.rmax ) THEN
305  iscale = 1
306  sigma = rmax / anrm
307  END IF
308  IF( iscale.EQ.1 )
309  $ CALL clascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
310 *
311 * Call CHETRD_2STAGE to reduce Hermitian matrix to tridiagonal form.
312 *
313  inde = 1
314  indtau = 1
315  indhous = indtau + n
316  indwrk = indhous + lhtrd
317  llwork = lwork - indwrk + 1
318 *
319  CALL chetrd_2stage( jobz, uplo, n, a, lda, w, rwork( inde ),
320  $ work( indtau ), work( indhous ), lhtrd,
321  $ work( indwrk ), llwork, iinfo )
322 *
323 * For eigenvalues only, call SSTERF. For eigenvectors, first call
324 * CUNGTR to generate the unitary matrix, then call CSTEQR.
325 *
326  IF( .NOT.wantz ) THEN
327  CALL ssterf( n, w, rwork( inde ), info )
328  ELSE
329  CALL cungtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
330  $ llwork, iinfo )
331  indwrk = inde + n
332  CALL csteqr( jobz, n, w, rwork( inde ), a, lda,
333  $ rwork( indwrk ), info )
334  END IF
335 *
336 * If matrix was scaled, then rescale eigenvalues appropriately.
337 *
338  IF( iscale.EQ.1 ) THEN
339  IF( info.EQ.0 ) THEN
340  imax = n
341  ELSE
342  imax = info - 1
343  END IF
344  CALL sscal( imax, one / sigma, w, 1 )
345  END IF
346 *
347 * Set WORK(1) to optimal complex workspace size.
348 *
349  work( 1 ) = lwmin
350 *
351  RETURN
352 *
353 * End of CHEEV_2STAGE
354 *
355  END
subroutine cungtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
CUNGTR
Definition: cungtr.f:125
subroutine csteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
CSTEQR
Definition: csteqr.f:134
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:145
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cheev_2stage(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO)
CHEEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE mat...
Definition: cheev_2stage.f:191
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
subroutine chetrd_2stage(VECT, UPLO, N, A, LDA, D, E, TAU, HOUS2, LHOUS2, WORK, LWORK, INFO)
CHETRD_2STAGE
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88