LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
cstein.f
Go to the documentation of this file.
1*> \brief \b CSTEIN
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cstein.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cstein.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cstein.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
22* IWORK, IFAIL, INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INFO, LDZ, M, N
26* ..
27* .. Array Arguments ..
28* INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
29* \$ IWORK( * )
30* REAL D( * ), E( * ), W( * ), WORK( * )
31* COMPLEX Z( LDZ, * )
32* ..
33*
34*
35*> \par Purpose:
36* =============
37*>
38*> \verbatim
39*>
40*> CSTEIN computes the eigenvectors of a real symmetric tridiagonal
41*> matrix T corresponding to specified eigenvalues, using inverse
42*> iteration.
43*>
44*> The maximum number of iterations allowed for each eigenvector is
45*> specified by an internal parameter MAXITS (currently set to 5).
46*>
47*> Although the eigenvectors are real, they are stored in a complex
48*> array, which may be passed to CUNMTR or CUPMTR for back
49*> transformation to the eigenvectors of a complex Hermitian matrix
50*> which was reduced to tridiagonal form.
51*>
52*> \endverbatim
53*
54* Arguments:
55* ==========
56*
57*> \param[in] N
58*> \verbatim
59*> N is INTEGER
60*> The order of the matrix. N >= 0.
61*> \endverbatim
62*>
63*> \param[in] D
64*> \verbatim
65*> D is REAL array, dimension (N)
66*> The n diagonal elements of the tridiagonal matrix T.
67*> \endverbatim
68*>
69*> \param[in] E
70*> \verbatim
71*> E is REAL array, dimension (N-1)
72*> The (n-1) subdiagonal elements of the tridiagonal matrix
73*> T, stored in elements 1 to N-1.
74*> \endverbatim
75*>
76*> \param[in] M
77*> \verbatim
78*> M is INTEGER
79*> The number of eigenvectors to be found. 0 <= M <= N.
80*> \endverbatim
81*>
82*> \param[in] W
83*> \verbatim
84*> W is REAL array, dimension (N)
85*> The first M elements of W contain the eigenvalues for
86*> which eigenvectors are to be computed. The eigenvalues
87*> should be grouped by split-off block and ordered from
88*> smallest to largest within the block. ( The output array
89*> W from SSTEBZ with ORDER = 'B' is expected here. )
90*> \endverbatim
91*>
92*> \param[in] IBLOCK
93*> \verbatim
94*> IBLOCK is INTEGER array, dimension (N)
95*> The submatrix indices associated with the corresponding
96*> eigenvalues in W; IBLOCK(i)=1 if eigenvalue W(i) belongs to
97*> the first submatrix from the top, =2 if W(i) belongs to
98*> the second submatrix, etc. ( The output array IBLOCK
99*> from SSTEBZ is expected here. )
100*> \endverbatim
101*>
102*> \param[in] ISPLIT
103*> \verbatim
104*> ISPLIT is INTEGER array, dimension (N)
105*> The splitting points, at which T breaks up into submatrices.
106*> The first submatrix consists of rows/columns 1 to
107*> ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1
108*> through ISPLIT( 2 ), etc.
109*> ( The output array ISPLIT from SSTEBZ is expected here. )
110*> \endverbatim
111*>
112*> \param[out] Z
113*> \verbatim
114*> Z is COMPLEX array, dimension (LDZ, M)
115*> The computed eigenvectors. The eigenvector associated
116*> with the eigenvalue W(i) is stored in the i-th column of
117*> Z. Any vector which fails to converge is set to its current
118*> iterate after MAXITS iterations.
119*> The imaginary parts of the eigenvectors are set to zero.
120*> \endverbatim
121*>
122*> \param[in] LDZ
123*> \verbatim
124*> LDZ is INTEGER
125*> The leading dimension of the array Z. LDZ >= max(1,N).
126*> \endverbatim
127*>
128*> \param[out] WORK
129*> \verbatim
130*> WORK is REAL array, dimension (5*N)
131*> \endverbatim
132*>
133*> \param[out] IWORK
134*> \verbatim
135*> IWORK is INTEGER array, dimension (N)
136*> \endverbatim
137*>
138*> \param[out] IFAIL
139*> \verbatim
140*> IFAIL is INTEGER array, dimension (M)
141*> On normal exit, all elements of IFAIL are zero.
142*> If one or more eigenvectors fail to converge after
143*> MAXITS iterations, then their indices are stored in
144*> array IFAIL.
145*> \endverbatim
146*>
147*> \param[out] INFO
148*> \verbatim
149*> INFO is INTEGER
150*> = 0: successful exit
151*> < 0: if INFO = -i, the i-th argument had an illegal value
152*> > 0: if INFO = i, then i eigenvectors failed to converge
153*> in MAXITS iterations. Their indices are stored in
154*> array IFAIL.
155*> \endverbatim
156*
157*> \par Internal Parameters:
158* =========================
159*>
160*> \verbatim
161*> MAXITS INTEGER, default = 5
162*> The maximum number of iterations performed.
163*>
164*> EXTRA INTEGER, default = 2
165*> The number of iterations performed after norm growth
166*> criterion is satisfied, should be at least 1.
167*> \endverbatim
168*
169* Authors:
170* ========
171*
172*> \author Univ. of Tennessee
173*> \author Univ. of California Berkeley
174*> \author Univ. of Colorado Denver
175*> \author NAG Ltd.
176*
177*> \ingroup stein
178*
179* =====================================================================
180 SUBROUTINE cstein( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
181 \$ IWORK, IFAIL, INFO )
182*
183* -- LAPACK computational routine --
184* -- LAPACK is a software package provided by Univ. of Tennessee, --
185* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
186*
187* .. Scalar Arguments ..
188 INTEGER INFO, LDZ, M, N
189* ..
190* .. Array Arguments ..
191 INTEGER IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
192 \$ iwork( * )
193 REAL D( * ), E( * ), W( * ), WORK( * )
194 COMPLEX Z( LDZ, * )
195* ..
196*
197* =====================================================================
198*
199* .. Parameters ..
200 COMPLEX CZERO, CONE
201 parameter( czero = ( 0.0e+0, 0.0e+0 ),
202 \$ cone = ( 1.0e+0, 0.0e+0 ) )
203 REAL ZERO, ONE, TEN, ODM3, ODM1
204 parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
205 \$ odm3 = 1.0e-3, odm1 = 1.0e-1 )
206 INTEGER MAXITS, EXTRA
207 parameter( maxits = 5, extra = 2 )
208* ..
209* .. Local Scalars ..
210 INTEGER B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
211 \$ indrv2, indrv3, indrv4, indrv5, its, j, j1,
212 \$ jblk, jmax, jr, nblk, nrmchk
213 REAL CTR, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
214 \$ scl, sep, stpcrt, tol, xj, xjm
215* ..
216* .. Local Arrays ..
217 INTEGER ISEED( 4 )
218* ..
219* .. External Functions ..
220 INTEGER ISAMAX
221 REAL SLAMCH, SNRM2
222 EXTERNAL isamax, slamch, snrm2
223* ..
224* .. External Subroutines ..
225 EXTERNAL scopy, slagtf, slagts, slarnv, sscal, xerbla
226* ..
227* .. Intrinsic Functions ..
228 INTRINSIC abs, cmplx, max, real, sqrt
229* ..
230* .. Executable Statements ..
231*
232* Test the input parameters.
233*
234 info = 0
235 DO 10 i = 1, m
236 ifail( i ) = 0
237 10 CONTINUE
238*
239 IF( n.LT.0 ) THEN
240 info = -1
241 ELSE IF( m.LT.0 .OR. m.GT.n ) THEN
242 info = -4
243 ELSE IF( ldz.LT.max( 1, n ) ) THEN
244 info = -9
245 ELSE
246 DO 20 j = 2, m
247 IF( iblock( j ).LT.iblock( j-1 ) ) THEN
248 info = -6
249 GO TO 30
250 END IF
251 IF( iblock( j ).EQ.iblock( j-1 ) .AND. w( j ).LT.w( j-1 ) )
252 \$ THEN
253 info = -5
254 GO TO 30
255 END IF
256 20 CONTINUE
257 30 CONTINUE
258 END IF
259*
260 IF( info.NE.0 ) THEN
261 CALL xerbla( 'CSTEIN', -info )
262 RETURN
263 END IF
264*
265* Quick return if possible
266*
267 IF( n.EQ.0 .OR. m.EQ.0 ) THEN
268 RETURN
269 ELSE IF( n.EQ.1 ) THEN
270 z( 1, 1 ) = cone
271 RETURN
272 END IF
273*
274* Get machine constants.
275*
276 eps = slamch( 'Precision' )
277*
278* Initialize seed for random number generator SLARNV.
279*
280 DO 40 i = 1, 4
281 iseed( i ) = 1
282 40 CONTINUE
283*
284* Initialize pointers.
285*
286 indrv1 = 0
287 indrv2 = indrv1 + n
288 indrv3 = indrv2 + n
289 indrv4 = indrv3 + n
290 indrv5 = indrv4 + n
291*
292* Compute eigenvectors of matrix blocks.
293*
294 j1 = 1
295 DO 180 nblk = 1, iblock( m )
296*
297* Find starting and ending indices of block nblk.
298*
299 IF( nblk.EQ.1 ) THEN
300 b1 = 1
301 ELSE
302 b1 = isplit( nblk-1 ) + 1
303 END IF
304 bn = isplit( nblk )
305 blksiz = bn - b1 + 1
306 IF( blksiz.EQ.1 )
307 \$ GO TO 60
308 gpind = j1
309*
310* Compute reorthogonalization criterion and stopping criterion.
311*
312 onenrm = abs( d( b1 ) ) + abs( e( b1 ) )
313 onenrm = max( onenrm, abs( d( bn ) )+abs( e( bn-1 ) ) )
314 DO 50 i = b1 + 1, bn - 1
315 onenrm = max( onenrm, abs( d( i ) )+abs( e( i-1 ) )+
316 \$ abs( e( i ) ) )
317 50 CONTINUE
318 ortol = odm3*onenrm
319*
320 stpcrt = sqrt( odm1 / blksiz )
321*
322* Loop through eigenvalues of block nblk.
323*
324 60 CONTINUE
325 jblk = 0
326 DO 170 j = j1, m
327 IF( iblock( j ).NE.nblk ) THEN
328 j1 = j
329 GO TO 180
330 END IF
331 jblk = jblk + 1
332 xj = w( j )
333*
334* Skip all the work if the block size is one.
335*
336 IF( blksiz.EQ.1 ) THEN
337 work( indrv1+1 ) = one
338 GO TO 140
339 END IF
340*
341* If eigenvalues j and j-1 are too close, add a relatively
342* small perturbation.
343*
344 IF( jblk.GT.1 ) THEN
345 eps1 = abs( eps*xj )
346 pertol = ten*eps1
347 sep = xj - xjm
348 IF( sep.LT.pertol )
349 \$ xj = xjm + pertol
350 END IF
351*
352 its = 0
353 nrmchk = 0
354*
355* Get random starting vector.
356*
357 CALL slarnv( 2, iseed, blksiz, work( indrv1+1 ) )
358*
359* Copy the matrix T so it won't be destroyed in factorization.
360*
361 CALL scopy( blksiz, d( b1 ), 1, work( indrv4+1 ), 1 )
362 CALL scopy( blksiz-1, e( b1 ), 1, work( indrv2+2 ), 1 )
363 CALL scopy( blksiz-1, e( b1 ), 1, work( indrv3+1 ), 1 )
364*
365* Compute LU factors with partial pivoting ( PT = LU )
366*
367 tol = zero
368 CALL slagtf( blksiz, work( indrv4+1 ), xj, 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 = isamax( blksiz, work( indrv1+1 ), 1 )
382 scl = blksiz*onenrm*max( eps,
383 \$ abs( work( indrv4+blksiz ) ) ) /
384 \$ abs( work( indrv1+jmax ) )
385 CALL sscal( blksiz, scl, work( indrv1+1 ), 1 )
386*
387* Solve the system LU = Pb.
388*
389 CALL slagts( -1, blksiz, work( indrv4+1 ), work( indrv2+2 ),
390 \$ work( indrv3+1 ), work( indrv5+1 ), iwork,
391 \$ work( indrv1+1 ), tol, iinfo )
392*
393* Reorthogonalize by modified Gram-Schmidt if eigenvalues are
394* close enough.
395*
396 IF( jblk.EQ.1 )
397 \$ GO TO 110
398 IF( abs( xj-xjm ).GT.ortol )
399 \$ gpind = j
400 IF( gpind.NE.j ) THEN
401 DO 100 i = gpind, j - 1
402 ctr = zero
403 DO 80 jr = 1, blksiz
404 ctr = ctr + work( indrv1+jr )*
405 \$ real( z( b1-1+jr, i ) )
406 80 CONTINUE
407 DO 90 jr = 1, blksiz
408 work( indrv1+jr ) = work( indrv1+jr ) -
409 \$ ctr*real( z( b1-1+jr, i ) )
410 90 CONTINUE
411 100 CONTINUE
412 END IF
413*
414* Check the infinity norm of the iterate.
415*
416 110 CONTINUE
417 jmax = isamax( blksiz, work( indrv1+1 ), 1 )
418 nrm = abs( work( indrv1+jmax ) )
419*
420* Continue for additional iterations after norm reaches
421* stopping criterion.
422*
423 IF( nrm.LT.stpcrt )
424 \$ GO TO 70
425 nrmchk = nrmchk + 1
426 IF( nrmchk.LT.extra+1 )
427 \$ GO TO 70
428*
429 GO TO 130
430*
431* If stopping criterion was not satisfied, update info and
432* store eigenvector number in array ifail.
433*
434 120 CONTINUE
435 info = info + 1
436 ifail( info ) = j
437*
438* Accept iterate as jth eigenvector.
439*
440 130 CONTINUE
441 scl = one / snrm2( blksiz, work( indrv1+1 ), 1 )
442 jmax = isamax( blksiz, work( indrv1+1 ), 1 )
443 IF( work( indrv1+jmax ).LT.zero )
444 \$ scl = -scl
445 CALL sscal( blksiz, scl, work( indrv1+1 ), 1 )
446 140 CONTINUE
447 DO 150 i = 1, n
448 z( i, j ) = czero
449 150 CONTINUE
450 DO 160 i = 1, blksiz
451 z( b1+i-1, j ) = cmplx( work( indrv1+i ), zero )
452 160 CONTINUE
453*
454* Save the shift to check eigenvalue spacing at next
455* iteration.
456*
457 xjm = xj
458*
459 170 CONTINUE
460 180 CONTINUE
461*
462 RETURN
463*
464* End of CSTEIN
465*
466 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
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
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
subroutine slarnv(idist, iseed, n, x)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition slarnv.f:97
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine cstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
CSTEIN
Definition cstein.f:182