LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlasd8()

subroutine dlasd8 ( integer  ICOMPQ,
integer  K,
double precision, dimension( * )  D,
double precision, dimension( * )  Z,
double precision, dimension( * )  VF,
double precision, dimension( * )  VL,
double precision, dimension( * )  DIFL,
double precision, dimension( lddifr, * )  DIFR,
integer  LDDIFR,
double precision, dimension( * )  DSIGMA,
double precision, dimension( * )  WORK,
integer  INFO 
)

DLASD8 finds the square roots of the roots of the secular equation, and stores, for each element in D, the distance to its two nearest poles. Used by sbdsdc.

Download DLASD8 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DLASD8 finds the square roots of the roots of the secular equation,
 as defined by the values in DSIGMA and Z. It makes the appropriate
 calls to DLASD4, and stores, for each  element in D, the distance
 to its two nearest poles (elements in DSIGMA). It also updates
 the arrays VF and VL, the first and last components of all the
 right singular vectors of the original bidiagonal matrix.

 DLASD8 is called from DLASD6.
Parameters
[in]ICOMPQ
          ICOMPQ is INTEGER
          Specifies whether singular vectors are to be computed in
          factored form in the calling routine:
          = 0: Compute singular values only.
          = 1: Compute singular vectors in factored form as well.
[in]K
          K is INTEGER
          The number of terms in the rational function to be solved
          by DLASD4.  K >= 1.
[out]D
          D is DOUBLE PRECISION array, dimension ( K )
          On output, D contains the updated singular values.
[in,out]Z
          Z is DOUBLE PRECISION array, dimension ( K )
          On entry, the first K elements of this array contain the
          components of the deflation-adjusted updating row vector.
          On exit, Z is updated.
[in,out]VF
          VF is DOUBLE PRECISION array, dimension ( K )
          On entry, VF contains  information passed through DBEDE8.
          On exit, VF contains the first K components of the first
          components of all right singular vectors of the bidiagonal
          matrix.
[in,out]VL
          VL is DOUBLE PRECISION array, dimension ( K )
          On entry, VL contains  information passed through DBEDE8.
          On exit, VL contains the first K components of the last
          components of all right singular vectors of the bidiagonal
          matrix.
[out]DIFL
          DIFL is DOUBLE PRECISION array, dimension ( K )
          On exit, DIFL(I) = D(I) - DSIGMA(I).
[out]DIFR
          DIFR is DOUBLE PRECISION array,
                   dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
                   dimension ( K ) if ICOMPQ = 0.
          On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
          defined and will not be referenced.

          If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
          normalizing factors for the right singular vector matrix.
[in]LDDIFR
          LDDIFR is INTEGER
          The leading dimension of DIFR, must be at least K.
[in,out]DSIGMA
          DSIGMA is DOUBLE PRECISION array, dimension ( K )
          On entry, the first K elements of this array contain the old
          roots of the deflated updating problem.  These are the poles
          of the secular equation.
          On exit, the elements of DSIGMA may be very slightly altered
          in value.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (3*K)
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  if INFO = 1, a singular value did not converge
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 164 of file dlasd8.f.

166*
167* -- LAPACK auxiliary routine --
168* -- LAPACK is a software package provided by Univ. of Tennessee, --
169* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
170*
171* .. Scalar Arguments ..
172 INTEGER ICOMPQ, INFO, K, LDDIFR
173* ..
174* .. Array Arguments ..
175 DOUBLE PRECISION D( * ), DIFL( * ), DIFR( LDDIFR, * ),
176 $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
177 $ Z( * )
178* ..
179*
180* =====================================================================
181*
182* .. Parameters ..
183 DOUBLE PRECISION ONE
184 parameter( one = 1.0d+0 )
185* ..
186* .. Local Scalars ..
187 INTEGER I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
188 DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
189* ..
190* .. External Subroutines ..
191 EXTERNAL dcopy, dlascl, dlasd4, dlaset, xerbla
192* ..
193* .. External Functions ..
194 DOUBLE PRECISION DDOT, DLAMC3, DNRM2
195 EXTERNAL ddot, dlamc3, dnrm2
196* ..
197* .. Intrinsic Functions ..
198 INTRINSIC abs, sign, sqrt
199* ..
200* .. Executable Statements ..
201*
202* Test the input parameters.
203*
204 info = 0
205*
206 IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
207 info = -1
208 ELSE IF( k.LT.1 ) THEN
209 info = -2
210 ELSE IF( lddifr.LT.k ) THEN
211 info = -9
212 END IF
213 IF( info.NE.0 ) THEN
214 CALL xerbla( 'DLASD8', -info )
215 RETURN
216 END IF
217*
218* Quick return if possible
219*
220 IF( k.EQ.1 ) THEN
221 d( 1 ) = abs( z( 1 ) )
222 difl( 1 ) = d( 1 )
223 IF( icompq.EQ.1 ) THEN
224 difl( 2 ) = one
225 difr( 1, 2 ) = one
226 END IF
227 RETURN
228 END IF
229*
230* Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
231* be computed with high relative accuracy (barring over/underflow).
232* This is a problem on machines without a guard digit in
233* add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
234* The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
235* which on any of these machines zeros out the bottommost
236* bit of DSIGMA(I) if it is 1; this makes the subsequent
237* subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
238* occurs. On binary machines with a guard digit (almost all
239* machines) it does not change DSIGMA(I) at all. On hexadecimal
240* and decimal machines with a guard digit, it slightly
241* changes the bottommost bits of DSIGMA(I). It does not account
242* for hexadecimal or decimal machines without guard digits
243* (we know of none). We use a subroutine call to compute
244* 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
245* this code.
246*
247 DO 10 i = 1, k
248 dsigma( i ) = dlamc3( dsigma( i ), dsigma( i ) ) - dsigma( i )
249 10 CONTINUE
250*
251* Book keeping.
252*
253 iwk1 = 1
254 iwk2 = iwk1 + k
255 iwk3 = iwk2 + k
256 iwk2i = iwk2 - 1
257 iwk3i = iwk3 - 1
258*
259* Normalize Z.
260*
261 rho = dnrm2( k, z, 1 )
262 CALL dlascl( 'G', 0, 0, rho, one, k, 1, z, k, info )
263 rho = rho*rho
264*
265* Initialize WORK(IWK3).
266*
267 CALL dlaset( 'A', k, 1, one, one, work( iwk3 ), k )
268*
269* Compute the updated singular values, the arrays DIFL, DIFR,
270* and the updated Z.
271*
272 DO 40 j = 1, k
273 CALL dlasd4( k, j, dsigma, z, work( iwk1 ), rho, d( j ),
274 $ work( iwk2 ), info )
275*
276* If the root finder fails, report the convergence failure.
277*
278 IF( info.NE.0 ) THEN
279 RETURN
280 END IF
281 work( iwk3i+j ) = work( iwk3i+j )*work( j )*work( iwk2i+j )
282 difl( j ) = -work( j )
283 difr( j, 1 ) = -work( j+1 )
284 DO 20 i = 1, j - 1
285 work( iwk3i+i ) = work( iwk3i+i )*work( i )*
286 $ work( iwk2i+i ) / ( dsigma( i )-
287 $ dsigma( j ) ) / ( dsigma( i )+
288 $ dsigma( j ) )
289 20 CONTINUE
290 DO 30 i = j + 1, k
291 work( iwk3i+i ) = work( iwk3i+i )*work( i )*
292 $ work( iwk2i+i ) / ( dsigma( i )-
293 $ dsigma( j ) ) / ( dsigma( i )+
294 $ dsigma( j ) )
295 30 CONTINUE
296 40 CONTINUE
297*
298* Compute updated Z.
299*
300 DO 50 i = 1, k
301 z( i ) = sign( sqrt( abs( work( iwk3i+i ) ) ), z( i ) )
302 50 CONTINUE
303*
304* Update VF and VL.
305*
306 DO 80 j = 1, k
307 diflj = difl( j )
308 dj = d( j )
309 dsigj = -dsigma( j )
310 IF( j.LT.k ) THEN
311 difrj = -difr( j, 1 )
312 dsigjp = -dsigma( j+1 )
313 END IF
314 work( j ) = -z( j ) / diflj / ( dsigma( j )+dj )
315 DO 60 i = 1, j - 1
316 work( i ) = z( i ) / ( dlamc3( dsigma( i ), dsigj )-diflj )
317 $ / ( dsigma( i )+dj )
318 60 CONTINUE
319 DO 70 i = j + 1, k
320 work( i ) = z( i ) / ( dlamc3( dsigma( i ), dsigjp )+difrj )
321 $ / ( dsigma( i )+dj )
322 70 CONTINUE
323 temp = dnrm2( k, work, 1 )
324 work( iwk2i+j ) = ddot( k, work, 1, vf, 1 ) / temp
325 work( iwk3i+j ) = ddot( k, work, 1, vl, 1 ) / temp
326 IF( icompq.EQ.1 ) THEN
327 difr( j, 2 ) = temp
328 END IF
329 80 CONTINUE
330*
331 CALL dcopy( k, work( iwk2 ), 1, vf, 1 )
332 CALL dcopy( k, work( iwk3 ), 1, vl, 1 )
333*
334 RETURN
335*
336* End of DLASD8
337*
double precision function dlamc3(A, B)
DLAMC3
Definition: dlamch.f:169
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine dlasd4(N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO)
DLASD4 computes the square root of the i-th updated eigenvalue of a positive symmetric rank-one modif...
Definition: dlasd4.f:153
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:82
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition: dnrm2.f90:89
Here is the call graph for this function:
Here is the caller graph for this function: