MB03ND

Number of singular values of a bidiagonal matrix less than a bound

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

Purpose

  To find the number of singular values of the bidiagonal matrix

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

  which are less than or equal to a given bound THETA.

  This routine is intended to be called only by other SLICOT
  routines.

Specification
      INTEGER FUNCTION MB03ND( N, THETA, Q2, E2, PIVMIN, INFO )
C     .. Scalar Arguments ..
      INTEGER           INFO, N
      DOUBLE PRECISION  PIVMIN, THETA
C     .. Array Arguments ..
      DOUBLE PRECISION  E2(*), Q2(*)

Arguments

Input/Output Parameters

  N       (input) INTEGER
          The order of the bidiagonal matrix J.  N >= 0.

  THETA   (input) DOUBLE PRECISION
          Given bound.
          Note: If THETA < 0.0 on entry, then MB03ND is set to 0
                as the singular values of J are non-negative.

  Q2      (input) DOUBLE PRECISION array, dimension (N)
          This array must contain the squares of the diagonal
          elements q(1),q(2),...,q(N) of the bidiagonal matrix J.
          That is, Q2(i) = J(i,i)**2 for i = 1,2,...,N.

  E2      (input) DOUBLE PRECISION array, dimension (N-1)
          This array must contain the squares of the superdiagonal
          elements e(1),e(2),...,e(N-1) of the bidiagonal matrix J.
          That is, E2(k) = J(k,k+1)**2 for k = 1,2,...,N-1.

  PIVMIN  (input) DOUBLE PRECISION
          The minimum absolute value of a "pivot" in the Sturm
          sequence loop.
          PIVMIN >= max( max( |q(i)|, |e(k)| )**2*sf_min, sf_min ),
          where i = 1,2,...,N, k = 1,2,...,N-1, and sf_min is at
          least the smallest number that can divide one without
          overflow (see LAPACK Library routine DLAMCH).
          Note that this condition is not checked by the routine.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit;
          < 0:  if INFO = -i, the i-th argument had an illegal
                value.

Method
  The computation of the number of singular values s(i) of J which
  are less than or equal to THETA is based on applying Sylvester's
  Law of Inertia, or equivalently, Sturm sequences [1,p.52] to the
  unreduced symmetric tridiagonal matrices associated with J as
  follows. Let T be the following 2N-by-2N symmetric matrix
  associated with J:

            | 0   J'|
       T =  |       |.
            | J   0 |

  (The eigenvalues of T are given by s(1),s(2),...,s(N),-s(1),-s(2),
  ...,-s(N)). Then, by permuting the rows and columns of T into the
  order 1, N+1, 2, N+2, ..., N, 2N it follows that T is orthogonally
  similar to the tridiagonal matrix T" with zeros on its diagonal
  and q(1), e(1), q(2), e(2), ..., e(N-1), q(N) on its offdiagonals
  [3,4]. If q(1),q(2),...,q(N) and e(1),e(2),...,e(N-1) are nonzero,
  Sylvester's Law of Inertia may be applied directly to T".
  Otherwise, T" is block diagonal and each diagonal block (which is
  then unreduced) must be analysed separately by applying
  Sylvester's Law of Inertia.

References
  [1] Parlett, B.N.
      The Symmetric Eigenvalue Problem.
      Prentice Hall, Englewood Cliffs, New Jersey, 1980.

  [2] Demmel, J. and Kahan, W.
      Computing Small Singular Values of Bidiagonal Matrices with
      Guaranteed High Relative Accuracy.
      Technical Report, Courant Inst., New York, March 1988.

  [3] Van Huffel, S. and Vandewalle, J.
      The Partial Total Least-Squares Algorithm.
      J. Comput. and Appl. Math., 21, pp. 333-341, 1988.

  [4] Golub, G.H. and Kahan, W.
      Calculating the Singular Values and Pseudo-inverse of a
      Matrix.
      SIAM J. Numer. Anal., Ser. B, 2, pp. 205-224, 1965.

  [5] Demmel, J.W., Dhillon, I. and Ren, H.
      On the Correctness of Parallel Bisection in Floating Point.
      Computer Science Division Technical Report UCB//CSD-94-805,
      University of California, Berkeley, CA 94720, March 1994.

Numerical Aspects
  The singular values s(i) could also be obtained with the use of
  the symmetric tridiagonal matrix T = J'J, whose eigenvalues are
  the squared singular values of J [4,p.213]. However, the method
  actually used by the routine is more accurate and equally
  efficient (see [2]).

  To avoid overflow, matrix J should be 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).

  With respect to accuracy the following condition holds (see [2]):

  If the established value is denoted by p, then at least p
  singular values of J are less than or equal to
  THETA/(1 - (3 x N - 1.5) x EPS) and no more than p singular values
  are less than or equal to
  THETA x (1 - (6 x N-2) x EPS)/(1 - (3 x N - 1.5) x EPS).

Further Comments
  None
Example

Program Text

*     MB03ND EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX
      PARAMETER        ( NMAX = 20 )
*     .. Local Scalars ..
      DOUBLE PRECISION PIVMIN, SAFMIN, THETA
      INTEGER          I, INFO, N, NUMSV
*     .. Local Arrays ..
      DOUBLE PRECISION E(NMAX-1), E2(NMAX-1), Q(NMAX), Q2(NMAX)
*     .. External Functions ..
      DOUBLE PRECISION DLAMCH
      EXTERNAL         DLAMCH
*     .. External Functions ..
      INTEGER          MB03ND
      EXTERNAL         MB03ND
*     .. Intrinsic Functions ..
      INTRINSIC        MAX
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) N, THETA
      IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99993 ) N
      ELSE
         READ ( NIN, FMT = * ) ( Q(I), I = 1,N )
         READ ( NIN, FMT = * ) ( E(I), I = 1,N-1 )
*        Print out the bidiagonal matrix J.
         WRITE ( NOUT, FMT = 99997 )
         DO 20 I = 1, N - 1
            WRITE ( NOUT, FMT = 99996 ) I, I, Q(I), I, (I+1), E(I)
   20    CONTINUE
         WRITE ( NOUT, FMT = 99995 ) N, N, Q(N)
*        Compute Q**2, E**2, and PIVMIN.
         Q2(N) = Q(N)**2
         PIVMIN = Q2(N)
         DO 40 I = 1, N - 1
            Q2(I) = Q(I)**2
            E2(I) = E(I)**2
            PIVMIN = MAX( PIVMIN, Q2(I), E2(I) )
   40    CONTINUE
         SAFMIN = DLAMCH( 'Safe minimum' )
         PIVMIN = MAX( PIVMIN*SAFMIN, SAFMIN )
*        Compute the number of singular values of J < =  THETA.
         NUMSV = MB03ND( N, THETA, Q2, E2, PIVMIN, INFO )
*
         IF ( INFO.NE.0 ) THEN
            WRITE ( NOUT, FMT = 99998 ) INFO
         ELSE
            WRITE ( NOUT, FMT = 99994 ) NUMSV, THETA
         END IF
      END IF
      STOP
*
99999 FORMAT (' MB03ND EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB03ND = ',I2)
99997 FORMAT (' The 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 (/' N is out of range.',/' N = ',I5)
      END
Program Data
 MB03ND EXAMPLE PROGRAM DATA
   5     5.0     0.0     0.0
   1.0  2.0  3.0  4.0  5.0
   2.0  3.0  4.0  5.0
Program Results
 MB03ND EXAMPLE PROGRAM RESULTS

 The Bidiagonal Matrix J is

 (1,1) =  1.0000   (1,2) =  2.0000
 (2,2) =  2.0000   (2,3) =  3.0000
 (3,3) =  3.0000   (3,4) =  4.0000
 (4,4) =  4.0000   (4,5) =  5.0000
 (5,5) =  5.0000

 J has  3 singular values < =   5.0000

Return to index