 LAPACK  3.10.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.

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.
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 276 of file dlasd7.f.

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