Documentation
Rate This Document
Findability
Accuracy
Completeness
Readability
SubmitClose

?gelss

Compute the minimum-norm solution to a linear least squares problem using SVD.

Interface Definition

C interface:

sgelss_(const int *m, const int *n, const int *nrhs, float *a, const int *lda, float *b, const int *ldb, float *s, const float *rcond, int *rank, float *work, const int *lwork, float *rwork, int *info);

dgelss_( const int *m, const int *n, const int *nrhs, double *a, const int *lda, double *b, const int *ldb, double *s, const double *rcond, int *rank, double *work, const int *lwork, double *rwork, int *info);

cgelss_(const int *m, const int *n, const int *nrhs, float _Complex *a, const int *lda, float _Complex *b, const int *ldb, float *s, const float *rcond, int *rank, float _Complex *work, const int *lwork, float *rwork, int *info);

zgelss_(const int *m, const int *n, const int *nrhs, double _Complex *a, const int *lda, double _Complex *b, const int *ldb, double *s, const double *rcond, int *rank, double _Complex *work, const int *lwork, double *rwork, int *info);

Fortran interface:

SGELSS(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO)

DGELSS(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO)

CGELSS(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO)

ZGELSS(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO)

Parameters

Parameters

Type

Description

Input/Output

m

Integer

Number of rows in matrix A. m ≥ 0.

Input

n

Integer

Number of columns in matrix A, n ≥ 0.

Input

nrhs

Integer

Number of right-hand side columns, nrhs ≥ 0.

Input

a
  • A single-precision floating-point array for sgelss
  • A double-precision floating-point array for dgelss
  • A complex single-precision array for cgelss
  • A complex double-precision array for zgelss

Array with a size of lda*n.

  • On entry: an m*n matrix.
  • On exit: The first min(m, n) rows of A are overwritten by its right singular vectors, stored row-wise.

Input/Output

lda

Integer

Leading dimension of matrix A. lda ≥ max(1, m).

Input

b
  • A single-precision floating-point array for sgelss
  • A double-precision floating-point array for dgelss
  • A complex single-precision array for cgelss
  • A complex double-precision array for zgelss

Right-hand side matrix, with a size of ldb*nrhs.

  • On entry: right-hand side matrix B with a size of m*nrhs.
  • On exit: B is overwritten by the N*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 modulus of elements n+1:m in that column.

Input/Output

ldb

Integer

Leading dimension of matrix b.

Input

s
  • A single-precision floating-point array for sgelss
  • A double-precision floating-point array for dgelss
  • A complex single-precision array for cgelss
  • A complex double-precision array for zgelss

Singular values of A in decreasing order.

Condition number of A in the 2-norm = S(1)/S(min(m, n))

Output

rcond
  • A single-precision floating-point number for sgelss
  • A double-precision floating-point number for dgelss
  • A single-precision floating-point number for cgelss
  • A double-precision floating-point number for zgelss

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.

Input

rank

Integer

Effective rank of A, that is, the number of singular values which are greater than RCOND*S(1).

Output

work
  • A single-precision floating-point array for sgelss
  • A double-precision floating-point array for dgelss
  • A complex single-precision array for cgelss
  • A complex double-precision array for zgelss

On exit, if INFO=0, work(0) returns the optimal lwork.

Output

lwork

Integer

Size of work.

Input

rwork (only available for the complex type)
  • A single-precision floating-point array for cgelss
  • A double-precision floating-point array for zgelss

The size must be at least max(1, 5*min(m, n)).

Output

info

Integer
  • 0: The execution is successful.
  • Smaller than 0: If info = -i, the i-th parameter has an illegal value.
  • Greater than 0: The algorithm for computing the SVD failed to converge. If info=i, i indicates the number of diagonal elements (of an intermediate bidiagonal form) which did not converge to zero.

Output

Dependency

#include "klapack.h"

Example

C interface:

const int m = 4;
const int n = 4;
const int nrhs = 2;
const int lda = 4;
const int ldb = 4;

double a[] = {72.1673, 66.1857, 64.7644, 28.0199, 91.4151,
              6.5180,  62.8483, 72.4323, 46.5760, 8.6928,
              28.9821, 42.1828, 18.6437, 99.8612, 35.6972,
              67.9812, 5.0880,  85.5035, 79.2945, 54.5920,
              28.6869, 49.7512, 7.5186,  28.6929, 84.6041};

double b[] = {9.4532, 1.5204, 2.2127, 0.9891, 7.1778,
              6.8955, 7.2465, 3.5019, 8.2268, 3.5287};

double *s = (double*)malloc(n * sizeof(double));
double rcond = -1.0;
int rank;

double qwork;
int lwork = -1;
int info = 0;

dgelss_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, &qwork, &lwork, &info);
if (info != 0) {
    printf("Error, info = %d\n", info);
    return info;
}

lwork = (int)qwork;
double *work = (double*)malloc(lwork * sizeof(double));

dgelss_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, work, &lwork, &info);
if (info != 0) {
    printf("Error, info = %d\n", info);
    return info;
}

/* output */
* a
* -0.552769       -0.578412       -0.312159       -0.512292
* -0.073442       0.602089        0.251939        -0.754070
* -0.827844       0.298713        0.251050        0.403012
* -0.061050       0.462277        -0.880941       0.080724
* b
* 0.153366        0.083995
* -0.538653       0.096517
* 1.080726        -0.174382
* -0.145340       0.022261
* rank
* 4

Fortran interface:

CHARACTER :: tarns = "N"
PARAMETER (n = 4) 
PARAMETER (m = 4)
PARAMETER (lda = 4) 
PARAMETER (ldb = 4)  
PARAMETER (nrhs = 2) 

INTEGER :: info = 0 
REAL(8) :: a(lda, n) 
REAL(8) :: b(ldb, nrhs)
REAL(8) :: s(n)
REAL(8), ALLOCATABLE :: work(:)
REAL(8) :: rcond = -1.0
REAL(8) :: qwork(1)
INTEGER :: lwork = -1
INTEGER :: rank
  
DATA a / 72.1673, 66.1857, 64.7644, 28.0199, 91.4151,
&         6.5180,  62.8483, 72.4323, 46.5760, 8.6928,
&        28.9821, 42.1828, 18.6437, 99.8612, 35.6972,
&        67.9812, 5.0880,  85.5035, 79.2945, 54.5920,
&        28.6869, 49.7512, 7.5186,  28.6929, 84.6041 / 
DATA b / 9.4532, 1.5204, 2.2127, 0.9891, 7.1778,
&        6.8955, 7.2465, 3.5019, 8.2268, 3.5287 /
EXTERNAL DGELSS
CALL DGELSS(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, qwork, lwork, info)
IF (info.NE.0) THEN 
    CALL EXIT(1) 
END IF
lwork = INT(qwork(1))  
ALLOCATE(work(lwork))
CALL DGELSS(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, info)
IF (info.NE.0) THEN 
    CALL EXIT(1) 
END IF
DEALLOCATE(work);
* 
* Output: 
* a
* -0.552769       -0.578412       -0.312159       -0.512292
* -0.073442       0.602089        0.251939        -0.754070
* -0.827844       0.298713        0.251050        0.403012
* -0.061050       0.462277        -0.880941       0.080724
* b
* 0.153366        0.083995
* -0.538653       0.096517
* 1.080726        -0.174382
* -0.145340       0.022261
* rank
* 4