LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 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.

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.```
Date
September 2012
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.4.2) --
284 * -- LAPACK is a software package provided by Univ. of Tennessee, --
285 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
286 * September 2012
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 tolerence
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 *
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:65

Here is the call graph for this function:

Here is the caller graph for this function: