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

◆ sstein()

subroutine sstein ( integer  n,
real, dimension( * )  d,
real, dimension( * )  e,
integer  m,
real, dimension( * )  w,
integer, dimension( * )  iblock,
integer, dimension( * )  isplit,
real, dimension( ldz, * )  z,
integer  ldz,
real, dimension( * )  work,
integer, dimension( * )  iwork,
integer, dimension( * )  ifail,
integer  info 
)

SSTEIN

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

Purpose:
 SSTEIN computes the eigenvectors of a real symmetric tridiagonal
 matrix T corresponding to specified eigenvalues, using inverse
 iteration.

 The maximum number of iterations allowed for each eigenvector is
 specified by an internal parameter MAXITS (currently set to 5).
Parameters
[in]N
          N is INTEGER
          The order of the matrix.  N >= 0.
[in]D
          D is REAL array, dimension (N)
          The n diagonal elements of the tridiagonal matrix T.
[in]E
          E is REAL array, dimension (N-1)
          The (n-1) subdiagonal elements of the tridiagonal matrix
          T, in elements 1 to N-1.
[in]M
          M is INTEGER
          The number of eigenvectors to be found.  0 <= M <= N.
[in]W
          W is REAL array, dimension (N)
          The first M elements of W contain the eigenvalues for
          which eigenvectors are to be computed.  The eigenvalues
          should be grouped by split-off block and ordered from
          smallest to largest within the block.  ( The output array
          W from SSTEBZ with ORDER = 'B' is expected here. )
[in]IBLOCK
          IBLOCK is INTEGER array, dimension (N)
          The submatrix indices associated with the corresponding
          eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
          the first submatrix from the top, =2 if W(i) belongs to
          the second submatrix, etc.  ( The output array IBLOCK
          from SSTEBZ is expected here. )
[in]ISPLIT
          ISPLIT is INTEGER array, dimension (N)
          The splitting points, at which T breaks up into submatrices.
          The first submatrix consists of rows/columns 1 to
          ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
          through ISPLIT( 2 ), etc.
          ( The output array ISPLIT from SSTEBZ is expected here. )
[out]Z
          Z is REAL array, dimension (LDZ, M)
          The computed eigenvectors.  The eigenvector associated
          with the eigenvalue W(i) is stored in the i-th column of
          Z.  Any vector which fails to converge is set to its current
          iterate after MAXITS iterations.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= max(1,N).
[out]WORK
          WORK is REAL array, dimension (5*N)
[out]IWORK
          IWORK is INTEGER array, dimension (N)
[out]IFAIL
          IFAIL is INTEGER array, dimension (M)
          On normal exit, all elements of IFAIL are zero.
          If one or more eigenvectors fail to converge after
          MAXITS iterations, then their indices are stored in
          array IFAIL.
[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
               in MAXITS iterations.  Their indices are stored in
               array IFAIL.
Internal Parameters:
  MAXITS  INTEGER, default = 5
          The maximum number of iterations performed.

  EXTRA   INTEGER, default = 2
          The number of iterations performed after norm growth
          criterion is satisfied, should be at least 1.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 172 of file sstein.f.

174*
175* -- LAPACK computational routine --
176* -- LAPACK is a software package provided by Univ. of Tennessee, --
177* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178*
179* .. Scalar Arguments ..
180 INTEGER INFO, LDZ, M, N
181* ..
182* .. Array Arguments ..
183 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
184 $ IWORK( * )
185 REAL D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
186* ..
187*
188* =====================================================================
189*
190* .. Parameters ..
191 REAL ZERO, ONE, TEN, ODM3, ODM1
192 parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
193 $ odm3 = 1.0e-3, odm1 = 1.0e-1 )
194 INTEGER MAXITS, EXTRA
195 parameter( maxits = 5, extra = 2 )
196* ..
197* .. Local Scalars ..
198 INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
199 $ INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1,
200 $ JBLK, JMAX, NBLK, NRMCHK
201 REAL CTR, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
202 $ SCL, SEP, STPCRT, TOL, XJ, XJM
203* ..
204* .. Local Arrays ..
205 INTEGER ISEED( 4 )
206* ..
207* .. External Functions ..
208 INTEGER ISAMAX
209 REAL SDOT, SLAMCH, SNRM2
210 EXTERNAL isamax, sdot, slamch, snrm2
211* ..
212* .. External Subroutines ..
213 EXTERNAL saxpy, scopy, slagtf, slagts, slarnv, sscal,
214 $ xerbla
215* ..
216* .. Intrinsic Functions ..
217 INTRINSIC abs, max, sqrt
218* ..
219* .. Executable Statements ..
220*
221* Test the input parameters.
222*
223 info = 0
224 DO 10 i = 1, m
225 ifail( i ) = 0
226 10 CONTINUE
227*
228 IF( n.LT.0 ) THEN
229 info = -1
230 ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
231 info = -4
232 ELSE IF( ldz.LT.max( 1, n ) ) THEN
233 info = -9
234 ELSE
235 DO 20 j = 2, m
236 IF( iblock( j ).LT.iblock( j-1 ) ) THEN
237 info = -6
238 GO TO 30
239 END IF
240 IF( iblock( j ).EQ.iblock( j-1 ) .AND. w( j ).LT.w( j-1 ) )
241 $ THEN
242 info = -5
243 GO TO 30
244 END IF
245 20 CONTINUE
246 30 CONTINUE
247 END IF
248*
249 IF( info.NE.0 ) THEN
250 CALL xerbla( 'SSTEIN', -info )
251 RETURN
252 END IF
253*
254* Quick return if possible
255*
256 IF( n.EQ.0 .OR. m.EQ.0 ) THEN
257 RETURN
258 ELSE IF( n.EQ.1 ) THEN
259 z( 1, 1 ) = one
260 RETURN
261 END IF
262*
263* Get machine constants.
264*
265 eps = slamch( 'Precision' )
266*
267* Initialize seed for random number generator SLARNV.
268*
269 DO 40 i = 1, 4
270 iseed( i ) = 1
271 40 CONTINUE
272*
273* Initialize pointers.
274*
275 indrv1 = 0
276 indrv2 = indrv1 + n
277 indrv3 = indrv2 + n
278 indrv4 = indrv3 + n
279 indrv5 = indrv4 + n
280*
281* Compute eigenvectors of matrix blocks.
282*
283 j1 = 1
284 DO 160 nblk = 1, iblock( m )
285*
286* Find starting and ending indices of block nblk.
287*
288 IF( nblk.EQ.1 ) THEN
289 b1 = 1
290 ELSE
291 b1 = isplit( nblk-1 ) + 1
292 END IF
293 bn = isplit( nblk )
294 blksiz = bn - b1 + 1
295 IF( blksiz.EQ.1 )
296 $ GO TO 60
297 gpind = j1
298*
299* Compute reorthogonalization criterion and stopping criterion.
300*
301 onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
302 onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
303 DO 50 i = b1 + 1, bn - 1
304 onenrm = max( onenrm, abs( d( i ) )+abs( e( i-1 ) )+
305 $ abs( e( i ) ) )
306 50 CONTINUE
307 ortol = odm3*onenrm
308*
309 stpcrt = sqrt( odm1 / blksiz )
310*
311* Loop through eigenvalues of block nblk.
312*
313 60 CONTINUE
314 jblk = 0
315 DO 150 j = j1, m
316 IF( iblock( j ).NE.nblk ) THEN
317 j1 = j
318 GO TO 160
319 END IF
320 jblk = jblk + 1
321 xj = w( j )
322*
323* Skip all the work if the block size is one.
324*
325 IF( blksiz.EQ.1 ) THEN
326 work( indrv1+1 ) = one
327 GO TO 120
328 END IF
329*
330* If eigenvalues j and j-1 are too close, add a relatively
331* small perturbation.
332*
333 IF( jblk.GT.1 ) THEN
334 eps1 = abs( eps*xj )
335 pertol = ten*eps1
336 sep = xj - xjm
337 IF( sep.LT.pertol )
338 $ xj = xjm + pertol
339 END IF
340*
341 its = 0
342 nrmchk = 0
343*
344* Get random starting vector.
345*
346 CALL slarnv( 2, iseed, blksiz, work( indrv1+1 ) )
347*
348* Copy the matrix T so it won't be destroyed in factorization.
349*
350 CALL scopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
351 CALL scopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
352 CALL scopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
353*
354* Compute LU factors with partial pivoting ( PT = LU )
355*
356 tol = zero
357 CALL slagtf( blksiz, work( indrv4+1 ), xj, work( indrv2+2 ),
358 $ work( indrv3+1 ), tol, work( indrv5+1 ), iwork,
359 $ iinfo )
360*
361* Update iteration count.
362*
363 70 CONTINUE
364 its = its + 1
365 IF( its.GT.maxits )
366 $ GO TO 100
367*
368* Normalize and scale the righthand side vector Pb.
369*
370 jmax = isamax( blksiz, work( indrv1+1 ), 1 )
371 scl = blksiz*onenrm*max( eps,
372 $ abs( work( indrv4+blksiz ) ) ) /
373 $ abs( work( indrv1+jmax ) )
374 CALL sscal( blksiz, scl, work( indrv1+1 ), 1 )
375*
376* Solve the system LU = Pb.
377*
378 CALL slagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
379 $ work( indrv3+1 ), work( indrv5+1 ), iwork,
380 $ work( indrv1+1 ), tol, iinfo )
381*
382* Reorthogonalize by modified Gram-Schmidt if eigenvalues are
383* close enough.
384*
385 IF( jblk.EQ.1 )
386 $ GO TO 90
387 IF( abs( xj-xjm ).GT.ortol )
388 $ gpind = j
389 IF( gpind.NE.j ) THEN
390 DO 80 i = gpind, j - 1
391 ctr = -sdot( blksiz, work( indrv1+1 ), 1, z( b1, i ),
392 $ 1 )
393 CALL saxpy( blksiz, ctr, z( b1, i ), 1,
394 $ work( indrv1+1 ), 1 )
395 80 CONTINUE
396 END IF
397*
398* Check the infinity norm of the iterate.
399*
400 90 CONTINUE
401 jmax = isamax( blksiz, work( indrv1+1 ), 1 )
402 nrm = abs( work( indrv1+jmax ) )
403*
404* Continue for additional iterations after norm reaches
405* stopping criterion.
406*
407 IF( nrm.LT.stpcrt )
408 $ GO TO 70
409 nrmchk = nrmchk + 1
410 IF( nrmchk.LT.extra+1 )
411 $ GO TO 70
412*
413 GO TO 110
414*
415* If stopping criterion was not satisfied, update info and
416* store eigenvector number in array ifail.
417*
418 100 CONTINUE
419 info = info + 1
420 ifail( info ) = j
421*
422* Accept iterate as jth eigenvector.
423*
424 110 CONTINUE
425 scl = one / snrm2( blksiz, work( indrv1+1 ), 1 )
426 jmax = isamax( blksiz, work( indrv1+1 ), 1 )
427 IF( work( indrv1+jmax ).LT.zero )
428 $ scl = -scl
429 CALL sscal( blksiz, scl, work( indrv1+1 ), 1 )
430 120 CONTINUE
431 DO 130 i = 1, n
432 z( i, j ) = zero
433 130 CONTINUE
434 DO 140 i = 1, blksiz
435 z( b1+i-1, j ) = work( indrv1+i )
436 140 CONTINUE
437*
438* Save the shift to check eigenvalue spacing at next
439* iteration.
440*
441 xjm = xj
442*
443 150 CONTINUE
444 160 CONTINUE
445*
446 RETURN
447*
448* End of SSTEIN
449*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
saxpy
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sdot
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
slagtf
subroutine slagtf(n, a, lambda, b, c, tol, d, in, info)
SLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix,...
Definition slagtf.f:156
slagts
subroutine slagts(job, n, a, b, c, d, in, y, tol, info)
SLAGTS solves the system of equations (T-λI)x = y or (T-λI)^Tx = y, where T is a general tridiagonal ...
Definition slagts.f:163
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slarnv
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition slarnv.f:97
snrm2
real(wp) function snrm2(n, x, incx)
SNRM2
Definition snrm2.f90:89
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: