ENH: Improving the travis dashboard name
[platform/upstream/lapack.git] / SRC / zstein.f
1 *> \brief \b ZSTEIN
2 *
3 *  =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 *            http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZSTEIN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstein.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstein.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstein.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE ZSTEIN( 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 *       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
31 *       COMPLEX*16         Z( LDZ, * )
32 *       ..
33 *
34 *
35 *> \par Purpose:
36 *  =============
37 *>
38 *> \verbatim
39 *>
40 *> ZSTEIN 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 ZUNMTR or ZUPMTR 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DSTEBZ 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 DSTEBZ 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 DSTEBZ is expected here. )
110 *> \endverbatim
111 *>
112 *> \param[out] Z
113 *> \verbatim
114 *>          Z is COMPLEX*16 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 DOUBLE PRECISION 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 *> \date November 2015
178 *
179 *> \ingroup complex16OTHERcomputational
180 *
181 *  =====================================================================
182       SUBROUTINE ZSTEIN( N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK,
183      $                   IWORK, IFAIL, INFO )
184 *
185 *  -- LAPACK computational routine (version 3.6.0) --
186 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
187 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
188 *     November 2015
189 *
190 *     .. Scalar Arguments ..
191       INTEGER            INFO, LDZ, M, N
192 *     ..
193 *     .. Array Arguments ..
194       INTEGER            IBLOCK( * ), IFAIL( * ), ISPLIT( * ),
195      $                   IWORK( * )
196       DOUBLE PRECISION   D( * ), E( * ), W( * ), WORK( * )
197       COMPLEX*16         Z( LDZ, * )
198 *     ..
199 *
200 * =====================================================================
201 *
202 *     .. Parameters ..
203       COMPLEX*16         CZERO, CONE
204       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
205      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
206       DOUBLE PRECISION   ZERO, ONE, TEN, ODM3, ODM1
207       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1,
208      $                   ODM3 = 1.0D-3, ODM1 = 1.0D-1 )
209       INTEGER            MAXITS, EXTRA
210       PARAMETER          ( MAXITS = 5, EXTRA = 2 )
211 *     ..
212 *     .. Local Scalars ..
213       INTEGER            B1, BLKSIZ, BN, GPIND, I, IINFO, INDRV1,
214      $                   INDRV2, INDRV3, INDRV4, INDRV5, ITS, J, J1,
215      $                   JBLK, JMAX, JR, NBLK, NRMCHK
216       DOUBLE PRECISION   DTPCRT, EPS, EPS1, NRM, ONENRM, ORTOL, PERTOL,
217      $                   SCL, SEP, TOL, XJ, XJM, ZTR
218 *     ..
219 *     .. Local Arrays ..
220       INTEGER            ISEED( 4 )
221 *     ..
222 *     .. External Functions ..
223       INTEGER            IDAMAX
224       DOUBLE PRECISION   DASUM, DLAMCH, DNRM2
225       EXTERNAL           IDAMAX, DASUM, DLAMCH, DNRM2
226 *     ..
227 *     .. External Subroutines ..
228       EXTERNAL           DCOPY, DLAGTF, DLAGTS, DLARNV, DSCAL, XERBLA
229 *     ..
230 *     .. Intrinsic Functions ..
231       INTRINSIC          ABS, DBLE, DCMPLX, MAX, SQRT
232 *     ..
233 *     .. Executable Statements ..
234 *
235 *     Test the input parameters.
236 *
237       INFO = 0
238       DO 10 I = 1, M
239          IFAIL( I ) = 0
240    10 CONTINUE
241 *
242       IF( N.LT.0 ) THEN
243          INFO = -1
244       ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
245          INFO = -4
246       ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN
247          INFO = -9
248       ELSE
249          DO 20 J = 2, M
250             IF( IBLOCK( J ).LT.IBLOCK( J-1 ) ) THEN
251                INFO = -6
252                GO TO 30
253             END IF
254             IF( IBLOCK( J ).EQ.IBLOCK( J-1 ) .AND. W( J ).LT.W( J-1 ) )
255      $           THEN
256                INFO = -5
257                GO TO 30
258             END IF
259    20    CONTINUE
260    30    CONTINUE
261       END IF
262 *
263       IF( INFO.NE.0 ) THEN
264          CALL XERBLA( 'ZSTEIN', -INFO )
265          RETURN
266       END IF
267 *
268 *     Quick return if possible
269 *
270       IF( N.EQ.0 .OR. M.EQ.0 ) THEN
271          RETURN
272       ELSE IF( N.EQ.1 ) THEN
273          Z( 1, 1 ) = CONE
274          RETURN
275       END IF
276 *
277 *     Get machine constants.
278 *
279       EPS = DLAMCH( 'Precision' )
280 *
281 *     Initialize seed for random number generator DLARNV.
282 *
283       DO 40 I = 1, 4
284          ISEED( I ) = 1
285    40 CONTINUE
286 *
287 *     Initialize pointers.
288 *
289       INDRV1 = 0
290       INDRV2 = INDRV1 + N
291       INDRV3 = INDRV2 + N
292       INDRV4 = INDRV3 + N
293       INDRV5 = INDRV4 + N
294 *
295 *     Compute eigenvectors of matrix blocks.
296 *
297       J1 = 1
298       DO 180 NBLK = 1, IBLOCK( M )
299 *
300 *        Find starting and ending indices of block nblk.
301 *
302          IF( NBLK.EQ.1 ) THEN
303             B1 = 1
304          ELSE
305             B1 = ISPLIT( NBLK-1 ) + 1
306          END IF
307          BN = ISPLIT( NBLK )
308          BLKSIZ = BN - B1 + 1
309          IF( BLKSIZ.EQ.1 )
310      $      GO TO 60
311          GPIND = J1
312 *
313 *        Compute reorthogonalization criterion and stopping criterion.
314 *
315          ONENRM = ABS( D( B1 ) ) + ABS( E( B1 ) )
316          ONENRM = MAX( ONENRM, ABS( D( BN ) )+ABS( E( BN-1 ) ) )
317          DO 50 I = B1 + 1, BN - 1
318             ONENRM = MAX( ONENRM, ABS( D( I ) )+ABS( E( I-1 ) )+
319      $               ABS( E( I ) ) )
320    50    CONTINUE
321          ORTOL = ODM3*ONENRM
322 *
323          DTPCRT = SQRT( ODM1 / BLKSIZ )
324 *
325 *        Loop through eigenvalues of block nblk.
326 *
327    60    CONTINUE
328          JBLK = 0
329          DO 170 J = J1, M
330             IF( IBLOCK( J ).NE.NBLK ) THEN
331                J1 = J
332                GO TO 180
333             END IF
334             JBLK = JBLK + 1
335             XJ = W( J )
336 *
337 *           Skip all the work if the block size is one.
338 *
339             IF( BLKSIZ.EQ.1 ) THEN
340                WORK( INDRV1+1 ) = ONE
341                GO TO 140
342             END IF
343 *
344 *           If eigenvalues j and j-1 are too close, add a relatively
345 *           small perturbation.
346 *
347             IF( JBLK.GT.1 ) THEN
348                EPS1 = ABS( EPS*XJ )
349                PERTOL = TEN*EPS1
350                SEP = XJ - XJM
351                IF( SEP.LT.PERTOL )
352      $            XJ = XJM + PERTOL
353             END IF
354 *
355             ITS = 0
356             NRMCHK = 0
357 *
358 *           Get random starting vector.
359 *
360             CALL DLARNV( 2, ISEED, BLKSIZ, WORK( INDRV1+1 ) )
361 *
362 *           Copy the matrix T so it won't be destroyed in factorization.
363 *
364             CALL DCOPY( BLKSIZ, D( B1 ), 1, WORK( INDRV4+1 ), 1 )
365             CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV2+2 ), 1 )
366             CALL DCOPY( BLKSIZ-1, E( B1 ), 1, WORK( INDRV3+1 ), 1 )
367 *
368 *           Compute LU factors with partial pivoting  ( PT = LU )
369 *
370             TOL = ZERO
371             CALL DLAGTF( BLKSIZ, WORK( INDRV4+1 ), XJ, WORK( INDRV2+2 ),
372      $                   WORK( INDRV3+1 ), TOL, WORK( INDRV5+1 ), IWORK,
373      $                   IINFO )
374 *
375 *           Update iteration count.
376 *
377    70       CONTINUE
378             ITS = ITS + 1
379             IF( ITS.GT.MAXITS )
380      $         GO TO 120
381 *
382 *           Normalize and scale the righthand side vector Pb.
383 *
384             JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
385             SCL = BLKSIZ*ONENRM*MAX( EPS,
386      $            ABS( WORK( INDRV4+BLKSIZ ) ) ) /
387      $            ABS( WORK( INDRV1+JMAX ) )
388             CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
389 *
390 *           Solve the system LU = Pb.
391 *
392             CALL DLAGTS( -1, BLKSIZ, WORK( INDRV4+1 ), WORK( INDRV2+2 ),
393      $                   WORK( INDRV3+1 ), WORK( INDRV5+1 ), IWORK,
394      $                   WORK( INDRV1+1 ), TOL, IINFO )
395 *
396 *           Reorthogonalize by modified Gram-Schmidt if eigenvalues are
397 *           close enough.
398 *
399             IF( JBLK.EQ.1 )
400      $         GO TO 110
401             IF( ABS( XJ-XJM ).GT.ORTOL )
402      $         GPIND = J
403             IF( GPIND.NE.J ) THEN
404                DO 100 I = GPIND, J - 1
405                   ZTR = ZERO
406                   DO 80 JR = 1, BLKSIZ
407                      ZTR = ZTR + WORK( INDRV1+JR )*
408      $                     DBLE( Z( B1-1+JR, I ) )
409    80             CONTINUE
410                   DO 90 JR = 1, BLKSIZ
411                      WORK( INDRV1+JR ) = WORK( INDRV1+JR ) -
412      $                                   ZTR*DBLE( Z( B1-1+JR, I ) )
413    90             CONTINUE
414   100          CONTINUE
415             END IF
416 *
417 *           Check the infinity norm of the iterate.
418 *
419   110       CONTINUE
420             JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
421             NRM = ABS( WORK( INDRV1+JMAX ) )
422 *
423 *           Continue for additional iterations after norm reaches
424 *           stopping criterion.
425 *
426             IF( NRM.LT.DTPCRT )
427      $         GO TO 70
428             NRMCHK = NRMCHK + 1
429             IF( NRMCHK.LT.EXTRA+1 )
430      $         GO TO 70
431 *
432             GO TO 130
433 *
434 *           If stopping criterion was not satisfied, update info and
435 *           store eigenvector number in array ifail.
436 *
437   120       CONTINUE
438             INFO = INFO + 1
439             IFAIL( INFO ) = J
440 *
441 *           Accept iterate as jth eigenvector.
442 *
443   130       CONTINUE
444             SCL = ONE / DNRM2( BLKSIZ, WORK( INDRV1+1 ), 1 )
445             JMAX = IDAMAX( BLKSIZ, WORK( INDRV1+1 ), 1 )
446             IF( WORK( INDRV1+JMAX ).LT.ZERO )
447      $         SCL = -SCL
448             CALL DSCAL( BLKSIZ, SCL, WORK( INDRV1+1 ), 1 )
449   140       CONTINUE
450             DO 150 I = 1, N
451                Z( I, J ) = CZERO
452   150       CONTINUE
453             DO 160 I = 1, BLKSIZ
454                Z( B1+I-1, J ) = DCMPLX( WORK( INDRV1+I ), ZERO )
455   160       CONTINUE
456 *
457 *           Save the shift to check eigenvalue spacing at next
458 *           iteration.
459 *
460             XJM = XJ
461 *
462   170    CONTINUE
463   180 CONTINUE
464 *
465       RETURN
466 *
467 *     End of ZSTEIN
468 *
469       END