文档
我要评分
获取效率
正确性
完整性
易理解
提交关闭

?gelsd

Compute the minimum-norm solution to a linear least squares problem using a divide and conquer method.

Interface Definition

C interface:

sgelsd_(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 *iwork, int *info);

dgelsd_(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 *iwork, int *info);

cgelsd_(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 *iwork, int *info);

zgelsd_(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 *iwork, int *info);

Fortran interface:

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

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

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

ZGELSD(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, IWORK, 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 sgelsd
  • A double-precision floating-point array for dgelsd
  • A complex single-precision array for cgelsd
  • A complex double-precision array for zgelsd

Array with a size of lda*n.

  • On entry: an m*n matrix.
  • On exit: A is destroyed.

Input/Output

lda

Integer

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

Input

b
  • A single-precision floating-point array for sgelsd
  • A double-precision floating-point array for dgelsd
  • A complex single-precision array for cgelsd
  • A complex double-precision array for zgelsd

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 sgelsd
  • A double-precision floating-point array for dgelsd
  • A complex single-precision array for cgelsd
  • A complex double-precision array for zgelsd

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 sgelsd
  • A double-precision floating-point number for dgelsd
  • A single-precision floating-point number for cgelsd
  • A double-precision floating-point number for zgelsd

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 sgelsd
  • A double-precision floating-point array for dgelsd
  • A complex single-precision array for cgelsd
  • A complex double-precision array for zgelsd

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 cgelsd
  • A double-precision floating-point array for zgelsd

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

Output

iwork

Integer array

The dimension is max(1, liwork). liwork=max(1, 3*minmn*nlvl + 11*minmn), minmn=min(m, n). On exit, if info=0, iwork(0) returns the optimal liwork.

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));
int *iwork = (int*)malloc((3 * n * n + 11 * n) * sizeof(int));
double rcond = -1.0;
int rank;

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

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

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

dgelsd_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, work, &lwork, iwork, &info);
if (info != 0) {
    printf("Error, info = %d\n", info);
    return info;
}
/* output */
* a
* -120.698833     159.246523      0.216251        0.376255
* 0.343169        -73.504609      -11.813549      -0.661117
* 0.335800        0.090502        84.406748       25.275449
* 0.145282        -0.806922       -0.918078       -3.175151
* 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)
INTEGER :: iwork(3*n*n+11*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 DGELSD
CALL DGELSD(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, qwork, lwork, iwork, info)
IF (info.NE.0) THEN 
    CALL EXIT(1) 
END IF
lwork = INT(qwork(1))  
ALLOCATE(work(lwork))
CALL DGELSD(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, qwork, lwork, iwork, info)
IF (info.NE.0) THEN 
    CALL EXIT(1) 
END IF
DEALLOCATE(work);
* 
* Output: 
* a
* -120.698833     159.246523      0.216251        0.376255
* 0.343169        -73.504609      -11.813549      -0.661117
* 0.335800        0.090502        84.406748       25.275449
* 0.145282        -0.806922       -0.918078       -3.175151
* b
* 0.153366        0.083995
* -0.538653       0.096517
* 1.080726        -0.174382
* -0.145340       0.022261
* rank
* 4