SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sstein2.f
Go to the documentation of this file.
1*
2*
3 SUBROUTINE sstein2( 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 REAL ORFAC
14* ..
15* .. Array Arguments ..
16 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
17 $ iwork( * )
18 REAL D( * ), E( * ), W( * ), WORK( * ), Z( LDZ, * )
19* ..
20*
21* Purpose
22* =======
23*
24* SSTEIN2 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) REAL array, dimension (N)
38* The n diagonal elements of the tridiagonal matrix T.
39*
40* E (input) REAL 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) REAL 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 SSTEBZ 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 SSTEBZ 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 SSTEBZ is expected here. )
67*
68* ORFAC (input) REAL
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) REAL 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) REAL 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 REAL ZERO, ONE, TEN, ODM1
114 parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
115 $ odm1 = 1.0e-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 REAL 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 ISAMAX
131 REAL SASUM, SDOT, SLAMCH, SNRM2
132 EXTERNAL isamax, sasum, sdot, slamch, snrm2
133* ..
134* .. External Subroutines ..
135 EXTERNAL saxpy, scopy, slagtf, slagts, slarnv, sscal,
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( 'SSTEIN2', -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 = slamch( 'Precision' )
190*
191* Initialize seed for random number generator SLARNV.
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 slarnv( 2, iseed, blksiz, work( indrv1+1 ) )
271*
272* Copy the matrix T so it won't be destroyed in factorization.
273*
274 CALL scopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
275 CALL scopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
276 CALL scopy( 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 slagtf( 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 $ sasum( blksiz, work( indrv1+1 ), 1 )
297 CALL sscal( blksiz, scl, work( indrv1+1 ), 1 )
298*
299* Solve the system LU = Pb.
300*
301 CALL slagts( -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 = -sdot( blksiz, work( indrv1+1 ), 1, z( b1, i ),
316 $ 1 )
317 CALL saxpy( 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 = isamax( 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 / snrm2( blksiz, work( indrv1+1 ), 1 )
350 jmax = isamax( blksiz, work( indrv1+1 ), 1 )
351 IF( work( indrv1+jmax ).LT.zero )
352 $ scl = -scl
353 CALL sscal( 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 SSTEIN2
373*
374 END
max
#define max(A, B)
Definition pcgemr.c:180
sstein2
subroutine sstein2(n, d, e, m, w, iblock, isplit, orfac, z, ldz, work, iwork, ifail, info)
Definition sstein2.f:5