1      SUBROUTINE DLAR1V( N, B1, BN, SIGMA, D, L, LD, LLD, GERSCH, Z,
2     $                   ZTZ, MINGMA, R, ISUPPZ, WORK )
3*
4*  -- LAPACK auxiliary routine (version 3.0) --
5*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6*     Courant Institute, Argonne National Lab, and Rice University
7*     June 30, 1999
8*
9*     .. Scalar Arguments ..
10      INTEGER            B1, BN, N, R
11      DOUBLE PRECISION   MINGMA, SIGMA, ZTZ
12*     ..
13*     .. Array Arguments ..
14      INTEGER            ISUPPZ( * )
15      DOUBLE PRECISION   D( * ), GERSCH( * ), L( * ), LD( * ), LLD( * ),
16     $                   WORK( * ), Z( * )
17*     ..
18*
19*  Purpose
20*  =======
21*
22*  DLAR1V computes the (scaled) r-th column of the inverse of
23*  the sumbmatrix in rows B1 through BN of the tridiagonal matrix
24*  L D L^T - sigma I. The following steps accomplish this computation :
25*  (a) Stationary qd transform,  L D L^T - sigma I = L(+) D(+) L(+)^T,
26*  (b) Progressive qd transform, L D L^T - sigma I = U(-) D(-) U(-)^T,
27*  (c) Computation of the diagonal elements of the inverse of
28*      L D L^T - sigma I by combining the above transforms, and choosing
29*      r as the index where the diagonal of the inverse is (one of the)
30*      largest in magnitude.
31*  (d) Computation of the (scaled) r-th column of the inverse using the
32*      twisted factorization obtained by combining the top part of the
33*      the stationary and the bottom part of the progressive transform.
34*
35*  Arguments
36*  =========
37*
38*  N        (input) INTEGER
39*           The order of the matrix L D L^T.
40*
41*  B1       (input) INTEGER
42*           First index of the submatrix of L D L^T.
43*
44*  BN       (input) INTEGER
45*           Last index of the submatrix of L D L^T.
46*
47*  SIGMA    (input) DOUBLE PRECISION
48*           The shift. Initially, when R = 0, SIGMA should be a good
49*           approximation to an eigenvalue of L D L^T.
50*
51*  L        (input) DOUBLE PRECISION array, dimension (N-1)
52*           The (n-1) subdiagonal elements of the unit bidiagonal matrix
53*           L, in elements 1 to N-1.
54*
55*  D        (input) DOUBLE PRECISION array, dimension (N)
56*           The n diagonal elements of the diagonal matrix D.
57*
58*  LD       (input) DOUBLE PRECISION array, dimension (N-1)
59*           The n-1 elements L(i)*D(i).
60*
61*  LLD      (input) DOUBLE PRECISION array, dimension (N-1)
62*           The n-1 elements L(i)*L(i)*D(i).
63*
64*  GERSCH   (input) DOUBLE PRECISION array, dimension (2*N)
65*           The n Gerschgorin intervals. These are used to restrict
66*           the initial search for R, when R is input as 0.
67*
68*  Z        (output) DOUBLE PRECISION array, dimension (N)
69*           The (scaled) r-th column of the inverse. Z(R) is returned
70*           to be 1.
71*
72*  ZTZ      (output) DOUBLE PRECISION
73*           The square of the norm of Z.
74*
75*  MINGMA   (output) DOUBLE PRECISION
76*           The reciprocal of the largest (in magnitude) diagonal
77*           element of the inverse of L D L^T - sigma I.
78*
79*  R        (input/output) INTEGER
80*           Initially, R should be input to be 0 and is then output as
81*           the index where the diagonal element of the inverse is
82*           largest in magnitude. In later iterations, this same value
83*           of R should be input.
84*
85*  ISUPPZ   (output) INTEGER array, dimension (2)
86*           The support of the vector in Z, i.e., the vector Z is
87*           nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
88*
89*  WORK     (workspace) DOUBLE PRECISION array, dimension (4*N)
90*
91*  Further Details
92*  ===============
93*
94*  Based on contributions by
95*     Inderjit Dhillon, IBM Almaden, USA
96*     Osni Marques, LBNL/NERSC, USA
97*
98*  =====================================================================
99*
100*     .. Parameters ..
101      INTEGER            BLKSIZ
102      PARAMETER          ( BLKSIZ = 32 )
103      DOUBLE PRECISION   ZERO, ONE
104      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
105*     ..
106*     .. Local Scalars ..
107      LOGICAL            SAWNAN
108      INTEGER            FROM, I, INDP, INDS, INDUMN, J, R1, R2, TO
109      DOUBLE PRECISION   DMINUS, DPLUS, EPS, S, TMP
110*     ..
111*     .. External Functions ..
112      DOUBLE PRECISION   DLAMCH
113      EXTERNAL           DLAMCH
114*     ..
115*     .. Intrinsic Functions ..
116      INTRINSIC          ABS, MAX, MIN
117*     ..
118*     .. Executable Statements ..
119*
120      EPS = DLAMCH( 'Precision' )
121      IF( R.EQ.0 ) THEN
122*
123*        Eliminate the top and bottom indices from the possible values
124*        of R where the desired eigenvector is largest in magnitude.
125*
126         R1 = B1
127         DO 10 I = B1, BN
128            IF( SIGMA.GE.GERSCH( 2*I-1 ) .OR. SIGMA.LE.GERSCH( 2*I ) )
129     $           THEN
130               R1 = I
131               GO TO 20
132            END IF
133   10    CONTINUE
134   20    CONTINUE
135         R2 = BN
136         DO 30 I = BN, B1, -1
137            IF( SIGMA.GE.GERSCH( 2*I-1 ) .OR. SIGMA.LE.GERSCH( 2*I ) )
138     $           THEN
139               R2 = I
140               GO TO 40
141            END IF
142   30    CONTINUE
143   40    CONTINUE
144      ELSE
145         R1 = R
146         R2 = R
147      END IF
148*
149      INDUMN = N
150      INDS = 2*N + 1
151      INDP = 3*N + 1
152      SAWNAN = .FALSE.
153*
154*     Compute the stationary transform (using the differential form)
155*     untill the index R2
156*
157      IF( B1.EQ.1 ) THEN
158         WORK( INDS ) = ZERO
159      ELSE
160         WORK( INDS ) = LLD( B1-1 )
161      END IF
162      S = WORK( INDS ) - SIGMA
163      DO 50 I = B1, R2 - 1
164         DPLUS = D( I ) + S
165         WORK( I ) = LD( I ) / DPLUS
166         WORK( INDS+I ) = S*WORK( I )*L( I )
167         S = WORK( INDS+I ) - SIGMA
168   50 CONTINUE
169*
170      IF( .NOT.( S.GT.ZERO .OR. S.LT.ONE ) ) THEN
171*
172*        Run a slower version of the above loop if a NaN is detected
173*
174         SAWNAN = .TRUE.
175         J = B1 + 1
176   60    CONTINUE
177         IF( WORK( INDS+J ).GT.ZERO .OR. WORK( INDS+J ).LT.ONE ) THEN
178            J = J + 1
179            GO TO 60
180         END IF
181         WORK( INDS+J ) = LLD( J )
182         S = WORK( INDS+J ) - SIGMA
183         DO 70 I = J + 1, R2 - 1
184            DPLUS = D( I ) + S
185            WORK( I ) = LD( I ) / DPLUS
186            IF( WORK( I ).EQ.ZERO ) THEN
187               WORK( INDS+I ) = LLD( I )
188            ELSE
189               WORK( INDS+I ) = S*WORK( I )*L( I )
190            END IF
191            S = WORK( INDS+I ) - SIGMA
192   70    CONTINUE
193      END IF
194      WORK( INDP+BN-1 ) = D( BN ) - SIGMA
195      DO 80 I = BN - 1, R1, -1
196         DMINUS = LLD( I ) + WORK( INDP+I )
197         TMP = D( I ) / DMINUS
198         WORK( INDUMN+I ) = L( I )*TMP
199         WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - SIGMA
200   80 CONTINUE
201      TMP = WORK( INDP+R1-1 )
202      IF( .NOT.( TMP.GT.ZERO .OR. TMP.LT.ONE ) ) THEN
203*
204*        Run a slower version of the above loop if a NaN is detected
205*
206         SAWNAN = .TRUE.
207         J = BN - 3
208   90    CONTINUE
209         IF( WORK( INDP+J ).GT.ZERO .OR. WORK( INDP+J ).LT.ONE ) THEN
210            J = J - 1
211            GO TO 90
212         END IF
213         WORK( INDP+J ) = D( J+1 ) - SIGMA
214         DO 100 I = J, R1, -1
215            DMINUS = LLD( I ) + WORK( INDP+I )
216            TMP = D( I ) / DMINUS
217            WORK( INDUMN+I ) = L( I )*TMP
218            IF( TMP.EQ.ZERO ) THEN
219               WORK( INDP+I-1 ) = D( I ) - SIGMA
220            ELSE
221               WORK( INDP+I-1 ) = WORK( INDP+I )*TMP - SIGMA
222            END IF
223  100    CONTINUE
224      END IF
225*
226*     Find the index (from R1 to R2) of the largest (in magnitude)
227*     diagonal element of the inverse
228*
229      MINGMA = WORK( INDS+R1-1 ) + WORK( INDP+R1-1 )
230      IF( MINGMA.EQ.ZERO )
231     $   MINGMA = EPS*WORK( INDS+R1-1 )
232      R = R1
233      DO 110 I = R1, R2 - 1
234         TMP = WORK( INDS+I ) + WORK( INDP+I )
235         IF( TMP.EQ.ZERO )
236     $      TMP = EPS*WORK( INDS+I )
237         IF( ABS( TMP ).LT.ABS( MINGMA ) ) THEN
238            MINGMA = TMP
239            R = I + 1
240         END IF
241  110 CONTINUE
242*
243*     Compute the (scaled) r-th column of the inverse
244*
245      ISUPPZ( 1 ) = B1
246      ISUPPZ( 2 ) = BN
247      Z( R ) = ONE
248      ZTZ = ONE
249      IF( .NOT.SAWNAN ) THEN
250         FROM = R - 1
251         TO = MAX( R-BLKSIZ, B1 )
252  120    CONTINUE
253         IF( FROM.GE.B1 ) THEN
254            DO 130 I = FROM, TO, -1
255               Z( I ) = -( WORK( I )*Z( I+1 ) )
256               ZTZ = ZTZ + Z( I )*Z( I )
257  130       CONTINUE
258            IF( ABS( Z( TO ) ).LE.EPS .AND. ABS( Z( TO+1 ) ).LE.EPS )
259     $           THEN
260               ISUPPZ( 1 ) = TO + 2
261            ELSE
262               FROM = TO - 1
263               TO = MAX( TO-BLKSIZ, B1 )
264               GO TO 120
265            END IF
266         END IF
267         FROM = R + 1
268         TO = MIN( R+BLKSIZ, BN )
269  140    CONTINUE
270         IF( FROM.LE.BN ) THEN
271            DO 150 I = FROM, TO
272               Z( I ) = -( WORK( INDUMN+I-1 )*Z( I-1 ) )
273               ZTZ = ZTZ + Z( I )*Z( I )
274  150       CONTINUE
275            IF( ABS( Z( TO ) ).LE.EPS .AND. ABS( Z( TO-1 ) ).LE.EPS )
276     $           THEN
277               ISUPPZ( 2 ) = TO - 2
278            ELSE
279               FROM = TO + 1
280               TO = MIN( TO+BLKSIZ, BN )
281               GO TO 140
282            END IF
283         END IF
284      ELSE
285         DO 160 I = R - 1, B1, -1
286            IF( Z( I+1 ).EQ.ZERO ) THEN
287               Z( I ) = -( LD( I+1 ) / LD( I ) )*Z( I+2 )
288            ELSE IF( ABS( Z( I+1 ) ).LE.EPS .AND. ABS( Z( I+2 ) ).LE.
289     $               EPS ) THEN
290               ISUPPZ( 1 ) = I + 3
291               GO TO 170
292            ELSE
293               Z( I ) = -( WORK( I )*Z( I+1 ) )
294            END IF
295            ZTZ = ZTZ + Z( I )*Z( I )
296  160    CONTINUE
297  170    CONTINUE
298         DO 180 I = R, BN - 1
299            IF( Z( I ).EQ.ZERO ) THEN
300               Z( I+1 ) = -( LD( I-1 ) / LD( I ) )*Z( I-1 )
301            ELSE IF( ABS( Z( I ) ).LE.EPS .AND. ABS( Z( I-1 ) ).LE.EPS )
302     $                THEN
303               ISUPPZ( 2 ) = I - 2
304               GO TO 190
305            ELSE
306               Z( I+1 ) = -( WORK( INDUMN+I )*Z( I ) )
307            END IF
308            ZTZ = ZTZ + Z( I+1 )*Z( I+1 )
309  180    CONTINUE
310  190    CONTINUE
311      END IF
312      DO 200 I = B1, ISUPPZ( 1 ) - 3
313         Z( I ) = ZERO
314  200 CONTINUE
315      DO 210 I = ISUPPZ( 2 ) + 3, BN
316         Z( I ) = ZERO
317  210 CONTINUE
318*
319      RETURN
320*
321*     End of DLAR1V
322*
323      END
324