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

◆ dstein()

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

DSTEIN

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

Purpose:
 DSTEIN 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 DOUBLE PRECISION array, dimension (N)
          The n diagonal elements of the tridiagonal matrix T.
[in]E
          E is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DSTEBZ 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 DSTEBZ 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 DSTEBZ is expected here. )
[out]Z
          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 dstein.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 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
186* ..
187*
188* =====================================================================
189*
190* .. Parameters ..
191 DOUBLE PRECISION ZERO, ONE, TEN, ODM3, ODM1
192 parameter( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
193 $ odm3 = 1.0d-3, odm1 = 1.0d-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 DOUBLE PRECISION DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
202 $ SCL, SEP, TOL, XJ, XJM, ZTR
203* ..
204* .. Local Arrays ..
205 INTEGER ISEED( 4 )
206* ..
207* .. External Functions ..
208 INTEGER IDAMAX
209 DOUBLE PRECISION DDOT, DLAMCH, DNRM2
210 EXTERNAL idamax, ddot, dlamch, dnrm2
211* ..
212* .. External Subroutines ..
213 EXTERNAL daxpy, dcopy, dlagtf, dlagts, dlarnv, dscal,
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( 'DSTEIN', -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 = dlamch( 'Precision' )
266*
267* Initialize seed for random number generator DLARNV.
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 dtpcrt = 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 dlarnv( 2, iseed, blksiz, work( indrv1+1 ) )
347*
348* Copy the matrix T so it won't be destroyed in factorization.
349*
350 CALL dcopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
351 CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
352 CALL dcopy( 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 dlagtf( 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 = idamax( blksiz, work( indrv1+1 ), 1 )
371 scl = blksiz*onenrm*max( eps,
372 $ abs( work( indrv4+blksiz ) ) ) /
373 $ abs( work( indrv1+jmax ) )
374 CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
375*
376* Solve the system LU = Pb.
377*
378 CALL dlagts( -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 ztr = -ddot( blksiz, work( indrv1+1 ), 1, z( b1, i ),
392 $ 1 )
393 CALL daxpy( blksiz, ztr, 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 = idamax( 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.dtpcrt )
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 / dnrm2( blksiz, work( indrv1+1 ), 1 )
426 jmax = idamax( blksiz, work( indrv1+1 ), 1 )
427 IF( work( indrv1+jmax ).LT.zero )
428 $ scl = -scl
429 CALL dscal( 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 DSTEIN
449*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
ddot
double precision function ddot(n, dx, incx, dy, incy)
DDOT
Definition ddot.f:82
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dlagtf
subroutine dlagtf(n, a, lambda, b, c, tol, d, in, info)
DLAGTF computes an LU factorization of a matrix T-λI, where T is a general tridiagonal matrix,...
Definition dlagtf.f:156
dlagts
subroutine dlagts(job, n, a, b, c, d, in, y, tol, info)
DLAGTS solves the system of equations (T-λI)x = y or (T-λI)^Tx = y, where T is a general tridiagonal ...
Definition dlagts.f:163
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlarnv
subroutine dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition dlarnv.f:97
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: