LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsbev_2stage.f
Go to the documentation of this file.
1*> \brief <b> DSBEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER 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 DSBEV_2STAGE + dependencies
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dsbev_2stage.f">
13*> [TGZ]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dsbev_2stage.f">
15*> [ZIP]</a>
16*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dsbev_2stage.f">
17*> [TXT]</a>
18*> \endhtmlonly
19*
20* Definition:
21* ===========
22*
23* SUBROUTINE DSBEV_2STAGE( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
24* WORK, LWORK, INFO )
25*
26* IMPLICIT NONE
27*
28* .. Scalar Arguments ..
29* CHARACTER JOBZ, UPLO
30* INTEGER INFO, KD, LDAB, LDZ, N, LWORK
31* ..
32* .. Array Arguments ..
33* DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> DSBEV_2STAGE computes all the eigenvalues and, optionally, eigenvectors of
43*> a real symmetric band 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] KD
72*> \verbatim
73*> KD is INTEGER
74*> The number of superdiagonals of the matrix A if UPLO = 'U',
75*> or the number of subdiagonals if UPLO = 'L'. KD >= 0.
76*> \endverbatim
77*>
78*> \param[in,out] AB
79*> \verbatim
80*> AB is DOUBLE PRECISION array, dimension (LDAB, N)
81*> On entry, the upper or lower triangle of the symmetric band
82*> matrix A, stored in the first KD+1 rows of the array. The
83*> j-th column of A is stored in the j-th column of the array AB
84*> as follows:
85*> if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
86*> if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
87*>
88*> On exit, AB is overwritten by values generated during the
89*> reduction to tridiagonal form. If UPLO = 'U', the first
90*> superdiagonal and the diagonal of the tridiagonal matrix T
91*> are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
92*> the diagonal and first subdiagonal of T are returned in the
93*> first two rows of AB.
94*> \endverbatim
95*>
96*> \param[in] LDAB
97*> \verbatim
98*> LDAB is INTEGER
99*> The leading dimension of the array AB. LDAB >= KD + 1.
100*> \endverbatim
101*>
102*> \param[out] W
103*> \verbatim
104*> W is DOUBLE PRECISION array, dimension (N)
105*> If INFO = 0, the eigenvalues in ascending order.
106*> \endverbatim
107*>
108*> \param[out] Z
109*> \verbatim
110*> Z is DOUBLE PRECISION array, dimension (LDZ, N)
111*> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
112*> eigenvectors of the matrix A, with the i-th column of Z
113*> holding the eigenvector associated with W(i).
114*> If JOBZ = 'N', then Z is not referenced.
115*> \endverbatim
116*>
117*> \param[in] LDZ
118*> \verbatim
119*> LDZ is INTEGER
120*> The leading dimension of the array Z. LDZ >= 1, and if
121*> JOBZ = 'V', LDZ >= max(1,N).
122*> \endverbatim
123*>
124*> \param[out] WORK
125*> \verbatim
126*> WORK is DOUBLE PRECISION array, dimension LWORK
127*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
128*> \endverbatim
129*>
130*> \param[in] LWORK
131*> \verbatim
132*> LWORK is INTEGER
133*> The length of the array WORK. LWORK >= 1, when N <= 1;
134*> otherwise
135*> If JOBZ = 'N' and N > 1, LWORK must be queried.
136*> LWORK = MAX(1, dimension) where
137*> dimension = (2KD+1)*N + KD*NTHREADS + N
138*> where KD is the size of the band.
139*> NTHREADS is the number of threads used when
140*> openMP compilation is enabled, otherwise =1.
141*> If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available.
142*>
143*> If LWORK = -1, then a workspace query is assumed; the routine
144*> only calculates the optimal size of the WORK array, returns
145*> this value as the first entry of the WORK array, and no error
146*> message related to LWORK is issued by XERBLA.
147*> \endverbatim
148*>
149*> \param[out] INFO
150*> \verbatim
151*> INFO is INTEGER
152*> = 0: successful exit
153*> < 0: if INFO = -i, the i-th argument had an illegal value
154*> > 0: if INFO = i, the algorithm failed to converge; i
155*> off-diagonal elements of an intermediate tridiagonal
156*> form did not converge to zero.
157*> \endverbatim
158*
159* Authors:
160* ========
161*
162*> \author Univ. of Tennessee
163*> \author Univ. of California Berkeley
164*> \author Univ. of Colorado Denver
165*> \author NAG Ltd.
166*
167*> \ingroup hbev_2stage
168*
169*> \par Further Details:
170* =====================
171*>
172*> \verbatim
173*>
174*> All details about the 2stage techniques are available in:
175*>
176*> Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
177*> Parallel reduction to condensed forms for symmetric eigenvalue problems
178*> using aggregated fine-grained and memory-aware kernels. In Proceedings
179*> of 2011 International Conference for High Performance Computing,
180*> Networking, Storage and Analysis (SC '11), New York, NY, USA,
181*> Article 8 , 11 pages.
182*> http://doi.acm.org/10.1145/2063384.2063394
183*>
184*> A. Haidar, J. Kurzak, P. Luszczek, 2013.
185*> An improved parallel singular value algorithm and its implementation
186*> for multicore hardware, In Proceedings of 2013 International Conference
187*> for High Performance Computing, Networking, Storage and Analysis (SC '13).
188*> Denver, Colorado, USA, 2013.
189*> Article 90, 12 pages.
190*> http://doi.acm.org/10.1145/2503210.2503292
191*>
192*> A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
193*> A novel hybrid CPU-GPU generalized eigensolver for electronic structure
194*> calculations based on fine-grained memory aware tasks.
195*> International Journal of High Performance Computing Applications.
196*> Volume 28 Issue 2, Pages 196-209, May 2014.
197*> http://hpc.sagepub.com/content/28/2/196
198*>
199*> \endverbatim
200*
201* =====================================================================
202 SUBROUTINE dsbev_2stage( JOBZ, UPLO, N, KD, AB, LDAB, W, Z, LDZ,
203 $ WORK, LWORK, INFO )
204*
205 IMPLICIT NONE
206*
207* -- LAPACK driver routine --
208* -- LAPACK is a software package provided by Univ. of Tennessee, --
209* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
210*
211* .. Scalar Arguments ..
212 CHARACTER JOBZ, UPLO
213 INTEGER INFO, KD, LDAB, LDZ, N, LWORK
214* ..
215* .. Array Arguments ..
216 DOUBLE PRECISION AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
217* ..
218*
219* =====================================================================
220*
221* .. Parameters ..
222 DOUBLE PRECISION ZERO, ONE
223 parameter( zero = 0.0d0, one = 1.0d0 )
224* ..
225* .. Local Scalars ..
226 LOGICAL LOWER, WANTZ, LQUERY
227 INTEGER IINFO, IMAX, INDE, INDWRK, ISCALE,
228 $ llwork, lwmin, lhtrd, lwtrd, ib, indhous
229 DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
230 $ smlnum
231* ..
232* .. External Functions ..
233 LOGICAL LSAME
234 INTEGER ILAENV2STAGE
235 DOUBLE PRECISION DLAMCH, DLANSB
236 EXTERNAL lsame, dlamch, dlansb, ilaenv2stage
237* ..
238* .. External Subroutines ..
239 EXTERNAL dlascl, dscal, dsteqr, dsterf, xerbla,
241* ..
242* .. Intrinsic Functions ..
243 INTRINSIC sqrt
244* ..
245* .. Executable Statements ..
246*
247* Test the input parameters.
248*
249 wantz = lsame( jobz, 'V' )
250 lower = lsame( uplo, 'L' )
251 lquery = ( lwork.EQ.-1 )
252*
253 info = 0
254 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
255 info = -1
256 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
257 info = -2
258 ELSE IF( n.LT.0 ) THEN
259 info = -3
260 ELSE IF( kd.LT.0 ) THEN
261 info = -4
262 ELSE IF( ldab.LT.kd+1 ) THEN
263 info = -6
264 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
265 info = -9
266 END IF
267*
268 IF( info.EQ.0 ) THEN
269 IF( n.LE.1 ) THEN
270 lwmin = 1
271 work( 1 ) = lwmin
272 ELSE
273 ib = ilaenv2stage( 2, 'DSYTRD_SB2ST', jobz,
274 $ n, kd, -1, -1 )
275 lhtrd = ilaenv2stage( 3, 'DSYTRD_SB2ST', jobz,
276 $ n, kd, ib, -1 )
277 lwtrd = ilaenv2stage( 4, 'DSYTRD_SB2ST', jobz,
278 $ n, kd, ib, -1 )
279 lwmin = n + lhtrd + lwtrd
280 work( 1 ) = lwmin
281 ENDIF
282*
283 IF( lwork.LT.lwmin .AND. .NOT.lquery )
284 $ info = -11
285 END IF
286*
287 IF( info.NE.0 ) THEN
288 CALL xerbla( 'DSBEV_2STAGE ', -info )
289 RETURN
290 ELSE IF( lquery ) THEN
291 RETURN
292 END IF
293*
294* Quick return if possible
295*
296 IF( n.EQ.0 )
297 $ RETURN
298*
299 IF( n.EQ.1 ) THEN
300 IF( lower ) THEN
301 w( 1 ) = ab( 1, 1 )
302 ELSE
303 w( 1 ) = ab( kd+1, 1 )
304 END IF
305 IF( wantz )
306 $ z( 1, 1 ) = one
307 RETURN
308 END IF
309*
310* Get machine constants.
311*
312 safmin = dlamch( 'Safe minimum' )
313 eps = dlamch( 'Precision' )
314 smlnum = safmin / eps
315 bignum = one / smlnum
316 rmin = sqrt( smlnum )
317 rmax = sqrt( bignum )
318*
319* Scale matrix to allowable range, if necessary.
320*
321 anrm = dlansb( 'M', uplo, n, kd, ab, ldab, work )
322 iscale = 0
323 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
324 iscale = 1
325 sigma = rmin / anrm
326 ELSE IF( anrm.GT.rmax ) THEN
327 iscale = 1
328 sigma = rmax / anrm
329 END IF
330 IF( iscale.EQ.1 ) THEN
331 IF( lower ) THEN
332 CALL dlascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
333 ELSE
334 CALL dlascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
335 END IF
336 END IF
337*
338* Call DSYTRD_SB2ST to reduce symmetric band matrix to tridiagonal form.
339*
340 inde = 1
341 indhous = inde + n
342 indwrk = indhous + lhtrd
343 llwork = lwork - indwrk + 1
344*
345 CALL dsytrd_sb2st( "N", jobz, uplo, n, kd, ab, ldab, w,
346 $ work( inde ), work( indhous ), lhtrd,
347 $ work( indwrk ), llwork, iinfo )
348*
349* For eigenvalues only, call DSTERF. For eigenvectors, call SSTEQR.
350*
351 IF( .NOT.wantz ) THEN
352 CALL dsterf( n, w, work( inde ), info )
353 ELSE
354 CALL dsteqr( jobz, n, w, work( inde ), z, ldz, work( indwrk ),
355 $ info )
356 END IF
357*
358* If matrix was scaled, then rescale eigenvalues appropriately.
359*
360 IF( iscale.EQ.1 ) THEN
361 IF( info.EQ.0 ) THEN
362 imax = n
363 ELSE
364 imax = info - 1
365 END IF
366 CALL dscal( imax, one / sigma, w, 1 )
367 END IF
368*
369* Set WORK(1) to optimal workspace size.
370*
371 work( 1 ) = lwmin
372*
373 RETURN
374*
375* End of DSBEV_2STAGE
376*
377 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dsbev_2stage(jobz, uplo, n, kd, ab, ldab, w, z, ldz, work, lwork, info)
DSBEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER m...
subroutine dsytrd_sb2st(stage1, vect, uplo, n, kd, ab, ldab, d, e, hous, lhous, work, lwork, info)
DSYTRD_SB2ST reduces a real symmetric band matrix A to real symmetric tridiagonal form T
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 dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:131
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86