LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages

◆ zstein()

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

ZSTEIN

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

Purpose:
!>
!> ZSTEIN 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).
!>
!> Although the eigenvectors are real, they are stored in a complex
!> array, which may be passed to ZUNMTR or ZUPMTR for back
!> transformation to the eigenvectors of a complex Hermitian matrix
!> which was reduced to tridiagonal form.
!>
!>
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, stored 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 COMPLEX*16 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.
!>          The imaginary parts of the eigenvectors are set to zero.
!>
[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 178 of file zstein.f.

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