使用对双对角矩阵的QR算法的SVD分解来求解奇异值或左右奇异向量。
C Interface:
void sbdsqr_(const char *uplo, const int *n, const int *ncvt, const int *nru, const int *ncc, float *d, float *e, float *vt, const int *ldvt, float *u, const int *ldu,
float *c, const int *ldc, float *work, int *info)
void dbdsqr_(const char *uplo, const int *n, const int *ncvt, const int *nru, const int *ncc, double *d, double *e, double *vt, const int *ldvt, double *u, const int *ldu,
double *c, const int *ldc, double *work, int *info)
void cbdsqr_(const char *uplo, const int *n, const int *ncvt, const int *nru, const int *ncc, float *d, float *e, float _Complex *vt, const int *ldvt, float _Complex *u, const int *ldu,
float _Complex *c, const int *ldc, float *work, int *info)
void zbdsqr_(const char *uplo, const int *n, const int *ncvt, const int *nru, const int *ncc, double *d, double *e, double _Complex *vt, const int *ldvt, double _Complex *u,
const int *ldu, double _Complex *c, const int *ldc, double *work, int *info)
Fortran Interface:
SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
参数 |
类型 |
说明 |
输入/输出 |
---|---|---|---|
uplo |
字符型 |
|
输入 |
n |
整数型 |
矩阵B的阶数,n>=0。 |
输入 |
nru |
整数型 |
矩阵U的行数,nru≥0。 |
输入 |
ncvt |
整数型 |
矩阵VT的列数,ncvt≥0。 |
输入 |
ncc |
整数型 |
矩阵C的列数,NCC≥0。 |
输入 |
d |
|
维度为n。
|
输入,输出 |
e |
|
维度为n-1。
|
输入 输出 |
vt |
|
维度未ldvt*ncvt,nvct=0时无意义。
|
输入,输出 |
ldvt |
整数型 |
vt的前导维度。
|
输入 |
u |
|
维度为LDU*N,NRU=0时无意义。
|
输入,输出 |
ldu |
整数型 |
矩阵U 的前导维度,ldu>=max(1,NRU)。 |
输入 |
c |
|
维度为 LDC*NCC,NCC=0时无意义。
|
输入,输出 |
ldc |
整数型 |
C的前导维度。
|
输入 |
rwork |
|
维度为4*N。 |
输出 |
info |
整数型 |
|
输出 |
#include "klapack.h"
C Interface
const char uplo = 'l'; const int n = 5; const int ncvt = 5; const int nru = 5; const int ncc = 5; const int ldu = 5; const int ldvt = 5; const int ldc = 5; double d[] = {-151.409804,-96.158378,30.604367,-49.958984,-65.000480}; double e[] = {187.838419,-61.829661,47.241856,4.133032}; double u[] = {0.562230,-0.822915,-0.081874,-0.001979,0.000007, 0.504441,0.270668,0.734742,0.363253,-0.021916, -0.297298,-0.191539,-0.135822,0.802484,-0.460975, -0.148954,-0.098365,-0.044364,0.423397,0.887076, -0.564680,-0.450761,0.658054,-0.211655,-0.010871}; double c[] = {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 vt[] = {-0.943444,-0.291675,0.133196,0.063327,0.055560, 0.331373,-0.816196,0.395311,0.189552,0.178367, -0.010230,0.452864,0.497591,0.289640,0.680678, 0.000389,-0.208325,-0.619098,-0.271648,0.706772, -0.000002,0.016272,0.441712,-0.895754,0.047431}; double *work = (double*)malloc(4 * (n-1) * sizeof(double)); int info = 0; dbdsqr_(&uplo, &n, &ncvt, &nru, &ncc, d, e, vt, &ldvt, u, &ldu, c, &ldc, work, &info); if (info != 0) { printf("Error, info = %d\n", info); return info; } * output: * d: * 254.071515 87.547400 67.835074 64.370735 14.897476 * e: * 0.000000 0.000000 0.000000 0.000000 * u: * -0.074378 -0.669531 -0.639450 -0.366581 0.054025 0.159978 -0.508435 0.027239 0.845380 -0.022155 * -0.082621 0.347676 -0.436891 0.259353 0.783682 -0.684158 -0.337057 0.510908 -0.079527 0.388547 * -0.702835 0.242353 -0.372103 0.278142 -0.481107 * vt: * 0.792096 0.561274 -0.189353 -0.143425 0.033614 -0.587068 0.711959 -0.120375 -0.230494 0.284338 * 0.154739 -0.190565 0.546605 -0.459085 0.655902 -0.063174 -0.042355 0.090008 -0.778081 -0.617013 * 0.000527 0.374135 0.801734 0.331865 -0.327278 * c: * -19.247910 110.078454 -62.583336 72.822482 -34.891161 -54.076617 90.246359 9.555710 17.064964 5.701693 * -20.142043 75.228107 44.453325 64.653600 -44.635543 26.877086 126.100144 5.668549 67.580452 -1.791622 * -25.484162 42.029730 -35.053810 77.698370 -40.669875
Fortran Interface
CHARACTER :: uplo = "L" PARAMETER (n = 5) PARAMETER (ncvt = 5) PARAMETER (nru = 5) PARAMETER (ncc = 5) PARAMETER (ldu = 5) PARAMETER (ldc = 5) PARAMETER (ldvt = 5) INTEGER :: info = 0 REAL(8) :: d(n) REAL(8) :: e(n-1) REAL(8) :: u(ldu, n) REAL(8) :: c(ldc, ncc) REAL(8) :: vt(ldvt, ncvt) REAL(8) :: work(4*(n-1)) DATA d / -151.409804,-96.158378,30.604367,-49.958984,-65.000480 / DATA e / 187.838419,-61.829661,47.241856,4.133032 / DATA u / 0.562230,-0.822915,-0.081874,-0.001979,0.000007, 0.504441,0.270668,0.734742,0.363253,-0.021916, -0.297298,-0.191539,-0.135822,0.802484,-0.460975, -0.148954,-0.098365,-0.044364,0.423397,0.887076, -0.564680,-0.450761,0.658054,-0.211655,-0.010871 / DATA c / 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 vt / -0.943444,-0.291675,0.133196,0.063327,0.055560, 0.331373,-0.816196,0.395311,0.189552,0.178367, -0.010230,0.452864,0.497591,0.289640,0.680678, 0.000389,-0.208325,-0.619098,-0.271648,0.706772, -0.000002,0.016272,0.441712,-0.895754,0.047431 / EXTERNAL DBDSQR CALL DBDSQR(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info); * output: * d: * 254.071515 87.547400 67.835074 64.370735 14.897476 * e: * 0.000000 0.000000 0.000000 0.000000 * u: * -0.074378 -0.669531 -0.639450 -0.366581 0.054025 0.159978 -0.508435 0.027239 0.845380 -0.022155 * -0.082621 0.347676 -0.436891 0.259353 0.783682 -0.684158 -0.337057 0.510908 -0.079527 0.388547 * -0.702835 0.242353 -0.372103 0.278142 -0.481107 * vt: * 0.792096 0.561274 -0.189353 -0.143425 0.033614 -0.587068 0.711959 -0.120375 -0.230494 0.284338 * 0.154739 -0.190565 0.546605 -0.459085 0.655902 -0.063174 -0.042355 0.090008 -0.778081 -0.617013 * 0.000527 0.374135 0.801734 0.331865 -0.327278 * c: * -19.247910 110.078454 -62.583336 72.822482 -34.891161 -54.076617 90.246359 9.555710 17.064964 5.701693 * -20.142043 75.228107 44.453325 64.653600 -44.635543 26.877086 126.100144 5.668549 67.580452 -1.791622 * -25.484162 42.029730 -35.053810 77.698370 -40.669875