LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsyevd_2stage.f
Go to the documentation of this file.
1*> \brief <b> DSYEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices</b>
2*
3* @precisions fortran d -> s
4*
5* =========== DOCUMENTATION ===========
6*
7* Online html documentation available at
8* http://www.netlib.org/lapack/explore-html/
9*
10*> \htmlonly
11*> Download DSYEVD_2STAGE + dependencies
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsyevd_2stage.f">
13*> [TGZ]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsyevd_2stage.f">
15*> [ZIP]</a>
16*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsyevd_2stage.f">
17*> [TXT]</a>
18*> \endhtmlonly
19*
20* Definition:
21* ===========
22*
23* SUBROUTINE DSYEVD_2STAGE( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK,
24* IWORK, LIWORK, INFO )
25*
26* IMPLICIT NONE
27*
28* .. Scalar Arguments ..
29* CHARACTER JOBZ, UPLO
30* INTEGER INFO, LDA, LIWORK, LWORK, N
31* ..
32* .. Array Arguments ..
33* INTEGER IWORK( * )
34* DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
35* ..
36*
37*
38*> \par Purpose:
39* =============
40*>
41*> \verbatim
42*>
43*> DSYEVD_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
44*> real symmetric matrix A using the 2stage technique for
45*> the reduction to tridiagonal. If eigenvectors are desired, it uses a
46*> divide and conquer algorithm.
47*>
48*> \endverbatim
49*
50* Arguments:
51* ==========
52*
53*> \param[in] JOBZ
54*> \verbatim
55*> JOBZ is CHARACTER*1
56*> = 'N': Compute eigenvalues only;
57*> = 'V': Compute eigenvalues and eigenvectors.
58*> Not available in this release.
59*> \endverbatim
60*>
61*> \param[in] UPLO
62*> \verbatim
63*> UPLO is CHARACTER*1
64*> = 'U': Upper triangle of A is stored;
65*> = 'L': Lower triangle of A is stored.
66*> \endverbatim
67*>
68*> \param[in] N
69*> \verbatim
70*> N is INTEGER
71*> The order of the matrix A. N >= 0.
72*> \endverbatim
73*>
74*> \param[in,out] A
75*> \verbatim
76*> A is DOUBLE PRECISION array, dimension (LDA, N)
77*> On entry, the symmetric matrix A. If UPLO = 'U', the
78*> leading N-by-N upper triangular part of A contains the
79*> upper triangular part of the matrix A. If UPLO = 'L',
80*> the leading N-by-N lower triangular part of A contains
81*> the lower triangular part of the matrix A.
82*> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
83*> orthonormal eigenvectors of the matrix A.
84*> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
85*> or the upper triangle (if UPLO='U') of A, including the
86*> diagonal, is destroyed.
87*> \endverbatim
88*>
89*> \param[in] LDA
90*> \verbatim
91*> LDA is INTEGER
92*> The leading dimension of the array A. LDA >= max(1,N).
93*> \endverbatim
94*>
95*> \param[out] W
96*> \verbatim
97*> W is DOUBLE PRECISION array, dimension (N)
98*> If INFO = 0, the eigenvalues in ascending order.
99*> \endverbatim
100*>
101*> \param[out] WORK
102*> \verbatim
103*> WORK is DOUBLE PRECISION array,
104*> dimension (LWORK)
105*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
106*> \endverbatim
107*>
108*> \param[in] LWORK
109*> \verbatim
110*> LWORK is INTEGER
111*> The dimension of the array WORK.
112*> If N <= 1, LWORK must be at least 1.
113*> If JOBZ = 'N' and N > 1, LWORK must be queried.
114*> LWORK = MAX(1, dimension) where
115*> dimension = max(stage1,stage2) + (KD+1)*N + 2*N+1
116*> = N*KD + N*max(KD+1,FACTOPTNB)
117*> + max(2*KD*KD, KD*NTHREADS)
118*> + (KD+1)*N + 2*N+1
119*> where KD is the blocking size of the reduction,
120*> FACTOPTNB is the blocking used by the QR or LQ
121*> algorithm, usually FACTOPTNB=128 is a good choice
122*> NTHREADS is the number of threads used when
123*> openMP compilation is enabled, otherwise =1.
124*> If JOBZ = 'V' and N > 1, LWORK must be at least
125*> 1 + 6*N + 2*N**2.
126*>
127*> If LWORK = -1, then a workspace query is assumed; the routine
128*> only calculates the optimal sizes of the WORK and IWORK
129*> arrays, returns these values as the first entries of the WORK
130*> and IWORK arrays, and no error message related to LWORK or
131*> LIWORK is issued by XERBLA.
132*> \endverbatim
133*>
134*> \param[out] IWORK
135*> \verbatim
136*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
137*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
138*> \endverbatim
139*>
140*> \param[in] LIWORK
141*> \verbatim
142*> LIWORK is INTEGER
143*> The dimension of the array IWORK.
144*> If N <= 1, LIWORK must be at least 1.
145*> If JOBZ = 'N' and N > 1, LIWORK must be at least 1.
146*> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
147*>
148*> If LIWORK = -1, then a workspace query is assumed; the
149*> routine only calculates the optimal sizes of the WORK and
150*> IWORK arrays, returns these values as the first entries of
151*> the WORK and IWORK arrays, and no error message related to
152*> LWORK or LIWORK is issued by XERBLA.
153*> \endverbatim
154*>
155*> \param[out] INFO
156*> \verbatim
157*> INFO is INTEGER
158*> = 0: successful exit
159*> < 0: if INFO = -i, the i-th argument had an illegal value
160*> > 0: if INFO = i and JOBZ = 'N', then the algorithm failed
161*> to converge; i off-diagonal elements of an intermediate
162*> tridiagonal form did not converge to zero;
163*> if INFO = i and JOBZ = 'V', then the algorithm failed
164*> to compute an eigenvalue while working on the submatrix
165*> lying in rows and columns INFO/(N+1) through
166*> mod(INFO,N+1).
167*> \endverbatim
168*
169* Authors:
170* ========
171*
172*> \author Univ. of Tennessee
173*> \author Univ. of California Berkeley
174*> \author Univ. of Colorado Denver
175*> \author NAG Ltd.
176*
177*> \ingroup heevd_2stage
178*
179*> \par Contributors:
180* ==================
181*>
182*> Jeff Rutter, Computer Science Division, University of California
183*> at Berkeley, USA \n
184*> Modified by Francoise Tisseur, University of Tennessee \n
185*> Modified description of INFO. Sven, 16 Feb 05. \n
186*> \par Further Details:
187* =====================
188*>
189*> \verbatim
190*>
191*> All details about the 2stage techniques are available in:
192*>
193*> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
194*> Parallel reduction to condensed forms for symmetric eigenvalue problems
195*> using aggregated fine-grained and memory-aware kernels. In Proceedings
196*> of 2011 International Conference for High Performance Computing,
197*> Networking, Storage and Analysis (SC '11), New York, NY, USA,
198*> Article 8 , 11 pages.
199*> http://doi.acm.org/10.1145/2063384.2063394
200*>
201*> A. Haidar, J. Kurzak, P. Luszczek, 2013.
202*> An improved parallel singular value algorithm and its implementation
203*> for multicore hardware, In Proceedings of 2013 International Conference
204*> for High Performance Computing, Networking, Storage and Analysis (SC '13).
205*> Denver, Colorado, USA, 2013.
206*> Article 90, 12 pages.
207*> http://doi.acm.org/10.1145/2503210.2503292
208*>
209*> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
210*> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
211*> calculations based on fine-grained memory aware tasks.
212*> International Journal of High Performance Computing Applications.
213*> Volume 28 Issue 2, Pages 196-209, May 2014.
214*> http://hpc.sagepub.com/content/28/2/196
215*>
216*> \endverbatim
217*
218* =====================================================================
219 SUBROUTINE dsyevd_2stage( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK,
220 $ IWORK, LIWORK, INFO )
221*
222 IMPLICIT NONE
223*
224* -- LAPACK driver routine --
225* -- LAPACK is a software package provided by Univ. of Tennessee, --
226* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
227*
228* .. Scalar Arguments ..
229 CHARACTER JOBZ, UPLO
230 INTEGER INFO, LDA, LIWORK, LWORK, N
231* ..
232* .. Array Arguments ..
233 INTEGER IWORK( * )
234 DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * )
235* ..
236*
237* =====================================================================
238*
239* .. Parameters ..
240 DOUBLE PRECISION ZERO, ONE
241 parameter( zero = 0.0d+0, one = 1.0d+0 )
242* ..
243* .. Local Scalars ..
244*
245 LOGICAL LOWER, LQUERY, WANTZ
246 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
247 $ liwmin, llwork, llwrk2, lwmin,
248 $ lhtrd, lwtrd, kd, ib, indhous
249 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
250 $ smlnum
251* ..
252* .. External Functions ..
253 LOGICAL LSAME
254 INTEGER ILAENV2STAGE
255 DOUBLE PRECISION DLAMCH, DLANSY
256 EXTERNAL lsame, dlamch, dlansy, ilaenv2stage
257* ..
258* .. External Subroutines ..
259 EXTERNAL dlacpy, dlascl, dormtr, dscal, dstedc, dsterf,
261* ..
262* .. Intrinsic Functions ..
263 INTRINSIC max, sqrt
264* ..
265* .. Executable Statements ..
266*
267* Test the input parameters.
268*
269 wantz = lsame( jobz, 'V' )
270 lower = lsame( uplo, 'L' )
271 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
272*
273 info = 0
274 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
275 info = -1
276 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
277 info = -2
278 ELSE IF( n.LT.0 ) THEN
279 info = -3
280 ELSE IF( lda.LT.max( 1, n ) ) THEN
281 info = -5
282 END IF
283*
284 IF( info.EQ.0 ) THEN
285 IF( n.LE.1 ) THEN
286 liwmin = 1
287 lwmin = 1
288 ELSE
289 kd = ilaenv2stage( 1, 'DSYTRD_2STAGE', jobz,
290 $ n, -1, -1, -1 )
291 ib = ilaenv2stage( 2, 'DSYTRD_2STAGE', jobz,
292 $ n, kd, -1, -1 )
293 lhtrd = ilaenv2stage( 3, 'DSYTRD_2STAGE', jobz,
294 $ n, kd, ib, -1 )
295 lwtrd = ilaenv2stage( 4, 'DSYTRD_2STAGE', jobz,
296 $ n, kd, ib, -1 )
297 IF( wantz ) THEN
298 liwmin = 3 + 5*n
299 lwmin = 1 + 6*n + 2*n**2
300 ELSE
301 liwmin = 1
302 lwmin = 2*n + 1 + lhtrd + lwtrd
303 END IF
304 END IF
305 work( 1 ) = lwmin
306 iwork( 1 ) = liwmin
307*
308 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
309 info = -8
310 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
311 info = -10
312 END IF
313 END IF
314*
315 IF( info.NE.0 ) THEN
316 CALL xerbla( 'DSYEVD_2STAGE', -info )
317 RETURN
318 ELSE IF( lquery ) THEN
319 RETURN
320 END IF
321*
322* Quick return if possible
323*
324 IF( n.EQ.0 )
325 $ RETURN
326*
327 IF( n.EQ.1 ) THEN
328 w( 1 ) = a( 1, 1 )
329 IF( wantz )
330 $ a( 1, 1 ) = one
331 RETURN
332 END IF
333*
334* Get machine constants.
335*
336 safmin = dlamch( 'Safe minimum' )
337 eps = dlamch( 'Precision' )
338 smlnum = safmin / eps
339 bignum = one / smlnum
340 rmin = sqrt( smlnum )
341 rmax = sqrt( bignum )
342*
343* Scale matrix to allowable range, if necessary.
344*
345 anrm = dlansy( 'M', uplo, n, a, lda, work )
346 iscale = 0
347 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
348 iscale = 1
349 sigma = rmin / anrm
350 ELSE IF( anrm.GT.rmax ) THEN
351 iscale = 1
352 sigma = rmax / anrm
353 END IF
354 IF( iscale.EQ.1 )
355 $ CALL dlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
356*
357* Call DSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
358*
359 inde = 1
360 indtau = inde + n
361 indhous = indtau + n
362 indwrk = indhous + lhtrd
363 llwork = lwork - indwrk + 1
364 indwk2 = indwrk + n*n
365 llwrk2 = lwork - indwk2 + 1
366*
367 CALL dsytrd_2stage( jobz, uplo, n, a, lda, w, work( inde ),
368 $ work( indtau ), work( indhous ), lhtrd,
369 $ work( indwrk ), llwork, iinfo )
370*
371* For eigenvalues only, call DSTERF. For eigenvectors, first call
372* DSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
373* tridiagonal matrix, then call DORMTR to multiply it by the
374* Householder transformations stored in A.
375*
376 IF( .NOT.wantz ) THEN
377 CALL dsterf( n, w, work( inde ), info )
378 ELSE
379* Not available in this release, and argument checking should not
380* let it getting here
381 RETURN
382 CALL dstedc( 'I', n, w, work( inde ), work( indwrk ), n,
383 $ work( indwk2 ), llwrk2, iwork, liwork, info )
384 CALL dormtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
385 $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
386 CALL dlacpy( 'A', n, n, work( indwrk ), n, a, lda )
387 END IF
388*
389* If matrix was scaled, then rescale eigenvalues appropriately.
390*
391 IF( iscale.EQ.1 )
392 $ CALL dscal( n, one / sigma, w, 1 )
393*
394 work( 1 ) = lwmin
395 iwork( 1 ) = liwmin
396*
397 RETURN
398*
399* End of DSYEVD_2STAGE
400*
401 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dsyevd_2stage(jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)
DSYEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY mat...
subroutine dsytrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
DSYTRD_2STAGE
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
DSTEDC
Definition dstedc.f:182
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
subroutine dormtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
DORMTR
Definition dormtr.f:171