LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ dlasd7()

subroutine dlasd7 ( integer  ICOMPQ,
integer  NL,
integer  NR,
integer  SQRE,
integer  K,
double precision, dimension( * )  D,
double precision, dimension( * )  Z,
double precision, dimension( * )  ZW,
double precision, dimension( * )  VF,
double precision, dimension( * )  VFW,
double precision, dimension( * )  VL,
double precision, dimension( * )  VLW,
double precision  ALPHA,
double precision  BETA,
double precision, dimension( * )  DSIGMA,
integer, dimension( * )  IDX,
integer, dimension( * )  IDXP,
integer, dimension( * )  IDXQ,
integer, dimension( * )  PERM,
integer  GIVPTR,
integer, dimension( ldgcol, * )  GIVCOL,
integer  LDGCOL,
double precision, dimension( ldgnum, * )  GIVNUM,
integer  LDGNUM,
double precision  C,
double precision  S,
integer  INFO 
)

DLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem. Used by sbdsdc.

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

Purpose:
 DLASD7 merges the two sets of singular values together into a single
 sorted set. Then it tries to deflate the size of the problem. There
 are two ways in which deflation can occur:  when two or more singular
 values are close together or if there is a tiny entry in the Z
 vector. For each such occurrence the order of the related
 secular equation problem is reduced by one.

 DLASD7 is called from DLASD6.
Parameters
[in]ICOMPQ
          ICOMPQ is INTEGER
          Specifies whether singular vectors are to be computed
          in compact form, as follows:
          = 0: Compute singular values only.
          = 1: Compute singular vectors of upper
               bidiagonal matrix in compact form.
[in]NL
          NL is INTEGER
         The row dimension of the upper block. NL >= 1.
[in]NR
          NR is INTEGER
         The row dimension of the lower block. NR >= 1.
[in]SQRE
          SQRE is INTEGER
         = 0: the lower block is an NR-by-NR square matrix.
         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.

         The bidiagonal matrix has
         N = NL + NR + 1 rows and
         M = N + SQRE >= N columns.
[out]K
          K is INTEGER
         Contains the dimension of the non-deflated matrix, this is
         the order of the related secular equation. 1 <= K <=N.
[in,out]D
          D is DOUBLE PRECISION array, dimension ( N )
         On entry D contains the singular values of the two submatrices
         to be combined. On exit D contains the trailing (N-K) updated
         singular values (those which were deflated) sorted into
         increasing order.
[out]Z
          Z is DOUBLE PRECISION array, dimension ( M )
         On exit Z contains the updating row vector in the secular
         equation.
[out]ZW
          ZW is DOUBLE PRECISION array, dimension ( M )
         Workspace for Z.
[in,out]VF
          VF is DOUBLE PRECISION array, dimension ( M )
         On entry, VF(1:NL+1) contains the first components of all
         right singular vectors of the upper block; and VF(NL+2:M)
         contains the first components of all right singular vectors
         of the lower block. On exit, VF contains the first components
         of all right singular vectors of the bidiagonal matrix.
[out]VFW
          VFW is DOUBLE PRECISION array, dimension ( M )
         Workspace for VF.
[in,out]VL
          VL is DOUBLE PRECISION array, dimension ( M )
         On entry, VL(1:NL+1) contains the  last components of all
         right singular vectors of the upper block; and VL(NL+2:M)
         contains the last components of all right singular vectors
         of the lower block. On exit, VL contains the last components
         of all right singular vectors of the bidiagonal matrix.
[out]VLW
          VLW is DOUBLE PRECISION array, dimension ( M )
         Workspace for VL.
[in]ALPHA
          ALPHA is DOUBLE PRECISION
         Contains the diagonal element associated with the added row.
[in]BETA
          BETA is DOUBLE PRECISION
         Contains the off-diagonal element associated with the added
         row.
[out]DSIGMA
          DSIGMA is DOUBLE PRECISION array, dimension ( N )
         Contains a copy of the diagonal elements (K-1 singular values
         and one zero) in the secular equation.
[out]IDX
          IDX is INTEGER array, dimension ( N )
         This will contain the permutation used to sort the contents of
         D into ascending order.
[out]IDXP
          IDXP is INTEGER array, dimension ( N )
         This will contain the permutation used to place deflated
         values of D at the end of the array. On output IDXP(2:K)
         points to the nondeflated D-values and IDXP(K+1:N)
         points to the deflated singular values.
[in]IDXQ
          IDXQ is INTEGER array, dimension ( N )
         This contains the permutation which separately sorts the two
         sub-problems in D into ascending order.  Note that entries in
         the first half of this permutation must first be moved one
         position backward; and entries in the second half
         must first have NL+1 added to their values.
[out]PERM
          PERM is INTEGER array, dimension ( N )
         The permutations (from deflation and sorting) to be applied
         to each singular block. Not referenced if ICOMPQ = 0.
[out]GIVPTR
          GIVPTR is INTEGER
         The number of Givens rotations which took place in this
         subproblem. Not referenced if ICOMPQ = 0.
[out]GIVCOL
          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
         Each pair of numbers indicates a pair of columns to take place
         in a Givens rotation. Not referenced if ICOMPQ = 0.
[in]LDGCOL
          LDGCOL is INTEGER
         The leading dimension of GIVCOL, must be at least N.
[out]GIVNUM
          GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
         Each number indicates the C or S value to be used in the
         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
[in]LDGNUM
          LDGNUM is INTEGER
         The leading dimension of GIVNUM, must be at least N.
[out]C
          C is DOUBLE PRECISION
         C contains garbage if SQRE =0 and the C-value of a Givens
         rotation related to the right null space if SQRE = 1.
[out]S
          S is DOUBLE PRECISION
         S contains garbage if SQRE =0 and the S-value of a Givens
         rotation related to the right null space if SQRE = 1.
[out]INFO
          INFO is INTEGER
         = 0:  successful exit.
         < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
December 2016
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 282 of file dlasd7.f.

282 *
283 * -- LAPACK auxiliary routine (version 3.7.0) --
284 * -- LAPACK is a software package provided by Univ. of Tennessee, --
285 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
286 * December 2016
287 *
288 * .. Scalar Arguments ..
289  INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
290  $ NR, SQRE
291  DOUBLE PRECISION ALPHA, BETA, C, S
292 * ..
293 * .. Array Arguments ..
294  INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
295  $ IDXQ( * ), PERM( * )
296  DOUBLE PRECISION D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
297  $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
298  $ ZW( * )
299 * ..
300 *
301 * =====================================================================
302 *
303 * .. Parameters ..
304  DOUBLE PRECISION ZERO, ONE, TWO, EIGHT
305  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
306  $ eight = 8.0d+0 )
307 * ..
308 * .. Local Scalars ..
309 *
310  INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
311  $ NLP1, NLP2
312  DOUBLE PRECISION EPS, HLFTOL, TAU, TOL, Z1
313 * ..
314 * .. External Subroutines ..
315  EXTERNAL dcopy, dlamrg, drot, xerbla
316 * ..
317 * .. External Functions ..
318  DOUBLE PRECISION DLAMCH, DLAPY2
319  EXTERNAL dlamch, dlapy2
320 * ..
321 * .. Intrinsic Functions ..
322  INTRINSIC abs, max
323 * ..
324 * .. Executable Statements ..
325 *
326 * Test the input parameters.
327 *
328  info = 0
329  n = nl + nr + 1
330  m = n + sqre
331 *
332  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
333  info = -1
334  ELSE IF( nl.LT.1 ) THEN
335  info = -2
336  ELSE IF( nr.LT.1 ) THEN
337  info = -3
338  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
339  info = -4
340  ELSE IF( ldgcol.LT.n ) THEN
341  info = -22
342  ELSE IF( ldgnum.LT.n ) THEN
343  info = -24
344  END IF
345  IF( info.NE.0 ) THEN
346  CALL xerbla( 'DLASD7', -info )
347  RETURN
348  END IF
349 *
350  nlp1 = nl + 1
351  nlp2 = nl + 2
352  IF( icompq.EQ.1 ) THEN
353  givptr = 0
354  END IF
355 *
356 * Generate the first part of the vector Z and move the singular
357 * values in the first part of D one position backward.
358 *
359  z1 = alpha*vl( nlp1 )
360  vl( nlp1 ) = zero
361  tau = vf( nlp1 )
362  DO 10 i = nl, 1, -1
363  z( i+1 ) = alpha*vl( i )
364  vl( i ) = zero
365  vf( i+1 ) = vf( i )
366  d( i+1 ) = d( i )
367  idxq( i+1 ) = idxq( i ) + 1
368  10 CONTINUE
369  vf( 1 ) = tau
370 *
371 * Generate the second part of the vector Z.
372 *
373  DO 20 i = nlp2, m
374  z( i ) = beta*vf( i )
375  vf( i ) = zero
376  20 CONTINUE
377 *
378 * Sort the singular values into increasing order
379 *
380  DO 30 i = nlp2, n
381  idxq( i ) = idxq( i ) + nlp1
382  30 CONTINUE
383 *
384 * DSIGMA, IDXC, IDXC, and ZW are used as storage space.
385 *
386  DO 40 i = 2, n
387  dsigma( i ) = d( idxq( i ) )
388  zw( i ) = z( idxq( i ) )
389  vfw( i ) = vf( idxq( i ) )
390  vlw( i ) = vl( idxq( i ) )
391  40 CONTINUE
392 *
393  CALL dlamrg( nl, nr, dsigma( 2 ), 1, 1, idx( 2 ) )
394 *
395  DO 50 i = 2, n
396  idxi = 1 + idx( i )
397  d( i ) = dsigma( idxi )
398  z( i ) = zw( idxi )
399  vf( i ) = vfw( idxi )
400  vl( i ) = vlw( idxi )
401  50 CONTINUE
402 *
403 * Calculate the allowable deflation tolerance
404 *
405  eps = dlamch( 'Epsilon' )
406  tol = max( abs( alpha ), abs( beta ) )
407  tol = eight*eight*eps*max( abs( d( n ) ), tol )
408 *
409 * There are 2 kinds of deflation -- first a value in the z-vector
410 * is small, second two (or more) singular values are very close
411 * together (their difference is small).
412 *
413 * If the value in the z-vector is small, we simply permute the
414 * array so that the corresponding singular value is moved to the
415 * end.
416 *
417 * If two values in the D-vector are close, we perform a two-sided
418 * rotation designed to make one of the corresponding z-vector
419 * entries zero, and then permute the array so that the deflated
420 * singular value is moved to the end.
421 *
422 * If there are multiple singular values then the problem deflates.
423 * Here the number of equal singular values are found. As each equal
424 * singular value is found, an elementary reflector is computed to
425 * rotate the corresponding singular subspace so that the
426 * corresponding components of Z are zero in this new basis.
427 *
428  k = 1
429  k2 = n + 1
430  DO 60 j = 2, n
431  IF( abs( z( j ) ).LE.tol ) THEN
432 *
433 * Deflate due to small z component.
434 *
435  k2 = k2 - 1
436  idxp( k2 ) = j
437  IF( j.EQ.n )
438  $ GO TO 100
439  ELSE
440  jprev = j
441  GO TO 70
442  END IF
443  60 CONTINUE
444  70 CONTINUE
445  j = jprev
446  80 CONTINUE
447  j = j + 1
448  IF( j.GT.n )
449  $ GO TO 90
450  IF( abs( z( j ) ).LE.tol ) THEN
451 *
452 * Deflate due to small z component.
453 *
454  k2 = k2 - 1
455  idxp( k2 ) = j
456  ELSE
457 *
458 * Check if singular values are close enough to allow deflation.
459 *
460  IF( abs( d( j )-d( jprev ) ).LE.tol ) THEN
461 *
462 * Deflation is possible.
463 *
464  s = z( jprev )
465  c = z( j )
466 *
467 * Find sqrt(a**2+b**2) without overflow or
468 * destructive underflow.
469 *
470  tau = dlapy2( c, s )
471  z( j ) = tau
472  z( jprev ) = zero
473  c = c / tau
474  s = -s / tau
475 *
476 * Record the appropriate Givens rotation
477 *
478  IF( icompq.EQ.1 ) THEN
479  givptr = givptr + 1
480  idxjp = idxq( idx( jprev )+1 )
481  idxj = idxq( idx( j )+1 )
482  IF( idxjp.LE.nlp1 ) THEN
483  idxjp = idxjp - 1
484  END IF
485  IF( idxj.LE.nlp1 ) THEN
486  idxj = idxj - 1
487  END IF
488  givcol( givptr, 2 ) = idxjp
489  givcol( givptr, 1 ) = idxj
490  givnum( givptr, 2 ) = c
491  givnum( givptr, 1 ) = s
492  END IF
493  CALL drot( 1, vf( jprev ), 1, vf( j ), 1, c, s )
494  CALL drot( 1, vl( jprev ), 1, vl( j ), 1, c, s )
495  k2 = k2 - 1
496  idxp( k2 ) = jprev
497  jprev = j
498  ELSE
499  k = k + 1
500  zw( k ) = z( jprev )
501  dsigma( k ) = d( jprev )
502  idxp( k ) = jprev
503  jprev = j
504  END IF
505  END IF
506  GO TO 80
507  90 CONTINUE
508 *
509 * Record the last singular value.
510 *
511  k = k + 1
512  zw( k ) = z( jprev )
513  dsigma( k ) = d( jprev )
514  idxp( k ) = jprev
515 *
516  100 CONTINUE
517 *
518 * Sort the singular values into DSIGMA. The singular values which
519 * were not deflated go into the first K slots of DSIGMA, except
520 * that DSIGMA(1) is treated separately.
521 *
522  DO 110 j = 2, n
523  jp = idxp( j )
524  dsigma( j ) = d( jp )
525  vfw( j ) = vf( jp )
526  vlw( j ) = vl( jp )
527  110 CONTINUE
528  IF( icompq.EQ.1 ) THEN
529  DO 120 j = 2, n
530  jp = idxp( j )
531  perm( j ) = idxq( idx( jp )+1 )
532  IF( perm( j ).LE.nlp1 ) THEN
533  perm( j ) = perm( j ) - 1
534  END IF
535  120 CONTINUE
536  END IF
537 *
538 * The deflated singular values go back into the last N - K slots of
539 * D.
540 *
541  CALL dcopy( n-k, dsigma( k+1 ), 1, d( k+1 ), 1 )
542 *
543 * Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
544 * VL(M).
545 *
546  dsigma( 1 ) = zero
547  hlftol = tol / two
548  IF( abs( dsigma( 2 ) ).LE.hlftol )
549  $ dsigma( 2 ) = hlftol
550  IF( m.GT.n ) THEN
551  z( 1 ) = dlapy2( z1, z( m ) )
552  IF( z( 1 ).LE.tol ) THEN
553  c = one
554  s = zero
555  z( 1 ) = tol
556  ELSE
557  c = z1 / z( 1 )
558  s = -z( m ) / z( 1 )
559  END IF
560  CALL drot( 1, vf( m ), 1, vf( 1 ), 1, c, s )
561  CALL drot( 1, vl( m ), 1, vl( 1 ), 1, c, s )
562  ELSE
563  IF( abs( z1 ).LE.tol ) THEN
564  z( 1 ) = tol
565  ELSE
566  z( 1 ) = z1
567  END IF
568  END IF
569 *
570 * Restore Z, VF, and VL.
571 *
572  CALL dcopy( k-1, zw( 2 ), 1, z( 2 ), 1 )
573  CALL dcopy( n-1, vfw( 2 ), 1, vf( 2 ), 1 )
574  CALL dcopy( n-1, vlw( 2 ), 1, vl( 2 ), 1 )
575 *
576  RETURN
577 *
578 * End of DLASD7
579 *
Here is the call graph for this function:
Here is the caller graph for this function:
dlamrg
subroutine dlamrg(N1, N2, A, DTRD1, DTRD2, INDEX)
DLAMRG creates a permutation list to merge the entries of two independently sorted sets into a single...
Definition: dlamrg.f:101
dlapy2
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:65
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:94
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:70