OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dstein2.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine dstein2 (n, d, e, m, w, iblock, isplit, orfac, z, ldz, work, iwork, ifail, info)

Function/Subroutine Documentation

◆ dstein2()

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

Definition at line 3 of file dstein2.f.

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*
dlarnv
subroutine dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition dlarnv.f:97
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 ma...
Definition dlagts.f:161
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
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
ddot
double precision function ddot(n, dx, incx, dy, incy)
DDOT
Definition ddot.f:82
dasum
double precision function dasum(n, dx, incx)
DASUM
Definition dasum.f:71
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
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
max
#define max(a, b)
Definition macros.h:21