NAG FL Interface
f08mbf (dbdsvdx)

1 Purpose

f08mbf computes all or selected singular values and, optionally, the corresponding left and right singular vectors of a real n by n (upper or lower) bidiagonal matrix B.

2 Specification

Fortran Interface
Subroutine f08mbf ( uplo, jobz, range, n, d, e, vl, vu, il, iu, ns, s, z, ldz, work, iwork, info)
Integer, Intent (In) :: n, il, iu, ldz
Integer, Intent (Out) :: ns, iwork(12*n), info
Real (Kind=nag_wp), Intent (In) :: d(n), e(n-1), vl, vu
Real (Kind=nag_wp), Intent (Inout) :: z(ldz,*)
Real (Kind=nag_wp), Intent (Out) :: s(n), work(14*n)
Character (1), Intent (In) :: uplo, jobz, range
C Header Interface
#include <nag.h>
void  f08mbf_ (const char *uplo, const char *jobz, const char *range, const Integer *n, const double d[], const double e[], const double *vl, const double *vu, const Integer *il, const Integer *iu, Integer *ns, double s[], double z[], const Integer *ldz, double work[], Integer iwork[], Integer *info, const Charlen length_uplo, const Charlen length_jobz, const Charlen length_range)
The routine may be called by the names f08mbf, nagf_lapackeig_dbdsvdx or its LAPACK name dbdsvdx.

3 Description

f08mbf computes the singular value decomposition (SVD) of a real n by n (upper or lower) bidiagonal matrix B as
B = USVT ,  
where S is a diagonal matrix with non-negative diagonal elements (the singular values of B), and U and VT are orthogonal matrices. The columns of U and V are the left and right singular vectors of B, respectively.
Given an upper bidiagonal matrix B with diagonal d= d1 d2 dn and superdiagonal e= e1 e2 eN-1 , f08mbf computes the singular value decomposition of B through the eigenvalues and eigenvectors of the n×2 by n×2 tridiagonal matrix
TGK= 0 d1 d1 0 e1 e1 0 d2 d2 . . . . .  
If s,u,v is a singular triplet of B with u=v=1, then s,q1 and -s,q2, q1=q2=1, are eigenpairs of TGK, with q1= v1,u1,v2,u2,,vn,un / 2, and q2= -v1,u1,-v2,u2,,-vn,un / 2.
Given a TGK matrix, you can either
  1. (i)compute -s,-v and change signs so that the singular values (and corresponding vectors) are already in descending order (as in f08kbf) or
  2. (ii)compute s,v and reorder the values (and corresponding vectors).
f08mbf implements (i) by calling f08jbf (bisection plus inverse iteration, to be replaced with a version of the Multiple Relative Robust Representation algorithm. (See Williams and Lang (2013).)
Alternative to computing all singular values of B, a selected set can be computed. The set is either those singular values lying in a given interval, σvl,vu, or those whose index (counting from largest to smallest in magnitude) lies in a given range 1il,,iun. In these cases, the corresponding left and right singular vectors can optionally be computed.

4 References

Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia https://www.netlib.org/lapack/lug
Williams P and Lang B (2013) A framework for the MR3 Algorithm: theory and implementation SIAM J. Sci. Comput. 35 740–766

5 Arguments

1: uplo Character(1) Input
On entry: indicates whether B is upper or lower bidiagonal.
uplo='U'
B is upper bidiagonal.
uplo='L'
B is lower bidiagonal.
Constraint: uplo='U' or 'L'.
2: jobz Character(1) Input
On entry: indicates whether singular vectors are computed.
jobz='N'
Only singular values are computed.
jobz='V'
Singular values and singular vectors are computed.
Constraint: jobz='N' or 'V'.
3: range Character(1) Input
On entry: indicates which singular values should be returned.
range='A'
All singular values will be found.
range='V'
All singular values in the half-open interval vl,vu will be found.
range='I'
The ilth through iuth singular values will be found.
Constraint: range='A', 'V' or 'I'.
4: n Integer Input
On entry: n, the order of the bidiagonal matrix B.
Constraint: n0.
5: dn Real (Kind=nag_wp) array Input
On entry: the diagonal elements d of the bidiagonal matrix B.
6: en-1 Real (Kind=nag_wp) array Input
On entry: the n-1 off-diagonal elements e of the bidiagonal matrix B.
7: vl Real (Kind=nag_wp) Input
On entry: if range='V', the lower bound of the interval to be searched for singular values.
If range='A' or 'I', vl is not referenced.
Constraint: if range='V', 0.0vl.
8: vu Real (Kind=nag_wp) Input
On entry: if range='V', the upper bound of the interval to be searched for singular values.
If range='A' or 'I', vu is not referenced.
Constraint: if range='V', vl<vu.
9: il Integer Input
10: iu Integer Input
On entry: if range='I', il and iu specify the indices (in ascending order) of the smallest and largest singular values to be returned, respectively.
If range='A' or 'V', il and iu are not referenced.
Constraints:
  • if range='I' and n=0, il=1 and iu=0;
  • if range='I' and n>0, 1 il iu n .
11: ns Integer Output
On exit: the total number of singular values found. 0nsn.
If range='A', ns=n.
If range='I', ns=iu-il+1.
12: sn Real (Kind=nag_wp) array Output
On exit: the first ns elements contain the selected singular values in ascending order.
13: zldz* Real (Kind=nag_wp) array Output
Note: the second dimension of the array z must be at least minn,ns+1 if jobz='V', and at least 1 otherwise.
On exit: if jobz='V', then if info=0 the first ns columns of z contain the singular vectors of the matrix B corresponding to the selected singular values, with U in rows 1 to n and V in rows n+1 to n×2, i.e.,
Z= U V .  
If jobz='N', then z is not referenced.
Note: the user must ensure that at least K=minn,ns+1 columns are supplied in the array Z. If range='V', the exact value of ns is not known in advance and an upper bound of at least n must be used.
14: ldz Integer Input
On entry: the first dimension of the array z as declared in the (sub)program from which f08mbf is called.
Constraints:
  • if jobz='V', ldz max2,n×2 ;
  • otherwise ldz1.
15: work14×n Real (Kind=nag_wp) array Workspace
16: iwork12×n Integer array Workspace
On exit: if jobz='V', then
  • if info=0, the first ns elements of iwork are zero;
  • if info>0, iwork contains the indices of the eigenvectors that failed to converge in f08jbf.
If jobz='N', iwork is not referenced.
17: info Integer Output
On exit: info=0 unless the routine detects an error (see Section 6).

6 Error Indicators and Warnings

info<0
If info=-i, argument i had an illegal value. An explanatory message is output, and execution of the program is terminated.
info>0
The algorithm failed to converge; value eigenvectors of the associated eigenproblem did not converge. Their indices are stored in array iwork.

7 Accuracy

Each computed singular value of B is accurate to nearly full relative precision, no matter how tiny the singular value. The ith computed singular value, s^i, satisfies the bound
s^i-si pnεsi  
where ε is the machine precision and pn is a modest function of n.
For bounds on the computed singular vectors, see Section 4.9.1 of Anderson et al. (1999). See also f08flf.

8 Parallelism and Performance

f08mbf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
f08mbf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

If only singular values are required, the total number of floating-point operations is approximately proportional to n2. When singular vectors are required the number of operations is bounded above by approximately the same number of operations as f08mef, but for large matrices f08mbf is usually much faster.
There is no complex analogue of f08mbf.

10 Example

This example computes the singular value decomposition of the upper bidiagonal matrix
B = 3.62 1.26 0.00 0.00 0.00 -2.41 -1.53 0.00 0.00 0.00 1.92 1.19 0.00 0.00 0.00 -1.43 .  

10.1 Program Text

Program Text (f08mbfe.f90)

10.2 Program Data

Program Data (f08mbfe.d)

10.3 Program Results

Program Results (f08mbfe.r)