OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
real
LAPACK » Symmetric Matrix » Eigenvalue

Functions

subroutine ssyev (jobz, uplo, n, a, lda, w, work, lwork, info)
  SSYEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
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 matrices
subroutine ssyevd (jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)
  SSYEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine ssyevd_2stage (jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)
  SSYEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine ssyevr (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, isuppz, work, lwork, iwork, liwork, info)
  SSYEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine ssyevr_2stage (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, isuppz, work, lwork, iwork, liwork, info)
  SSYEVR_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine ssyevx (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, iwork, ifail, info)
  SSYEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine ssyevx_2stage (jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, iwork, ifail, info)
  SSYEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
subroutine ssygv (itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, info)
 SSYGV
subroutine ssygv_2stage (itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, info)
 SSYGV_2STAGE
subroutine ssygvd (itype, jobz, uplo, n, a, lda, b, ldb, w, work, lwork, iwork, liwork, info)
 SSYGVD
subroutine ssygvx (itype, jobz, range, uplo, n, a, lda, b, ldb, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, iwork, ifail, info)
 SSYGVX

Detailed Description

This is the group of real eigenvalue driver functions for SY matrices

Function Documentation

◆ ssyev()

subroutine ssyev ( character jobz,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( * ) w,
real, dimension( * ) work,
integer lwork,
integer info )

SSYEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
!>
!> SSYEV computes all eigenvalues and, optionally, eigenvectors of a
!> real symmetric matrix A.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[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,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
!>          orthonormal eigenvectors of the matrix A.
!>          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
!>          or the upper triangle (if UPLO='U') of A, including the
!>          diagonal, is destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The length of the array WORK.  LWORK >= max(1,3*N-1).
!>          For optimal efficiency, LWORK >= (NB+2)*N,
!>          where NB is the blocksize for SSYTRD returned by ILAENV.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK 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.

Definition at line 131 of file ssyev.f.

132*
133* -- LAPACK driver routine --
134* -- LAPACK is a software package provided by Univ. of Tennessee, --
135* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
136*
137* .. Scalar Arguments ..
138 CHARACTER JOBZ, UPLO
139 INTEGER INFO, LDA, LWORK, N
140* ..
141* .. Array Arguments ..
142 REAL A( LDA, * ), W( * ), WORK( * )
143* ..
144*
145* =====================================================================
146*
147* .. Parameters ..
148 REAL ZERO, ONE
149 parameter( zero = 0.0e0, one = 1.0e0 )
150* ..
151* .. Local Scalars ..
152 LOGICAL LOWER, LQUERY, WANTZ
153 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
154 $ LLWORK, LWKOPT, NB
155 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
156 $ SMLNUM
157* ..
158* .. External Functions ..
159 LOGICAL LSAME
160 INTEGER ILAENV
161 REAL SLAMCH, SLANSY
162 EXTERNAL ilaenv, lsame, slamch, slansy
163* ..
164* .. External Subroutines ..
165 EXTERNAL slascl, sorgtr, sscal, ssteqr, ssterf, ssytrd,
166 $ xerbla
167* ..
168* .. Intrinsic Functions ..
169 INTRINSIC max, sqrt
170* ..
171* .. Executable Statements ..
172*
173* Test the input parameters.
174*
175 wantz = lsame( jobz, 'V' )
176 lower = lsame( uplo, 'L' )
177 lquery = ( lwork.EQ.-1 )
178*
179 info = 0
180 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
181 info = -1
182 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
183 info = -2
184 ELSE IF( n.LT.0 ) THEN
185 info = -3
186 ELSE IF( lda.LT.max( 1, n ) ) THEN
187 info = -5
188 END IF
189*
190 IF( info.EQ.0 ) THEN
191 nb = ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 )
192 lwkopt = max( 1, ( nb+2 )*n )
193 work( 1 ) = lwkopt
194*
195 IF( lwork.LT.max( 1, 3*n-1 ) .AND. .NOT.lquery )
196 $ info = -8
197 END IF
198*
199 IF( info.NE.0 ) THEN
200 CALL xerbla( 'SSYEV ', -info )
201 RETURN
202 ELSE IF( lquery ) THEN
203 RETURN
204 END IF
205*
206* Quick return if possible
207*
208 IF( n.EQ.0 ) THEN
209 RETURN
210 END IF
211*
212 IF( n.EQ.1 ) THEN
213 w( 1 ) = a( 1, 1 )
214 work( 1 ) = 2
215 IF( wantz )
216 $ a( 1, 1 ) = one
217 RETURN
218 END IF
219*
220* Get machine constants.
221*
222 safmin = slamch( 'Safe minimum' )
223 eps = slamch( 'Precision' )
224 smlnum = safmin / eps
225 bignum = one / smlnum
226 rmin = sqrt( smlnum )
227 rmax = sqrt( bignum )
228*
229* Scale matrix to allowable range, if necessary.
230*
231 anrm = slansy( 'M', uplo, n, a, lda, work )
232 iscale = 0
233 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
234 iscale = 1
235 sigma = rmin / anrm
236 ELSE IF( anrm.GT.rmax ) THEN
237 iscale = 1
238 sigma = rmax / anrm
239 END IF
240 IF( iscale.EQ.1 )
241 $ CALL slascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
242*
243* Call SSYTRD to reduce symmetric matrix to tridiagonal form.
244*
245 inde = 1
246 indtau = inde + n
247 indwrk = indtau + n
248 llwork = lwork - indwrk + 1
249 CALL ssytrd( uplo, n, a, lda, w, work( inde ), work( indtau ),
250 $ work( indwrk ), llwork, iinfo )
251*
252* For eigenvalues only, call SSTERF. For eigenvectors, first call
253* SORGTR to generate the orthogonal matrix, then call SSTEQR.
254*
255 IF( .NOT.wantz ) THEN
256 CALL ssterf( n, w, work( inde ), info )
257 ELSE
258 CALL sorgtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
259 $ llwork, iinfo )
260 CALL ssteqr( jobz, n, w, work( inde ), a, lda, work( indtau ),
261 $ info )
262 END IF
263*
264* If matrix was scaled, then rescale eigenvalues appropriately.
265*
266 IF( iscale.EQ.1 ) THEN
267 IF( info.EQ.0 ) THEN
268 imax = n
269 ELSE
270 imax = info - 1
271 END IF
272 CALL sscal( imax, one / sigma, w, 1 )
273 END IF
274*
275* Set WORK(1) to optimal workspace size.
276*
277 work( 1 ) = lwkopt
278*
279 RETURN
280*
281* End of SSYEV
282*
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
slascl
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
ssteqr
subroutine ssteqr(compz, n, d, e, z, ldz, work, info)
SSTEQR
Definition ssteqr.f:131
ssterf
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
sorgtr
subroutine sorgtr(uplo, n, a, lda, tau, work, lwork, info)
SORGTR
Definition sorgtr.f:123
slansy
real function slansy(norm, uplo, n, a, lda, work)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition slansy.f:122
ssytrd
subroutine ssytrd(uplo, n, a, lda, d, e, tau, work, lwork, info)
SSYTRD
Definition ssytrd.f:192
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
max
#define max(a, b)
Definition macros.h:21

◆ ssyev_2stage()

subroutine ssyev_2stage ( character jobz,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( * ) w,
real, dimension( * ) work,
integer lwork,
integer info )

SSYEV_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
!>
!> SSYEV_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
!> real symmetric matrix A using the 2stage technique for
!> the reduction to tridiagonal.
!>
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,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
!>          orthonormal eigenvectors of the matrix A.
!>          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
!>          or the upper triangle (if UPLO='U') of A, including the
!>          diagonal, is destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[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 = max(stage1,stage2) + (KD+1)*N + 2*N
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 2*N
!>                                   where KD is the blocking size of the reduction,
!>                                   FACTOPTNB is the blocking used by the QR or LQ
!>                                   algorithm, usually FACTOPTNB=128 is a good choice
!>                                   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 size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK 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 181 of file ssyev_2stage.f.

183*
184 IMPLICIT NONE
185*
186* -- LAPACK driver routine --
187* -- LAPACK is a software package provided by Univ. of Tennessee, --
188* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189*
190* .. Scalar Arguments ..
191 CHARACTER JOBZ, UPLO
192 INTEGER INFO, LDA, LWORK, N
193* ..
194* .. Array Arguments ..
195 REAL A( LDA, * ), W( * ), WORK( * )
196* ..
197*
198* =====================================================================
199*
200* .. Parameters ..
201 REAL ZERO, ONE
202 parameter( zero = 0.0e0, one = 1.0e0 )
203* ..
204* .. Local Scalars ..
205 LOGICAL LOWER, LQUERY, WANTZ
206 INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
207 $ LLWORK, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
208 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
209 $ SMLNUM
210* ..
211* .. External Functions ..
212 LOGICAL LSAME
213 INTEGER ILAENV2STAGE
214 REAL SLAMCH, SLANSY
215 EXTERNAL lsame, slamch, slansy, ilaenv2stage
216* ..
217* .. External Subroutines ..
218 EXTERNAL slascl, sorgtr, sscal, ssteqr, ssterf,
219 $ xerbla, ssytrd_2stage
220* ..
221* .. Intrinsic Functions ..
222 INTRINSIC max, sqrt
223* ..
224* .. Executable Statements ..
225*
226* Test the input parameters.
227*
228 wantz = lsame( jobz, 'V' )
229 lower = lsame( uplo, 'L' )
230 lquery = ( lwork.EQ.-1 )
231*
232 info = 0
233 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
234 info = -1
235 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
236 info = -2
237 ELSE IF( n.LT.0 ) THEN
238 info = -3
239 ELSE IF( lda.LT.max( 1, n ) ) THEN
240 info = -5
241 END IF
242*
243 IF( info.EQ.0 ) THEN
244 kd = ilaenv2stage( 1, 'SSYTRD_2STAGE', jobz, n, -1, -1, -1 )
245 ib = ilaenv2stage( 2, 'SSYTRD_2STAGE', jobz, n, kd, -1, -1 )
246 lhtrd = ilaenv2stage( 3, 'SSYTRD_2STAGE', jobz, n, kd, ib, -1 )
247 lwtrd = ilaenv2stage( 4, 'SSYTRD_2STAGE', jobz, n, kd, ib, -1 )
248 lwmin = 2*n + lhtrd + lwtrd
249 work( 1 ) = lwmin
250*
251 IF( lwork.LT.lwmin .AND. .NOT.lquery )
252 $ info = -8
253 END IF
254*
255 IF( info.NE.0 ) THEN
256 CALL xerbla( 'SSYEV_2STAGE ', -info )
257 RETURN
258 ELSE IF( lquery ) THEN
259 RETURN
260 END IF
261*
262* Quick return if possible
263*
264 IF( n.EQ.0 ) THEN
265 RETURN
266 END IF
267*
268 IF( n.EQ.1 ) THEN
269 w( 1 ) = a( 1, 1 )
270 work( 1 ) = 2
271 IF( wantz )
272 $ a( 1, 1 ) = one
273 RETURN
274 END IF
275*
276* Get machine constants.
277*
278 safmin = slamch( 'Safe minimum' )
279 eps = slamch( 'Precision' )
280 smlnum = safmin / eps
281 bignum = one / smlnum
282 rmin = sqrt( smlnum )
283 rmax = sqrt( bignum )
284*
285* Scale matrix to allowable range, if necessary.
286*
287 anrm = slansy( 'M', uplo, n, a, lda, work )
288 iscale = 0
289 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
290 iscale = 1
291 sigma = rmin / anrm
292 ELSE IF( anrm.GT.rmax ) THEN
293 iscale = 1
294 sigma = rmax / anrm
295 END IF
296 IF( iscale.EQ.1 )
297 $ CALL slascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
298*
299* Call SSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
300*
301 inde = 1
302 indtau = inde + n
303 indhous = indtau + n
304 indwrk = indhous + lhtrd
305 llwork = lwork - indwrk + 1
306*
307 CALL ssytrd_2stage( jobz, uplo, n, a, lda, w, work( inde ),
308 $ work( indtau ), work( indhous ), lhtrd,
309 $ work( indwrk ), llwork, iinfo )
310*
311* For eigenvalues only, call SSTERF. For eigenvectors, first call
312* SORGTR to generate the orthogonal matrix, then call SSTEQR.
313*
314 IF( .NOT.wantz ) THEN
315 CALL ssterf( n, w, work( inde ), info )
316 ELSE
317* Not available in this release, and argument checking should not
318* let it getting here
319 RETURN
320 CALL sorgtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
321 $ llwork, iinfo )
322 CALL ssteqr( jobz, n, w, work( inde ), a, lda, work( indtau ),
323 $ info )
324 END IF
325*
326* If matrix was scaled, then rescale eigenvalues appropriately.
327*
328 IF( iscale.EQ.1 ) THEN
329 IF( info.EQ.0 ) THEN
330 imax = n
331 ELSE
332 imax = info - 1
333 END IF
334 CALL sscal( imax, one / sigma, w, 1 )
335 END IF
336*
337* Set WORK(1) to optimal workspace size.
338*
339 work( 1 ) = lwmin
340*
341 RETURN
342*
343* End of SSYEV_2STAGE
344*
ilaenv2stage
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
Definition ilaenv2stage.f:149
ssytrd_2stage
subroutine ssytrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
SSYTRD_2STAGE
Definition ssytrd_2stage.f:224

◆ ssyevd()

subroutine ssyevd ( character jobz,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( * ) w,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

SSYEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
!>
!> SSYEVD computes all eigenvalues and, optionally, eigenvectors of a
!> real symmetric matrix A. If eigenvectors are desired, it uses a
!> divide and conquer algorithm.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
!> Because of large use of BLAS of level 3, SSYEVD needs N**2 more
!> workspace than SSYEVX.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[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,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
!>          orthonormal eigenvectors of the matrix A.
!>          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
!>          or the upper triangle (if UPLO='U') of A, including the
!>          diagonal, is destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[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 dimension of the array WORK.
!>          If N <= 1,               LWORK must be at least 1.
!>          If JOBZ = 'N' and N > 1, LWORK must be at least 2*N+1.
!>          If JOBZ = 'V' and N > 1, LWORK must be at least
!>                                                1 + 6*N + 2*N**2.
!>
!>          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 N <= 1,                LIWORK must be at least 1.
!>          If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
!>          If JOBZ  = 'V' and N > 1, 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 and JOBZ = 'N', then the algorithm failed
!>                to converge; i off-diagonal elements of an intermediate
!>                tridiagonal form did not converge to zero;
!>                if INFO = i and JOBZ = 'V', then the algorithm failed
!>                to compute an eigenvalue while working on the submatrix
!>                lying in rows and columns INFO/(N+1) through
!>                mod(INFO,N+1).
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee
Modified description of INFO. Sven, 16 Feb 05.

Definition at line 181 of file ssyevd.f.

183*
184* -- LAPACK driver routine --
185* -- LAPACK is a software package provided by Univ. of Tennessee, --
186* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187*
188* .. Scalar Arguments ..
189 CHARACTER JOBZ, UPLO
190 INTEGER INFO, LDA, LIWORK, LWORK, N
191* ..
192* .. Array Arguments ..
193 INTEGER IWORK( * )
194 REAL A( LDA, * ), W( * ), WORK( * )
195* ..
196*
197* =====================================================================
198*
199* .. Parameters ..
200 REAL ZERO, ONE
201 parameter( zero = 0.0e+0, one = 1.0e+0 )
202* ..
203* .. Local Scalars ..
204*
205 LOGICAL LOWER, LQUERY, WANTZ
206 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
207 $ LIOPT, LIWMIN, LLWORK, LLWRK2, LOPT, LWMIN
208 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
209 $ SMLNUM
210* ..
211* .. External Functions ..
212 LOGICAL LSAME
213 INTEGER ILAENV
214 REAL SLAMCH, SLANSY
215 EXTERNAL ilaenv, lsame, slamch, slansy
216* ..
217* .. External Subroutines ..
218 EXTERNAL slacpy, slascl, sormtr, sscal, sstedc, ssterf,
219 $ ssytrd, xerbla
220* ..
221* .. Intrinsic Functions ..
222 INTRINSIC max, sqrt
223* ..
224* .. Executable Statements ..
225*
226* Test the input parameters.
227*
228 wantz = lsame( jobz, 'V' )
229 lower = lsame( uplo, 'L' )
230 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
231*
232 info = 0
233 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
234 info = -1
235 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
236 info = -2
237 ELSE IF( n.LT.0 ) THEN
238 info = -3
239 ELSE IF( lda.LT.max( 1, n ) ) THEN
240 info = -5
241 END IF
242*
243 IF( info.EQ.0 ) THEN
244 IF( n.LE.1 ) THEN
245 liwmin = 1
246 lwmin = 1
247 lopt = lwmin
248 liopt = liwmin
249 ELSE
250 IF( wantz ) THEN
251 liwmin = 3 + 5*n
252 lwmin = 1 + 6*n + 2*n**2
253 ELSE
254 liwmin = 1
255 lwmin = 2*n + 1
256 END IF
257 lopt = max( lwmin, 2*n +
258 $ ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 ) )
259 liopt = liwmin
260 END IF
261 work( 1 ) = lopt
262 iwork( 1 ) = liopt
263*
264 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
265 info = -8
266 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
267 info = -10
268 END IF
269 END IF
270*
271 IF( info.NE.0 ) THEN
272 CALL xerbla( 'SSYEVD', -info )
273 RETURN
274 ELSE IF( lquery ) THEN
275 RETURN
276 END IF
277*
278* Quick return if possible
279*
280 IF( n.EQ.0 )
281 $ RETURN
282*
283 IF( n.EQ.1 ) THEN
284 w( 1 ) = a( 1, 1 )
285 IF( wantz )
286 $ a( 1, 1 ) = one
287 RETURN
288 END IF
289*
290* Get machine constants.
291*
292 safmin = slamch( 'Safe minimum' )
293 eps = slamch( 'Precision' )
294 smlnum = safmin / eps
295 bignum = one / smlnum
296 rmin = sqrt( smlnum )
297 rmax = sqrt( bignum )
298*
299* Scale matrix to allowable range, if necessary.
300*
301 anrm = slansy( 'M', uplo, n, a, lda, work )
302 iscale = 0
303 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
304 iscale = 1
305 sigma = rmin / anrm
306 ELSE IF( anrm.GT.rmax ) THEN
307 iscale = 1
308 sigma = rmax / anrm
309 END IF
310 IF( iscale.EQ.1 )
311 $ CALL slascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
312*
313* Call SSYTRD to reduce symmetric matrix to tridiagonal form.
314*
315 inde = 1
316 indtau = inde + n
317 indwrk = indtau + n
318 llwork = lwork - indwrk + 1
319 indwk2 = indwrk + n*n
320 llwrk2 = lwork - indwk2 + 1
321*
322 CALL ssytrd( uplo, n, a, lda, w, work( inde ), work( indtau ),
323 $ work( indwrk ), llwork, iinfo )
324*
325* For eigenvalues only, call SSTERF. For eigenvectors, first call
326* SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
327* tridiagonal matrix, then call SORMTR to multiply it by the
328* Householder transformations stored in A.
329*
330 IF( .NOT.wantz ) THEN
331 CALL ssterf( n, w, work( inde ), info )
332 ELSE
333 CALL sstedc( 'I', n, w, work( inde ), work( indwrk ), n,
334 $ work( indwk2 ), llwrk2, iwork, liwork, info )
335 CALL sormtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
336 $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
337 CALL slacpy( 'A', n, n, work( indwrk ), n, a, lda )
338 END IF
339*
340* If matrix was scaled, then rescale eigenvalues appropriately.
341*
342 IF( iscale.EQ.1 )
343 $ CALL sscal( n, one / sigma, w, 1 )
344*
345 work( 1 ) = lopt
346 iwork( 1 ) = liopt
347*
348 RETURN
349*
350* End of SSYEVD
351*
slacpy
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
sstedc
subroutine sstedc(compz, n, d, e, z, ldz, work, lwork, iwork, liwork, info)
SSTEDC
Definition sstedc.f:188
sormtr
subroutine sormtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
SORMTR
Definition sormtr.f:172

◆ ssyevd_2stage()

subroutine ssyevd_2stage ( character jobz,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( * ) w,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

SSYEVD_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
!>
!> SSYEVD_2STAGE computes all eigenvalues and, optionally, eigenvectors of a
!> real symmetric matrix A using the 2stage technique for
!> the reduction to tridiagonal. If eigenvectors are desired, it uses a
!> divide and conquer algorithm.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
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,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
!>          orthonormal eigenvectors of the matrix A.
!>          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
!>          or the upper triangle (if UPLO='U') of A, including the
!>          diagonal, is destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[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 dimension of the array WORK.
!>          If N <= 1,               LWORK must be at least 1.
!>          If JOBZ = 'N' and N > 1, LWORK must be queried.
!>                                   LWORK = MAX(1, dimension) where
!>                                   dimension = max(stage1,stage2) + (KD+1)*N + 2*N+1
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 2*N+1
!>                                   where KD is the blocking size of the reduction,
!>                                   FACTOPTNB is the blocking used by the QR or LQ
!>                                   algorithm, usually FACTOPTNB=128 is a good choice
!>                                   NTHREADS is the number of threads used when
!>                                   openMP compilation is enabled, otherwise =1.
!>          If JOBZ = 'V' and N > 1, LWORK must be at least
!>                                                1 + 6*N + 2*N**2.
!>
!>          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 N <= 1,                LIWORK must be at least 1.
!>          If JOBZ  = 'N' and N > 1, LIWORK must be at least 1.
!>          If JOBZ  = 'V' and N > 1, 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 and JOBZ = 'N', then the algorithm failed
!>                to converge; i off-diagonal elements of an intermediate
!>                tridiagonal form did not converge to zero;
!>                if INFO = i and JOBZ = 'V', then the algorithm failed
!>                to compute an eigenvalue while working on the submatrix
!>                lying in rows and columns INFO/(N+1) through
!>                mod(INFO,N+1).
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA
Modified by Francoise Tisseur, University of Tennessee
Modified description of INFO. Sven, 16 Feb 05.
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 225 of file ssyevd_2stage.f.

227*
228 IMPLICIT NONE
229*
230* -- LAPACK driver routine --
231* -- LAPACK is a software package provided by Univ. of Tennessee, --
232* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
233*
234* .. Scalar Arguments ..
235 CHARACTER JOBZ, UPLO
236 INTEGER INFO, LDA, LIWORK, LWORK, N
237* ..
238* .. Array Arguments ..
239 INTEGER IWORK( * )
240 REAL A( LDA, * ), W( * ), WORK( * )
241* ..
242*
243* =====================================================================
244*
245* .. Parameters ..
246 REAL ZERO, ONE
247 parameter( zero = 0.0e+0, one = 1.0e+0 )
248* ..
249* .. Local Scalars ..
250*
251 LOGICAL LOWER, LQUERY, WANTZ
252 INTEGER IINFO, INDE, INDTAU, INDWK2, INDWRK, ISCALE,
253 $ LIWMIN, LLWORK, LLWRK2, LWMIN,
254 $ LHTRD, LWTRD, KD, IB, INDHOUS
255 REAL ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
256 $ SMLNUM
257* ..
258* .. External Functions ..
259 LOGICAL LSAME
260 INTEGER ILAENV2STAGE
261 REAL SLAMCH, SLANSY
262 EXTERNAL lsame, slamch, slansy, ilaenv2stage
263* ..
264* .. External Subroutines ..
265 EXTERNAL slacpy, slascl, sormtr, sscal, sstedc, ssterf,
266 $ ssytrd_2stage, xerbla
267* ..
268* .. Intrinsic Functions ..
269 INTRINSIC max, sqrt
270* ..
271* .. Executable Statements ..
272*
273* Test the input parameters.
274*
275 wantz = lsame( jobz, 'V' )
276 lower = lsame( uplo, 'L' )
277 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
278*
279 info = 0
280 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
281 info = -1
282 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
283 info = -2
284 ELSE IF( n.LT.0 ) THEN
285 info = -3
286 ELSE IF( lda.LT.max( 1, n ) ) THEN
287 info = -5
288 END IF
289*
290 IF( info.EQ.0 ) THEN
291 IF( n.LE.1 ) THEN
292 liwmin = 1
293 lwmin = 1
294 ELSE
295 kd = ilaenv2stage( 1, 'SSYTRD_2STAGE', jobz,
296 $ n, -1, -1, -1 )
297 ib = ilaenv2stage( 2, 'SSYTRD_2STAGE', jobz,
298 $ n, kd, -1, -1 )
299 lhtrd = ilaenv2stage( 3, 'SSYTRD_2STAGE', jobz,
300 $ n, kd, ib, -1 )
301 lwtrd = ilaenv2stage( 4, 'SSYTRD_2STAGE', jobz,
302 $ n, kd, ib, -1 )
303 IF( wantz ) THEN
304 liwmin = 3 + 5*n
305 lwmin = 1 + 6*n + 2*n**2
306 ELSE
307 liwmin = 1
308 lwmin = 2*n + 1 + lhtrd + lwtrd
309 END IF
310 END IF
311 work( 1 ) = lwmin
312 iwork( 1 ) = liwmin
313*
314 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
315 info = -8
316 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
317 info = -10
318 END IF
319 END IF
320*
321 IF( info.NE.0 ) THEN
322 CALL xerbla( 'SSYEVD_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 ) = a( 1, 1 )
335 IF( wantz )
336 $ a( 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 = slansy( 'M', uplo, n, a, lda, 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 )
361 $ CALL slascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
362*
363* Call SSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
364*
365 inde = 1
366 indtau = inde + n
367 indhous = indtau + n
368 indwrk = indhous + lhtrd
369 llwork = lwork - indwrk + 1
370 indwk2 = indwrk + n*n
371 llwrk2 = lwork - indwk2 + 1
372*
373 CALL ssytrd_2stage( jobz, uplo, n, a, lda, w, work( inde ),
374 $ work( indtau ), work( indhous ), lhtrd,
375 $ work( indwrk ), llwork, iinfo )
376*
377* For eigenvalues only, call SSTERF. For eigenvectors, first call
378* SSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
379* tridiagonal matrix, then call SORMTR to multiply it by the
380* Householder transformations stored in A.
381*
382 IF( .NOT.wantz ) THEN
383 CALL ssterf( n, w, work( inde ), info )
384 ELSE
385* Not available in this release, and argument checking should not
386* let it getting here
387 RETURN
388 CALL sstedc( 'I', n, w, work( inde ), work( indwrk ), n,
389 $ work( indwk2 ), llwrk2, iwork, liwork, info )
390 CALL sormtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
391 $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
392 CALL slacpy( 'A', n, n, work( indwrk ), n, a, lda )
393 END IF
394*
395* If matrix was scaled, then rescale eigenvalues appropriately.
396*
397 IF( iscale.EQ.1 )
398 $ CALL sscal( n, one / sigma, w, 1 )
399*
400 work( 1 ) = lwmin
401 iwork( 1 ) = liwmin
402*
403 RETURN
404*
405* End of SSYEVD_2STAGE
406*

◆ ssyevr()

subroutine ssyevr ( character jobz,
character range,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real vl,
real vu,
integer il,
integer iu,
real abstol,
integer m,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
integer, dimension( * ) isuppz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

SSYEVR computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
!>
!> SSYEVR computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric matrix A.  Eigenvalues and eigenvectors can be
!> selected by specifying either a range of values or a range of
!> indices for the desired eigenvalues.
!>
!> SSYEVR first reduces the matrix A to tridiagonal form T with a call
!> to SSYTRD.  Then, whenever possible, SSYEVR calls SSTEMR to compute
!> the eigenspectrum using Relatively Robust Representations.  SSTEMR
!> computes eigenvalues by the dqds algorithm, while orthogonal
!> eigenvectors are computed from various  L D L^T representations
!> (also known as Relatively Robust Representations). Gram-Schmidt
!> orthogonalization is avoided as far as possible. More specifically,
!> the various steps of the algorithm are as follows.
!>
!> For each unreduced block (submatrix) of T,
!>    (a) Compute T - sigma I  = L D L^T, so that L and D
!>        define all the wanted eigenvalues to high relative accuracy.
!>        This means that small relative changes in the entries of D and L
!>        cause only small relative changes in the eigenvalues and
!>        eigenvectors. The standard (unfactored) representation of the
!>        tridiagonal matrix T does not have this property in general.
!>    (b) Compute the eigenvalues to suitable accuracy.
!>        If the eigenvectors are desired, the algorithm attains full
!>        accuracy of the computed eigenvalues only right before
!>        the corresponding vectors have to be computed, see steps c) and d).
!>    (c) For each cluster of close eigenvalues, select a new
!>        shift close to the cluster, find a new factorization, and refine
!>        the shifted eigenvalues to suitable accuracy.
!>    (d) For each eigenvalue with a large enough relative separation compute
!>        the corresponding eigenvector by forming a rank revealing twisted
!>        factorization. Go back to (c) for any clusters that remain.
!>
!> The desired accuracy of the output can be specified by the input
!> parameter ABSTOL.
!>
!> For more details, see SSTEMR's documentation and:
!> - Inderjit S. Dhillon and Beresford N. Parlett: 
!>   Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
!> - Inderjit Dhillon and Beresford Parlett:  SIAM Journal on Matrix Analysis and Applications, Vol. 25,
!>   2004.  Also LAPACK Working Note 154.
!> - Inderjit Dhillon: ,
!>   Computer Science Division Technical Report No. UCB/CSD-97-971,
!>   UC Berkeley, May 1997.
!>
!>
!> Note 1 : SSYEVR calls SSTEMR when the full spectrum is requested
!> on machines which conform to the ieee-754 floating point standard.
!> SSYEVR calls SSTEBZ and SSTEIN on non-ieee machines and
!> when partial spectrum requests are made.
!>
!> Normal execution of SSTEMR may create NaNs and infinities and
!> hence may abort due to a floating point exception in environments
!> which do not handle NaNs and infinities in the ieee standard default
!> manner.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found.
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>          For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and
!>          SSTEIN are called
!>
[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,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>          On exit, the lower triangle (if UPLO='L') or the upper
!>          triangle (if UPLO='U') of A, including the diagonal, is
!>          destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]VL
!>          VL is REAL
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is REAL
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is REAL
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing A to tridiagonal form.
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
!>          If high relative accuracy is important, set ABSTOL to
!>          SLAMCH( 'Safe minimum' ).  Doing so will guarantee that
!>          eigenvalues are computed to high relative accuracy when
!>          possible in future releases.  The current code does not
!>          make any guarantees about high relative accuracy, but
!>          future releases will. See J. Barlow and J. Demmel,
!>          , LAPACK Working Note #7, for a discussion
!>          of which matrices define their eigenvalues to high relative
!>          accuracy.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>          Supplying N columns is always safe.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]ISUPPZ
!>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
!>          The support of the eigenvectors in Z, i.e., the indices
!>          indicating the nonzero elements in Z. The i-th eigenvector
!>          is nonzero only in elements ISUPPZ( 2*i-1 ) through
!>          ISUPPZ( 2*i ). This is an output of SSTEMR (tridiagonal
!>          matrix). The support of the eigenvectors of A is typically
!>          1:N because of the orthogonal transformations applied by SORMTR.
!>          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  LWORK >= max(1,26*N).
!>          For optimal efficiency, LWORK >= (NB+6)*N,
!>          where NB is the max of the blocksize for SSYTRD and SORMTR
!>          returned by ILAENV.
!>
!>          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 LWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.  LIWORK >= max(1,10*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:  Internal error
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Inderjit Dhillon, IBM Almaden, USA
Osni Marques, LBNL/NERSC, USA
Ken Stanley, Computer Science Division, University of California at Berkeley, USA
Jason Riedy, Computer Science Division, University of California at Berkeley, USA

Definition at line 333 of file ssyevr.f.

336*
337* -- LAPACK driver routine --
338* -- LAPACK is a software package provided by Univ. of Tennessee, --
339* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
340*
341* .. Scalar Arguments ..
342 CHARACTER JOBZ, RANGE, UPLO
343 INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, N
344 REAL ABSTOL, VL, VU
345* ..
346* .. Array Arguments ..
347 INTEGER ISUPPZ( * ), IWORK( * )
348 REAL A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
349* ..
350*
351* =====================================================================
352*
353* .. Parameters ..
354 REAL ZERO, ONE, TWO
355 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
356* ..
357* .. Local Scalars ..
358 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
359 $ WANTZ, TRYRAC
360 CHARACTER ORDER
361 INTEGER I, IEEEOK, IINFO, IMAX, INDD, INDDD, INDE,
362 $ INDEE, INDIBL, INDIFL, INDISP, INDIWO, INDTAU,
363 $ INDWK, INDWKN, ISCALE, J, JJ, LIWMIN,
364 $ LLWORK, LLWRKN, LWKOPT, LWMIN, NB, NSPLIT
365 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
366 $ SIGMA, SMLNUM, TMP1, VLL, VUU
367* ..
368* .. External Functions ..
369 LOGICAL LSAME
370 INTEGER ILAENV
371 REAL SLAMCH, SLANSY
372 EXTERNAL lsame, ilaenv, slamch, slansy
373* ..
374* .. External Subroutines ..
375 EXTERNAL scopy, sormtr, sscal, sstebz, sstemr, sstein,
376 $ ssterf, sswap, ssytrd, xerbla
377* ..
378* .. Intrinsic Functions ..
379 INTRINSIC max, min, sqrt
380* ..
381* .. Executable Statements ..
382*
383* Test the input parameters.
384*
385 ieeeok = ilaenv( 10, 'SSYEVR', 'N', 1, 2, 3, 4 )
386*
387 lower = lsame( uplo, 'L' )
388 wantz = lsame( jobz, 'V' )
389 alleig = lsame( range, 'A' )
390 valeig = lsame( range, 'V' )
391 indeig = lsame( range, 'I' )
392*
393 lquery = ( ( lwork.EQ.-1 ) .OR. ( liwork.EQ.-1 ) )
394*
395 lwmin = max( 1, 26*n )
396 liwmin = max( 1, 10*n )
397*
398 info = 0
399 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
400 info = -1
401 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
402 info = -2
403 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
404 info = -3
405 ELSE IF( n.LT.0 ) THEN
406 info = -4
407 ELSE IF( lda.LT.max( 1, n ) ) THEN
408 info = -6
409 ELSE
410 IF( valeig ) THEN
411 IF( n.GT.0 .AND. vu.LE.vl )
412 $ info = -8
413 ELSE IF( indeig ) THEN
414 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
415 info = -9
416 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
417 info = -10
418 END IF
419 END IF
420 END IF
421 IF( info.EQ.0 ) THEN
422 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
423 info = -15
424 END IF
425 END IF
426*
427 IF( info.EQ.0 ) THEN
428 nb = ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 )
429 nb = max( nb, ilaenv( 1, 'SORMTR', uplo, n, -1, -1, -1 ) )
430 lwkopt = max( ( nb+1 )*n, lwmin )
431 work( 1 ) = lwkopt
432 iwork( 1 ) = liwmin
433*
434 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
435 info = -18
436 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
437 info = -20
438 END IF
439 END IF
440*
441 IF( info.NE.0 ) THEN
442 CALL xerbla( 'SSYEVR', -info )
443 RETURN
444 ELSE IF( lquery ) THEN
445 RETURN
446 END IF
447*
448* Quick return if possible
449*
450 m = 0
451 IF( n.EQ.0 ) THEN
452 work( 1 ) = 1
453 RETURN
454 END IF
455*
456 IF( n.EQ.1 ) THEN
457 work( 1 ) = 26
458 IF( alleig .OR. indeig ) THEN
459 m = 1
460 w( 1 ) = a( 1, 1 )
461 ELSE
462 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
463 m = 1
464 w( 1 ) = a( 1, 1 )
465 END IF
466 END IF
467 IF( wantz ) THEN
468 z( 1, 1 ) = one
469 isuppz( 1 ) = 1
470 isuppz( 2 ) = 1
471 END IF
472 RETURN
473 END IF
474*
475* Get machine constants.
476*
477 safmin = slamch( 'Safe minimum' )
478 eps = slamch( 'Precision' )
479 smlnum = safmin / eps
480 bignum = one / smlnum
481 rmin = sqrt( smlnum )
482 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
483*
484* Scale matrix to allowable range, if necessary.
485*
486 iscale = 0
487 abstll = abstol
488 IF (valeig) THEN
489 vll = vl
490 vuu = vu
491 END IF
492 anrm = slansy( 'M', uplo, n, a, lda, work )
493 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
494 iscale = 1
495 sigma = rmin / anrm
496 ELSE IF( anrm.GT.rmax ) THEN
497 iscale = 1
498 sigma = rmax / anrm
499 END IF
500 IF( iscale.EQ.1 ) THEN
501 IF( lower ) THEN
502 DO 10 j = 1, n
503 CALL sscal( n-j+1, sigma, a( j, j ), 1 )
504 10 CONTINUE
505 ELSE
506 DO 20 j = 1, n
507 CALL sscal( j, sigma, a( 1, j ), 1 )
508 20 CONTINUE
509 END IF
510 IF( abstol.GT.0 )
511 $ abstll = abstol*sigma
512 IF( valeig ) THEN
513 vll = vl*sigma
514 vuu = vu*sigma
515 END IF
516 END IF
517
518* Initialize indices into workspaces. Note: The IWORK indices are
519* used only if SSTERF or SSTEMR fail.
520
521* WORK(INDTAU:INDTAU+N-1) stores the scalar factors of the
522* elementary reflectors used in SSYTRD.
523 indtau = 1
524* WORK(INDD:INDD+N-1) stores the tridiagonal's diagonal entries.
525 indd = indtau + n
526* WORK(INDE:INDE+N-1) stores the off-diagonal entries of the
527* tridiagonal matrix from SSYTRD.
528 inde = indd + n
529* WORK(INDDD:INDDD+N-1) is a copy of the diagonal entries over
530* -written by SSTEMR (the SSTERF path copies the diagonal to W).
531 inddd = inde + n
532* WORK(INDEE:INDEE+N-1) is a copy of the off-diagonal entries over
533* -written while computing the eigenvalues in SSTERF and SSTEMR.
534 indee = inddd + n
535* INDWK is the starting offset of the left-over workspace, and
536* LLWORK is the remaining workspace size.
537 indwk = indee + n
538 llwork = lwork - indwk + 1
539
540* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in SSTEBZ and
541* stores the block indices of each of the M<=N eigenvalues.
542 indibl = 1
543* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in SSTEBZ and
544* stores the starting and finishing indices of each block.
545 indisp = indibl + n
546* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
547* that corresponding to eigenvectors that fail to converge in
548* SSTEIN. This information is discarded; if any fail, the driver
549* returns INFO > 0.
550 indifl = indisp + n
551* INDIWO is the offset of the remaining integer workspace.
552 indiwo = indifl + n
553
554*
555* Call SSYTRD to reduce symmetric matrix to tridiagonal form.
556*
557 CALL ssytrd( uplo, n, a, lda, work( indd ), work( inde ),
558 $ work( indtau ), work( indwk ), llwork, iinfo )
559*
560* If all eigenvalues are desired
561* then call SSTERF or SSTEMR and SORMTR.
562*
563 test = .false.
564 IF( indeig ) THEN
565 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
566 test = .true.
567 END IF
568 END IF
569 IF( ( alleig.OR.test ) .AND. ( ieeeok.EQ.1 ) ) THEN
570 IF( .NOT.wantz ) THEN
571 CALL scopy( n, work( indd ), 1, w, 1 )
572 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
573 CALL ssterf( n, w, work( indee ), info )
574 ELSE
575 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
576 CALL scopy( n, work( indd ), 1, work( inddd ), 1 )
577*
578 IF (abstol .LE. two*n*eps) THEN
579 tryrac = .true.
580 ELSE
581 tryrac = .false.
582 END IF
583 CALL sstemr( jobz, 'A', n, work( inddd ), work( indee ),
584 $ vl, vu, il, iu, m, w, z, ldz, n, isuppz,
585 $ tryrac, work( indwk ), lwork, iwork, liwork,
586 $ info )
587*
588*
589*
590* Apply orthogonal matrix used in reduction to tridiagonal
591* form to eigenvectors returned by SSTEMR.
592*
593 IF( wantz .AND. info.EQ.0 ) THEN
594 indwkn = inde
595 llwrkn = lwork - indwkn + 1
596 CALL sormtr( 'L', uplo, 'N', n, m, a, lda,
597 $ work( indtau ), z, ldz, work( indwkn ),
598 $ llwrkn, iinfo )
599 END IF
600 END IF
601*
602*
603 IF( info.EQ.0 ) THEN
604* Everything worked. Skip SSTEBZ/SSTEIN. IWORK(:) are
605* undefined.
606 m = n
607 GO TO 30
608 END IF
609 info = 0
610 END IF
611*
612* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
613* Also call SSTEBZ and SSTEIN if SSTEMR fails.
614*
615 IF( wantz ) THEN
616 order = 'B'
617 ELSE
618 order = 'E'
619 END IF
620
621 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
622 $ work( indd ), work( inde ), m, nsplit, w,
623 $ iwork( indibl ), iwork( indisp ), work( indwk ),
624 $ iwork( indiwo ), info )
625*
626 IF( wantz ) THEN
627 CALL sstein( n, work( indd ), work( inde ), m, w,
628 $ iwork( indibl ), iwork( indisp ), z, ldz,
629 $ work( indwk ), iwork( indiwo ), iwork( indifl ),
630 $ info )
631*
632* Apply orthogonal matrix used in reduction to tridiagonal
633* form to eigenvectors returned by SSTEIN.
634*
635 indwkn = inde
636 llwrkn = lwork - indwkn + 1
637 CALL sormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
638 $ ldz, work( indwkn ), llwrkn, iinfo )
639 END IF
640*
641* If matrix was scaled, then rescale eigenvalues appropriately.
642*
643* Jump here if SSTEMR/SSTEIN succeeded.
644 30 CONTINUE
645 IF( iscale.EQ.1 ) THEN
646 IF( info.EQ.0 ) THEN
647 imax = m
648 ELSE
649 imax = info - 1
650 END IF
651 CALL sscal( imax, one / sigma, w, 1 )
652 END IF
653*
654* If eigenvalues are not in order, then sort them, along with
655* eigenvectors. Note: We do not sort the IFAIL portion of IWORK.
656* It may not be initialized (if SSTEMR/SSTEIN succeeded), and we do
657* not return this detailed information to the user.
658*
659 IF( wantz ) THEN
660 DO 50 j = 1, m - 1
661 i = 0
662 tmp1 = w( j )
663 DO 40 jj = j + 1, m
664 IF( w( jj ).LT.tmp1 ) THEN
665 i = jj
666 tmp1 = w( jj )
667 END IF
668 40 CONTINUE
669*
670 IF( i.NE.0 ) THEN
671 w( i ) = w( j )
672 w( j ) = tmp1
673 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
674 END IF
675 50 CONTINUE
676 END IF
677*
678* Set WORK(1) to optimal workspace size.
679*
680 work( 1 ) = lwkopt
681 iwork( 1 ) = liwmin
682*
683 RETURN
684*
685* End of SSYEVR
686*
sstebz
subroutine sstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
SSTEBZ
Definition sstebz.f:273
sstemr
subroutine sstemr(jobz, range, n, d, e, vl, vu, il, iu, m, w, z, ldz, nzc, isuppz, tryrac, work, lwork, iwork, liwork, info)
SSTEMR
Definition sstemr.f:321
sstein
subroutine sstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
SSTEIN
Definition sstein.f:174
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
min
#define min(a, b)
Definition macros.h:20

◆ ssyevr_2stage()

subroutine ssyevr_2stage ( character jobz,
character range,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real vl,
real vu,
integer il,
integer iu,
real abstol,
integer m,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
integer, dimension( * ) isuppz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

SSYEVR_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
!>
!> SSYEVR_2STAGE computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric matrix A using the 2stage technique for
!> the reduction to tridiagonal.  Eigenvalues and eigenvectors can be
!> selected by specifying either a range of values or a range of
!> indices for the desired eigenvalues.
!>
!> SSYEVR_2STAGE first reduces the matrix A to tridiagonal form T with a call
!> to SSYTRD.  Then, whenever possible, SSYEVR_2STAGE calls SSTEMR to compute
!> the eigenspectrum using Relatively Robust Representations.  SSTEMR
!> computes eigenvalues by the dqds algorithm, while orthogonal
!> eigenvectors are computed from various  L D L^T representations
!> (also known as Relatively Robust Representations). Gram-Schmidt
!> orthogonalization is avoided as far as possible. More specifically,
!> the various steps of the algorithm are as follows.
!>
!> For each unreduced block (submatrix) of T,
!>    (a) Compute T - sigma I  = L D L^T, so that L and D
!>        define all the wanted eigenvalues to high relative accuracy.
!>        This means that small relative changes in the entries of D and L
!>        cause only small relative changes in the eigenvalues and
!>        eigenvectors. The standard (unfactored) representation of the
!>        tridiagonal matrix T does not have this property in general.
!>    (b) Compute the eigenvalues to suitable accuracy.
!>        If the eigenvectors are desired, the algorithm attains full
!>        accuracy of the computed eigenvalues only right before
!>        the corresponding vectors have to be computed, see steps c) and d).
!>    (c) For each cluster of close eigenvalues, select a new
!>        shift close to the cluster, find a new factorization, and refine
!>        the shifted eigenvalues to suitable accuracy.
!>    (d) For each eigenvalue with a large enough relative separation compute
!>        the corresponding eigenvector by forming a rank revealing twisted
!>        factorization. Go back to (c) for any clusters that remain.
!>
!> The desired accuracy of the output can be specified by the input
!> parameter ABSTOL.
!>
!> For more details, see SSTEMR's documentation and:
!> - Inderjit S. Dhillon and Beresford N. Parlett: 
!>   Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
!> - Inderjit Dhillon and Beresford Parlett:  SIAM Journal on Matrix Analysis and Applications, Vol. 25,
!>   2004.  Also LAPACK Working Note 154.
!> - Inderjit Dhillon: ,
!>   Computer Science Division Technical Report No. UCB/CSD-97-971,
!>   UC Berkeley, May 1997.
!>
!>
!> Note 1 : SSYEVR_2STAGE calls SSTEMR when the full spectrum is requested
!> on machines which conform to the ieee-754 floating point standard.
!> SSYEVR_2STAGE calls SSTEBZ and SSTEIN on non-ieee machines and
!> when partial spectrum requests are made.
!>
!> Normal execution of SSTEMR may create NaNs and infinities and
!> hence may abort due to a floating point exception in environments
!> which do not handle NaNs and infinities in the ieee standard default
!> manner.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>                  Not available in this release.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found.
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>          For RANGE = 'V' or 'I' and IU - IL < N - 1, SSTEBZ and
!>          SSTEIN are called
!>
[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,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>          On exit, the lower triangle (if UPLO='L') or the upper
!>          triangle (if UPLO='U') of A, including the diagonal, is
!>          destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]VL
!>          VL is REAL
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is REAL
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is REAL
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing A to tridiagonal form.
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
!>          If high relative accuracy is important, set ABSTOL to
!>          SLAMCH( 'Safe minimum' ).  Doing so will guarantee that
!>          eigenvalues are computed to high relative accuracy when
!>          possible in future releases.  The current code does not
!>          make any guarantees about high relative accuracy, but
!>          future releases will. See J. Barlow and J. Demmel,
!>          , LAPACK Working Note #7, for a discussion
!>          of which matrices define their eigenvalues to high relative
!>          accuracy.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          The first M elements contain the selected eigenvalues in
!>          ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>          Supplying N columns is always safe.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]ISUPPZ
!>          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
!>          The support of the eigenvectors in Z, i.e., the indices
!>          indicating the nonzero elements in Z. The i-th eigenvector
!>          is nonzero only in elements ISUPPZ( 2*i-1 ) through
!>          ISUPPZ( 2*i ). This is an output of SSTEMR (tridiagonal
!>          matrix). The support of the eigenvectors of A is typically 
!>          1:N because of the orthogonal transformations applied by SORMTR.
!>          Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.  
!>          If JOBZ = 'N' and N > 1, LWORK must be queried.
!>                                   LWORK = MAX(1, 26*N, dimension) where
!>                                   dimension = max(stage1,stage2) + (KD+1)*N + 5*N
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 5*N
!>                                   where KD is the blocking size of the reduction,
!>                                   FACTOPTNB is the blocking used by the QR or LQ
!>                                   algorithm, usually FACTOPTNB=128 is a good choice
!>                                   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 size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
!>          On exit, if INFO = 0, IWORK(1) returns the optimal LWORK.
!>
[in]LIWORK
!>          LIWORK is INTEGER
!>          The dimension of the array IWORK.  LIWORK >= max(1,10*N).
!>
!>          If LIWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal size of the IWORK array,
!>          returns this value as the first entry of the IWORK array, and
!>          no error message related to 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:  Internal error
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Inderjit Dhillon, IBM Almaden, USA
Osni Marques, LBNL/NERSC, USA
Ken Stanley, Computer Science Division, University of California at Berkeley, USA
Jason Riedy, Computer Science Division, University of California at Berkeley, USA
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 378 of file ssyevr_2stage.f.

381*
382 IMPLICIT NONE
383*
384* -- LAPACK driver routine --
385* -- LAPACK is a software package provided by Univ. of Tennessee, --
386* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
387*
388* .. Scalar Arguments ..
389 CHARACTER JOBZ, RANGE, UPLO
390 INTEGER IL, INFO, IU, LDA, LDZ, LIWORK, LWORK, M, N
391 REAL ABSTOL, VL, VU
392* ..
393* .. Array Arguments ..
394 INTEGER ISUPPZ( * ), IWORK( * )
395 REAL A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
396* ..
397*
398* =====================================================================
399*
400* .. Parameters ..
401 REAL ZERO, ONE, TWO
402 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
403* ..
404* .. Local Scalars ..
405 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, VALEIG, WANTZ,
406 $ TRYRAC, TEST
407 CHARACTER ORDER
408 INTEGER I, IEEEOK, IINFO, IMAX, INDD, INDDD, INDE,
409 $ INDEE, INDIBL, INDIFL, INDISP, INDIWO, INDTAU,
410 $ INDWK, INDWKN, ISCALE, J, JJ, LIWMIN,
411 $ LLWORK, LLWRKN, LWMIN, NSPLIT,
412 $ LHTRD, LWTRD, KD, IB, INDHOUS
413 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
414 $ SIGMA, SMLNUM, TMP1, VLL, VUU
415* ..
416* .. External Functions ..
417 LOGICAL LSAME
418 INTEGER ILAENV, ILAENV2STAGE
419 REAL SLAMCH, SLANSY
420 EXTERNAL lsame, slamch, slansy, ilaenv, ilaenv2stage
421* ..
422* .. External Subroutines ..
423 EXTERNAL scopy, sormtr, sscal, sstebz, sstemr, sstein,
424 $ ssterf, sswap, ssytrd_2stage, xerbla
425* ..
426* .. Intrinsic Functions ..
427 INTRINSIC max, min, sqrt
428* ..
429* .. Executable Statements ..
430*
431* Test the input parameters.
432*
433 ieeeok = ilaenv( 10, 'SSYEVR', 'N', 1, 2, 3, 4 )
434*
435 lower = lsame( uplo, 'L' )
436 wantz = lsame( jobz, 'V' )
437 alleig = lsame( range, 'A' )
438 valeig = lsame( range, 'V' )
439 indeig = lsame( range, 'I' )
440*
441 lquery = ( ( lwork.EQ.-1 ) .OR. ( liwork.EQ.-1 ) )
442*
443 kd = ilaenv2stage( 1, 'SSYTRD_2STAGE', jobz, n, -1, -1, -1 )
444 ib = ilaenv2stage( 2, 'SSYTRD_2STAGE', jobz, n, kd, -1, -1 )
445 lhtrd = ilaenv2stage( 3, 'SSYTRD_2STAGE', jobz, n, kd, ib, -1 )
446 lwtrd = ilaenv2stage( 4, 'SSYTRD_2STAGE', jobz, n, kd, ib, -1 )
447 lwmin = max( 26*n, 5*n + lhtrd + lwtrd )
448 liwmin = max( 1, 10*n )
449*
450 info = 0
451 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
452 info = -1
453 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
454 info = -2
455 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
456 info = -3
457 ELSE IF( n.LT.0 ) THEN
458 info = -4
459 ELSE IF( lda.LT.max( 1, n ) ) THEN
460 info = -6
461 ELSE
462 IF( valeig ) THEN
463 IF( n.GT.0 .AND. vu.LE.vl )
464 $ info = -8
465 ELSE IF( indeig ) THEN
466 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
467 info = -9
468 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
469 info = -10
470 END IF
471 END IF
472 END IF
473 IF( info.EQ.0 ) THEN
474 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
475 info = -15
476 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
477 info = -18
478 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
479 info = -20
480 END IF
481 END IF
482*
483 IF( info.EQ.0 ) THEN
484* NB = ILAENV( 1, 'SSYTRD', UPLO, N, -1, -1, -1 )
485* NB = MAX( NB, ILAENV( 1, 'SORMTR', UPLO, N, -1, -1, -1 ) )
486* LWKOPT = MAX( ( NB+1 )*N, LWMIN )
487 work( 1 ) = lwmin
488 iwork( 1 ) = liwmin
489 END IF
490*
491 IF( info.NE.0 ) THEN
492 CALL xerbla( 'SSYEVR_2STAGE', -info )
493 RETURN
494 ELSE IF( lquery ) THEN
495 RETURN
496 END IF
497*
498* Quick return if possible
499*
500 m = 0
501 IF( n.EQ.0 ) THEN
502 work( 1 ) = 1
503 RETURN
504 END IF
505*
506 IF( n.EQ.1 ) THEN
507 work( 1 ) = 26
508 IF( alleig .OR. indeig ) THEN
509 m = 1
510 w( 1 ) = a( 1, 1 )
511 ELSE
512 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
513 m = 1
514 w( 1 ) = a( 1, 1 )
515 END IF
516 END IF
517 IF( wantz ) THEN
518 z( 1, 1 ) = one
519 isuppz( 1 ) = 1
520 isuppz( 2 ) = 1
521 END IF
522 RETURN
523 END IF
524*
525* Get machine constants.
526*
527 safmin = slamch( 'Safe minimum' )
528 eps = slamch( 'Precision' )
529 smlnum = safmin / eps
530 bignum = one / smlnum
531 rmin = sqrt( smlnum )
532 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
533*
534* Scale matrix to allowable range, if necessary.
535*
536 iscale = 0
537 abstll = abstol
538 IF (valeig) THEN
539 vll = vl
540 vuu = vu
541 END IF
542 anrm = slansy( 'M', uplo, n, a, lda, work )
543 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
544 iscale = 1
545 sigma = rmin / anrm
546 ELSE IF( anrm.GT.rmax ) THEN
547 iscale = 1
548 sigma = rmax / anrm
549 END IF
550 IF( iscale.EQ.1 ) THEN
551 IF( lower ) THEN
552 DO 10 j = 1, n
553 CALL sscal( n-j+1, sigma, a( j, j ), 1 )
554 10 CONTINUE
555 ELSE
556 DO 20 j = 1, n
557 CALL sscal( j, sigma, a( 1, j ), 1 )
558 20 CONTINUE
559 END IF
560 IF( abstol.GT.0 )
561 $ abstll = abstol*sigma
562 IF( valeig ) THEN
563 vll = vl*sigma
564 vuu = vu*sigma
565 END IF
566 END IF
567
568* Initialize indices into workspaces. Note: The IWORK indices are
569* used only if SSTERF or SSTEMR fail.
570
571* WORK(INDTAU:INDTAU+N-1) stores the scalar factors of the
572* elementary reflectors used in SSYTRD.
573 indtau = 1
574* WORK(INDD:INDD+N-1) stores the tridiagonal's diagonal entries.
575 indd = indtau + n
576* WORK(INDE:INDE+N-1) stores the off-diagonal entries of the
577* tridiagonal matrix from SSYTRD.
578 inde = indd + n
579* WORK(INDDD:INDDD+N-1) is a copy of the diagonal entries over
580* -written by SSTEMR (the SSTERF path copies the diagonal to W).
581 inddd = inde + n
582* WORK(INDEE:INDEE+N-1) is a copy of the off-diagonal entries over
583* -written while computing the eigenvalues in SSTERF and SSTEMR.
584 indee = inddd + n
585* INDHOUS is the starting offset Householder storage of stage 2
586 indhous = indee + n
587* INDWK is the starting offset of the left-over workspace, and
588* LLWORK is the remaining workspace size.
589 indwk = indhous + lhtrd
590 llwork = lwork - indwk + 1
591
592
593* IWORK(INDIBL:INDIBL+M-1) corresponds to IBLOCK in SSTEBZ and
594* stores the block indices of each of the M<=N eigenvalues.
595 indibl = 1
596* IWORK(INDISP:INDISP+NSPLIT-1) corresponds to ISPLIT in SSTEBZ and
597* stores the starting and finishing indices of each block.
598 indisp = indibl + n
599* IWORK(INDIFL:INDIFL+N-1) stores the indices of eigenvectors
600* that corresponding to eigenvectors that fail to converge in
601* SSTEIN. This information is discarded; if any fail, the driver
602* returns INFO > 0.
603 indifl = indisp + n
604* INDIWO is the offset of the remaining integer workspace.
605 indiwo = indifl + n
606
607*
608* Call SSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
609*
610*
611 CALL ssytrd_2stage( jobz, uplo, n, a, lda, work( indd ),
612 $ work( inde ), work( indtau ), work( indhous ),
613 $ lhtrd, work( indwk ), llwork, iinfo )
614*
615* If all eigenvalues are desired
616* then call SSTERF or SSTEMR and SORMTR.
617*
618 test = .false.
619 IF( indeig ) THEN
620 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
621 test = .true.
622 END IF
623 END IF
624 IF( ( alleig.OR.test ) .AND. ( ieeeok.EQ.1 ) ) THEN
625 IF( .NOT.wantz ) THEN
626 CALL scopy( n, work( indd ), 1, w, 1 )
627 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
628 CALL ssterf( n, w, work( indee ), info )
629 ELSE
630 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
631 CALL scopy( n, work( indd ), 1, work( inddd ), 1 )
632*
633 IF (abstol .LE. two*n*eps) THEN
634 tryrac = .true.
635 ELSE
636 tryrac = .false.
637 END IF
638 CALL sstemr( jobz, 'A', n, work( inddd ), work( indee ),
639 $ vl, vu, il, iu, m, w, z, ldz, n, isuppz,
640 $ tryrac, work( indwk ), lwork, iwork, liwork,
641 $ info )
642*
643*
644*
645* Apply orthogonal matrix used in reduction to tridiagonal
646* form to eigenvectors returned by SSTEMR.
647*
648 IF( wantz .AND. info.EQ.0 ) THEN
649 indwkn = inde
650 llwrkn = lwork - indwkn + 1
651 CALL sormtr( 'L', uplo, 'N', n, m, a, lda,
652 $ work( indtau ), z, ldz, work( indwkn ),
653 $ llwrkn, iinfo )
654 END IF
655 END IF
656*
657*
658 IF( info.EQ.0 ) THEN
659* Everything worked. Skip SSTEBZ/SSTEIN. IWORK(:) are
660* undefined.
661 m = n
662 GO TO 30
663 END IF
664 info = 0
665 END IF
666*
667* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
668* Also call SSTEBZ and SSTEIN if SSTEMR fails.
669*
670 IF( wantz ) THEN
671 order = 'B'
672 ELSE
673 order = 'E'
674 END IF
675
676 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
677 $ work( indd ), work( inde ), m, nsplit, w,
678 $ iwork( indibl ), iwork( indisp ), work( indwk ),
679 $ iwork( indiwo ), info )
680*
681 IF( wantz ) THEN
682 CALL sstein( n, work( indd ), work( inde ), m, w,
683 $ iwork( indibl ), iwork( indisp ), z, ldz,
684 $ work( indwk ), iwork( indiwo ), iwork( indifl ),
685 $ info )
686*
687* Apply orthogonal matrix used in reduction to tridiagonal
688* form to eigenvectors returned by SSTEIN.
689*
690 indwkn = inde
691 llwrkn = lwork - indwkn + 1
692 CALL sormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
693 $ ldz, work( indwkn ), llwrkn, iinfo )
694 END IF
695*
696* If matrix was scaled, then rescale eigenvalues appropriately.
697*
698* Jump here if SSTEMR/SSTEIN succeeded.
699 30 CONTINUE
700 IF( iscale.EQ.1 ) THEN
701 IF( info.EQ.0 ) THEN
702 imax = m
703 ELSE
704 imax = info - 1
705 END IF
706 CALL sscal( imax, one / sigma, w, 1 )
707 END IF
708*
709* If eigenvalues are not in order, then sort them, along with
710* eigenvectors. Note: We do not sort the IFAIL portion of IWORK.
711* It may not be initialized (if SSTEMR/SSTEIN succeeded), and we do
712* not return this detailed information to the user.
713*
714 IF( wantz ) THEN
715 DO 50 j = 1, m - 1
716 i = 0
717 tmp1 = w( j )
718 DO 40 jj = j + 1, m
719 IF( w( jj ).LT.tmp1 ) THEN
720 i = jj
721 tmp1 = w( jj )
722 END IF
723 40 CONTINUE
724*
725 IF( i.NE.0 ) THEN
726 w( i ) = w( j )
727 w( j ) = tmp1
728 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
729 END IF
730 50 CONTINUE
731 END IF
732*
733* Set WORK(1) to optimal workspace size.
734*
735 work( 1 ) = lwmin
736 iwork( 1 ) = liwmin
737*
738 RETURN
739*
740* End of SSYEVR_2STAGE
741*

◆ ssyevx()

subroutine ssyevx ( character jobz,
character range,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real vl,
real vu,
integer il,
integer iu,
real abstol,
integer m,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

SSYEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
!>
!> SSYEVX computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric matrix A.  Eigenvalues and eigenvectors can be
!> selected by specifying either a range of values or a range of indices
!> for the desired eigenvalues.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found.
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>
[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,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>          On exit, the lower triangle (if UPLO='L') or the upper
!>          triangle (if UPLO='U') of A, including the diagonal, is
!>          destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]VL
!>          VL is REAL
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is REAL
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is REAL
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing A to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*SLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          On normal exit, the first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If an eigenvector fails to converge, then that column of Z
!>          contains the latest approximation to the eigenvector, and the
!>          index of the eigenvector is returned in IFAIL.
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[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 (MAX(1,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 8*N.
!>          For optimal efficiency, LWORK >= (NB+3)*N,
!>          where NB is the max of the blocksize for SSYTRD and SORMTR
!>          returned by ILAENV.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, then i eigenvectors failed to converge.
!>                Their indices are stored in array IFAIL.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 250 of file ssyevx.f.

253*
254* -- LAPACK driver routine --
255* -- LAPACK is a software package provided by Univ. of Tennessee, --
256* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
257*
258* .. Scalar Arguments ..
259 CHARACTER JOBZ, RANGE, UPLO
260 INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
261 REAL ABSTOL, VL, VU
262* ..
263* .. Array Arguments ..
264 INTEGER IFAIL( * ), IWORK( * )
265 REAL A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
266* ..
267*
268* =====================================================================
269*
270* .. Parameters ..
271 REAL ZERO, ONE
272 parameter( zero = 0.0e+0, one = 1.0e+0 )
273* ..
274* .. Local Scalars ..
275 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
276 $ WANTZ
277 CHARACTER ORDER
278 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
279 $ INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
280 $ ITMP1, J, JJ, LLWORK, LLWRKN, LWKMIN,
281 $ LWKOPT, NB, NSPLIT
282 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
283 $ SIGMA, SMLNUM, TMP1, VLL, VUU
284* ..
285* .. External Functions ..
286 LOGICAL LSAME
287 INTEGER ILAENV
288 REAL SLAMCH, SLANSY
289 EXTERNAL lsame, ilaenv, slamch, slansy
290* ..
291* .. External Subroutines ..
292 EXTERNAL scopy, slacpy, sorgtr, sormtr, sscal, sstebz,
293 $ sstein, ssteqr, ssterf, sswap, ssytrd, xerbla
294* ..
295* .. Intrinsic Functions ..
296 INTRINSIC max, min, sqrt
297* ..
298* .. Executable Statements ..
299*
300* Test the input parameters.
301*
302 lower = lsame( uplo, 'L' )
303 wantz = lsame( jobz, 'V' )
304 alleig = lsame( range, 'A' )
305 valeig = lsame( range, 'V' )
306 indeig = lsame( range, 'I' )
307 lquery = ( lwork.EQ.-1 )
308*
309 info = 0
310 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
311 info = -1
312 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
313 info = -2
314 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
315 info = -3
316 ELSE IF( n.LT.0 ) THEN
317 info = -4
318 ELSE IF( lda.LT.max( 1, n ) ) THEN
319 info = -6
320 ELSE
321 IF( valeig ) THEN
322 IF( n.GT.0 .AND. vu.LE.vl )
323 $ info = -8
324 ELSE IF( indeig ) THEN
325 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
326 info = -9
327 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
328 info = -10
329 END IF
330 END IF
331 END IF
332 IF( info.EQ.0 ) THEN
333 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
334 info = -15
335 END IF
336 END IF
337*
338 IF( info.EQ.0 ) THEN
339 IF( n.LE.1 ) THEN
340 lwkmin = 1
341 work( 1 ) = lwkmin
342 ELSE
343 lwkmin = 8*n
344 nb = ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 )
345 nb = max( nb, ilaenv( 1, 'SORMTR', uplo, n, -1, -1, -1 ) )
346 lwkopt = max( lwkmin, ( nb + 3 )*n )
347 work( 1 ) = lwkopt
348 END IF
349*
350 IF( lwork.LT.lwkmin .AND. .NOT.lquery )
351 $ info = -17
352 END IF
353*
354 IF( info.NE.0 ) THEN
355 CALL xerbla( 'SSYEVX', -info )
356 RETURN
357 ELSE IF( lquery ) THEN
358 RETURN
359 END IF
360*
361* Quick return if possible
362*
363 m = 0
364 IF( n.EQ.0 ) THEN
365 RETURN
366 END IF
367*
368 IF( n.EQ.1 ) THEN
369 IF( alleig .OR. indeig ) THEN
370 m = 1
371 w( 1 ) = a( 1, 1 )
372 ELSE
373 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
374 m = 1
375 w( 1 ) = a( 1, 1 )
376 END IF
377 END IF
378 IF( wantz )
379 $ z( 1, 1 ) = one
380 RETURN
381 END IF
382*
383* Get machine constants.
384*
385 safmin = slamch( 'Safe minimum' )
386 eps = slamch( 'Precision' )
387 smlnum = safmin / eps
388 bignum = one / smlnum
389 rmin = sqrt( smlnum )
390 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
391*
392* Scale matrix to allowable range, if necessary.
393*
394 iscale = 0
395 abstll = abstol
396 IF( valeig ) THEN
397 vll = vl
398 vuu = vu
399 END IF
400 anrm = slansy( 'M', uplo, n, a, lda, work )
401 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
402 iscale = 1
403 sigma = rmin / anrm
404 ELSE IF( anrm.GT.rmax ) THEN
405 iscale = 1
406 sigma = rmax / anrm
407 END IF
408 IF( iscale.EQ.1 ) THEN
409 IF( lower ) THEN
410 DO 10 j = 1, n
411 CALL sscal( n-j+1, sigma, a( j, j ), 1 )
412 10 CONTINUE
413 ELSE
414 DO 20 j = 1, n
415 CALL sscal( j, sigma, a( 1, j ), 1 )
416 20 CONTINUE
417 END IF
418 IF( abstol.GT.0 )
419 $ abstll = abstol*sigma
420 IF( valeig ) THEN
421 vll = vl*sigma
422 vuu = vu*sigma
423 END IF
424 END IF
425*
426* Call SSYTRD to reduce symmetric matrix to tridiagonal form.
427*
428 indtau = 1
429 inde = indtau + n
430 indd = inde + n
431 indwrk = indd + n
432 llwork = lwork - indwrk + 1
433 CALL ssytrd( uplo, n, a, lda, work( indd ), work( inde ),
434 $ work( indtau ), work( indwrk ), llwork, iinfo )
435*
436* If all eigenvalues are desired and ABSTOL is less than or equal to
437* zero, then call SSTERF or SORGTR and SSTEQR. If this fails for
438* some eigenvalue, then try SSTEBZ.
439*
440 test = .false.
441 IF( indeig ) THEN
442 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
443 test = .true.
444 END IF
445 END IF
446 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
447 CALL scopy( n, work( indd ), 1, w, 1 )
448 indee = indwrk + 2*n
449 IF( .NOT.wantz ) THEN
450 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
451 CALL ssterf( n, w, work( indee ), info )
452 ELSE
453 CALL slacpy( 'A', n, n, a, lda, z, ldz )
454 CALL sorgtr( uplo, n, z, ldz, work( indtau ),
455 $ work( indwrk ), llwork, iinfo )
456 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
457 CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
458 $ work( indwrk ), info )
459 IF( info.EQ.0 ) THEN
460 DO 30 i = 1, n
461 ifail( i ) = 0
462 30 CONTINUE
463 END IF
464 END IF
465 IF( info.EQ.0 ) THEN
466 m = n
467 GO TO 40
468 END IF
469 info = 0
470 END IF
471*
472* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
473*
474 IF( wantz ) THEN
475 order = 'B'
476 ELSE
477 order = 'E'
478 END IF
479 indibl = 1
480 indisp = indibl + n
481 indiwo = indisp + n
482 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
483 $ work( indd ), work( inde ), m, nsplit, w,
484 $ iwork( indibl ), iwork( indisp ), work( indwrk ),
485 $ iwork( indiwo ), info )
486*
487 IF( wantz ) THEN
488 CALL sstein( n, work( indd ), work( inde ), m, w,
489 $ iwork( indibl ), iwork( indisp ), z, ldz,
490 $ work( indwrk ), iwork( indiwo ), ifail, info )
491*
492* Apply orthogonal matrix used in reduction to tridiagonal
493* form to eigenvectors returned by SSTEIN.
494*
495 indwkn = inde
496 llwrkn = lwork - indwkn + 1
497 CALL sormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
498 $ ldz, work( indwkn ), llwrkn, iinfo )
499 END IF
500*
501* If matrix was scaled, then rescale eigenvalues appropriately.
502*
503 40 CONTINUE
504 IF( iscale.EQ.1 ) THEN
505 IF( info.EQ.0 ) THEN
506 imax = m
507 ELSE
508 imax = info - 1
509 END IF
510 CALL sscal( imax, one / sigma, w, 1 )
511 END IF
512*
513* If eigenvalues are not in order, then sort them, along with
514* eigenvectors.
515*
516 IF( wantz ) THEN
517 DO 60 j = 1, m - 1
518 i = 0
519 tmp1 = w( j )
520 DO 50 jj = j + 1, m
521 IF( w( jj ).LT.tmp1 ) THEN
522 i = jj
523 tmp1 = w( jj )
524 END IF
525 50 CONTINUE
526*
527 IF( i.NE.0 ) THEN
528 itmp1 = iwork( indibl+i-1 )
529 w( i ) = w( j )
530 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
531 w( j ) = tmp1
532 iwork( indibl+j-1 ) = itmp1
533 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
534 IF( info.NE.0 ) THEN
535 itmp1 = ifail( i )
536 ifail( i ) = ifail( j )
537 ifail( j ) = itmp1
538 END IF
539 END IF
540 60 CONTINUE
541 END IF
542*
543* Set WORK(1) to optimal workspace size.
544*
545 work( 1 ) = lwkopt
546*
547 RETURN
548*
549* End of SSYEVX
550*

◆ ssyevx_2stage()

subroutine ssyevx_2stage ( character jobz,
character range,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real vl,
real vu,
integer il,
integer iu,
real abstol,
integer m,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

SSYEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
!>
!> SSYEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
!> of a real symmetric matrix A using the 2stage technique for
!> the reduction to tridiagonal.  Eigenvalues and eigenvectors can be
!> selected by specifying either a range of values or a range of indices
!> for the desired eigenvalues.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>                  Not available in this release.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found.
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>
[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,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>          On exit, the lower triangle (if UPLO='L') or the upper
!>          triangle (if UPLO='U') of A, including the diagonal, is
!>          destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in]VL
!>          VL is REAL
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is REAL
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is REAL
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing A to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*SLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          On normal exit, the first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If an eigenvector fails to converge, then that column of Z
!>          contains the latest approximation to the eigenvector, and the
!>          index of the eigenvector is returned in IFAIL.
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[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 (MAX(1,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, 8*N, dimension) where
!>                                   dimension = max(stage1,stage2) + (KD+1)*N + 3*N
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 3*N
!>                                   where KD is the blocking size of the reduction,
!>                                   FACTOPTNB is the blocking used by the QR or LQ
!>                                   algorithm, usually FACTOPTNB=128 is a good choice
!>                                   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 size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, then i eigenvectors failed to converge.
!>                Their indices are stored in array IFAIL.
!>
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 297 of file ssyevx_2stage.f.

300*
301 IMPLICIT NONE
302*
303* -- LAPACK driver routine --
304* -- LAPACK is a software package provided by Univ. of Tennessee, --
305* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
306*
307* .. Scalar Arguments ..
308 CHARACTER JOBZ, RANGE, UPLO
309 INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
310 REAL ABSTOL, VL, VU
311* ..
312* .. Array Arguments ..
313 INTEGER IFAIL( * ), IWORK( * )
314 REAL A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * )
315* ..
316*
317* =====================================================================
318*
319* .. Parameters ..
320 REAL ZERO, ONE
321 parameter( zero = 0.0e+0, one = 1.0e+0 )
322* ..
323* .. Local Scalars ..
324 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
325 $ WANTZ
326 CHARACTER ORDER
327 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
328 $ INDISP, INDIWO, INDTAU, INDWKN, INDWRK, ISCALE,
329 $ ITMP1, J, JJ, LLWORK, LLWRKN,
330 $ NSPLIT, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
331 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
332 $ SIGMA, SMLNUM, TMP1, VLL, VUU
333* ..
334* .. External Functions ..
335 LOGICAL LSAME
336 INTEGER ILAENV2STAGE
337 REAL SLAMCH, SLANSY
338 EXTERNAL lsame, slamch, slansy, ilaenv2stage
339* ..
340* .. External Subroutines ..
341 EXTERNAL scopy, slacpy, sorgtr, sormtr, sscal, sstebz,
342 $ sstein, ssteqr, ssterf, sswap, xerbla,
343 $ ssytrd_2stage
344* ..
345* .. Intrinsic Functions ..
346 INTRINSIC max, min, sqrt
347* ..
348* .. Executable Statements ..
349*
350* Test the input parameters.
351*
352 lower = lsame( uplo, 'L' )
353 wantz = lsame( jobz, 'V' )
354 alleig = lsame( range, 'A' )
355 valeig = lsame( range, 'V' )
356 indeig = lsame( range, 'I' )
357 lquery = ( lwork.EQ.-1 )
358*
359 info = 0
360 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
361 info = -1
362 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
363 info = -2
364 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
365 info = -3
366 ELSE IF( n.LT.0 ) THEN
367 info = -4
368 ELSE IF( lda.LT.max( 1, n ) ) THEN
369 info = -6
370 ELSE
371 IF( valeig ) THEN
372 IF( n.GT.0 .AND. vu.LE.vl )
373 $ info = -8
374 ELSE IF( indeig ) THEN
375 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
376 info = -9
377 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
378 info = -10
379 END IF
380 END IF
381 END IF
382 IF( info.EQ.0 ) THEN
383 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
384 info = -15
385 END IF
386 END IF
387*
388 IF( info.EQ.0 ) THEN
389 IF( n.LE.1 ) THEN
390 lwmin = 1
391 work( 1 ) = lwmin
392 ELSE
393 kd = ilaenv2stage( 1, 'SSYTRD_2STAGE', jobz,
394 $ n, -1, -1, -1 )
395 ib = ilaenv2stage( 2, 'SSYTRD_2STAGE', jobz,
396 $ n, kd, -1, -1 )
397 lhtrd = ilaenv2stage( 3, 'SSYTRD_2STAGE', jobz,
398 $ n, kd, ib, -1 )
399 lwtrd = ilaenv2stage( 4, 'SSYTRD_2STAGE', jobz,
400 $ n, kd, ib, -1 )
401 lwmin = max( 8*n, 3*n + lhtrd + lwtrd )
402 work( 1 ) = lwmin
403 END IF
404*
405 IF( lwork.LT.lwmin .AND. .NOT.lquery )
406 $ info = -17
407 END IF
408*
409 IF( info.NE.0 ) THEN
410 CALL xerbla( 'SSYEVX_2STAGE', -info )
411 RETURN
412 ELSE IF( lquery ) THEN
413 RETURN
414 END IF
415*
416* Quick return if possible
417*
418 m = 0
419 IF( n.EQ.0 ) THEN
420 RETURN
421 END IF
422*
423 IF( n.EQ.1 ) THEN
424 IF( alleig .OR. indeig ) THEN
425 m = 1
426 w( 1 ) = a( 1, 1 )
427 ELSE
428 IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
429 m = 1
430 w( 1 ) = a( 1, 1 )
431 END IF
432 END IF
433 IF( wantz )
434 $ z( 1, 1 ) = one
435 RETURN
436 END IF
437*
438* Get machine constants.
439*
440 safmin = slamch( 'Safe minimum' )
441 eps = slamch( 'Precision' )
442 smlnum = safmin / eps
443 bignum = one / smlnum
444 rmin = sqrt( smlnum )
445 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
446*
447* Scale matrix to allowable range, if necessary.
448*
449 iscale = 0
450 abstll = abstol
451 IF( valeig ) THEN
452 vll = vl
453 vuu = vu
454 END IF
455 anrm = slansy( 'M', uplo, n, a, lda, work )
456 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
457 iscale = 1
458 sigma = rmin / anrm
459 ELSE IF( anrm.GT.rmax ) THEN
460 iscale = 1
461 sigma = rmax / anrm
462 END IF
463 IF( iscale.EQ.1 ) THEN
464 IF( lower ) THEN
465 DO 10 j = 1, n
466 CALL sscal( n-j+1, sigma, a( j, j ), 1 )
467 10 CONTINUE
468 ELSE
469 DO 20 j = 1, n
470 CALL sscal( j, sigma, a( 1, j ), 1 )
471 20 CONTINUE
472 END IF
473 IF( abstol.GT.0 )
474 $ abstll = abstol*sigma
475 IF( valeig ) THEN
476 vll = vl*sigma
477 vuu = vu*sigma
478 END IF
479 END IF
480*
481* Call SSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
482*
483 indtau = 1
484 inde = indtau + n
485 indd = inde + n
486 indhous = indd + n
487 indwrk = indhous + lhtrd
488 llwork = lwork - indwrk + 1
489*
490 CALL ssytrd_2stage( jobz, uplo, n, a, lda, work( indd ),
491 $ work( inde ), work( indtau ), work( indhous ),
492 $ lhtrd, work( indwrk ), llwork, iinfo )
493*
494* If all eigenvalues are desired and ABSTOL is less than or equal to
495* zero, then call SSTERF or SORGTR and SSTEQR. If this fails for
496* some eigenvalue, then try SSTEBZ.
497*
498 test = .false.
499 IF( indeig ) THEN
500 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
501 test = .true.
502 END IF
503 END IF
504 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
505 CALL scopy( n, work( indd ), 1, w, 1 )
506 indee = indwrk + 2*n
507 IF( .NOT.wantz ) THEN
508 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
509 CALL ssterf( n, w, work( indee ), info )
510 ELSE
511 CALL slacpy( 'A', n, n, a, lda, z, ldz )
512 CALL sorgtr( uplo, n, z, ldz, work( indtau ),
513 $ work( indwrk ), llwork, iinfo )
514 CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
515 CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
516 $ work( indwrk ), info )
517 IF( info.EQ.0 ) THEN
518 DO 30 i = 1, n
519 ifail( i ) = 0
520 30 CONTINUE
521 END IF
522 END IF
523 IF( info.EQ.0 ) THEN
524 m = n
525 GO TO 40
526 END IF
527 info = 0
528 END IF
529*
530* Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
531*
532 IF( wantz ) THEN
533 order = 'B'
534 ELSE
535 order = 'E'
536 END IF
537 indibl = 1
538 indisp = indibl + n
539 indiwo = indisp + n
540 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
541 $ work( indd ), work( inde ), m, nsplit, w,
542 $ iwork( indibl ), iwork( indisp ), work( indwrk ),
543 $ iwork( indiwo ), info )
544*
545 IF( wantz ) THEN
546 CALL sstein( n, work( indd ), work( inde ), m, w,
547 $ iwork( indibl ), iwork( indisp ), z, ldz,
548 $ work( indwrk ), iwork( indiwo ), ifail, info )
549*
550* Apply orthogonal matrix used in reduction to tridiagonal
551* form to eigenvectors returned by SSTEIN.
552*
553 indwkn = inde
554 llwrkn = lwork - indwkn + 1
555 CALL sormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
556 $ ldz, work( indwkn ), llwrkn, iinfo )
557 END IF
558*
559* If matrix was scaled, then rescale eigenvalues appropriately.
560*
561 40 CONTINUE
562 IF( iscale.EQ.1 ) THEN
563 IF( info.EQ.0 ) THEN
564 imax = m
565 ELSE
566 imax = info - 1
567 END IF
568 CALL sscal( imax, one / sigma, w, 1 )
569 END IF
570*
571* If eigenvalues are not in order, then sort them, along with
572* eigenvectors.
573*
574 IF( wantz ) THEN
575 DO 60 j = 1, m - 1
576 i = 0
577 tmp1 = w( j )
578 DO 50 jj = j + 1, m
579 IF( w( jj ).LT.tmp1 ) THEN
580 i = jj
581 tmp1 = w( jj )
582 END IF
583 50 CONTINUE
584*
585 IF( i.NE.0 ) THEN
586 itmp1 = iwork( indibl+i-1 )
587 w( i ) = w( j )
588 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
589 w( j ) = tmp1
590 iwork( indibl+j-1 ) = itmp1
591 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
592 IF( info.NE.0 ) THEN
593 itmp1 = ifail( i )
594 ifail( i ) = ifail( j )
595 ifail( j ) = itmp1
596 END IF
597 END IF
598 60 CONTINUE
599 END IF
600*
601* Set WORK(1) to optimal workspace size.
602*
603 work( 1 ) = lwmin
604*
605 RETURN
606*
607* End of SSYEVX_2STAGE
608*

◆ ssygv()

subroutine ssygv ( integer itype,
character jobz,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) w,
real, dimension( * ) work,
integer lwork,
integer info )

SSYGV

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

Purpose:
!>
!> SSYGV computes all the eigenvalues, and optionally, the eigenvectors
!> of a real generalized symmetric-definite eigenproblem, of the form
!> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
!> Here A and B are assumed to be symmetric and B is also
!> positive definite.
!>
Parameters
[in]ITYPE
!>          ITYPE is INTEGER
!>          Specifies the problem type to be solved:
!>          = 1:  A*x = (lambda)*B*x
!>          = 2:  A*B*x = (lambda)*x
!>          = 3:  B*A*x = (lambda)*x
!>
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangles of A and B are stored;
!>          = 'L':  Lower triangles of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>
!>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
!>          matrix Z of eigenvectors.  The eigenvectors are normalized
!>          as follows:
!>          if ITYPE = 1 or 2, Z**T*B*Z = I;
!>          if ITYPE = 3, Z**T*inv(B)*Z = I.
!>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
!>          or the lower triangle (if UPLO='L') of A, including the
!>          diagonal, is destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the symmetric positive definite matrix B.
!>          If UPLO = 'U', the leading N-by-N upper triangular part of B
!>          contains the upper triangular part of the matrix B.
!>          If UPLO = 'L', the leading N-by-N lower triangular part of B
!>          contains the lower triangular part of the matrix B.
!>
!>          On exit, if INFO <= N, the part of B containing the matrix is
!>          overwritten by the triangular factor U or L from the Cholesky
!>          factorization B = U**T*U or B = L*L**T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The length of the array WORK.  LWORK >= max(1,3*N-1).
!>          For optimal efficiency, LWORK >= (NB+2)*N,
!>          where NB is the blocksize for SSYTRD returned by ILAENV.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK 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:  SPOTRF or SSYEV returned an error code:
!>             <= N:  if INFO = i, SSYEV failed to converge;
!>                    i off-diagonal elements of an intermediate
!>                    tridiagonal form did not converge to zero;
!>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
!>                    minor of order i of B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 173 of file ssygv.f.

175*
176* -- LAPACK driver routine --
177* -- LAPACK is a software package provided by Univ. of Tennessee, --
178* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
179*
180* .. Scalar Arguments ..
181 CHARACTER JOBZ, UPLO
182 INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
183* ..
184* .. Array Arguments ..
185 REAL A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
186* ..
187*
188* =====================================================================
189*
190* .. Parameters ..
191 REAL ONE
192 parameter( one = 1.0e+0 )
193* ..
194* .. Local Scalars ..
195 LOGICAL LQUERY, UPPER, WANTZ
196 CHARACTER TRANS
197 INTEGER LWKMIN, LWKOPT, NB, NEIG
198* ..
199* .. External Functions ..
200 LOGICAL LSAME
201 INTEGER ILAENV
202 EXTERNAL ilaenv, lsame
203* ..
204* .. External Subroutines ..
205 EXTERNAL spotrf, ssyev, ssygst, strmm, strsm, xerbla
206* ..
207* .. Intrinsic Functions ..
208 INTRINSIC max
209* ..
210* .. Executable Statements ..
211*
212* Test the input parameters.
213*
214 wantz = lsame( jobz, 'V' )
215 upper = lsame( uplo, 'U' )
216 lquery = ( lwork.EQ.-1 )
217*
218 info = 0
219 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
220 info = -1
221 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
222 info = -2
223 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
224 info = -3
225 ELSE IF( n.LT.0 ) THEN
226 info = -4
227 ELSE IF( lda.LT.max( 1, n ) ) THEN
228 info = -6
229 ELSE IF( ldb.LT.max( 1, n ) ) THEN
230 info = -8
231 END IF
232*
233 IF( info.EQ.0 ) THEN
234 lwkmin = max( 1, 3*n - 1 )
235 nb = ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 )
236 lwkopt = max( lwkmin, ( nb + 2 )*n )
237 work( 1 ) = lwkopt
238*
239 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
240 info = -11
241 END IF
242 END IF
243*
244 IF( info.NE.0 ) THEN
245 CALL xerbla( 'SSYGV ', -info )
246 RETURN
247 ELSE IF( lquery ) THEN
248 RETURN
249 END IF
250*
251* Quick return if possible
252*
253 IF( n.EQ.0 )
254 $ RETURN
255*
256* Form a Cholesky factorization of B.
257*
258 CALL spotrf( uplo, n, b, ldb, info )
259 IF( info.NE.0 ) THEN
260 info = n + info
261 RETURN
262 END IF
263*
264* Transform problem to standard eigenvalue problem and solve.
265*
266 CALL ssygst( itype, uplo, n, a, lda, b, ldb, info )
267 CALL ssyev( jobz, uplo, n, a, lda, w, work, lwork, info )
268*
269 IF( wantz ) THEN
270*
271* Backtransform eigenvectors to the original problem.
272*
273 neig = n
274 IF( info.GT.0 )
275 $ neig = info - 1
276 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
277*
278* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
279* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
280*
281 IF( upper ) THEN
282 trans = 'N'
283 ELSE
284 trans = 'T'
285 END IF
286*
287 CALL strsm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
288 $ b, ldb, a, lda )
289*
290 ELSE IF( itype.EQ.3 ) THEN
291*
292* For B*A*x=(lambda)*x;
293* backtransform eigenvectors: x = L*y or U**T*y
294*
295 IF( upper ) THEN
296 trans = 'T'
297 ELSE
298 trans = 'N'
299 END IF
300*
301 CALL strmm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
302 $ b, ldb, a, lda )
303 END IF
304 END IF
305*
306 work( 1 ) = lwkopt
307 RETURN
308*
309* End of SSYGV
310*
spotrf
subroutine spotrf(uplo, n, a, lda, info)
SPOTRF
Definition spotrf.f:107
ssygst
subroutine ssygst(itype, uplo, n, a, lda, b, ldb, info)
SSYGST
Definition ssygst.f:127
ssyev
subroutine ssyev(jobz, uplo, n, a, lda, w, work, lwork, info)
SSYEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
Definition ssyev.f:132
strsm
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181
strmm
subroutine strmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRMM
Definition strmm.f:177

◆ ssygv_2stage()

subroutine ssygv_2stage ( integer itype,
character jobz,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) w,
real, dimension( * ) work,
integer lwork,
integer info )

SSYGV_2STAGE

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

Purpose:
!>
!> SSYGV_2STAGE computes all the eigenvalues, and optionally, the eigenvectors
!> of a real generalized symmetric-definite eigenproblem, of the form
!> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.
!> Here A and B are assumed to be symmetric and B is also
!> positive definite.
!> This routine use the 2stage technique for the reduction to tridiagonal
!> which showed higher performance on recent architecture and for large
!> sizes N>2000.
!>
Parameters
[in]ITYPE
!>          ITYPE is INTEGER
!>          Specifies the problem type to be solved:
!>          = 1:  A*x = (lambda)*B*x
!>          = 2:  A*B*x = (lambda)*x
!>          = 3:  B*A*x = (lambda)*x
!>
[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 triangles of A and B are stored;
!>          = 'L':  Lower triangles of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>
!>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
!>          matrix Z of eigenvectors.  The eigenvectors are normalized
!>          as follows:
!>          if ITYPE = 1 or 2, Z**T*B*Z = I;
!>          if ITYPE = 3, Z**T*inv(B)*Z = I.
!>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
!>          or the lower triangle (if UPLO='L') of A, including the
!>          diagonal, is destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the symmetric positive definite matrix B.
!>          If UPLO = 'U', the leading N-by-N upper triangular part of B
!>          contains the upper triangular part of the matrix B.
!>          If UPLO = 'L', the leading N-by-N lower triangular part of B
!>          contains the lower triangular part of the matrix B.
!>
!>          On exit, if INFO <= N, the part of B containing the matrix is
!>          overwritten by the triangular factor U or L from the Cholesky
!>          factorization B = U**T*U or B = L*L**T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,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 = max(stage1,stage2) + (KD+1)*N + 2*N
!>                                             = N*KD + N*max(KD+1,FACTOPTNB) 
!>                                               + max(2*KD*KD, KD*NTHREADS) 
!>                                               + (KD+1)*N + 2*N
!>                                   where KD is the blocking size of the reduction,
!>                                   FACTOPTNB is the blocking used by the QR or LQ
!>                                   algorithm, usually FACTOPTNB=128 is a good choice
!>                                   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 size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK 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:  SPOTRF or SSYEV returned an error code:
!>             <= N:  if INFO = i, SSYEV failed to converge;
!>                    i off-diagonal elements of an intermediate
!>                    tridiagonal form did not converge to zero;
!>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
!>                    minor of order i of B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
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 224 of file ssygv_2stage.f.

226*
227 IMPLICIT NONE
228*
229* -- LAPACK driver routine --
230* -- LAPACK is a software package provided by Univ. of Tennessee, --
231* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
232*
233* .. Scalar Arguments ..
234 CHARACTER JOBZ, UPLO
235 INTEGER INFO, ITYPE, LDA, LDB, LWORK, N
236* ..
237* .. Array Arguments ..
238 REAL A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
239* ..
240*
241* =====================================================================
242*
243* .. Parameters ..
244 REAL ONE
245 parameter( one = 1.0e+0 )
246* ..
247* .. Local Scalars ..
248 LOGICAL LQUERY, UPPER, WANTZ
249 CHARACTER TRANS
250 INTEGER NEIG, LWMIN, LHTRD, LWTRD, KD, IB
251* ..
252* .. External Functions ..
253 LOGICAL LSAME
254 INTEGER ILAENV2STAGE
255 EXTERNAL lsame, ilaenv2stage
256* ..
257* .. External Subroutines ..
258 EXTERNAL spotrf, ssygst, strmm, strsm, xerbla,
259 $ ssyev_2stage
260* ..
261* .. Intrinsic Functions ..
262 INTRINSIC max
263* ..
264* .. Executable Statements ..
265*
266* Test the input parameters.
267*
268 wantz = lsame( jobz, 'V' )
269 upper = lsame( uplo, 'U' )
270 lquery = ( lwork.EQ.-1 )
271*
272 info = 0
273 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
274 info = -1
275 ELSE IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
276 info = -2
277 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
278 info = -3
279 ELSE IF( n.LT.0 ) THEN
280 info = -4
281 ELSE IF( lda.LT.max( 1, n ) ) THEN
282 info = -6
283 ELSE IF( ldb.LT.max( 1, n ) ) THEN
284 info = -8
285 END IF
286*
287 IF( info.EQ.0 ) THEN
288 kd = ilaenv2stage( 1, 'SSYTRD_2STAGE', jobz, n, -1, -1, -1 )
289 ib = ilaenv2stage( 2, 'SSYTRD_2STAGE', jobz, n, kd, -1, -1 )
290 lhtrd = ilaenv2stage( 3, 'SSYTRD_2STAGE', jobz, n, kd, ib, -1 )
291 lwtrd = ilaenv2stage( 4, 'SSYTRD_2STAGE', jobz, n, kd, ib, -1 )
292 lwmin = 2*n + lhtrd + lwtrd
293 work( 1 ) = lwmin
294*
295 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
296 info = -11
297 END IF
298 END IF
299*
300 IF( info.NE.0 ) THEN
301 CALL xerbla( 'SSYGV_2STAGE ', -info )
302 RETURN
303 ELSE IF( lquery ) THEN
304 RETURN
305 END IF
306*
307* Quick return if possible
308*
309 IF( n.EQ.0 )
310 $ RETURN
311*
312* Form a Cholesky factorization of B.
313*
314 CALL spotrf( uplo, n, b, ldb, info )
315 IF( info.NE.0 ) THEN
316 info = n + info
317 RETURN
318 END IF
319*
320* Transform problem to standard eigenvalue problem and solve.
321*
322 CALL ssygst( itype, uplo, n, a, lda, b, ldb, info )
323 CALL ssyev_2stage( jobz, uplo, n, a, lda, w, work, lwork, info )
324*
325 IF( wantz ) THEN
326*
327* Backtransform eigenvectors to the original problem.
328*
329 neig = n
330 IF( info.GT.0 )
331 $ neig = info - 1
332 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
333*
334* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
335* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
336*
337 IF( upper ) THEN
338 trans = 'N'
339 ELSE
340 trans = 'T'
341 END IF
342*
343 CALL strsm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
344 $ b, ldb, a, lda )
345*
346 ELSE IF( itype.EQ.3 ) THEN
347*
348* For B*A*x=(lambda)*x;
349* backtransform eigenvectors: x = L*y or U**T*y
350*
351 IF( upper ) THEN
352 trans = 'T'
353 ELSE
354 trans = 'N'
355 END IF
356*
357 CALL strmm( 'Left', uplo, trans, 'Non-unit', n, neig, one,
358 $ b, ldb, a, lda )
359 END IF
360 END IF
361*
362 work( 1 ) = lwmin
363 RETURN
364*
365* End of SSYGV_2STAGE
366*
ssyev_2stage
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 matr...
Definition ssyev_2stage.f:183

◆ ssygvd()

subroutine ssygvd ( integer itype,
character jobz,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( * ) w,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

SSYGVD

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

Purpose:
!>
!> SSYGVD computes all the eigenvalues, and optionally, the eigenvectors
!> of a real generalized symmetric-definite eigenproblem, of the form
!> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A and
!> B are assumed to be symmetric and B is also positive definite.
!> If eigenvectors are desired, it uses a divide and conquer algorithm.
!>
!> The divide and conquer algorithm makes very mild assumptions about
!> floating point arithmetic. It will work on machines with a guard
!> digit in add/subtract, or on those binary machines without guard
!> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
!> Cray-2. It could conceivably fail on hexadecimal or decimal machines
!> without guard digits, but we know of none.
!>
Parameters
[in]ITYPE
!>          ITYPE is INTEGER
!>          Specifies the problem type to be solved:
!>          = 1:  A*x = (lambda)*B*x
!>          = 2:  A*B*x = (lambda)*x
!>          = 3:  B*A*x = (lambda)*x
!>
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangles of A and B are stored;
!>          = 'L':  Lower triangles of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrices A and B.  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>
!>          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
!>          matrix Z of eigenvectors.  The eigenvectors are normalized
!>          as follows:
!>          if ITYPE = 1 or 2, Z**T*B*Z = I;
!>          if ITYPE = 3, Z**T*inv(B)*Z = I.
!>          If JOBZ = 'N', then on exit the upper triangle (if UPLO='U')
!>          or the lower triangle (if UPLO='L') of A, including the
!>          diagonal, is destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the symmetric matrix B.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of B contains the
!>          upper triangular part of the matrix B.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of B contains
!>          the lower triangular part of the matrix B.
!>
!>          On exit, if INFO <= N, the part of B containing the matrix is
!>          overwritten by the triangular factor U or L from the Cholesky
!>          factorization B = U**T*U or B = L*L**T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the eigenvalues in ascending order.
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          If N <= 1,               LWORK >= 1.
!>          If JOBZ = 'N' and N > 1, LWORK >= 2*N+1.
!>          If JOBZ = 'V' and N > 1, LWORK >= 1 + 6*N + 2*N**2.
!>
!>          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 N <= 1,                LIWORK >= 1.
!>          If JOBZ  = 'N' and N > 1, LIWORK >= 1.
!>          If JOBZ  = 'V' and N > 1, LIWORK >= 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:  SPOTRF or SSYEVD returned an error code:
!>             <= N:  if INFO = i and JOBZ = 'N', then the algorithm
!>                    failed to converge; i off-diagonal elements of an
!>                    intermediate tridiagonal form did not converge to
!>                    zero;
!>                    if INFO = i and JOBZ = 'V', then the algorithm
!>                    failed to compute an eigenvalue while working on
!>                    the submatrix lying in rows and columns INFO/(N+1)
!>                    through mod(INFO,N+1);
!>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
!>                    minor of order i of B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
!>
!>  Modified so that no backsubstitution is performed if SSYEVD fails to
!>  converge (NEIG in old code could be greater than N causing out of
!>  bounds reference to A - reported by Ralf Meyer).  Also corrected the
!>  description of INFO and the test on ITYPE. Sven, 16 Feb 05.
!>
Contributors:
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Definition at line 225 of file ssygvd.f.

227*
228* -- LAPACK driver routine --
229* -- LAPACK is a software package provided by Univ. of Tennessee, --
230* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231*
232* .. Scalar Arguments ..
233 CHARACTER JOBZ, UPLO
234 INTEGER INFO, ITYPE, LDA, LDB, LIWORK, LWORK, N
235* ..
236* .. Array Arguments ..
237 INTEGER IWORK( * )
238 REAL A( LDA, * ), B( LDB, * ), W( * ), WORK( * )
239* ..
240*
241* =====================================================================
242*
243* .. Parameters ..
244 REAL ONE
245 parameter( one = 1.0e+0 )
246* ..
247* .. Local Scalars ..
248 LOGICAL LQUERY, UPPER, WANTZ
249 CHARACTER TRANS
250 INTEGER LIOPT, LIWMIN, LOPT, LWMIN
251* ..
252* .. External Functions ..
253 LOGICAL LSAME
254 EXTERNAL lsame
255* ..
256* .. External Subroutines ..
257 EXTERNAL spotrf, ssyevd, ssygst, strmm, strsm, xerbla
258* ..
259* .. Intrinsic Functions ..
260 INTRINSIC max, real
261* ..
262* .. Executable Statements ..
263*
264* Test the input parameters.
265*
266 wantz = lsame( jobz, 'V' )
267 upper = lsame( uplo, 'U' )
268 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
269*
270 info = 0
271 IF( n.LE.1 ) THEN
272 liwmin = 1
273 lwmin = 1
274 ELSE IF( wantz ) THEN
275 liwmin = 3 + 5*n
276 lwmin = 1 + 6*n + 2*n**2
277 ELSE
278 liwmin = 1
279 lwmin = 2*n + 1
280 END IF
281 lopt = lwmin
282 liopt = liwmin
283 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
284 info = -1
285 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
286 info = -2
287 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
288 info = -3
289 ELSE IF( n.LT.0 ) THEN
290 info = -4
291 ELSE IF( lda.LT.max( 1, n ) ) THEN
292 info = -6
293 ELSE IF( ldb.LT.max( 1, n ) ) THEN
294 info = -8
295 END IF
296*
297 IF( info.EQ.0 ) THEN
298 work( 1 ) = lopt
299 iwork( 1 ) = liopt
300*
301 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
302 info = -11
303 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
304 info = -13
305 END IF
306 END IF
307*
308 IF( info.NE.0 ) THEN
309 CALL xerbla( 'SSYGVD', -info )
310 RETURN
311 ELSE IF( lquery ) THEN
312 RETURN
313 END IF
314*
315* Quick return if possible
316*
317 IF( n.EQ.0 )
318 $ RETURN
319*
320* Form a Cholesky factorization of B.
321*
322 CALL spotrf( uplo, n, b, ldb, info )
323 IF( info.NE.0 ) THEN
324 info = n + info
325 RETURN
326 END IF
327*
328* Transform problem to standard eigenvalue problem and solve.
329*
330 CALL ssygst( itype, uplo, n, a, lda, b, ldb, info )
331 CALL ssyevd( jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork,
332 $ info )
333 lopt = max( real( lopt ), real( work( 1 ) ) )
334 liopt = max( real( liopt ), real( iwork( 1 ) ) )
335*
336 IF( wantz .AND. info.EQ.0 ) THEN
337*
338* Backtransform eigenvectors to the original problem.
339*
340 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
341*
342* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
343* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
344*
345 IF( upper ) THEN
346 trans = 'N'
347 ELSE
348 trans = 'T'
349 END IF
350*
351 CALL strsm( 'Left', uplo, trans, 'Non-unit', n, n, one,
352 $ b, ldb, a, lda )
353*
354 ELSE IF( itype.EQ.3 ) THEN
355*
356* For B*A*x=(lambda)*x;
357* backtransform eigenvectors: x = L*y or U**T*y
358*
359 IF( upper ) THEN
360 trans = 'T'
361 ELSE
362 trans = 'N'
363 END IF
364*
365 CALL strmm( 'Left', uplo, trans, 'Non-unit', n, n, one,
366 $ b, ldb, a, lda )
367 END IF
368 END IF
369*
370 work( 1 ) = lopt
371 iwork( 1 ) = liopt
372*
373 RETURN
374*
375* End of SSYGVD
376*
ssyevd
subroutine ssyevd(jobz, uplo, n, a, lda, w, work, lwork, iwork, liwork, info)
SSYEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
Definition ssyevd.f:183

◆ ssygvx()

subroutine ssygvx ( integer itype,
character jobz,
character range,
character uplo,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real vl,
real vu,
integer il,
integer iu,
real abstol,
integer m,
real, dimension( * ) w,
real, dimension( ldz, * ) z,
integer ldz,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

SSYGVX

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

Purpose:
!>
!> SSYGVX computes selected eigenvalues, and optionally, eigenvectors
!> of a real generalized symmetric-definite eigenproblem, of the form
!> A*x=(lambda)*B*x,  A*Bx=(lambda)*x,  or B*A*x=(lambda)*x.  Here A
!> and B are assumed to be symmetric and B is also positive definite.
!> Eigenvalues and eigenvectors can be selected by specifying either a
!> range of values or a range of indices for the desired eigenvalues.
!>
Parameters
[in]ITYPE
!>          ITYPE is INTEGER
!>          Specifies the problem type to be solved:
!>          = 1:  A*x = (lambda)*B*x
!>          = 2:  A*B*x = (lambda)*x
!>          = 3:  B*A*x = (lambda)*x
!>
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found.
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found.
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A and B are stored;
!>          = 'L':  Lower triangle of A and B are stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix pencil (A,B).  N >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA, N)
!>          On entry, the symmetric matrix A.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of A contains the
!>          upper triangular part of the matrix A.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of A contains
!>          the lower triangular part of the matrix A.
!>
!>          On exit, the lower triangle (if UPLO='L') or the upper
!>          triangle (if UPLO='U') of A, including the diagonal, is
!>          destroyed.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!>
[in,out]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, the symmetric matrix B.  If UPLO = 'U', the
!>          leading N-by-N upper triangular part of B contains the
!>          upper triangular part of the matrix B.  If UPLO = 'L',
!>          the leading N-by-N lower triangular part of B contains
!>          the lower triangular part of the matrix B.
!>
!>          On exit, if INFO <= N, the part of B containing the matrix is
!>          overwritten by the triangular factor U or L from the Cholesky
!>          factorization B = U**T*U or B = L*L**T.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B.  LDB >= max(1,N).
!>
[in]VL
!>          VL is REAL
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is REAL
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is REAL
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing C to tridiagonal form, where C is the symmetric
!>          matrix of the standard symmetric problem to which the
!>          generalized problem is transformed.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*DLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*SLAMCH('S').
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          On normal exit, the first M elements contain the selected
!>          eigenvalues in ascending order.
!>
[out]Z
!>          Z is REAL array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'N', then Z is not referenced.
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          The eigenvectors are normalized as follows:
!>          if ITYPE = 1 or 2, Z**T*B*Z = I;
!>          if ITYPE = 3, Z**T*inv(B)*Z = I.
!>
!>          If an eigenvector fails to converge, then that column of Z
!>          contains the latest approximation to the eigenvector, and the
!>          index of the eigenvector is returned in IFAIL.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[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 (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The length of the array WORK.  LWORK >= max(1,8*N).
!>          For optimal efficiency, LWORK >= (NB+3)*N,
!>          where NB is the blocksize for SSYTRD returned by ILAENV.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  SPOTRF or SSYEVX returned an error code:
!>             <= N:  if INFO = i, SSYEVX failed to converge;
!>                    i eigenvectors failed to converge.  Their indices
!>                    are stored in array IFAIL.
!>             > N:   if INFO = N + i, for 1 <= i <= N, then the leading
!>                    minor of order i of B is not positive definite.
!>                    The factorization of B could not be completed and
!>                    no eigenvalues or eigenvectors were computed.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA

Definition at line 294 of file ssygvx.f.

297*
298* -- LAPACK driver routine --
299* -- LAPACK is a software package provided by Univ. of Tennessee, --
300* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
301*
302* .. Scalar Arguments ..
303 CHARACTER JOBZ, RANGE, UPLO
304 INTEGER IL, INFO, ITYPE, IU, LDA, LDB, LDZ, LWORK, M, N
305 REAL ABSTOL, VL, VU
306* ..
307* .. Array Arguments ..
308 INTEGER IFAIL( * ), IWORK( * )
309 REAL A( LDA, * ), B( LDB, * ), W( * ), WORK( * ),
310 $ Z( LDZ, * )
311* ..
312*
313* =====================================================================
314*
315* .. Parameters ..
316 REAL ONE
317 parameter( one = 1.0e+0 )
318* ..
319* .. Local Scalars ..
320 LOGICAL ALLEIG, INDEIG, LQUERY, UPPER, VALEIG, WANTZ
321 CHARACTER TRANS
322 INTEGER LWKMIN, LWKOPT, NB
323* ..
324* .. External Functions ..
325 LOGICAL LSAME
326 INTEGER ILAENV
327 EXTERNAL ilaenv, lsame
328* ..
329* .. External Subroutines ..
330 EXTERNAL spotrf, ssyevx, ssygst, strmm, strsm, xerbla
331* ..
332* .. Intrinsic Functions ..
333 INTRINSIC max, min
334* ..
335* .. Executable Statements ..
336*
337* Test the input parameters.
338*
339 upper = lsame( uplo, 'U' )
340 wantz = lsame( jobz, 'V' )
341 alleig = lsame( range, 'A' )
342 valeig = lsame( range, 'V' )
343 indeig = lsame( range, 'I' )
344 lquery = ( lwork.EQ.-1 )
345*
346 info = 0
347 IF( itype.LT.1 .OR. itype.GT.3 ) THEN
348 info = -1
349 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
350 info = -2
351 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
352 info = -3
353 ELSE IF( .NOT.( upper .OR. lsame( uplo, 'L' ) ) ) THEN
354 info = -4
355 ELSE IF( n.LT.0 ) THEN
356 info = -5
357 ELSE IF( lda.LT.max( 1, n ) ) THEN
358 info = -7
359 ELSE IF( ldb.LT.max( 1, n ) ) THEN
360 info = -9
361 ELSE
362 IF( valeig ) THEN
363 IF( n.GT.0 .AND. vu.LE.vl )
364 $ info = -11
365 ELSE IF( indeig ) THEN
366 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
367 info = -12
368 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
369 info = -13
370 END IF
371 END IF
372 END IF
373 IF (info.EQ.0) THEN
374 IF (ldz.LT.1 .OR. (wantz .AND. ldz.LT.n)) THEN
375 info = -18
376 END IF
377 END IF
378*
379 IF( info.EQ.0 ) THEN
380 lwkmin = max( 1, 8*n )
381 nb = ilaenv( 1, 'SSYTRD', uplo, n, -1, -1, -1 )
382 lwkopt = max( lwkmin, ( nb + 3 )*n )
383 work( 1 ) = lwkopt
384*
385 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
386 info = -20
387 END IF
388 END IF
389*
390 IF( info.NE.0 ) THEN
391 CALL xerbla( 'SSYGVX', -info )
392 RETURN
393 ELSE IF( lquery ) THEN
394 RETURN
395 END IF
396*
397* Quick return if possible
398*
399 m = 0
400 IF( n.EQ.0 ) THEN
401 RETURN
402 END IF
403*
404* Form a Cholesky factorization of B.
405*
406 CALL spotrf( uplo, n, b, ldb, info )
407 IF( info.NE.0 ) THEN
408 info = n + info
409 RETURN
410 END IF
411*
412* Transform problem to standard eigenvalue problem and solve.
413*
414 CALL ssygst( itype, uplo, n, a, lda, b, ldb, info )
415 CALL ssyevx( jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol,
416 $ m, w, z, ldz, work, lwork, iwork, ifail, info )
417*
418 IF( wantz ) THEN
419*
420* Backtransform eigenvectors to the original problem.
421*
422 IF( info.GT.0 )
423 $ m = info - 1
424 IF( itype.EQ.1 .OR. itype.EQ.2 ) THEN
425*
426* For A*x=(lambda)*B*x and A*B*x=(lambda)*x;
427* backtransform eigenvectors: x = inv(L)**T*y or inv(U)*y
428*
429 IF( upper ) THEN
430 trans = 'N'
431 ELSE
432 trans = 'T'
433 END IF
434*
435 CALL strsm( 'Left', uplo, trans, 'Non-unit', n, m, one, b,
436 $ ldb, z, ldz )
437*
438 ELSE IF( itype.EQ.3 ) THEN
439*
440* For B*A*x=(lambda)*x;
441* backtransform eigenvectors: x = L*y or U**T*y
442*
443 IF( upper ) THEN
444 trans = 'T'
445 ELSE
446 trans = 'N'
447 END IF
448*
449 CALL strmm( 'Left', uplo, trans, 'Non-unit', n, m, one, b,
450 $ ldb, z, ldz )
451 END IF
452 END IF
453*
454* Set WORK(1) to optimal workspace size.
455*
456 work( 1 ) = lwkopt
457*
458 RETURN
459*
460* End of SSYGVX
461*
ssyevx
subroutine ssyevx(jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, work, lwork, iwork, ifail, info)
SSYEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices
Definition ssyevx.f:253