?gels

使用矩阵的QR或LQ分解来求解超定或欠定线性方程组。

接口定义

C Interface：

sgels_(const char *trans, const int *m, const int *n, const int *nrhs, float *a, const int *lda, float *b, const int *ldb, float *work, const int *lwork, int *info);

dgels_(const char *trans, const int *m, const int *n, const int *nrhs, double *a, const int *lda, double *b, const int *ldb, double *work, const int *lwork, int *info);

cgels_(const char *trans, const int *m, const int *n, const int *nrhs, float _Complex *a, const int *lda, float _Complex *b, const int *ldb, float _Complex *work, const int *lwork, int *info);

zgels_(const char *trans, const int *m, const int *n, const int *nrhs, double _Complex *a, const int *lda, double _Complex *b, const int *ldb, double _Complex *work, const int *lwork, int *info);

Fortran Interface：

SGELS(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO);

DGELS(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO);

CGELS(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO);

ZGELS(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO);

参数

参数	类型	说明	输入/输出
trans	字符型	'N'：矩阵A不转置。 'C'：矩阵A转置。	输入
m	整数型	矩阵A的行数，m≥0。	输入
n	整数型	矩阵A的列数，n≥0。	输入
nrhs	整数型	右端项列数，nrhs≥0。	输入
a	在sgels中为单精度浮点型数组。在dgels中为双精度浮点型数组。在cgels中为单精度复数型数组。在zgels中为双精度复数型数组。	大小为ldan。进入时：为mn的矩阵。退出时：若m≥n，则为GEQRF返回的QR分解结果。若m<n，则为GELQF返回的LQ分解。	输入，输出
lda	整数型	矩阵A的主维，lda≥max(1, m)	输入
b	在sgels中为单精度浮点型数组。在dgels中为双精度浮点型数组。在cgels中为单精度复数型数组。在zgels中为双精度复数型数组。	右端矩阵，大小为ldbnrhs。进入时：若trans='N'，则为mnrhs的矩阵。若trans='C'或'T'，则为n*nrhs的矩阵。退出时：若info=0 若trans='N'，且m≥n，则矩阵的第1到n行包含了最小二乘向量；同时，每列的残差平方和由该列的第n+1到m个元素给出。若trans='N'，且m<n, 则矩阵的第1到n行包含了最小范数解向量。若trans='C'或'T'，且m>=n，则矩阵的第1到n行包含了最小范数解向量。若trans='C'或'T'，且m<n，则矩阵的第1到n行包含了最小二乘向量；同时，每列的残差平方和由该列的第n+1到m个元素给出。	输入，输出
ldb	整数型	矩阵b的主维。	输入
work	在sgels中为单精度浮点型数组。在dgels中为双精度浮点型数组。在cgels中为单精度复数型数组。在zgels中为双精度复数型数组。	工作数组。Info=0时，work(0)返回最优lwork大小。	输出
lwork	整数型	work的大小。	输入
Info	整数型	状态值：等于0：表示成功。小于0：info=-i，表示第i个参数非法。大于0：算法出错。	输出

依赖

#include "klapack.h"

示例

C Interface：

const char trans = 'N';
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 qwork;
int lwork = -1;
int info = 0;

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

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

dgels_(&trans, &m, &n, &nrhs, a, &lda, b, &ldb, work, &lwork, &info);
if (info != 0) {
    printf("Error, info = %d\n", info);
    return info;
}
/* output */
* a
* -120.698833     -108.770377     -57.958881      -100.842591
* 0.343169        75.924679       38.025488       -19.044303
* 0.335800        0.031671        7.685677        67.409117
* 0.145282        -0.313887       -0.556039       33.759113
* b
* 0.153366        0.083995
* -0.538653       0.096517
* 1.080726        -0.174382
* -0.145340       0.022261

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), ALLOCATABLE :: work(:)
REAL(8) :: qwork(1)
INTEGER :: lwork = -1
  
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 DGELS
CALL DGELS(trans, m, n, nrhs, a, lda, b, ldb, qwork, lwork, info)
IF (info.NE.0) THEN 
    CALL EXIT(1) 
END IF
lwork = INT(qwork(1))  
ALLOCATE(work(lwork))
CALL DGELS(trans, m, n, nrhs, a, lda, b, ldb, work, lwork, info)
IF (info.NE.0) THEN 
    CALL EXIT(1) 
END IF
DEALLOCATE(work);
* 
* Output: 
* 
a
* -120.698833     -108.770377     -57.958881      -100.842591
* 0.343169        75.924679       38.025488       -19.044303
* 0.335800        0.031671        7.685677        67.409117
* 0.145282        -0.313887       -0.556039       33.759113
* b
* 0.153366        0.083995
* -0.538653       0.096517
* 1.080726        -0.174382
* -0.145340       0.022261

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