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