?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 |
|
Array with a size of lda*n.
|
Input/Output |
lda |
Integer |
Leading dimension of matrix A. lda ≥ max(1, m). |
Input |
b |
|
Right-hand side matrix, with a size of ldb*nrhs.
|
Input/Output |
ldb |
Integer |
Leading dimension of matrix b. |
Input |
s |
|
Singular values of A in decreasing order. Condition number of A in the 2-norm = S(1)/S(min(m, n)) |
Output |
rcond |
|
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 |
|
On exit, if INFO=0, work(0) returns the optimal lwork. |
Output |
lwork |
Integer |
Size of work. |
Input |
rwork (only for complex types) |
|
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 |
|
Output |
Dependency
#include "klapack.h"
Example
C interface:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 | 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 :: trans = "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