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

Go to the source code of this file.

Functions/Subroutines

subroutine slasd2 (NL, NR, SQRE, K, D, Z, ALPHA, BETA, U, LDU, VT, LDVT, DSIGMA, U2, LDU2, VT2, LDVT2, IDXP, IDX, IDXC, IDXQ, COLTYP, INFO)
 SLASD2 merges the two sets of singular values together into a single sorted set. Used by sbdsdc. More...
 

Function/Subroutine Documentation

subroutine slasd2 ( integer  NL,
integer  NR,
integer  SQRE,
integer  K,
real, dimension( * )  D,
real, dimension( * )  Z,
real  ALPHA,
real  BETA,
real, dimension( ldu, * )  U,
integer  LDU,
real, dimension( ldvt, * )  VT,
integer  LDVT,
real, dimension( * )  DSIGMA,
real, dimension( ldu2, * )  U2,
integer  LDU2,
real, dimension( ldvt2, * )  VT2,
integer  LDVT2,
integer, dimension( * )  IDXP,
integer, dimension( * )  IDX,
integer, dimension( * )  IDXC,
integer, dimension( * )  IDXQ,
integer, dimension( * )  COLTYP,
integer  INFO 
)

SLASD2 merges the two sets of singular values together into a single sorted set. Used by sbdsdc.

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

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

 SLASD2 is called from SLASD1.
Parameters
[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 (N)
         On exit Z contains the updating row vector in the secular
         equation.
[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.
[in,out]U
          U is REAL array, dimension (LDU,N)
         On entry U contains the left singular vectors of two
         submatrices in the two square blocks with corners at (1,1),
         (NL, NL), and (NL+2, NL+2), (N,N).
         On exit U contains the trailing (N-K) updated left singular
         vectors (those which were deflated) in its last N-K columns.
[in]LDU
          LDU is INTEGER
         The leading dimension of the array U.  LDU >= N.
[in,out]VT
          VT is REAL array, dimension (LDVT,M)
         On entry VT**T contains the right singular vectors of two
         submatrices in the two square blocks with corners at (1,1),
         (NL+1, NL+1), and (NL+2, NL+2), (M,M).
         On exit VT**T contains the trailing (N-K) updated right singular
         vectors (those which were deflated) in its last N-K columns.
         In case SQRE =1, the last row of VT spans the right null
         space.
[in]LDVT
          LDVT is INTEGER
         The leading dimension of the array VT.  LDVT >= M.
[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]U2
          U2 is REAL array, dimension (LDU2,N)
         Contains a copy of the first K-1 left singular vectors which
         will be used by SLASD3 in a matrix multiply (SGEMM) to solve
         for the new left singular vectors. U2 is arranged into four
         blocks. The first block contains a column with 1 at NL+1 and
         zero everywhere else; the second block contains non-zero
         entries only at and above NL; the third contains non-zero
         entries only below NL+1; and the fourth is dense.
[in]LDU2
          LDU2 is INTEGER
         The leading dimension of the array U2.  LDU2 >= N.
[out]VT2
          VT2 is REAL array, dimension (LDVT2,N)
         VT2**T contains a copy of the first K right singular vectors
         which will be used by SLASD3 in a matrix multiply (SGEMM) to
         solve for the new right singular vectors. VT2 is arranged into
         three blocks. The first block contains a row that corresponds
         to the special 0 diagonal element in SIGMA; the second block
         contains non-zeros only at and before NL +1; the third block
         contains non-zeros only at and after  NL +2.
[in]LDVT2
          LDVT2 is INTEGER
         The leading dimension of the array VT2.  LDVT2 >= M.
[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.
[out]IDX
          IDX is INTEGER array, dimension (N)
         This will contain the permutation used to sort the contents of
         D into ascending order.
[out]IDXC
          IDXC is INTEGER array, dimension (N)
         This will contain the permutation used to arrange the columns
         of the deflated U matrix into three groups:  the first group
         contains non-zero entries only at and above NL, the second
         contains non-zero entries only below NL+2, and the third is
         dense.
[in,out]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 hlaf 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]COLTYP
          COLTYP is INTEGER array, dimension (N)
         As workspace, this will contain a label which will indicate
         which of the following types a column in the U2 matrix or a
         row in the VT2 matrix is:
         1 : non-zero in the upper half only
         2 : non-zero in the lower half only
         3 : dense
         4 : deflated

         On exit, it is an array of dimension 4, with COLTYP(I) being
         the dimension of the I-th type columns.
[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 271 of file slasd2.f.

271 *
272 * -- LAPACK auxiliary routine (version 3.4.2) --
273 * -- LAPACK is a software package provided by Univ. of Tennessee, --
274 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
275 * September 2012
276 *
277 * .. Scalar Arguments ..
278  INTEGER info, k, ldu, ldu2, ldvt, ldvt2, nl, nr, sqre
279  REAL alpha, beta
280 * ..
281 * .. Array Arguments ..
282  INTEGER coltyp( * ), idx( * ), idxc( * ), idxp( * ),
283  $ idxq( * )
284  REAL d( * ), dsigma( * ), u( ldu, * ),
285  $ u2( ldu2, * ), vt( ldvt, * ), vt2( ldvt2, * ),
286  $ z( * )
287 * ..
288 *
289 * =====================================================================
290 *
291 * .. Parameters ..
292  REAL zero, one, two, eight
293  parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0,
294  $ eight = 8.0e+0 )
295 * ..
296 * .. Local Arrays ..
297  INTEGER ctot( 4 ), psm( 4 )
298 * ..
299 * .. Local Scalars ..
300  INTEGER ct, i, idxi, idxj, idxjp, j, jp, jprev, k2, m,
301  $ n, nlp1, nlp2
302  REAL c, eps, hlftol, s, tau, tol, z1
303 * ..
304 * .. External Functions ..
305  REAL slamch, slapy2
306  EXTERNAL slamch, slapy2
307 * ..
308 * .. External Subroutines ..
309  EXTERNAL scopy, slacpy, slamrg, slaset, srot, xerbla
310 * ..
311 * .. Intrinsic Functions ..
312  INTRINSIC abs, max
313 * ..
314 * .. Executable Statements ..
315 *
316 * Test the input parameters.
317 *
318  info = 0
319 *
320  IF( nl.LT.1 ) THEN
321  info = -1
322  ELSE IF( nr.LT.1 ) THEN
323  info = -2
324  ELSE IF( ( sqre.NE.1 ) .AND. ( sqre.NE.0 ) ) THEN
325  info = -3
326  END IF
327 *
328  n = nl + nr + 1
329  m = n + sqre
330 *
331  IF( ldu.LT.n ) THEN
332  info = -10
333  ELSE IF( ldvt.LT.m ) THEN
334  info = -12
335  ELSE IF( ldu2.LT.n ) THEN
336  info = -15
337  ELSE IF( ldvt2.LT.m ) THEN
338  info = -17
339  END IF
340  IF( info.NE.0 ) THEN
341  CALL xerbla( 'SLASD2', -info )
342  RETURN
343  END IF
344 *
345  nlp1 = nl + 1
346  nlp2 = nl + 2
347 *
348 * Generate the first part of the vector Z; and move the singular
349 * values in the first part of D one position backward.
350 *
351  z1 = alpha*vt( nlp1, nlp1 )
352  z( 1 ) = z1
353  DO 10 i = nl, 1, -1
354  z( i+1 ) = alpha*vt( i, nlp1 )
355  d( i+1 ) = d( i )
356  idxq( i+1 ) = idxq( i ) + 1
357  10 CONTINUE
358 *
359 * Generate the second part of the vector Z.
360 *
361  DO 20 i = nlp2, m
362  z( i ) = beta*vt( i, nlp2 )
363  20 CONTINUE
364 *
365 * Initialize some reference arrays.
366 *
367  DO 30 i = 2, nlp1
368  coltyp( i ) = 1
369  30 CONTINUE
370  DO 40 i = nlp2, n
371  coltyp( i ) = 2
372  40 CONTINUE
373 *
374 * Sort the singular values into increasing order
375 *
376  DO 50 i = nlp2, n
377  idxq( i ) = idxq( i ) + nlp1
378  50 CONTINUE
379 *
380 * DSIGMA, IDXC, IDXC, and the first column of U2
381 * are used as storage space.
382 *
383  DO 60 i = 2, n
384  dsigma( i ) = d( idxq( i ) )
385  u2( i, 1 ) = z( idxq( i ) )
386  idxc( i ) = coltyp( idxq( i ) )
387  60 CONTINUE
388 *
389  CALL slamrg( nl, nr, dsigma( 2 ), 1, 1, idx( 2 ) )
390 *
391  DO 70 i = 2, n
392  idxi = 1 + idx( i )
393  d( i ) = dsigma( idxi )
394  z( i ) = u2( idxi, 1 )
395  coltyp( i ) = idxc( idxi )
396  70 CONTINUE
397 *
398 * Calculate the allowable deflation tolerance
399 *
400  eps = slamch( 'Epsilon' )
401  tol = max( abs( alpha ), abs( beta ) )
402  tol = eight*eps*max( abs( d( n ) ), tol )
403 *
404 * There are 2 kinds of deflation -- first a value in the z-vector
405 * is small, second two (or more) singular values are very close
406 * together (their difference is small).
407 *
408 * If the value in the z-vector is small, we simply permute the
409 * array so that the corresponding singular value is moved to the
410 * end.
411 *
412 * If two values in the D-vector are close, we perform a two-sided
413 * rotation designed to make one of the corresponding z-vector
414 * entries zero, and then permute the array so that the deflated
415 * singular value is moved to the end.
416 *
417 * If there are multiple singular values then the problem deflates.
418 * Here the number of equal singular values are found. As each equal
419 * singular value is found, an elementary reflector is computed to
420 * rotate the corresponding singular subspace so that the
421 * corresponding components of Z are zero in this new basis.
422 *
423  k = 1
424  k2 = n + 1
425  DO 80 j = 2, n
426  IF( abs( z( j ) ).LE.tol ) THEN
427 *
428 * Deflate due to small z component.
429 *
430  k2 = k2 - 1
431  idxp( k2 ) = j
432  coltyp( j ) = 4
433  IF( j.EQ.n )
434  $ GO TO 120
435  ELSE
436  jprev = j
437  GO TO 90
438  END IF
439  80 CONTINUE
440  90 CONTINUE
441  j = jprev
442  100 CONTINUE
443  j = j + 1
444  IF( j.GT.n )
445  $ GO TO 110
446  IF( abs( z( j ) ).LE.tol ) THEN
447 *
448 * Deflate due to small z component.
449 *
450  k2 = k2 - 1
451  idxp( k2 ) = j
452  coltyp( j ) = 4
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  c = c / tau
469  s = -s / tau
470  z( j ) = tau
471  z( jprev ) = zero
472 *
473 * Apply back the Givens rotation to the left and right
474 * singular vector matrices.
475 *
476  idxjp = idxq( idx( jprev )+1 )
477  idxj = idxq( idx( j )+1 )
478  IF( idxjp.LE.nlp1 ) THEN
479  idxjp = idxjp - 1
480  END IF
481  IF( idxj.LE.nlp1 ) THEN
482  idxj = idxj - 1
483  END IF
484  CALL srot( n, u( 1, idxjp ), 1, u( 1, idxj ), 1, c, s )
485  CALL srot( m, vt( idxjp, 1 ), ldvt, vt( idxj, 1 ), ldvt, c,
486  $ s )
487  IF( coltyp( j ).NE.coltyp( jprev ) ) THEN
488  coltyp( j ) = 3
489  END IF
490  coltyp( jprev ) = 4
491  k2 = k2 - 1
492  idxp( k2 ) = jprev
493  jprev = j
494  ELSE
495  k = k + 1
496  u2( k, 1 ) = z( jprev )
497  dsigma( k ) = d( jprev )
498  idxp( k ) = jprev
499  jprev = j
500  END IF
501  END IF
502  GO TO 100
503  110 CONTINUE
504 *
505 * Record the last singular value.
506 *
507  k = k + 1
508  u2( k, 1 ) = z( jprev )
509  dsigma( k ) = d( jprev )
510  idxp( k ) = jprev
511 *
512  120 CONTINUE
513 *
514 * Count up the total number of the various types of columns, then
515 * form a permutation which positions the four column types into
516 * four groups of uniform structure (although one or more of these
517 * groups may be empty).
518 *
519  DO 130 j = 1, 4
520  ctot( j ) = 0
521  130 CONTINUE
522  DO 140 j = 2, n
523  ct = coltyp( j )
524  ctot( ct ) = ctot( ct ) + 1
525  140 CONTINUE
526 *
527 * PSM(*) = Position in SubMatrix (of types 1 through 4)
528 *
529  psm( 1 ) = 2
530  psm( 2 ) = 2 + ctot( 1 )
531  psm( 3 ) = psm( 2 ) + ctot( 2 )
532  psm( 4 ) = psm( 3 ) + ctot( 3 )
533 *
534 * Fill out the IDXC array so that the permutation which it induces
535 * will place all type-1 columns first, all type-2 columns next,
536 * then all type-3's, and finally all type-4's, starting from the
537 * second column. This applies similarly to the rows of VT.
538 *
539  DO 150 j = 2, n
540  jp = idxp( j )
541  ct = coltyp( jp )
542  idxc( psm( ct ) ) = j
543  psm( ct ) = psm( ct ) + 1
544  150 CONTINUE
545 *
546 * Sort the singular values and corresponding singular vectors into
547 * DSIGMA, U2, and VT2 respectively. The singular values/vectors
548 * which were not deflated go into the first K slots of DSIGMA, U2,
549 * and VT2 respectively, while those which were deflated go into the
550 * last N - K slots, except that the first column/row will be treated
551 * separately.
552 *
553  DO 160 j = 2, n
554  jp = idxp( j )
555  dsigma( j ) = d( jp )
556  idxj = idxq( idx( idxp( idxc( j ) ) )+1 )
557  IF( idxj.LE.nlp1 ) THEN
558  idxj = idxj - 1
559  END IF
560  CALL scopy( n, u( 1, idxj ), 1, u2( 1, j ), 1 )
561  CALL scopy( m, vt( idxj, 1 ), ldvt, vt2( j, 1 ), ldvt2 )
562  160 CONTINUE
563 *
564 * Determine DSIGMA(1), DSIGMA(2) and Z(1)
565 *
566  dsigma( 1 ) = zero
567  hlftol = tol / two
568  IF( abs( dsigma( 2 ) ).LE.hlftol )
569  $ dsigma( 2 ) = hlftol
570  IF( m.GT.n ) THEN
571  z( 1 ) = slapy2( z1, z( m ) )
572  IF( z( 1 ).LE.tol ) THEN
573  c = one
574  s = zero
575  z( 1 ) = tol
576  ELSE
577  c = z1 / z( 1 )
578  s = z( m ) / z( 1 )
579  END IF
580  ELSE
581  IF( abs( z1 ).LE.tol ) THEN
582  z( 1 ) = tol
583  ELSE
584  z( 1 ) = z1
585  END IF
586  END IF
587 *
588 * Move the rest of the updating row to Z.
589 *
590  CALL scopy( k-1, u2( 2, 1 ), 1, z( 2 ), 1 )
591 *
592 * Determine the first column of U2, the first row of VT2 and the
593 * last row of VT.
594 *
595  CALL slaset( 'A', n, 1, zero, zero, u2, ldu2 )
596  u2( nlp1, 1 ) = one
597  IF( m.GT.n ) THEN
598  DO 170 i = 1, nlp1
599  vt( m, i ) = -s*vt( nlp1, i )
600  vt2( 1, i ) = c*vt( nlp1, i )
601  170 CONTINUE
602  DO 180 i = nlp2, m
603  vt2( 1, i ) = s*vt( m, i )
604  vt( m, i ) = c*vt( m, i )
605  180 CONTINUE
606  ELSE
607  CALL scopy( m, vt( nlp1, 1 ), ldvt, vt2( 1, 1 ), ldvt2 )
608  END IF
609  IF( m.GT.n ) THEN
610  CALL scopy( m, vt( m, 1 ), ldvt, vt2( m, 1 ), ldvt2 )
611  END IF
612 *
613 * The deflated singular values and their corresponding vectors go
614 * into the back of D, U, and V respectively.
615 *
616  IF( n.GT.k ) THEN
617  CALL scopy( n-k, dsigma( k+1 ), 1, d( k+1 ), 1 )
618  CALL slacpy( 'A', n, n-k, u2( 1, k+1 ), ldu2, u( 1, k+1 ),
619  $ ldu )
620  CALL slacpy( 'A', n-k, m, vt2( k+1, 1 ), ldvt2, vt( k+1, 1 ),
621  $ ldvt )
622  END IF
623 *
624 * Copy CTOT into COLTYP for referencing in SLASD3.
625 *
626  DO 190 j = 1, 4
627  coltyp( j ) = ctot( j )
628  190 CONTINUE
629 *
630  RETURN
631 *
632 * End of SLASD2
633 *
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
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
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
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: