MB02UD

Minimum norm least squares solution of op(R) X = alpha B, or X op(R) = alpha B, with R upper triangular, using singular value decomposition

[Specification] [Arguments] [Method] [References] [Comments] [Example]

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

  None
Program Data
  None
Program 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