MB04YD

Partial diagonalization of a bidiagonal matrix

[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

  To partially diagonalize the bidiagonal matrix

            |q(1) e(1)  0    ...       0      |
            | 0   q(2) e(2)            .      |
        J = | .                        .      |                  (1)
            | .                  e(MIN(M,N)-1)|
            | 0   ...        ...  q(MIN(M,N)) |

  using QR or QL iterations in such a way that J is split into
  unreduced bidiagonal submatrices whose singular values are either
  all larger than a given bound or are all smaller than (or equal
  to) this bound. The left- and right-hand Givens rotations
  performed on J (corresponding to each QR or QL iteration step) may
  be optionally accumulated in the arrays U and V.

Specification
      SUBROUTINE MB04YD( JOBU, JOBV, M, N, RANK, THETA, Q, E, U, LDU, V,
     $                   LDV, INUL, TOL, RELTOL, DWORK, LDWORK, IWARN,
     $                   INFO )
C     .. Scalar Arguments ..
      CHARACTER         JOBU, JOBV
      INTEGER           INFO, IWARN, LDU, LDV, LDWORK, M, N, RANK
      DOUBLE PRECISION  RELTOL, THETA, TOL
C     .. Array Arguments ..
      LOGICAL           INUL(*)
      DOUBLE PRECISION  DWORK(*), E(*), Q(*), U(LDU,*), V(LDV,*)

Arguments

Mode Parameters

  JOBU    CHARACTER*1
          Indicates whether the user wishes to accumulate in a
          matrix U the left-hand Givens rotations, as follows:
          = 'N':  Do not form U;
          = 'I':  U is initialized to the M-by-MIN(M,N) submatrix of
                  the unit matrix and the left-hand Givens rotations
                  are accumulated in U;
          = 'U':  The given matrix U is updated by the left-hand
                  Givens rotations used in the calculation.

  JOBV    CHARACTER*1
          Indicates whether the user wishes to accumulate in a
          matrix V the right-hand Givens rotations, as follows:
          = 'N':  Do not form V;
          = 'I':  V is initialized to the N-by-MIN(M,N) submatrix of
                  the unit matrix and the right-hand Givens
                  rotations are accumulated in V;
          = 'U':  The given matrix V is updated by the right-hand
                  Givens rotations used in the calculation.

Input/Output Parameters
  M       (input) INTEGER
          The number of rows in matrix U.  M >= 0.

  N       (input) INTEGER
          The number of rows in matrix V.  N >= 0.

  RANK    (input/output) INTEGER
          On entry, if RANK < 0, then the rank of matrix J is
          computed by the routine as the number of singular values
          larger than THETA.
          Otherwise, RANK must specify the rank of matrix J.
          RANK <= MIN(M,N).
          On exit, if RANK < 0 on entry, then RANK contains the
          computed rank of J. That is, the number of singular
          values of J larger than THETA.
          Otherwise, the user-supplied value of RANK may be
          changed by the routine on exit if the RANK-th and the
          (RANK+1)-th singular values of J are considered to be
          equal. See also the parameter TOL.

  THETA   (input/output) DOUBLE PRECISION
          On entry, if RANK < 0, then THETA must specify an upper
          bound on the smallest singular values of J. THETA >= 0.0.
          Otherwise, THETA must specify an initial estimate (t say)
          for computing an upper bound such that precisely RANK
          singular values are greater than this bound.
          If THETA < 0.0, then t is computed by the routine.
          On exit, if RANK >= 0 on entry, then THETA contains the
          computed upper bound such that precisely RANK singular
          values of J are greater than THETA + TOL.
          Otherwise, THETA is unchanged.

  Q       (input/output) DOUBLE PRECISION array, dimension
          (MIN(M,N))
          On entry, this array must contain the diagonal elements
          q(1),q(2),...,q(MIN(M,N)) of the bidiagonal matrix J. That
          is, Q(i) = J(i,i) for i = 1,2,...,MIN(M,N).
          On exit, this array contains the leading diagonal of the
          transformed bidiagonal matrix J.

  E       (input/output) DOUBLE PRECISION array, dimension
          (MIN(M,N)-1)
          On entry, this array must contain the superdiagonal
          elements e(1),e(2),...,e(MIN(M,N)-1) of the bidiagonal
          matrix J. That is, E(k) = J(k,k+1) for k = 1,2,...,
          MIN(M,N)-1.
          On exit, this array contains the superdiagonal of the
          transformed bidiagonal matrix J.

  U       (input/output) DOUBLE PRECISION array, dimension (LDU,*)
          On entry, if JOBU = 'U', the leading M-by-MIN(M,N) part
          of this array must contain a left transformation matrix
          applied to the original matrix of the problem, and
          on exit, the leading M-by-MIN(M,N) part of this array
          contains the product of the input matrix U and the
          left-hand Givens rotations.
          On exit, if JOBU = 'I', then the leading M-by-MIN(M,N)
          part of this array contains the matrix of accumulated
          left-hand Givens rotations used.
          If JOBU = 'N', the array U is not referenced and can be
          supplied as a dummy array (i.e. set parameter LDU = 1 and
          declare this array to be U(1,1) in the calling program).

  LDU     INTEGER
          The leading dimension of array U. If JOBU = 'U' or
          JOBU = 'I', LDU >= MAX(1,M); if JOBU = 'N', LDU >= 1.

  V       (input/output) DOUBLE PRECISION array, dimension (LDV,*)
          On entry, if JOBV = 'U', the leading N-by-MIN(M,N) part
          of this array must contain a right transformation matrix
          applied to the original matrix of the problem, and
          on exit, the leading N-by-MIN(M,N) part of this array
          contains the product of the input matrix V and the
          right-hand Givens rotations.
          On exit, if JOBV = 'I', then the leading N-by-MIN(M,N)
          part of this array contains the matrix of accumulated
          right-hand Givens rotations used.
          If JOBV = 'N', the array V is not referenced and can be
          supplied as a dummy array (i.e. set parameter LDV = 1 and
          declare this array to be V(1,1) in the calling program).

  LDV     INTEGER
          The leading dimension of array V. If JOBV = 'U' or
          JOBV = 'I', LDV >= MAX(1,N); if JOBV = 'N', LDV >= 1.

  INUL    (input/output) LOGICAL array, dimension (MIN(M,N))
          On entry, the leading MIN(M,N) elements of this array must
          be set to .FALSE. unless the i-th columns of U (if JOBU =
          'U') and V (if JOBV = 'U') already contain a computed base
          vector of the desired singular subspace of the original
          matrix, in which case INUL(i) must be set to .TRUE.
          for 1 <= i <= MIN(M,N).
          On exit, the indices of the elements of this array with
          value .TRUE. indicate the indices of the diagonal entries
          of J which belong to those bidiagonal submatrices whose
          singular values are all less than or equal to THETA.

Tolerances
  TOL     DOUBLE PRECISION
          This parameter defines the multiplicity of singular values
          by considering all singular values within an interval of
          length TOL as coinciding. TOL is used in checking how many
          singular values are less than or equal to THETA. Also in
          computing an appropriate upper bound THETA by a bisection
          method, TOL is used as a stopping criterion defining the
          minimum (absolute) subinterval width. TOL is also taken
          as an absolute tolerance for negligible elements in the
          QR/QL iterations. If the user sets TOL to be less than or
          equal to 0, then the tolerance is taken as
          EPS * MAX(ABS(Q(i)), ABS(E(k))), where EPS is the
          machine precision (see LAPACK Library routine DLAMCH),
          i = 1,2,...,MIN(M,N) and k = 1,2,...,MIN(M,N)-1.

  RELTOL  DOUBLE PRECISION
          This parameter specifies the minimum relative width of an
          interval. When an interval is narrower than TOL, or than
          RELTOL times the larger (in magnitude) endpoint, then it
          is considered to be sufficiently small and bisection has
          converged. If the user sets RELTOL to be less than
          BASE * EPS, where BASE is machine radix and EPS is machine
          precision (see LAPACK Library routine DLAMCH), then the
          tolerance is taken as BASE * EPS.

Workspace
  DWORK   DOUBLE PRECISION array, dimension (LDWORK)

  LDWORK  INTEGER
          The length of the array DWORK.
          LDWORK >= MAX(1,6*MIN(M,N)-5), if JOBU = 'I' or 'U', or
                                            JOBV = 'I' or 'U';
          LDWORK >= MAX(1,4*MIN(M,N)-3), if JOBU = 'N' and
                                            JOBV = 'N'.

Warning Indicator
  IWARN   INTEGER
          = 0:  no warning;
          = 1:  if the rank of the bidiagonal matrix J (as specified
                by the user) has been lowered because a singular
                value of multiplicity larger than 1 was found.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value; this includes values like RANK > MIN(M,N), or
                THETA < 0.0 and RANK < 0;
          = 1:  if the maximum number of QR/QL iteration steps
                (30*MIN(M,N)) has been exceeded.

Method
  If the upper bound THETA is not specified by the user, then it is
  computed by the routine (using a bisection method) such that
  precisely (MIN(M,N) - RANK) singular values of J are less than or
  equal to THETA + TOL.

  The method used by the routine (see [1]) then proceeds as follows.

  The unreduced bidiagonal submatrices of J(j), where J(j) is the
  transformed bidiagonal matrix after the j-th iteration step, are
  classified into the following three classes:

  - C1 contains the bidiagonal submatrices with all singular values
    > THETA,
  - C2 contains the bidiagonal submatrices with all singular values
    <= THETA and
  - C3 contains the bidiagonal submatrices with singular values
    > THETA and also singular values <= THETA.

  If C3 is empty, then the partial diagonalization is complete, and
  RANK is the sum of the dimensions of the bidiagonal submatrices of
  C1.
  Otherwise, QR or QL iterations are performed on each bidiagonal
  submatrix of C3, until this bidiagonal submatrix has been split
  into two bidiagonal submatrices. These two submatrices are then
  classified and the iterations are restarted.
  If the upper left diagonal element of the bidiagonal submatrix is
  larger than its lower right diagonal element, then QR iterations
  are performed, else QL iterations are used. The shift is taken as
  the smallest diagonal element of the bidiagonal submatrix (in
  magnitude) unless its value exceeds THETA, in which case it is
  taken as zero.

References
  [1] Van Huffel, S., Vandewalle, J. and Haegemans, A.
      An efficient and reliable algorithm for computing the
      singular subspace of a matrix associated with its smallest
      singular values.
      J. Comput. and Appl. Math., 19, pp. 313-330, 1987.

Numerical Aspects
  The algorithm is backward stable.

  To avoid overflow, matrix J is scaled so that its largest element
  is no greater than  overflow**(1/2) * underflow**(1/4) in absolute
  value (and not much smaller than that, for maximal accuracy).

Further Comments
  None
Example

Program Text

*     MB04YD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO
      PARAMETER        ( ZERO = 0.0D0 )
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          MMAX, NMAX
      PARAMETER        ( MMAX = 20, NMAX = 20 )
      INTEGER          MNMIN
      PARAMETER        ( MNMIN = MIN( MMAX, NMAX ) )
      INTEGER          LDU, LDV
      PARAMETER        ( LDU = MMAX, LDV = NMAX )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = 6*MNMIN - 5 )
*     .. Local Scalars ..
      DOUBLE PRECISION RELTOL, THETA, TOL
      INTEGER          I, INFO, IWARN, J, M, MINMN, N, RANK, RANK1
      CHARACTER*1      JOBU, JOBV
      LOGICAL          LJOBUU, LJOBVU
*     .. Local Arrays ..
      DOUBLE PRECISION DWORK(LDWORK), E(MNMIN-1), Q(MNMIN),
     $                 U(LDU,MNMIN), V(LDV,MNMIN)
      LOGICAL          INUL(MNMIN)
*     .. External Functions ..
      LOGICAL          LSAME
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         MB04YD
*     .. Intrinsic Functions ..
      INTRINSIC        MIN
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) M, N, THETA, RANK, TOL, RELTOL, JOBU, JOBV
      MINMN = MIN( M, N )
      IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
         WRITE ( NOUT, FMT = 99988 ) M
      ELSE IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99987 ) N
      ELSE IF ( RANK.GT.MINMN ) THEN
         WRITE ( NOUT, FMT = 99986 ) RANK
      ELSE IF ( RANK.LT.0 .AND. THETA.LT.ZERO ) THEN
         WRITE ( NOUT, FMT = 99985 ) THETA
      ELSE
         READ ( NIN, FMT = * ) ( Q(I), I = 1,MINMN )
         READ ( NIN, FMT = * ) ( E(I), I = 1,MINMN-1 )
         RANK1 = RANK
         LJOBUU = LSAME( JOBU, 'U' )
         LJOBVU = LSAME( JOBV, 'U' )
         IF ( LJOBUU ) READ ( NIN, FMT = * )
     $                      ( ( U(I,J), J = 1,MINMN ), I = 1,M )
         IF ( LJOBVU ) READ ( NIN, FMT = * )
     $                      ( ( V(I,J), J = 1,MINMN ), I = 1,N )
*        Initialise the array INUL.
         DO 20 I = 1, MINMN
            INUL(I) = .FALSE.
   20    CONTINUE
         IF ( LJOBUU.OR.LJOBVU ) READ ( NIN, FMT = * )
     $                                ( INUL(I), I = 1,MINMN )
*        Compute the number of singular values of J > THETA.
         CALL MB04YD( JOBU, JOBV, M, N, RANK, THETA, Q, E, U, LDU, V,
     $                LDV, INUL, TOL, RELTOL, DWORK, LDWORK, IWARN,
     $                INFO )
*
         IF ( INFO.NE.0 ) THEN
            WRITE ( NOUT, FMT = 99998 ) INFO
         ELSE
            IF ( IWARN.NE.0 ) THEN
               WRITE ( NOUT, FMT = 99993 ) IWARN
               WRITE ( NOUT, FMT = 99984 ) RANK
            END IF
            WRITE ( NOUT, FMT = 99997 )
            DO 160 I = 1, MINMN - 1
               WRITE ( NOUT, FMT = 99996 ) I, I, Q(I), I, (I+1), E(I)
  160       CONTINUE
            WRITE ( NOUT, FMT = 99995 ) MINMN, MINMN, Q(MINMN)
            IF ( RANK1.LT.0 ) WRITE ( NOUT, FMT = 99994 ) RANK, THETA
            IF ( .NOT.LSAME( JOBV, 'N' ) ) THEN
               WRITE ( NOUT, FMT = 99992 )
               DO 180 I = 1, N
                  WRITE ( NOUT, FMT = 99991 ) ( V(I,J), J = 1,MINMN )
  180          CONTINUE
            END IF
            IF ( ( .NOT.LSAME( JOBU, 'N' ) ) .AND.
     $           ( .NOT.LSAME( JOBV, 'N' ) ) )
     $           WRITE ( NOUT, FMT = 99990 )
            IF ( .NOT.LSAME( JOBU, 'N' ) ) THEN
               WRITE ( NOUT, FMT = 99989 )
               DO 200 I = 1, M
                  WRITE ( NOUT, FMT = 99991 ) ( U(I,J), J = 1,MINMN )
  200          CONTINUE
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' MB04YD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB04YD = ',I2)
99997 FORMAT (' The transformed bidiagonal matrix J is',/)
99996 FORMAT (2(' (',I1,',',I1,') = ',F7.4,2X))
99995 FORMAT (' (',I1,',',I1,') = ',F7.4)
99994 FORMAT (/' J has ',I2,' singular values >',F7.4,/)
99993 FORMAT (' IWARN on exit from MB04YD = ',I2,/)
99992 FORMAT (' The product of the right-hand Givens rotation matrices',
     $       ' equals ')
99991 FORMAT (20(1X,F8.4))
99990 FORMAT (' ')
99989 FORMAT (' The product of the left-hand Givens rotation matrices ',
     $       'equals ')
99988 FORMAT (/' M is out of range.',/' M = ',I5)
99987 FORMAT (/' N is out of range.',/' N = ',I5)
99986 FORMAT (/' RANK is out of range.',/' RANK = ',I5)
99985 FORMAT (/' THETA must be at least zero.',/' THETA = ',F8.4)
99984 FORMAT (/' The computed rank of matrix J = ',I3,/)
      END
Program Data
 MB04YD EXAMPLE PROGRAM DATA
   5     5     2.0     -1     0.0     0.0     N     N
   1.0  2.0  3.0  4.0  5.0
   2.0  3.0  4.0  5.0
Program Results
 MB04YD EXAMPLE PROGRAM RESULTS

 The transformed bidiagonal matrix J is

 (1,1) =  0.4045   (1,2) =  0.0000
 (2,2) =  1.9839   (2,3) =  0.0000
 (3,3) =  3.4815   (3,4) =  0.0128
 (4,4) =  5.3723   (4,5) =  0.0273
 (5,5) =  7.9948

 J has  3 singular values > 2.0000


Return to index