LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ slasd7()

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

SLASD7 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 SLASD7 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 SLASD7 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.

 SLASD7 is called from SLASD6.
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 REAL 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 REAL array, dimension ( M )
         On exit Z contains the updating row vector in the secular
         equation.
[out]ZW
          ZW is REAL array, dimension ( M )
         Workspace for Z.
[in,out]VF
          VF is REAL 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 REAL array, dimension ( M )
         Workspace for VF.
[in,out]VL
          VL is REAL 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 REAL array, dimension ( M )
         Workspace for VL.
[in]ALPHA
          ALPHA is REAL
         Contains the diagonal element associated with the added row.
[in]BETA
          BETA is REAL
         Contains the off-diagonal element associated with the added
         row.
[out]DSIGMA
          DSIGMA is REAL 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 REAL 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 REAL
         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 REAL
         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.
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 276 of file slasd7.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  REAL ALPHA, BETA, C, S
289 * ..
290 * .. Array Arguments ..
291  INTEGER GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
292  $ IDXQ( * ), PERM( * )
293  REAL D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
294  $ VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
295  $ ZW( * )
296 * ..
297 *
298 * =====================================================================
299 *
300 * .. Parameters ..
301  REAL ZERO, ONE, TWO, EIGHT
302  parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0,
303  $ eight = 8.0e+0 )
304 * ..
305 * .. Local Scalars ..
306 *
307  INTEGER I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
308  $ NLP1, NLP2
309  REAL EPS, HLFTOL, TAU, TOL, Z1
310 * ..
311 * .. External Subroutines ..
312  EXTERNAL scopy, slamrg, srot, xerbla
313 * ..
314 * .. External Functions ..
315  REAL SLAMCH, SLAPY2
316  EXTERNAL slamch, slapy2
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( 'SLASD7', -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 slamrg( 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 = slamch( '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 = slapy2( 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 srot( 1, vf( jprev ), 1, vf( j ), 1, c, s )
491  CALL srot( 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 scopy( 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 ) = slapy2( 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 srot( 1, vf( m ), 1, vf( 1 ), 1, c, s )
558  CALL srot( 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 scopy( k-1, zw( 2 ), 1, z( 2 ), 1 )
570  CALL scopy( n-1, vfw( 2 ), 1, vf( 2 ), 1 )
571  CALL scopy( n-1, vlw( 2 ), 1, vl( 2 ), 1 )
572 *
573  RETURN
574 *
575 * End of SLASD7
576 *
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:63
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slamrg(N1, N2, A, STRD1, STRD2, INDEX)
SLAMRG creates a permutation list to merge the entries of two independently sorted sets into a single...
Definition: slamrg.f:99
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: