Contents
cgelsd - compute the minimum-norm solution to a real linear
least squares problem
SUBROUTINE CGELSD(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK,
LWORK, RWORK, IWORK, INFO)
COMPLEX A(LDA,*), B(LDB,*), WORK(*)
INTEGER M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
INTEGER IWORK(*)
REAL RCOND
REAL S(*), RWORK(*)
SUBROUTINE CGELSD_64(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
WORK, LWORK, RWORK, IWORK, INFO)
COMPLEX A(LDA,*), B(LDB,*), WORK(*)
INTEGER*8 M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
INTEGER*8 IWORK(*)
REAL RCOND
REAL S(*), RWORK(*)
F95 INTERFACE
SUBROUTINE GELSD([M], [N], [NRHS], A, [LDA], B, [LDB], S, RCOND,
RANK, [WORK], [LWORK], [RWORK], [IWORK], [INFO])
COMPLEX, DIMENSION(:) :: WORK
COMPLEX, DIMENSION(:,:) :: A, B
INTEGER :: M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
INTEGER, DIMENSION(:) :: IWORK
REAL :: RCOND
REAL, DIMENSION(:) :: S, RWORK
SUBROUTINE GELSD_64([M], [N], [NRHS], A, [LDA], B, [LDB], S, RCOND,
RANK, [WORK], [LWORK], [RWORK], [IWORK], [INFO])
COMPLEX, DIMENSION(:) :: WORK
COMPLEX, DIMENSION(:,:) :: A, B
INTEGER(8) :: M, N, NRHS, LDA, LDB, RANK, LWORK, INFO
INTEGER(8), DIMENSION(:) :: IWORK
REAL :: RCOND
REAL, DIMENSION(:) :: S, RWORK
C INTERFACE
#include <sunperf.h>
void cgelsd(int m, int n, int nrhs, complex *a, int lda,
complex *b, int ldb, float *s, float rcond, int
*rank, int *info);
void cgelsd_64(long m, long n, long nrhs, complex *a, long
lda, complex *b, long ldb, float *s, float rcond,
long *rank, long *info);
cgelsd 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 tranformations, 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 the matrix A. M >= 0.
N (input) The number of columns of the matrix 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,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. If M >= N, LWORK >= 2*N +
N*NRHS. If M < N, LWORK >= 2*M + M*NRHS. 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.
RWORK (workspace)
If M >= N, LRWORK >= 8*N + 2*N*SMLSIZ + 8*N*NLVL +
N*NRHS. If M < N, LRWORK >= 8*M + 2*M*SMLSIZ +
8*M*NLVL + M*NRHS. 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
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