Contents
dgelsd - compute the minimum-norm solution to a real linear
least squares problem
SUBROUTINE DGELSD(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK,
LWORK, IWORK, INFO)
INTEGER M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
INTEGER IWORK(*)
DOUBLE PRECISION RCOND
DOUBLE PRECISION A(LDA,*), B(LDB,*), S(*), WORK(*)
SUBROUTINE DGELSD_64(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
WORK, LWORK, IWORK, INFO)
INTEGER*8 M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
INTEGER*8 IWORK(*)
DOUBLE PRECISION RCOND
DOUBLE PRECISION A(LDA,*), B(LDB,*), S(*), WORK(*)
F95 INTERFACE
SUBROUTINE GELSD([M], [N], [NRHS], A, [LDA], B, [LDB], S, RCOND,
RANK, [WORK], [LWORK], [IWORK], [INFO])
INTEGER :: M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
INTEGER, DIMENSION(:) :: IWORK
REAL(8) :: RCOND
REAL(8), DIMENSION(:) :: S, WORK
REAL(8), DIMENSION(:,:) :: A, B
SUBROUTINE GELSD_64([M], [N], [NRHS], A, [LDA], B, [LDB], S, RCOND,
RANK, [WORK], [LWORK], [IWORK], [INFO])
INTEGER(8) :: M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
INTEGER(8), DIMENSION(:) :: IWORK
REAL(8) :: RCOND
REAL(8), DIMENSION(:) :: S, WORK
REAL(8), DIMENSION(:,:) :: A, B
C INTERFACE
#include <sunperf.h>
void dgelsd(int m, int n, int nrhs, double *a, int lda, dou-
ble *b, int ldb, double *s, double rcond, int
*rank, int *info);
void dgelsd_64(long m, long n, long nrhs, double *a, long
lda, double *b, long ldb, double *s, double rcond,
long *rank, long *info);
dgelsd computes the minimum-norm solution to a real linear
least squares problem:
minimize 2-norm(| b - A*x |)
using the singular value decomposition (SVD) of A. A is an
M-by-N matrix which may be rank-deficient.
Several right hand side vectors b and solution vectors x can
be handled in a single call; they are stored as the columns
of the M-by-NRHS right hand side matrix B and the N-by-NRHS
solution matrix X.
The problem is solved in three steps:
(1) Reduce the coefficient matrix A to bidiagonal form with
Householder transformations, reducing the original prob-
lem
into a "bidiagonal least squares problem" (BLS)
(2) Solve the BLS using a divide and conquer approach.
(3) Apply back all the Householder tranformations to solve
the original least squares problem.
The effective rank of A is determined by treating as zero
those singular values which are less than RCOND times the
largest singular value.
The divide and conquer algorithm makes very mild assumptions
about floating point arithmetic. It will work on machines
with a guard digit in add/subtract, or on those binary
machines without guard digits which subtract like the Cray
X-MP, Cray Y-MP, Cray C-90, or Cray-2. It could conceivably
fail on hexadecimal or decimal machines without guard
digits, but we know of none.
M (input) The number of rows of A. M >= 0.
N (input) The number of columns of A. N >= 0.
NRHS (input)
The number of right hand sides, i.e., the number
of columns of the matrices B and X. NRHS >= 0.
A (input/output)
On entry, the M-by-N matrix A. On exit, A has
been destroyed.
LDA (input)
The leading dimension of the array A. LDA >=
max(1,M).
B (input/output)
On entry, the M-by-NRHS right hand side matrix B.
On exit, B is overwritten by the N-by-NRHS solu-
tion matrix X. If m >= n and RANK = n, the resi-
dual sum-of-squares for the solution in the i-th
column is given by the sum of squares of elements
n+1:m in that column.
LDB (input)
The leading dimension of the array B. LDB >=
max(1,max(M,N)).
S (output)
The singular values of A in decreasing order. The
condition number of A in the 2-norm =
S(1)/S(min(m,n)).
RCOND (input)
RCOND is used to determine the effective rank of
A. Singular values S(i) <= RCOND*S(1) are treated
as zero. If RCOND < 0, machine precision is used
instead.
RANK (output)
The effective rank of A, i.e., the number of
singular values which are greater than RCOND*S(1).
WORK (workspace)
On exit, if INFO = 0, WORK(1) returns the optimal
LWORK.
LWORK (input)
The dimension of the array WORK. LWORK >= 1. The
exact minimum amount of workspace needed depends
on M, N and NRHS. As long as LWORK is at least
12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS *
(SMLSIZ+1)**2, if M is greater than or equal to N
or 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS +
(SMLSIZ+1)**2, if M is less than N, the code will
execute correctly. SMLSIZ is returned by ILAENV
and is equal to the maximum size of the subprob-
lems at the bottom of the computation tree (usu-
ally about 25), and NLVL = INT( LOG_2( MIN( M,N
)/(SMLSIZ+1) ) ) + 1 For good performance, LWORK
should generally be larger.
If LWORK = -1, then a workspace query is assumed;
the routine only calculates the optimal size of
the WORK array, returns this value as the first
entry of the WORK array, and no error message
related to LWORK is issued by XERBLA.
IWORK (workspace)
LIWORK >= 3 * MINMN * NLVL + 11 * MINMN, where
MINMN = MIN( M,N ).
INFO (output)
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an ille-
gal value.
> 0: the algorithm for computing the SVD failed
to converge; if INFO = i, i off-diagonal elements
of an intermediate bidiagonal form did not con-
verge to zero.
Based on contributions by
Ming Gu and Ren-Cang Li, Computer Science Division,
University of California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA