**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.

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,*)

**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.

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.

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.

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'.

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.

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.

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.

[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.

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).

None

**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

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

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