?gelsd

使用分治法对线性最小二乘问题求解其最小范数解。

接口定义

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)

参数

参数	类型	说明	输入/输出
m	整数型	矩阵A的行数，m≥0。	输入
n	整数型	矩阵A的列数，n≥0。	输入
nrhs	整数型	右端项列数，nrhs≥0。	输入
a	在sgelsd中为单精度浮点型数组。在dgelsd中为双精度浮点型数组。在cgelsd中为单精度复数型数组。在zgelsd中为双精度复数型数组。	大小为ldan。进入时：为mn的矩阵。退出时，A被销毁。	输入，输出
lda	整数型	矩阵A的主维，lda≥max(1, m)。	输入
b	在sgelsd中为单精度浮点型数组。在dgelsd中为双精度浮点型数组。在cgelsd中为单精度复数型数组。在zgelsd中为双精度复数型数组。	右端矩阵，大小为ldbnrhs。进入时：mnrhs大小的右端项矩阵B。退出时，B被N*NRHS大小的解矩阵X覆盖。如果m≥n且RANK=n，则第i列中解的剩余平方和由该列中元素n+1:m的平方和给出。	输入，输出
ldb	整数型	矩阵b的主维。	输入
s	在sgelsd中为单精度浮点型数组。在dgelsd中为双精度浮点型数组。在cgelsd中为单精度复数型数组。在zgelsd中为双精度复数型数组。	A的奇异值按降序排列。 A的条件数二范数= S(1)/S(min(m,n))。	输出
rcond	在sgelsd中为单精度浮点型数。在dgelsd中为双精度浮点型数。在cgelsd中为单精度浮点型数。在zgelsd中为双精度浮点型数。	RCOND用于确定A的有效秩。奇异值S(i)≤RCOND*S(1)被视为零。如果RCOND<0，则使用机器精度。	输入
rank	整数型	A的有效秩，即大于RCOND*S(1)的奇异值的数量。	输出
work	在sgelsd中为单精度浮点型数组。在dgelsd中为双精度浮点型数组。在cgelsd中为单精度复数型数组。在zgelsd中为双精度复数型数组。	退出时，如果INFO=0，则work(0)返回最佳lwork。	输出
lwork	整数型	work的大小。	输入
rwork（复数特有）	在cgelsd中为单精度浮点型数组。在zgelsd中为双精度浮点型数组。	大小至少为max(1, 5*min(m, n))。	输出
iwork	整数型数组	维度是max(1,liwork), liwork=max(1, 3minmnnlvl + 11*minmn), minmn=min(m,n)。退出时，如果info=0,iwork(0)返回最佳liwork。	输出
info	整数型	=0：表示成功。 <0：info=-i，表示第i个参数非法。 >0：算法计算SVD时无法收敛，若info=i，则中间双对角线形式的i个对角线元素没有收敛到零。	输出

依赖

#include "klapack.h"

示例

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

父主题： 最小二乘问题函数