?bdsqr

使用对双对角矩阵的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

字符型

  • 'U':B是上双矩阵。
  • 'L':B是下双矩阵。

输入

n

整数型

矩阵B的阶数,n>=0。

输入

nru

整数型

矩阵U的行数,nru≥0。

输入

ncvt

整数型

矩阵VT的列数,ncvt≥0。

输入

ncc

整数型

矩阵C的列数,NCC≥0。

输入

d

  • 在sbdsqr中为单精度浮点型数组
  • 在dbdsqr中为双精度浮点型数组
  • 在cbdsqr中为单精度浮点型数组
  • 在zbdsqr中为双精度浮点型数组

维度为n。

  • 进入时:为双对角矩阵B的n个对角元素
  • 退出时:
    • 若info =0,则B的奇异值为下降趋势

输入,输出

e

  • 在sbdsqr中为单精度浮点型数组
  • 在dbdsqr中为双精度浮点型数组
  • 在cbdsqr中为单精度浮点型数组
  • 在zbdsqr中为双精度浮点型数组

维度为n-1。

  • 进入时:为双对角矩阵B的n-1个非对角元素
  • 退出时:
    • 若info =0,则e被摧毁。
    • 若info >0, 则D和E会容纳与输入正交化相等的双对角矩阵的对角和超对角元素。

输入

输出

vt

  • 在sbdsqr中为单精度浮点型数组
  • 在dbdsqr中为双精度浮点型数组
  • 在cbdsqr中为单精度复数型数组
  • 在zbdsqr中为双精度复数型数组

维度未ldvt*ncvt,nvct=0时无意义。

  • 进入时:为N*NCVT大小的矩阵VT。
  • 退出时:被P^H*VT覆盖。

输入,输出

ldvt

整数型

vt的前导维度。

  • ncvt>0时,ldvt>=max(1,N)。
  • ncvt=0时,ldvt>-1。

输入

u

  • 在sbdsqr中为单精度浮点型数组
  • 在dbdsqr中为双精度浮点型数组
  • 在cbdsqr中为单精度复数型数组
  • 在zbdsqr中为双精度复数型数组

维度为LDU*N,NRU=0时无意义。

  • 进入时,NRU*N的矩阵。
  • 退出时,U被U*Q覆盖。

输入,输出

ldu

整数型

矩阵U 的前导维度,ldu>=max(1,NRU)。

输入

c

  • 在sbdsqr中为单精度浮点型数组
  • 在dbdsqr中为双精度浮点型数组
  • 在cbdsqr中为单精度复数型数组
  • 在zbdsqr中为双精度复数型数组

维度为 LDC*NCC,NCC=0时无意义。

  • 进入时,N*NCC的矩阵C。
  • 退出时,被Q^H*C覆盖。

输入,输出

ldc

整数型

C的前导维度。

  • 如果NCC>0,LDC>=max(1,N)。
  • 如果NCC=0,LDC>=1。

输入

rwork

  • 在sbdsqr中为单精度浮点型数组
  • 在dbdsqr中为双精度浮点型数组
  • 在cbdsqr中为单精度浮点型数组
  • 在zbdsqr中为双精度浮点型数组

维度为4*N。

输出

info

整数型

  • 等于0:表示成功。
  • 小于0:info=-i,表示第i个参数非法。
  • 大于0:计算出错。

输出

依赖

#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