?gesvd
一般矩阵SVD分解。
接口定义
C Interface：
sgesvd_(const char *jobu, const char *jobvt, const int *m, const int *n, float *a, const int *lda, float *s, float *u, const int *ldu, float *vt,
const int *ldvt, float *work, const int *lwork, float *rwork, int *info);
dgesvd_(const char *jobu, const char *jobvt, const int *m, const int *n, double *a, const int *lda, double *s, double *u, const int *ldu, double *vt,
const int *ldvt, double *work, const int *lwork, double *rwork, int *info);
cgesvd_(const char *jobu, const char *jobvt, const int *m, const int *n, float _Complex *a, const int *lda, float *s, float _Complex *u, const int *ldu, float _Complex *vt,
const int *ldvt, float _Complex *work, const int *lwork, float *rwork, int *info);
zgesvd_(const char *jobu, const char *jobvt, const int *m, const int *n, double _Complex *a, const int *lda, double *s, double _Complex *u, const int *ldu, double _Complex *vt,
const int *ldvt, double _Complex *work, const int *lwork, double *rwork, int *info);
Fortran Interface：
SGESVD(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO);
DGESVD(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO);
CGESVD(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO);
ZGESVD(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO);
参数
参数 |
类型 |
说明 |
输入/输出 |
|---|---|---|---|
jobu |
字符型 |
|
输入 |
jobvt |
字符型 |
|
输入 |
m |
整数型 |
矩阵a的行数，m≥0。 |
输入 |
n |
整数型 |
矩阵a的列数，n≥0。 |
输入 |
a |
|
大小为lda*n。 输入时：矩阵a为m*n的矩阵。 输出时：
|
输入，输出 |
lda |
|
矩阵a的主维，lda≥max(1,m)。 |
输入 |
s |
|
大小为min(m, n)。 输出为矩阵a的奇异值，且按s(i)>=s(i+1)的规则排序。 |
输出 |
u |
|
大小为ldu*ucol，其中当jobu='A'时，ucol=m; jobu='S'时，ucol=min(m,n)。 若jobu='A': 其为m*m的酉矩阵u。 若jobu='S': 其包含了U矩阵(左奇异值向量，列优先存储)的前min(m,n)列。 若jobu='N' 或 'O': 其不被使用。 |
输出 |
ldu |
整数型 |
矩阵u的主维，ldu>=1; 当jobu='S'或'A'时， ldu>=m。 |
输入 |
vt |
|
大小为ldvt*n。 若jobvt='A': 其为n*n的酉矩阵v**H。 若jobvt='S': 其包含矩阵v**H(右奇异值向量, 按行优先存储)的前min(m,n)行。 若jobvt='N'或'O': 其不被使用。 |
输出 |
ldvt |
整数型 |
矩阵vt的主维，ldvt>=1; 若jobvt='A', 则ldvt>=n; 若jobvt='S'，则ldvt>=min(m,n)。 |
输入 |
work |
|
大小为max(1, lwork)。 退出时，若info=0, work(1)返回最优lwork大小。 |
输出 |
lwork |
整数型 |
work数组的大小，lwork>=max(1, 2*min(m,n)+max(m,n)) 若lwork=-1, 则表明函数只计算最优的work大小 |
输入 |
rwork |
|
大小为5*min(m,n) 退出时，若info>0, rwork(1:min(m,n)-1)包含了上双对角矩阵B的未收敛超对角线元素，而其对角线元素保存在s中 |
输出 |
info |
整数型 |
|
输出 |
依赖
#include "klapack.h"
示例
C Interface：
const char jobu = 'A';
const char jobvt = 'A';
const int m = 5;
const int n = 5;
const int lda = 5;
const int ldu = 5;
const int ldvt = 5;
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 *s = (double*)malloc(m * sizeof(double));
double *u = (double*)malloc(ldu*m * sizeof(double));
double *vt = (double*)malloc(ldvt*n * sizeof(double));
double qwork;
int lwork = -1;
int info = 0;
dgesvd_(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &qwork, &lwork, &info);
if (info != 0) {
printf("Error, info = %d\n", info);
return info;
}
lwork = (int)qwork;
double *work = (double*)malloc(lwork * sizeof(double));
dgesvd_(&jobu, &jobvt, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &info);
if (info != 0) {
printf("Error, info = %d\n", info);
return info;
}
/* output */
* s
* 254.071516 87.547400 67.835073 64.370735 14.897476
* u
* -0.395825 -0.377874 -0.458327 -0.488090 -0.502234 -0.155273 -0.217110 0.348128 0.628260 -0.642534 * -0.312281 0.362918 -0.694714 0.528618 0.093312 -0.399962 0.761762 0.321232 -0.289476 -0.269743
* -0.749482 -0.313183 0.287969 0.061841 0.503430
* vt
* -0.562230 -0.504441 -0.297298 0.148954 -0.564680 -0.340949 0.391046 -0.060652 0.818829 0.238069 * -0.403926 0.372759 0.728616 -0.186514 -0.363629 -0.527964 0.375181 -0.468379 -0.520847 0.299721 * -0.354610 -0.559386 0.397083 -0.035531 0.634351
Fortran Interface：
CHARACTER :: jobu = "A"
CHARACTER :: jobvt = "A"
PARAMETER (m = 5)
PARAMETER (n = 5)
PARAMETER (lda = 5)
PARAMETER (ldu = 5)
PARAMETER (ldvt = 5)
INTEGER :: info = 0
REAL(8) :: a(lad, n)
REAL(8) :: s(m)
REAL(8) :: u(ldu, m)
REAL(8) :: vt(ldvt, n)
REAL(8), ALLOCATABLE :: work(:)
INTEGER :: lwork = -1
REAL(8 :: qwork
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/
EXTERNAL DGESVD
CALL DGESVD(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, qwork, lwork, info);
lwork = (int)qwork;
work = (double *)malloc(sizeof(double) * lwork);
CALL DGESVD(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info);
*
* Output:
* s
* 254.071516 87.547400 67.835073 64.370735 14.897476
* u
* -0.395825 -0.377874 -0.458327 -0.488090 -0.502234 -0.155273 -0.217110 0.348128 0.628260 -0.642534 * -0.312281 0.362918 -0.694714 0.528618 0.093312 -0.399962 0.761762 0.321232 -0.289476 -0.269743
* -0.749482 -0.313183 0.287969 0.061841 0.503430
* vt
* -0.562230 -0.504441 -0.297298 0.148954 -0.564680 -0.340949 0.391046 -0.060652 0.818829 0.238069 * -0.403926 0.372759 0.728616 -0.186514 -0.363629 -0.527964 0.375181 -0.468379 -0.520847 0.299721 * -0.354610 -0.559386 0.397083 -0.035531 0.634351