## 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
```