LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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 at least 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.
Date
November 2015
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 168 of file dlasd8.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: