LAPACK  3.9.1
LAPACK: Linear Algebra PACKage
zlar1v.f
Go to the documentation of this file.
1 *> \brief \b ZLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the tridiagonal matrix LDLT - λI.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZLAR1V + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlar1v.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlar1v.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlar1v.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLAR1V( N, B1, BN, LAMBDA, D, L, LD, LLD,
22 * PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
23 * R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
24 *
25 * .. Scalar Arguments ..
26 * LOGICAL WANTNC
27 * INTEGER B1, BN, N, NEGCNT, R
28 * DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
29 * $ RQCORR, ZTZ
30 * ..
31 * .. Array Arguments ..
32 * INTEGER ISUPPZ( * )
33 * DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
34 * $ WORK( * )
35 * COMPLEX*16 Z( * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> ZLAR1V computes the (scaled) r-th column of the inverse of
45 *> the sumbmatrix in rows B1 through BN of the tridiagonal matrix
46 *> L D L**T - sigma I. When sigma is close to an eigenvalue, the
47 *> computed vector is an accurate eigenvector. Usually, r corresponds
48 *> to the index where the eigenvector is largest in magnitude.
49 *> The following steps accomplish this computation :
50 *> (a) Stationary qd transform, L D L**T - sigma I = L(+) D(+) L(+)**T,
51 *> (b) Progressive qd transform, L D L**T - sigma I = U(-) D(-) U(-)**T,
52 *> (c) Computation of the diagonal elements of the inverse of
53 *> L D L**T - sigma I by combining the above transforms, and choosing
54 *> r as the index where the diagonal of the inverse is (one of the)
55 *> largest in magnitude.
56 *> (d) Computation of the (scaled) r-th column of the inverse using the
57 *> twisted factorization obtained by combining the top part of the
58 *> the stationary and the bottom part of the progressive transform.
59 *> \endverbatim
60 *
61 * Arguments:
62 * ==========
63 *
64 *> \param[in] N
65 *> \verbatim
66 *> N is INTEGER
67 *> The order of the matrix L D L**T.
68 *> \endverbatim
69 *>
70 *> \param[in] B1
71 *> \verbatim
72 *> B1 is INTEGER
73 *> First index of the submatrix of L D L**T.
74 *> \endverbatim
75 *>
76 *> \param[in] BN
77 *> \verbatim
78 *> BN is INTEGER
79 *> Last index of the submatrix of L D L**T.
80 *> \endverbatim
81 *>
82 *> \param[in] LAMBDA
83 *> \verbatim
84 *> LAMBDA is DOUBLE PRECISION
85 *> The shift. In order to compute an accurate eigenvector,
86 *> LAMBDA should be a good approximation to an eigenvalue
87 *> of L D L**T.
88 *> \endverbatim
89 *>
90 *> \param[in] L
91 *> \verbatim
92 *> L is DOUBLE PRECISION array, dimension (N-1)
93 *> The (n-1) subdiagonal elements of the unit bidiagonal matrix
94 *> L, in elements 1 to N-1.
95 *> \endverbatim
96 *>
97 *> \param[in] D
98 *> \verbatim
99 *> D is DOUBLE PRECISION array, dimension (N)
100 *> The n diagonal elements of the diagonal matrix D.
101 *> \endverbatim
102 *>
103 *> \param[in] LD
104 *> \verbatim
105 *> LD is DOUBLE PRECISION array, dimension (N-1)
106 *> The n-1 elements L(i)*D(i).
107 *> \endverbatim
108 *>
109 *> \param[in] LLD
110 *> \verbatim
111 *> LLD is DOUBLE PRECISION array, dimension (N-1)
112 *> The n-1 elements L(i)*L(i)*D(i).
113 *> \endverbatim
114 *>
115 *> \param[in] PIVMIN
116 *> \verbatim
117 *> PIVMIN is DOUBLE PRECISION
118 *> The minimum pivot in the Sturm sequence.
119 *> \endverbatim
120 *>
121 *> \param[in] GAPTOL
122 *> \verbatim
123 *> GAPTOL is DOUBLE PRECISION
124 *> Tolerance that indicates when eigenvector entries are negligible
125 *> w.r.t. their contribution to the residual.
126 *> \endverbatim
127 *>
128 *> \param[in,out] Z
129 *> \verbatim
130 *> Z is COMPLEX*16 array, dimension (N)
131 *> On input, all entries of Z must be set to 0.
132 *> On output, Z contains the (scaled) r-th column of the
133 *> inverse. The scaling is such that Z(R) equals 1.
134 *> \endverbatim
135 *>
136 *> \param[in] WANTNC
137 *> \verbatim
138 *> WANTNC is LOGICAL
139 *> Specifies whether NEGCNT has to be computed.
140 *> \endverbatim
141 *>
142 *> \param[out] NEGCNT
143 *> \verbatim
144 *> NEGCNT is INTEGER
145 *> If WANTNC is .TRUE. then NEGCNT = the number of pivots < pivmin
146 *> in the matrix factorization L D L**T, and NEGCNT = -1 otherwise.
147 *> \endverbatim
148 *>
149 *> \param[out] ZTZ
150 *> \verbatim
151 *> ZTZ is DOUBLE PRECISION
152 *> The square of the 2-norm of Z.
153 *> \endverbatim
154 *>
155 *> \param[out] MINGMA
156 *> \verbatim
157 *> MINGMA is DOUBLE PRECISION
158 *> The reciprocal of the largest (in magnitude) diagonal
159 *> element of the inverse of L D L**T - sigma I.
160 *> \endverbatim
161 *>
162 *> \param[in,out] R
163 *> \verbatim
164 *> R is INTEGER
165 *> The twist index for the twisted factorization used to
166 *> compute Z.
167 *> On input, 0 <= R <= N. If R is input as 0, R is set to
168 *> the index where (L D L**T - sigma I)^{-1} is largest
169 *> in magnitude. If 1 <= R <= N, R is unchanged.
170 *> On output, R contains the twist index used to compute Z.
171 *> Ideally, R designates the position of the maximum entry in the
172 *> eigenvector.
173 *> \endverbatim
174 *>
175 *> \param[out] ISUPPZ
176 *> \verbatim
177 *> ISUPPZ is INTEGER array, dimension (2)
178 *> The support of the vector in Z, i.e., the vector Z is
179 *> nonzero only in elements ISUPPZ(1) through ISUPPZ( 2 ).
180 *> \endverbatim
181 *>
182 *> \param[out] NRMINV
183 *> \verbatim
184 *> NRMINV is DOUBLE PRECISION
185 *> NRMINV = 1/SQRT( ZTZ )
186 *> \endverbatim
187 *>
188 *> \param[out] RESID
189 *> \verbatim
190 *> RESID is DOUBLE PRECISION
191 *> The residual of the FP vector.
192 *> RESID = ABS( MINGMA )/SQRT( ZTZ )
193 *> \endverbatim
194 *>
195 *> \param[out] RQCORR
196 *> \verbatim
197 *> RQCORR is DOUBLE PRECISION
198 *> The Rayleigh Quotient correction to LAMBDA.
199 *> RQCORR = MINGMA*TMP
200 *> \endverbatim
201 *>
202 *> \param[out] WORK
203 *> \verbatim
204 *> WORK is DOUBLE PRECISION array, dimension (4*N)
205 *> \endverbatim
206 *
207 * Authors:
208 * ========
209 *
210 *> \author Univ. of Tennessee
211 *> \author Univ. of California Berkeley
212 *> \author Univ. of Colorado Denver
213 *> \author NAG Ltd.
214 *
215 *> \ingroup complex16OTHERauxiliary
216 *
217 *> \par Contributors:
218 * ==================
219 *>
220 *> Beresford Parlett, University of California, Berkeley, USA \n
221 *> Jim Demmel, University of California, Berkeley, USA \n
222 *> Inderjit Dhillon, University of Texas, Austin, USA \n
223 *> Osni Marques, LBNL/NERSC, USA \n
224 *> Christof Voemel, University of California, Berkeley, USA
225 *
226 * =====================================================================
227  SUBROUTINE zlar1v( N, B1, BN, LAMBDA, D, L, LD, LLD,
228  $ PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA,
229  $ R, ISUPPZ, NRMINV, RESID, RQCORR, WORK )
230 *
231 * -- LAPACK auxiliary routine --
232 * -- LAPACK is a software package provided by Univ. of Tennessee, --
233 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
234 *
235 * .. Scalar Arguments ..
236  LOGICAL WANTNC
237  INTEGER B1, BN, N, NEGCNT, R
238  DOUBLE PRECISION GAPTOL, LAMBDA, MINGMA, NRMINV, PIVMIN, RESID,
239  $ rqcorr, ztz
240 * ..
241 * .. Array Arguments ..
242  INTEGER ISUPPZ( * )
243  DOUBLE PRECISION D( * ), L( * ), LD( * ), LLD( * ),
244  $ work( * )
245  COMPLEX*16 Z( * )
246 * ..
247 *
248 * =====================================================================
249 *
250 * .. Parameters ..
251  DOUBLE PRECISION ZERO, ONE
252  PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
253  COMPLEX*16 CONE
254  parameter( cone = ( 1.0d0, 0.0d0 ) )
255 
256 * ..
257 * .. Local Scalars ..
258  LOGICAL SAWNAN1, SAWNAN2
259  INTEGER I, INDLPL, INDP, INDS, INDUMN, NEG1, NEG2, R1,
260  $ r2
261  DOUBLE PRECISION DMINUS, DPLUS, EPS, S, TMP
262 * ..
263 * .. External Functions ..
264  LOGICAL DISNAN
265  DOUBLE PRECISION DLAMCH
266  EXTERNAL disnan, dlamch
267 * ..
268 * .. Intrinsic Functions ..
269  INTRINSIC abs, dble
270 * ..
271 * .. Executable Statements ..
272 *
273  eps = dlamch( 'Precision' )
274 
275 
276  IF( r.EQ.0 ) THEN
277  r1 = b1
278  r2 = bn
279  ELSE
280  r1 = r
281  r2 = r
282  END IF
283 
284 * Storage for LPLUS
285  indlpl = 0
286 * Storage for UMINUS
287  indumn = n
288  inds = 2*n + 1
289  indp = 3*n + 1
290 
291  IF( b1.EQ.1 ) THEN
292  work( inds ) = zero
293  ELSE
294  work( inds+b1-1 ) = lld( b1-1 )
295  END IF
296 
297 *
298 * Compute the stationary transform (using the differential form)
299 * until the index R2.
300 *
301  sawnan1 = .false.
302  neg1 = 0
303  s = work( inds+b1-1 ) - lambda
304  DO 50 i = b1, r1 - 1
305  dplus = d( i ) + s
306  work( indlpl+i ) = ld( i ) / dplus
307  IF(dplus.LT.zero) neg1 = neg1 + 1
308  work( inds+i ) = s*work( indlpl+i )*l( i )
309  s = work( inds+i ) - lambda
310  50 CONTINUE
311  sawnan1 = disnan( s )
312  IF( sawnan1 ) GOTO 60
313  DO 51 i = r1, r2 - 1
314  dplus = d( i ) + s
315  work( indlpl+i ) = ld( i ) / dplus
316  work( inds+i ) = s*work( indlpl+i )*l( i )
317  s = work( inds+i ) - lambda
318  51 CONTINUE
319  sawnan1 = disnan( s )
320 *
321  60 CONTINUE
322  IF( sawnan1 ) THEN
323 * Runs a slower version of the above loop if a NaN is detected
324  neg1 = 0
325  s = work( inds+b1-1 ) - lambda
326  DO 70 i = b1, r1 - 1
327  dplus = d( i ) + s
328  IF(abs(dplus).LT.pivmin) dplus = -pivmin
329  work( indlpl+i ) = ld( i ) / dplus
330  IF(dplus.LT.zero) neg1 = neg1 + 1
331  work( inds+i ) = s*work( indlpl+i )*l( i )
332  IF( work( indlpl+i ).EQ.zero )
333  $ work( inds+i ) = lld( i )
334  s = work( inds+i ) - lambda
335  70 CONTINUE
336  DO 71 i = r1, r2 - 1
337  dplus = d( i ) + s
338  IF(abs(dplus).LT.pivmin) dplus = -pivmin
339  work( indlpl+i ) = ld( i ) / dplus
340  work( inds+i ) = s*work( indlpl+i )*l( i )
341  IF( work( indlpl+i ).EQ.zero )
342  $ work( inds+i ) = lld( i )
343  s = work( inds+i ) - lambda
344  71 CONTINUE
345  END IF
346 *
347 * Compute the progressive transform (using the differential form)
348 * until the index R1
349 *
350  sawnan2 = .false.
351  neg2 = 0
352  work( indp+bn-1 ) = d( bn ) - lambda
353  DO 80 i = bn - 1, r1, -1
354  dminus = lld( i ) + work( indp+i )
355  tmp = d( i ) / dminus
356  IF(dminus.LT.zero) neg2 = neg2 + 1
357  work( indumn+i ) = l( i )*tmp
358  work( indp+i-1 ) = work( indp+i )*tmp - lambda
359  80 CONTINUE
360  tmp = work( indp+r1-1 )
361  sawnan2 = disnan( tmp )
362 
363  IF( sawnan2 ) THEN
364 * Runs a slower version of the above loop if a NaN is detected
365  neg2 = 0
366  DO 100 i = bn-1, r1, -1
367  dminus = lld( i ) + work( indp+i )
368  IF(abs(dminus).LT.pivmin) dminus = -pivmin
369  tmp = d( i ) / dminus
370  IF(dminus.LT.zero) neg2 = neg2 + 1
371  work( indumn+i ) = l( i )*tmp
372  work( indp+i-1 ) = work( indp+i )*tmp - lambda
373  IF( tmp.EQ.zero )
374  $ work( indp+i-1 ) = d( i ) - lambda
375  100 CONTINUE
376  END IF
377 *
378 * Find the index (from R1 to R2) of the largest (in magnitude)
379 * diagonal element of the inverse
380 *
381  mingma = work( inds+r1-1 ) + work( indp+r1-1 )
382  IF( mingma.LT.zero ) neg1 = neg1 + 1
383  IF( wantnc ) THEN
384  negcnt = neg1 + neg2
385  ELSE
386  negcnt = -1
387  ENDIF
388  IF( abs(mingma).EQ.zero )
389  $ mingma = eps*work( inds+r1-1 )
390  r = r1
391  DO 110 i = r1, r2 - 1
392  tmp = work( inds+i ) + work( indp+i )
393  IF( tmp.EQ.zero )
394  $ tmp = eps*work( inds+i )
395  IF( abs( tmp ).LE.abs( mingma ) ) THEN
396  mingma = tmp
397  r = i + 1
398  END IF
399  110 CONTINUE
400 *
401 * Compute the FP vector: solve N^T v = e_r
402 *
403  isuppz( 1 ) = b1
404  isuppz( 2 ) = bn
405  z( r ) = cone
406  ztz = one
407 *
408 * Compute the FP vector upwards from R
409 *
410  IF( .NOT.sawnan1 .AND. .NOT.sawnan2 ) THEN
411  DO 210 i = r-1, b1, -1
412  z( i ) = -( work( indlpl+i )*z( i+1 ) )
413  IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
414  $ THEN
415  z( i ) = zero
416  isuppz( 1 ) = i + 1
417  GOTO 220
418  ENDIF
419  ztz = ztz + dble( z( i )*z( i ) )
420  210 CONTINUE
421  220 CONTINUE
422  ELSE
423 * Run slower loop if NaN occurred.
424  DO 230 i = r - 1, b1, -1
425  IF( z( i+1 ).EQ.zero ) THEN
426  z( i ) = -( ld( i+1 ) / ld( i ) )*z( i+2 )
427  ELSE
428  z( i ) = -( work( indlpl+i )*z( i+1 ) )
429  END IF
430  IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
431  $ THEN
432  z( i ) = zero
433  isuppz( 1 ) = i + 1
434  GO TO 240
435  END IF
436  ztz = ztz + dble( z( i )*z( i ) )
437  230 CONTINUE
438  240 CONTINUE
439  ENDIF
440 
441 * Compute the FP vector downwards from R in blocks of size BLKSIZ
442  IF( .NOT.sawnan1 .AND. .NOT.sawnan2 ) THEN
443  DO 250 i = r, bn-1
444  z( i+1 ) = -( work( indumn+i )*z( i ) )
445  IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
446  $ THEN
447  z( i+1 ) = zero
448  isuppz( 2 ) = i
449  GO TO 260
450  END IF
451  ztz = ztz + dble( z( i+1 )*z( i+1 ) )
452  250 CONTINUE
453  260 CONTINUE
454  ELSE
455 * Run slower loop if NaN occurred.
456  DO 270 i = r, bn - 1
457  IF( z( i ).EQ.zero ) THEN
458  z( i+1 ) = -( ld( i-1 ) / ld( i ) )*z( i-1 )
459  ELSE
460  z( i+1 ) = -( work( indumn+i )*z( i ) )
461  END IF
462  IF( (abs(z(i))+abs(z(i+1)))* abs(ld(i)).LT.gaptol )
463  $ THEN
464  z( i+1 ) = zero
465  isuppz( 2 ) = i
466  GO TO 280
467  END IF
468  ztz = ztz + dble( z( i+1 )*z( i+1 ) )
469  270 CONTINUE
470  280 CONTINUE
471  END IF
472 *
473 * Compute quantities for convergence test
474 *
475  tmp = one / ztz
476  nrminv = sqrt( tmp )
477  resid = abs( mingma )*nrminv
478  rqcorr = mingma*tmp
479 *
480 *
481  RETURN
482 *
483 * End of ZLAR1V
484 *
485  END
subroutine zlar1v(N, B1, BN, LAMBDA, D, L, LD, LLD, PIVMIN, GAPTOL, Z, WANTNC, NEGCNT, ZTZ, MINGMA, R, ISUPPZ, NRMINV, RESID, RQCORR, WORK)
ZLAR1V computes the (scaled) r-th column of the inverse of the submatrix in rows b1 through bn of the...
Definition: zlar1v.f:230