LAPACK  3.10.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 *> \ingroup complexHEeigen
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 cheev_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  REAL RWORK( * ), W( * )
202  COMPLEX A( LDA, * ), WORK( * )
203 * ..
204 *
205 * =====================================================================
206 *
207 * .. Parameters ..
208  REAL ZERO, ONE
209  parameter( zero = 0.0e0, one = 1.0e0 )
210  COMPLEX CONE
211  parameter( cone = ( 1.0e0, 0.0e0 ) )
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  REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
218  $ smlnum
219 * ..
220 * .. External Functions ..
221  LOGICAL LSAME
222  INTEGER ILAENV2STAGE
223  REAL SLAMCH, CLANHE
224  EXTERNAL lsame, slamch, clanhe, ilaenv2stage
225 * ..
226 * .. External Subroutines ..
227  EXTERNAL sscal, ssterf, xerbla, clascl, csteqr,
229 * ..
230 * .. Intrinsic Functions ..
231  INTRINSIC real, 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, 'CHETRD_2STAGE', jobz, n, -1, -1, -1 )
254  ib = ilaenv2stage( 2, 'CHETRD_2STAGE', jobz, n, kd, -1, -1 )
255  lhtrd = ilaenv2stage( 3, 'CHETRD_2STAGE', jobz, n, kd, ib, -1 )
256  lwtrd = ilaenv2stage( 4, 'CHETRD_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( 'CHEEV_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 ) = real( 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 = slamch( 'Safe minimum' )
288  eps = slamch( '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 = clanhe( '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 clascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
307 *
308 * Call CHETRD_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 chetrd_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 SSTERF. For eigenvectors, first call
321 * CUNGTR to generate the unitary matrix, then call CSTEQR.
322 *
323  IF( .NOT.wantz ) THEN
324  CALL ssterf( n, w, rwork( inde ), info )
325  ELSE
326  CALL cungtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
327  $ llwork, iinfo )
328  indwrk = inde + n
329  CALL csteqr( 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 sscal( 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 CHEEV_2STAGE
351 *
352  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:86
subroutine chetrd_2stage(VECT, UPLO, N, A, LDA, D, E, TAU, HOUS2, LHOUS2, WORK, LWORK, INFO)
CHETRD_2STAGE
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 matr...
Definition: cheev_2stage.f:189
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:143
subroutine csteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
CSTEQR
Definition: csteqr.f:132
subroutine cungtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
CUNGTR
Definition: cungtr.f:123
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79