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
NoneExample
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.0Program 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