Purpose
To compute the minimum norm least squares solution of one of the following linear systems op(R)*X = alpha*B, (1) X*op(R) = alpha*B, (2) where alpha is a real scalar, op(R) is either R or its transpose, R', R is an L-by-L real upper triangular matrix, B is an M-by-N real matrix, and L = M for (1), or L = N for (2). Singular value decomposition, R = Q*S*P', is used, assuming that R is rank deficient.Specification
SUBROUTINE MB02UD( FACT, SIDE, TRANS, JOBP, M, N, ALPHA, RCOND, $ RANK, R, LDR, Q, LDQ, SV, B, LDB, RP, LDRP, $ DWORK, LDWORK, INFO ) C .. Scalar Arguments .. CHARACTER FACT, JOBP, SIDE, TRANS INTEGER INFO, LDB, LDQ, LDR, LDRP, LDWORK, M, N, RANK DOUBLE PRECISION ALPHA, RCOND C .. Array Arguments .. DOUBLE PRECISION B(LDB,*), DWORK(*), Q(LDQ,*), R(LDR,*), $ RP(LDRP,*), SV(*)Arguments
Mode Parameters
FACT CHARACTER*1 Specifies whether R has been previously factored or not, as follows: = 'F': R has been factored and its rank and singular value decomposition, R = Q*S*P', are available; = 'N': R has not been factored and its singular value decomposition, R = Q*S*P', should be computed. SIDE CHARACTER*1 Specifies whether op(R) appears on the left or right of X as follows: = 'L': Solve op(R)*X = alpha*B (op(R) is on the left); = 'R': Solve X*op(R) = alpha*B (op(R) is on the right). TRANS CHARACTER*1 Specifies the form of op(R) to be used as follows: = 'N': op(R) = R; = 'T': op(R) = R'; = 'C': op(R) = R'. JOBP CHARACTER*1 Specifies whether or not the pseudoinverse of R is to be computed or it is available as follows: = 'P': Compute pinv(R), if FACT = 'N', or use pinv(R), if FACT = 'F'; = 'N': Do not compute or use pinv(R).Input/Output Parameters
M (input) INTEGER The number of rows of the matrix B. M >= 0. N (input) INTEGER The number of columns of the matrix B. N >= 0. ALPHA (input) DOUBLE PRECISION The scalar alpha. When alpha is zero then B need not be set before entry. RCOND (input) DOUBLE PRECISION RCOND is used to determine the effective rank of R. Singular values of R satisfying Sv(i) <= RCOND*Sv(1) are treated as zero. If RCOND <= 0, then EPS is used instead, where EPS is the relative machine precision (see LAPACK Library routine DLAMCH). RCOND <= 1. RCOND is not used if FACT = 'F'. RANK (input or output) INTEGER The rank of matrix R. RANK is an input parameter when FACT = 'F', and an output parameter when FACT = 'N'. L >= RANK >= 0. R (input/output) DOUBLE PRECISION array, dimension (LDR,L) On entry, if FACT = 'F', the leading L-by-L part of this array must contain the L-by-L orthogonal matrix P' from singular value decomposition, R = Q*S*P', of the matrix R; if JOBP = 'P', the first RANK rows of P' are assumed to be scaled by inv(S(1:RANK,1:RANK)). On entry, if FACT = 'N', the leading L-by-L upper triangular part of this array must contain the upper triangular matrix R. On exit, if INFO = 0, the leading L-by-L part of this array contains the L-by-L orthogonal matrix P', with its first RANK rows scaled by inv(S(1:RANK,1:RANK)), when JOBP = 'P'. LDR INTEGER The leading dimension of array R. LDR >= MAX(1,L). Q (input or output) DOUBLE PRECISION array, dimension (LDQ,L) On entry, if FACT = 'F', the leading L-by-L part of this array must contain the L-by-L orthogonal matrix Q from singular value decomposition, R = Q*S*P', of the matrix R. If FACT = 'N', this array need not be set on entry, and on exit, if INFO = 0, the leading L-by-L part of this array contains the orthogonal matrix Q. LDQ INTEGER The leading dimension of array Q. LDQ >= MAX(1,L). SV (input or output) DOUBLE PRECISION array, dimension (L) On entry, if FACT = 'F', the first RANK entries of this array must contain the reciprocal of the largest RANK singular values of the matrix R, and the last L-RANK entries of this array must contain the remaining singular values of R sorted in descending order. If FACT = 'N', this array need not be set on input, and on exit, if INFO = 0, the first RANK entries of this array contain the reciprocal of the largest RANK singular values of the matrix R, and the last L-RANK entries of this array contain the remaining singular values of R sorted in descending order. B (input/output) DOUBLE PRECISION array, dimension (LDB,N) On entry, if ALPHA <> 0, the leading M-by-N part of this array must contain the matrix B. On exit, if INFO = 0 and RANK > 0, the leading M-by-N part of this array contains the M-by-N solution matrix X. LDB INTEGER The leading dimension of array B. LDB >= MAX(1,M). RP (input or output) DOUBLE PRECISION array, dimension (LDRP,L) On entry, if FACT = 'F', JOBP = 'P', and RANK > 0, the leading L-by-L part of this array must contain the L-by-L matrix pinv(R), the Moore-Penrose pseudoinverse of R. On exit, if FACT = 'N', JOBP = 'P', and RANK > 0, the leading L-by-L part of this array contains the L-by-L matrix pinv(R), the Moore-Penrose pseudoinverse of R. If JOBP = 'N', this array is not referenced. LDRP INTEGER The leading dimension of array RP. LDRP >= MAX(1,L), if JOBP = 'P'. LDRP >= 1, if JOBP = 'N'.Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK) On exit, if INFO = 0, DWORK(1) returns the optimal LDWORK; if INFO = i, 1 <= i <= L, then DWORK(2:L) contain the unconverged superdiagonal elements of an upper bidiagonal matrix D whose diagonal is in SV (not necessarily sorted). D satisfies R = Q*D*P', so it has the same singular values as R, and singular vectors related by Q and P'. LDWORK INTEGER The length of the array DWORK. LDWORK >= MAX(1,L), if FACT = 'F'; LDWORK >= MAX(1,5*L), if FACT = 'N'. For optimum performance LDWORK should be larger than MAX(1,L,M*N), if FACT = 'F'; MAX(1,5*L,M*N), if FACT = 'N'.Error Indicator
INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value; > 0: if INFO = i, i = 1:L, the SVD algorithm has failed to converge. In this case INFO specifies how many superdiagonals did not converge (see the description of DWORK); this failure is not likely to occur.Method
The L-by-L upper triangular matrix R is factored as R = Q*S*P', if FACT = 'N', using SLICOT Library routine MB03UD, where Q and P are L-by-L orthogonal matrices and S is an L-by-L diagonal matrix with non-negative diagonal elements, SV(1), SV(2), ..., SV(L), ordered decreasingly. Then, the effective rank of R is estimated, and matrix (or matrix-vector) products and scalings are used to compute X. If FACT = 'F', only matrix (or matrix-vector) products and scalings are performed.Further Comments
Option JOBP = 'P' should be used only if the pseudoinverse is really needed. Usually, it is possible to avoid the use of pseudoinverse, by computing least squares solutions. The routine uses BLAS 3 calculations if LDWORK >= M*N, and BLAS 2 calculations, otherwise. No advantage of any additional workspace larger than L is taken for matrix products, but the routine can be called repeatedly for chunks of columns of B, if LDWORK < M*N.Example
Program Text
NoneProgram Data
NoneProgram Results
None
Click here to get a compressed (gzip) tar file containing the source code of the routine, the example program, data, documentation, and related files.
Return to index