DGELSD(3S)DGELSD(3S)NAMEDGELSD - compute the minimum-norm solution to a real linear least squares
problem
SYNOPSIS
SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK,
LWORK, IWORK, INFO )
INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
DOUBLE PRECISION RCOND
INTEGER IWORK( * )
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
IMPLEMENTATION
These routines are part of the SCSL Scientific Library and can be loaded
using either the -lscs or the -lscs_mp option. The -lscs_mp option
directs the linker to use the multi-processor version of the library.
When linking to SCSL with -lscs or -lscs_mp, the default integer size is
4 bytes (32 bits). Another version of SCSL is available in which integers
are 8 bytes (64 bits). This version allows the user access to larger
memory sizes and helps when porting legacy Cray codes. It can be loaded
by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
only one of the two versions; 4-byte integer and 8-byte integer library
calls cannot be mixed.
PURPOSEDGELSD 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 problem
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
Page 1
DGELSD(3S)DGELSD(3S)
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.
ARGUMENTS
M (input) INTEGER
The number of rows of A. M >= 0.
N (input) INTEGER
The number of columns of A. N >= 0.
NRHS (input) INTEGER
The number of right hand sides, i.e., the number of columns of
the matrices B and X. NRHS >= 0.
A (input) DOUBLE PRECISION array, dimension (LDA,N)
On entry, the M-by-N matrix A. On exit, A has been destroyed.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
On entry, the M-by-NRHS right hand side matrix B. On exit, B is
overwritten by the N-by-NRHS solution matrix X. If m >= n and
RANK = n, the residual 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) INTEGER
The leading dimension of the array B. LDB >= max(1,max(M,N)).
S (output) DOUBLE PRECISION array, dimension (min(M,N))
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) DOUBLE PRECISION
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) INTEGER
The effective rank of A, i.e., the number of singular values
which are greater than RCOND*S(1).
WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
LWORK (input) INTEGER
The dimension of the array WORK. LWORK must be at least 1. The
exact minimum amount of workspace needed depends on M, N and
Page 2
DGELSD(3S)DGELSD(3S)
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 subproblems at
the bottom of the computation tree (usually about 25), and NLVL =
MAX( 0, 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) INTEGER array, dimension (LIWORK)
LIWORK >= 3 * MINMN * NLVL + 11 * MINMN, where MINMN = MIN( M,N
).
INFO (output) INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal 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 converge to zero.
FURTHER DETAILS
Based on contributions by
Ming Gu and Ren-Cang Li, Computer Science Division, University of
California at Berkeley, USA
Osni Marques, LBNL/NERSC, USA
SEE ALSOINTRO_LAPACK(3S), INTRO_SCSL(3S)
This man page is available only online.
Page 3
