LAPACK  3.8.0
LAPACK: Linear Algebra PACKage
ssyev_2stage.f
Go to the documentation of this file.
1 *> \brief <b> SSYEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices</b>
2 *
3 * @generated from dsyev_2stage.f, fortran d -> s, Sat Nov 5 23:55:51 2016
4 *
5 * =========== DOCUMENTATION ===========
6 *
7 * Online html documentation available at
8 * http://www.netlib.org/lapack/explore-html/
9 *
10 *> \htmlonly
11 *> Download SSYEV_2STAGE + dependencies
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ssyevd_2stage.f">
13 *> [TGZ]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ssyevd_2stage.f">
15 *> [ZIP]</a>
16 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ssyevd_2stage.f">
17 *> [TXT]</a>
18 *> \endhtmlonly
19 *
20 * Definition:
21 * ===========
22 *
23 * SUBROUTINE SSYEV_2STAGE( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK,
24 * 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 A( LDA, * ), W( * ), WORK( * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> SSYEV_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
43 *> real symmetric matrix A using the 2stage technique for
44 *> the reduction to tridiagonal.
45 *> \endverbatim
46 *
47 * Arguments:
48 * ==========
49 *
50 *> \param[in] JOBZ
51 *> \verbatim
52 *> JOBZ is CHARACTER*1
53 *> = 'N': Compute eigenvalues only;
54 *> = 'V': Compute eigenvalues and eigenvectors.
55 *> Not available in this release.
56 *> \endverbatim
57 *>
58 *> \param[in] UPLO
59 *> \verbatim
60 *> UPLO is CHARACTER*1
61 *> = 'U': Upper triangle of A is stored;
62 *> = 'L': Lower triangle of A is stored.
63 *> \endverbatim
64 *>
65 *> \param[in] N
66 *> \verbatim
67 *> N is INTEGER
68 *> The order of the matrix A. N >= 0.
69 *> \endverbatim
70 *>
71 *> \param[in,out] A
72 *> \verbatim
73 *> A is REAL array, dimension (LDA, N)
74 *> On entry, the symmetric matrix A. If UPLO = 'U', the
75 *> leading N-by-N upper triangular part of A contains the
76 *> upper triangular part of the matrix A. If UPLO = 'L',
77 *> the leading N-by-N lower triangular part of A contains
78 *> the lower triangular part of the matrix A.
79 *> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
80 *> orthonormal eigenvectors of the matrix A.
81 *> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
82 *> or the upper triangle (if UPLO='U') of A, including the
83 *> diagonal, is destroyed.
84 *> \endverbatim
85 *>
86 *> \param[in] LDA
87 *> \verbatim
88 *> LDA is INTEGER
89 *> The leading dimension of the array A. LDA >= max(1,N).
90 *> \endverbatim
91 *>
92 *> \param[out] W
93 *> \verbatim
94 *> W is REAL array, dimension (N)
95 *> If INFO = 0, the eigenvalues in ascending order.
96 *> \endverbatim
97 *>
98 *> \param[out] WORK
99 *> \verbatim
100 *> WORK is REAL array, dimension LWORK
101 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
102 *> \endverbatim
103 *>
104 *> \param[in] LWORK
105 *> \verbatim
106 *> LWORK is INTEGER
107 *> The length of the array WORK. LWORK >= 1, when N <= 1;
108 *> otherwise
109 *> If JOBZ = 'N' and N > 1, LWORK must be queried.
110 *> LWORK = MAX(1, dimension) where
111 *> dimension = max(stage1,stage2) + (KD+1)*N + 2*N
112 *> = N*KD + N*max(KD+1,FACTOPTNB)
113 *> + max(2*KD*KD, KD*NTHREADS)
114 *> + (KD+1)*N + 2*N
115 *> where KD is the blocking size of the reduction,
116 *> FACTOPTNB is the blocking used by the QR or LQ
117 *> algorithm, usually FACTOPTNB=128 is a good choice
118 *> NTHREADS is the number of threads used when
119 *> openMP compilation is enabled, otherwise =1.
120 *> If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available
121 *>
122 *> If LWORK = -1, then a workspace query is assumed; the routine
123 *> only calculates the optimal size of the WORK array, returns
124 *> this value as the first entry of the WORK array, and no error
125 *> message related to LWORK is issued by XERBLA.
126 *> \endverbatim
127 *>
128 *> \param[out] INFO
129 *> \verbatim
130 *> INFO is INTEGER
131 *> = 0: successful exit
132 *> < 0: if INFO = -i, the i-th argument had an illegal value
133 *> > 0: if INFO = i, the algorithm failed to converge; i
134 *> off-diagonal elements of an intermediate tridiagonal
135 *> form did not converge to zero.
136 *> \endverbatim
137 *
138 * Authors:
139 * ========
140 *
141 *> \author Univ. of Tennessee
142 *> \author Univ. of California Berkeley
143 *> \author Univ. of Colorado Denver
144 *> \author NAG Ltd.
145 *
146 *> \date November 2017
147 *
148 *> \ingroup realSYeigen
149 *
150 *> \par Further Details:
151 * =====================
152 *>
153 *> \verbatim
154 *>
155 *> All details about the 2stage techniques are available in:
156 *>
157 *> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
158 *> Parallel reduction to condensed forms for symmetric eigenvalue problems
159 *> using aggregated fine-grained and memory-aware kernels. In Proceedings
160 *> of 2011 International Conference for High Performance Computing,
161 *> Networking, Storage and Analysis (SC '11), New York, NY, USA,
162 *> Article 8 , 11 pages.
163 *> http://doi.acm.org/10.1145/2063384.2063394
164 *>
165 *> A. Haidar, J. Kurzak, P. Luszczek, 2013.
166 *> An improved parallel singular value algorithm and its implementation
167 *> for multicore hardware, In Proceedings of 2013 International Conference
168 *> for High Performance Computing, Networking, Storage and Analysis (SC '13).
169 *> Denver, Colorado, USA, 2013.
170 *> Article 90, 12 pages.
171 *> http://doi.acm.org/10.1145/2503210.2503292
172 *>
173 *> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
174 *> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
175 *> calculations based on fine-grained memory aware tasks.
176 *> International Journal of High Performance Computing Applications.
177 *> Volume 28 Issue 2, Pages 196-209, May 2014.
178 *> http://hpc.sagepub.com/content/28/2/196
179 *>
180 *> \endverbatim
181 *
182 * =====================================================================
183  SUBROUTINE ssyev_2stage( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK,
184  $ INFO )
185 *
186  IMPLICIT NONE
187 *
188 * -- LAPACK driver routine (version 3.8.0) --
189 * -- LAPACK is a software package provided by Univ. of Tennessee, --
190 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
191 * November 2017
192 *
193 * .. Scalar Arguments ..
194  CHARACTER JOBZ, UPLO
195  INTEGER INFO, LDA, LWORK, N
196 * ..
197 * .. Array Arguments ..
198  REAL A( lda, * ), W( * ), WORK( * )
199 * ..
200 *
201 * =====================================================================
202 *
203 * .. Parameters ..
204  REAL ZERO, ONE
205  parameter( zero = 0.0e0, one = 1.0e0 )
206 * ..
207 * .. Local Scalars ..
208  LOGICAL LOWER, LQUERY, WANTZ
209  INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
210  $ llwork, lwmin, lhtrd, lwtrd, kd, ib, indhous
211  REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
212  $ smlnum
213 * ..
214 * .. External Functions ..
215  LOGICAL LSAME
216  INTEGER ILAENV2STAGE
217  REAL SLAMCH, SLANSY
218  EXTERNAL lsame, slamch, slansy, ilaenv2stage
219 * ..
220 * .. External Subroutines ..
221  EXTERNAL slascl, sorgtr, sscal, ssteqr, ssterf,
223 * ..
224 * .. Intrinsic Functions ..
225  INTRINSIC max, sqrt
226 * ..
227 * .. Executable Statements ..
228 *
229 * Test the input parameters.
230 *
231  wantz = lsame( jobz, 'V' )
232  lower = lsame( uplo, 'L' )
233  lquery = ( lwork.EQ.-1 )
234 *
235  info = 0
236  IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
237  info = -1
238  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
239  info = -2
240  ELSE IF( n.LT.0 ) THEN
241  info = -3
242  ELSE IF( lda.LT.max( 1, n ) ) THEN
243  info = -5
244  END IF
245 *
246  IF( info.EQ.0 ) THEN
247  kd = ilaenv2stage( 1, 'SSYTRD_2STAGE', jobz, n, -1, -1, -1 )
248  ib = ilaenv2stage( 2, 'SSYTRD_2STAGE', jobz, n, kd, -1, -1 )
249  lhtrd = ilaenv2stage( 3, 'SSYTRD_2STAGE', jobz, n, kd, ib, -1 )
250  lwtrd = ilaenv2stage( 4, 'SSYTRD_2STAGE', jobz, n, kd, ib, -1 )
251  lwmin = 2*n + lhtrd + lwtrd
252  work( 1 ) = lwmin
253 *
254  IF( lwork.LT.lwmin .AND. .NOT.lquery )
255  $ info = -8
256  END IF
257 *
258  IF( info.NE.0 ) THEN
259  CALL xerbla( 'SSYEV_2STAGE ', -info )
260  RETURN
261  ELSE IF( lquery ) THEN
262  RETURN
263  END IF
264 *
265 * Quick return if possible
266 *
267  IF( n.EQ.0 ) THEN
268  RETURN
269  END IF
270 *
271  IF( n.EQ.1 ) THEN
272  w( 1 ) = a( 1, 1 )
273  work( 1 ) = 2
274  IF( wantz )
275  $ a( 1, 1 ) = one
276  RETURN
277  END IF
278 *
279 * Get machine constants.
280 *
281  safmin = slamch( 'Safe minimum' )
282  eps = slamch( 'Precision' )
283  smlnum = safmin / eps
284  bignum = one / smlnum
285  rmin = sqrt( smlnum )
286  rmax = sqrt( bignum )
287 *
288 * Scale matrix to allowable range, if necessary.
289 *
290  anrm = slansy( 'M', uplo, n, a, lda, work )
291  iscale = 0
292  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
293  iscale = 1
294  sigma = rmin / anrm
295  ELSE IF( anrm.GT.rmax ) THEN
296  iscale = 1
297  sigma = rmax / anrm
298  END IF
299  IF( iscale.EQ.1 )
300  $ CALL slascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
301 *
302 * Call SSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
303 *
304  inde = 1
305  indtau = inde + n
306  indhous = indtau + n
307  indwrk = indhous + lhtrd
308  llwork = lwork - indwrk + 1
309 *
310  CALL ssytrd_2stage( jobz, uplo, n, a, lda, w, work( inde ),
311  $ work( indtau ), work( indhous ), lhtrd,
312  $ work( indwrk ), llwork, iinfo )
313 *
314 * For eigenvalues only, call SSTERF. For eigenvectors, first call
315 * SORGTR to generate the orthogonal matrix, then call SSTEQR.
316 *
317  IF( .NOT.wantz ) THEN
318  CALL ssterf( n, w, work( inde ), info )
319  ELSE
320 * Not available in this release, and agrument checking should not
321 * let it getting here
322  RETURN
323  CALL sorgtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
324  $ llwork, iinfo )
325  CALL ssteqr( jobz, n, w, work( inde ), a, lda, work( indtau ),
326  $ info )
327  END IF
328 *
329 * If matrix was scaled, then rescale eigenvalues appropriately.
330 *
331  IF( iscale.EQ.1 ) THEN
332  IF( info.EQ.0 ) THEN
333  imax = n
334  ELSE
335  imax = info - 1
336  END IF
337  CALL sscal( imax, one / sigma, w, 1 )
338  END IF
339 *
340 * Set WORK(1) to optimal workspace size.
341 *
342  work( 1 ) = lwmin
343 *
344  RETURN
345 *
346 * End of SSYEV_2STAGE
347 *
348  END
subroutine ssytrd_2stage(VECT, UPLO, N, A, LDA, D, E, TAU, HOUS2, LHOUS2, WORK, LWORK, INFO)
SSYTRD_2STAGE
subroutine sorgtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
SORGTR
Definition: sorgtr.f:125
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:145
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine ssyev_2stage(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO)
SSYEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY mat...
Definition: ssyev_2stage.f:185
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81