LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ ssbevd_2stage()

subroutine ssbevd_2stage ( character  jobz,
character  uplo,
integer  n,
integer  kd,
real, dimension( ldab, * )  ab,
integer  ldab,
real, dimension( * )  w,
real, dimension( ldz, * )  z,
integer  ldz,
real, dimension( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
integer  liwork,
integer  info 
)

SSBEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

Download SSBEVD_2STAGE + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SSBEVD_2STAGE computes all the eigenvalues and, optionally, eigenvectors of
 a real symmetric band matrix A using the 2stage technique for
 the reduction to tridiagonal. If eigenvectors are desired, it uses
 a divide and conquer algorithm.
Parameters
[in]JOBZ
          JOBZ is CHARACTER*1
          = 'N':  Compute eigenvalues only;
          = 'V':  Compute eigenvalues and eigenvectors.
                  Not available in this release.
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  Upper triangle of A is stored;
          = 'L':  Lower triangle of A is stored.
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.
[in]KD
          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
[in,out]AB
          AB is REAL array, dimension (LDAB, N)
          On entry, the upper or lower triangle of the symmetric band
          matrix A, stored in the first KD+1 rows of the array.  The
          j-th column of A is stored in the j-th column of the array AB
          as follows:
          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).

          On exit, AB is overwritten by values generated during the
          reduction to tridiagonal form.  If UPLO = 'U', the first
          superdiagonal and the diagonal of the tridiagonal matrix T
          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
          the diagonal and first subdiagonal of T are returned in the
          first two rows of AB.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KD + 1.
[out]W
          W is REAL array, dimension (N)
          If INFO = 0, the eigenvalues in ascending order.
[out]Z
          Z is REAL array, dimension (LDZ, N)
          If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
          eigenvectors of the matrix A, with the i-th column of Z
          holding the eigenvector associated with W(i).
          If JOBZ = 'N', then Z is not referenced.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', LDZ >= max(1,N).
[out]WORK
          WORK is REAL array, dimension LWORK
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The length of the array WORK. LWORK >= 1, when N <= 1;
          otherwise
          If JOBZ = 'N' and N > 1, LWORK must be queried.
                                   LWORK = MAX(1, dimension) where
                                   dimension = (2KD+1)*N + KD*NTHREADS + N
                                   where KD is the size of the band.
                                   NTHREADS is the number of threads used when
                                   openMP compilation is enabled, otherwise =1.
          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal sizes of the WORK and IWORK
          arrays, returns these values as the first entries of the WORK
          and IWORK arrays, and no error message related to LWORK or
          LIWORK is issued by XERBLA.
[out]IWORK
          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
[in]LIWORK
          LIWORK is INTEGER
          The dimension of the array IWORK.
          If JOBZ  = 'N' or N <= 1, LIWORK must be at least 1.
          If JOBZ  = 'V' and N > 2, LIWORK must be at least 3 + 5*N.

          If LIWORK = -1, then a workspace query is assumed; the
          routine only calculates the optimal sizes of the WORK and
          IWORK arrays, returns these values as the first entries of
          the WORK and IWORK arrays, and no error message related to
          LWORK or LIWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, the algorithm failed to converge; i
                off-diagonal elements of an intermediate tridiagonal
                form did not converge to zero.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  All details about the 2stage techniques are available in:

  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
  Parallel reduction to condensed forms for symmetric eigenvalue problems
  using aggregated fine-grained and memory-aware kernels. In Proceedings
  of 2011 International Conference for High Performance Computing,
  Networking, Storage and Analysis (SC '11), New York, NY, USA,
  Article 8 , 11 pages.
  http://doi.acm.org/10.1145/2063384.2063394

  A. Haidar, J. Kurzak, P. Luszczek, 2013.
  An improved parallel singular value algorithm and its implementation
  for multicore hardware, In Proceedings of 2013 International Conference
  for High Performance Computing, Networking, Storage and Analysis (SC '13).
  Denver, Colorado, USA, 2013.
  Article 90, 12 pages.
  http://doi.acm.org/10.1145/2503210.2503292

  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
  A novel hybrid CPU-GPU generalized eigensolver for electronic structure
  calculations based on fine-grained memory aware tasks.
  International Journal of High Performance Computing Applications.
  Volume 28 Issue 2, Pages 196-209, May 2014.
  http://hpc.sagepub.com/content/28/2/196

Definition at line 226 of file ssbevd_2stage.f.

228*
229 IMPLICIT NONE
230*
231* -- LAPACK driver routine --
232* -- LAPACK is a software package provided by Univ. of Tennessee, --
233* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
234*
235* .. Scalar Arguments ..
236 CHARACTER JOBZ, UPLO
237 INTEGER INFO, KD, LDAB, LDZ, LIWORK, LWORK, N
238* ..
239* .. Array Arguments ..
240 INTEGER IWORK( * )
241 REAL AB( LDAB, * ), W( * ), WORK( * ), Z( LDZ, * )
242* ..
243*
244* =====================================================================
245*
246* .. Parameters ..
247 REAL ZERO, ONE
248 parameter( zero = 0.0e+0, one = 1.0e+0 )
249* ..
250* .. Local Scalars ..
251 LOGICAL LOWER, LQUERY, WANTZ
252 INTEGER IINFO, INDE, INDWK2, INDWRK, ISCALE, LIWMIN,
253 $ LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS,
254 $ LLWRK2
255 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
256 $ SMLNUM
257* ..
258* .. External Functions ..
259 LOGICAL LSAME
260 INTEGER ILAENV2STAGE
261 REAL SLAMCH, SLANSB, SROUNDUP_LWORK
262 EXTERNAL lsame, slamch, slansb, ilaenv2stage,
264* ..
265* .. External Subroutines ..
266 EXTERNAL sgemm, slacpy, slascl, sscal, sstedc,
268* ..
269* .. Intrinsic Functions ..
270 INTRINSIC sqrt
271* ..
272* .. Executable Statements ..
273*
274* Test the input parameters.
275*
276 wantz = lsame( jobz, 'V' )
277 lower = lsame( uplo, 'L' )
278 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
279*
280 info = 0
281 IF( n.LE.1 ) THEN
282 liwmin = 1
283 lwmin = 1
284 ELSE
285 ib = ilaenv2stage( 2, 'SSYTRD_SB2ST', jobz, n, kd, -1, -1 )
286 lhtrd = ilaenv2stage( 3, 'SSYTRD_SB2ST', jobz, n, kd, ib, -1 )
287 lwtrd = ilaenv2stage( 4, 'SSYTRD_SB2ST', jobz, n, kd, ib, -1 )
288 IF( wantz ) THEN
289 liwmin = 3 + 5*n
290 lwmin = 1 + 5*n + 2*n**2
291 ELSE
292 liwmin = 1
293 lwmin = max( 2*n, n+lhtrd+lwtrd )
294 END IF
295 END IF
296 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
297 info = -1
298 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
299 info = -2
300 ELSE IF( n.LT.0 ) THEN
301 info = -3
302 ELSE IF( kd.LT.0 ) THEN
303 info = -4
304 ELSE IF( ldab.LT.kd+1 ) THEN
305 info = -6
306 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
307 info = -9
308 END IF
309*
310 IF( info.EQ.0 ) THEN
311 work( 1 ) = sroundup_lwork(lwmin)
312 iwork( 1 ) = liwmin
313*
314 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
315 info = -11
316 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
317 info = -13
318 END IF
319 END IF
320*
321 IF( info.NE.0 ) THEN
322 CALL xerbla( 'SSBEVD_2STAGE', -info )
323 RETURN
324 ELSE IF( lquery ) THEN
325 RETURN
326 END IF
327*
328* Quick return if possible
329*
330 IF( n.EQ.0 )
331 $ RETURN
332*
333 IF( n.EQ.1 ) THEN
334 w( 1 ) = ab( 1, 1 )
335 IF( wantz )
336 $ z( 1, 1 ) = one
337 RETURN
338 END IF
339*
340* Get machine constants.
341*
342 safmin = slamch( 'Safe minimum' )
343 eps = slamch( 'Precision' )
344 smlnum = safmin / eps
345 bignum = one / smlnum
346 rmin = sqrt( smlnum )
347 rmax = sqrt( bignum )
348*
349* Scale matrix to allowable range, if necessary.
350*
351 anrm = slansb( 'M', uplo, n, kd, ab, ldab, work )
352 iscale = 0
353 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
354 iscale = 1
355 sigma = rmin / anrm
356 ELSE IF( anrm.GT.rmax ) THEN
357 iscale = 1
358 sigma = rmax / anrm
359 END IF
360 IF( iscale.EQ.1 ) THEN
361 IF( lower ) THEN
362 CALL slascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
363 ELSE
364 CALL slascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
365 END IF
366 END IF
367*
368* Call SSYTRD_SB2ST to reduce band symmetric matrix to tridiagonal form.
369*
370 inde = 1
371 indhous = inde + n
372 indwrk = indhous + lhtrd
373 llwork = lwork - indwrk + 1
374 indwk2 = indwrk + n*n
375 llwrk2 = lwork - indwk2 + 1
376*
377 CALL ssytrd_sb2st( "N", jobz, uplo, n, kd, ab, ldab, w,
378 $ work( inde ), work( indhous ), lhtrd,
379 $ work( indwrk ), llwork, iinfo )
380*
381* For eigenvalues only, call SSTERF. For eigenvectors, call SSTEDC.
382*
383 IF( .NOT.wantz ) THEN
384 CALL ssterf( n, w, work( inde ), info )
385 ELSE
386 CALL sstedc( 'I', n, w, work( inde ), work( indwrk ), n,
387 $ work( indwk2 ), llwrk2, iwork, liwork, info )
388 CALL sgemm( 'N', 'N', n, n, n, one, z, ldz, work( indwrk ), n,
389 $ zero, work( indwk2 ), n )
390 CALL slacpy( 'A', n, n, work( indwk2 ), n, z, ldz )
391 END IF
392*
393* If matrix was scaled, then rescale eigenvalues appropriately.
394*
395 IF( iscale.EQ.1 )
396 $ CALL sscal( n, one / sigma, w, 1 )
397*
398 work( 1 ) = sroundup_lwork(lwmin)
399 iwork( 1 ) = liwmin
400 RETURN
401*
402* End of SSBEVD_2STAGE
403*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine ssytrd_sb2st(stage1, vect, uplo, n, kd, ab, ldab, d, e, hous, lhous, work, lwork, info)
SSYTRD_SB2ST reduces a real symmetric band matrix A to real symmetric tridiagonal form T
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function slansb(norm, uplo, n, k, ab, ldab, work)
SLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansb.f:129
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:143
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
SSTEDC
Definition sstedc.f:182
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
Here is the call graph for this function:
Here is the caller graph for this function: