ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
dstein2.f
Go to the documentation of this file.
1 *
2 *
3  SUBROUTINE dstein2( N, D, E, M, W, IBLOCK, ISPLIT, ORFAC, Z, LDZ,
4  $ WORK, IWORK, IFAIL, INFO )
5 *
6 * -- ScaLAPACK routine (version 1.7) --
7 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
8 * and University of California, Berkeley.
9 * May 1, 1997
10 *
11 * .. Scalar Arguments ..
12  INTEGER INFO, LDZ, M, N
13  DOUBLE PRECISION ORFAC
14 * ..
15 * .. Array Arguments ..
16  INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
17  $ iwork( * )
18  DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DSTEIN2 computes the eigenvectors of a real symmetric tridiagonal
25 * matrix T corresponding to specified eigenvalues, using inverse
26 * iteration.
27 *
28 * The maximum number of iterations allowed for each eigenvector is
29 * specified by an internal parameter MAXITS (currently set to 5).
30 *
31 * Arguments
32 * =========
33 *
34 * N (input) INTEGER
35 * The order of the matrix. N >= 0.
36 *
37 * D (input) DOUBLE PRECISION array, dimension (N)
38 * The n diagonal elements of the tridiagonal matrix T.
39 *
40 * E (input) DOUBLE PRECISION array, dimension (N)
41 * The (n-1) subdiagonal elements of the tridiagonal matrix
42 * T, in elements 1 to N-1. E(N) need not be set.
43 *
44 * M (input) INTEGER
45 * The number of eigenvectors to be found. 0 <= M <= N.
46 *
47 * W (input) DOUBLE PRECISION array, dimension (N)
48 * The first M elements of W contain the eigenvalues for
49 * which eigenvectors are to be computed. The eigenvalues
50 * should be grouped by split-off block and ordered from
51 * smallest to largest within the block. ( The output array
52 * W from DSTEBZ with ORDER = 'B' is expected here. )
53 *
54 * IBLOCK (input) INTEGER array, dimension (N)
55 * The submatrix indices associated with the corresponding
56 * eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
57 * the first submatrix from the top, =2 if W(i) belongs to
58 * the second submatrix, etc. ( The output array IBLOCK
59 * from DSTEBZ is expected here. )
60 *
61 * ISPLIT (input) INTEGER array, dimension (N)
62 * The splitting points, at which T breaks up into submatrices.
63 * The first submatrix consists of rows/columns 1 to
64 * ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
65 * through ISPLIT( 2 ), etc.
66 * ( The output array ISPLIT from DSTEBZ is expected here. )
67 *
68 * ORFAC (input) DOUBLE PRECISION
69 * ORFAC specifies which eigenvectors should be
70 * orthogonalized. Eigenvectors that correspond to eigenvalues
71 * which are within ORFAC*||T|| of each other are to be
72 * orthogonalized.
73 *
74 * Z (output) DOUBLE PRECISION array, dimension (LDZ, M)
75 * The computed eigenvectors. The eigenvector associated
76 * with the eigenvalue W(i) is stored in the i-th column of
77 * Z. Any vector which fails to converge is set to its current
78 * iterate after MAXITS iterations.
79 *
80 * LDZ (input) INTEGER
81 * The leading dimension of the array Z. LDZ >= max(1,N).
82 *
83 * WORK (workspace) DOUBLE PRECISION array, dimension (5*N)
84 *
85 * IWORK (workspace) INTEGER array, dimension (N)
86 *
87 * IFAIL (output) INTEGER array, dimension (M)
88 * On normal exit, all elements of IFAIL are zero.
89 * If one or more eigenvectors fail to converge after
90 * MAXITS iterations, then their indices are stored in
91 * array IFAIL.
92 *
93 * INFO (output) INTEGER
94 * = 0: successful exit.
95 * < 0: if INFO = -i, the i-th argument had an illegal value
96 * > 0: if INFO = i, then i eigenvectors failed to converge
97 * in MAXITS iterations. Their indices are stored in
98 * array IFAIL.
99 *
100 * Internal Parameters
101 * ===================
102 *
103 * MAXITS INTEGER, default = 5
104 * The maximum number of iterations performed.
105 *
106 * EXTRA INTEGER, default = 2
107 * The number of iterations performed after norm growth
108 * criterion is satisfied, should be at least 1.
109 *
110 * =====================================================================
111 *
112 * .. Parameters ..
113  DOUBLE PRECISION ZERO, ONE, TEN, ODM1
114  parameter( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
115  $ odm1 = 1.0d-1 )
116  INTEGER MAXITS, EXTRA
117  parameter( maxits = 5, extra = 2 )
118 * ..
119 * .. Local Scalars ..
120  INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
121  $ indrv2, indrv3, indrv4, indrv5, its, j, j1,
122  $ jblk, jmax, nblk, nrmchk
123  DOUBLE PRECISION EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL, SCL,
124  $ sep, stpcrt, tol, xj, xjm, ztr
125 * ..
126 * .. Local Arrays ..
127  INTEGER ISEED( 4 )
128 * ..
129 * .. External Functions ..
130  INTEGER IDAMAX
131  DOUBLE PRECISION DASUM, DDOT, DLAMCH, DNRM2
132  EXTERNAL idamax, dasum, ddot, dlamch, dnrm2
133 * ..
134 * .. External Subroutines ..
135  EXTERNAL daxpy, dcopy, dlagtf, dlagts, dlarnv, dscal,
136  $ xerbla
137 * ..
138 * .. Intrinsic Functions ..
139  INTRINSIC abs, max, sqrt
140 * ..
141 * .. Executable Statements ..
142 *
143 * Test the input parameters.
144 *
145  info = 0
146  DO 10 i = 1, m
147  ifail( i ) = 0
148  10 CONTINUE
149 *
150  IF( n.LT.0 ) THEN
151  info = -1
152  ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
153  info = -4
154  ELSE IF( orfac.LT.zero ) THEN
155  info = -8
156  ELSE IF( ldz.LT.max( 1, n ) ) THEN
157  info = -10
158  ELSE
159  DO 20 j = 2, m
160  IF( iblock( j ).LT.iblock( j-1 ) ) THEN
161  info = -6
162  GO TO 30
163  END IF
164  IF( iblock( j ).EQ.iblock( j-1 ) .AND. w( j ).LT.w( j-1 ) )
165  $ THEN
166  info = -5
167  GO TO 30
168  END IF
169  20 CONTINUE
170  30 CONTINUE
171  END IF
172 *
173  IF( info.NE.0 ) THEN
174  CALL xerbla( 'DSTEIN2', -info )
175  RETURN
176  END IF
177 *
178 * Quick return if possible
179 *
180  IF( n.EQ.0 .OR. m.EQ.0 ) THEN
181  RETURN
182  ELSE IF( n.EQ.1 ) THEN
183  z( 1, 1 ) = one
184  RETURN
185  END IF
186 *
187 * Get machine constants.
188 *
189  eps = dlamch( 'Precision' )
190 *
191 * Initialize seed for random number generator DLARNV.
192 *
193  DO 40 i = 1, 4
194  iseed( i ) = 1
195  40 CONTINUE
196 *
197 * Initialize pointers.
198 *
199  indrv1 = 0
200  indrv2 = indrv1 + n
201  indrv3 = indrv2 + n
202  indrv4 = indrv3 + n
203  indrv5 = indrv4 + n
204 *
205 * Compute eigenvectors of matrix blocks.
206 *
207  j1 = 1
208  DO 160 nblk = 1, iblock( m )
209 *
210 * Find starting and ending indices of block nblk.
211 *
212  IF( nblk.EQ.1 ) THEN
213  b1 = 1
214  ELSE
215  b1 = isplit( nblk-1 ) + 1
216  END IF
217  bn = isplit( nblk )
218  blksiz = bn - b1 + 1
219  IF( blksiz.EQ.1 )
220  $ GO TO 60
221  gpind = j1
222 *
223 * Compute reorthogonalization criterion and stopping criterion.
224 *
225  onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
226  onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
227  DO 50 i = b1 + 1, bn - 1
228  onenrm = max( onenrm, abs( d( i ) )+abs( e( i-1 ) )+
229  $ abs( e( i ) ) )
230  50 CONTINUE
231  ortol = orfac*onenrm
232 *
233  stpcrt = sqrt( odm1 / blksiz )
234 *
235 * Loop through eigenvalues of block nblk.
236 *
237  60 CONTINUE
238  jblk = 0
239  DO 150 j = j1, m
240  IF( iblock( j ).NE.nblk ) THEN
241  j1 = j
242  GO TO 160
243  END IF
244  jblk = jblk + 1
245  xj = w( j )
246 *
247 * Skip all the work if the block size is one.
248 *
249  IF( blksiz.EQ.1 ) THEN
250  work( indrv1+1 ) = one
251  GO TO 120
252  END IF
253 *
254 * If eigenvalues j and j-1 are too close, add a relatively
255 * small perturbation.
256 *
257  IF( jblk.GT.1 ) THEN
258  eps1 = abs( eps*xj )
259  pertol = ten*eps1
260  sep = xj - xjm
261  IF( sep.LT.pertol )
262  $ xj = xjm + pertol
263  END IF
264 *
265  its = 0
266  nrmchk = 0
267 *
268 * Get random starting vector.
269 *
270  CALL dlarnv( 2, iseed, blksiz, work( indrv1+1 ) )
271 *
272 * Copy the matrix T so it won't be destroyed in factorization.
273 *
274  CALL dcopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
275  CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
276  CALL dcopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
277 *
278 * Compute LU factors with partial pivoting ( PT = LU )
279 *
280  tol = zero
281  CALL dlagtf( blksiz, work( indrv4+1 ), xj, work( indrv2+2 ),
282  $ work( indrv3+1 ), tol, work( indrv5+1 ), iwork,
283  $ iinfo )
284 *
285 * Update iteration count.
286 *
287  70 CONTINUE
288  its = its + 1
289  IF( its.GT.maxits )
290  $ GO TO 100
291 *
292 * Normalize and scale the righthand side vector Pb.
293 *
294  scl = blksiz*onenrm*max( eps,
295  $ abs( work( indrv4+blksiz ) ) ) /
296  $ dasum( blksiz, work( indrv1+1 ), 1 )
297  CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
298 *
299 * Solve the system LU = Pb.
300 *
301  CALL dlagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
302  $ work( indrv3+1 ), work( indrv5+1 ), iwork,
303  $ work( indrv1+1 ), tol, iinfo )
304 *
305 * Reorthogonalize by modified Gram-Schmidt if eigenvalues are
306 * close enough.
307 *
308  IF( jblk.EQ.1 )
309  $ GO TO 90
310  IF( abs( xj-xjm ).GT.ortol )
311  $ gpind = j
312 *
313  IF( gpind.NE.j ) THEN
314  DO 80 i = gpind, j - 1
315  ztr = -ddot( blksiz, work( indrv1+1 ), 1, z( b1, i ),
316  $ 1 )
317  CALL daxpy( blksiz, ztr, z( b1, i ), 1,
318  $ work( indrv1+1 ), 1 )
319  80 CONTINUE
320  END IF
321 *
322 * Check the infinity norm of the iterate.
323 *
324  90 CONTINUE
325  jmax = idamax( blksiz, work( indrv1+1 ), 1 )
326  nrm = abs( work( indrv1+jmax ) )
327 *
328 * Continue for additional iterations after norm reaches
329 * stopping criterion.
330 *
331  IF( nrm.LT.stpcrt )
332  $ GO TO 70
333  nrmchk = nrmchk + 1
334  IF( nrmchk.LT.extra+1 )
335  $ GO TO 70
336 *
337  GO TO 110
338 *
339 * If stopping criterion was not satisfied, update info and
340 * store eigenvector number in array ifail.
341 *
342  100 CONTINUE
343  info = info + 1
344  ifail( info ) = j
345 *
346 * Accept iterate as jth eigenvector.
347 *
348  110 CONTINUE
349  scl = one / dnrm2( blksiz, work( indrv1+1 ), 1 )
350  jmax = idamax( blksiz, work( indrv1+1 ), 1 )
351  IF( work( indrv1+jmax ).LT.zero )
352  $ scl = -scl
353  CALL dscal( blksiz, scl, work( indrv1+1 ), 1 )
354  120 CONTINUE
355  DO 130 i = 1, n
356  z( i, j ) = zero
357  130 CONTINUE
358  DO 140 i = 1, blksiz
359  z( b1+i-1, j ) = work( indrv1+i )
360  140 CONTINUE
361 *
362 * Save the shift to check eigenvalue spacing at next
363 * iteration.
364 *
365  xjm = xj
366 *
367  150 CONTINUE
368  160 CONTINUE
369 *
370  RETURN
371 *
372 * End of DSTEIN2
373 *
374  END
max
#define max(A, B)
Definition: pcgemr.c:180
dstein2
subroutine dstein2(N, D, E, M, W, IBLOCK, ISPLIT, ORFAC, Z, LDZ, WORK, IWORK, IFAIL, INFO)
Definition: dstein2.f:5