LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
slasd7.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine slasd7 (ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL, VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, 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. More...
 

Function/Subroutine Documentation

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

Definition at line 282 of file slasd7.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  REAL alpha, beta, c, s
292 * ..
293 * .. Array Arguments ..
294  INTEGER givcol( ldgcol, * ), idx( * ), idxp( * ),
295  $ idxq( * ), perm( * )
296  REAL d( * ), dsigma( * ), givnum( ldgnum, * ),
297  $ vf( * ), vfw( * ), vl( * ), vlw( * ), z( * ),
298  $ zw( * )
299 * ..
300 *
301 * =====================================================================
302 *
303 * .. Parameters ..
304  REAL zero, one, two, eight
305  parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0,
306  $ eight = 8.0e+0 )
307 * ..
308 * .. Local Scalars ..
309 *
310  INTEGER i, idxi, idxj, idxjp, j, jp, jprev, k2, m, n,
311  $ nlp1, nlp2
312  REAL eps, hlftol, tau, tol, z1
313 * ..
314 * .. External Subroutines ..
315  EXTERNAL scopy, slamrg, srot, xerbla
316 * ..
317 * .. External Functions ..
318  REAL slamch, slapy2
319  EXTERNAL slamch, slapy2
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( 'SLASD7', -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 slamrg( 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 = slamch( '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 = slapy2( 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 srot( 1, vf( jprev ), 1, vf( j ), 1, c, s )
494  CALL srot( 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 scopy( 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 ) = slapy2( 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 srot( 1, vf( m ), 1, vf( 1 ), 1, c, s )
561  CALL srot( 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 scopy( k-1, zw( 2 ), 1, z( 2 ), 1 )
573  CALL scopy( n-1, vfw( 2 ), 1, vf( 2 ), 1 )
574  CALL scopy( n-1, vlw( 2 ), 1, vl( 2 ), 1 )
575 *
576  RETURN
577 *
578 * End of SLASD7
579 *
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:101
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
real function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
Definition: slapy2.f:65
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:53

Here is the call graph for this function:

Here is the caller graph for this function: