Man page - gesvdq(3)
Packages contains this manual
- hptrd(3)
- potri(3)
- xerbla_array(3)
- ggsvd_driver_grp(3)
- hfrk(3)
- getsqr_comp_grp(3)
- laed6(3)
- gtrfs(3)
- lasdq(3)
- gglse(3)
- la_xisnan_la_isnan(3)
- unmr2(3)
- hetrs_aa(3)
- tpttr(3)
- gerz_comp_grp(3)
- potrf(3)
- hegv_driver(3)
- laqps(3)
- ggqr_comp_grp(3)
- ilalc(3)
- ung2r(3)
- heevd(3)
- pstf2(3)
- lacn2(3)
- ptrfs(3)
- ungrq(3)
- gelqf(3)
- ppsv_comp(3)
- blas2_full(3)
- gemlqt(3)
- unml2(3)
- tplqt(3)
- tpcon(3)
- getf2(3)
- ggbak(3)
- bdsvd_driver(3)
- lamch(3)
- gelq(3)
- gebal(3)
- laqr1(3)
- ptsvx(3)
- lahr2(3)
- larscl2(3)
- geqrt(3)
- larfb(3)
- gtsv_comp(3)
- gesvd_aux(3)
- hbevx_2stage(3)
- hbgvx(3)
- tprfs(3)
- params_grp(3)
- lahef(3)
- laqr_group(3)
- unmqr(3)
- tgsy2(3)
- tfsv_comp(3)
- ggls_driver_grp(3)
- geev(3)
- latrd(3)
- unbdb4(3)
- bbcsd(3)
- lange(3)
- gelq_comp3(3)
- gttrs(3)
- lasy2(3)
- hetf2_rook(3)
- gtsv(3)
- lalsd(3)
- lanhb(3)
- laqhb(3)
- hgeqz(3)
- gesvj(3)
- gsvj0(3)
- ungtsqr_row(3)
- gelq_comp1(3)
- gemmtr(3)
- pbequ(3)
- heev_driver(3)
- unhr_col(3)
- syconvf_rook(3)
- getc2(3)
- syconv(3)
- norm_grp(3)
- larrc(3)
- laqr4(3)
- posv_comp(3)
- geev_driver_grp(3)
- heev_comp(3)
- pfsv(3)
- trevc3(3)
- gesv_driver_grp(3)
- reflector_aux_grp(3)
- langt(3)
- lacrt(3)
- latdf(3)
- hetrs_aa_2stage(3)
- lamc1(3)
- hpev_driver(3)
- hegvd(3)
- pptri(3)
- geqrt3(3)
- gelqt3(3)
- lasd5(3)
- laeda(3)
- geqr(3)
- lamtsqr(3)
- heev(3)
- hpev_comp(3)
- larfg(3)
- blas2_grp(3)
- hesv_rook(3)
- laexc(3)
- hetrd(3)
- geesx(3)
- ppsvx(3)
- blas_top(3)
- gtts2(3)
- la_herpvgrw(3)
- hpevx(3)
- ggevx(3)
- lahqr(3)
- gelq_comp_grp(3)
- hesv_comp_v3(3)
- tplqt2(3)
- hpev(3)
- hbtrd(3)
- getrs(3)
- hecon_3(3)
- lasrt(3)
- lanhe(3)
- gesv_comp(3)
- gbequ(3)
- hetrf_rk(3)
- laqr3(3)
- heev_comp_grp(3)
- ungtsqr(3)
- ppcon(3)
- ggrq_comp_grp(3)
- larmm(3)
- ieeeck(3)
- geqrf(3)
- solve_aux_grp(3)
- herfs(3)
- posvx(3)
- posvxx(3)
- gges3(3)
- hbgvd(3)
- lantb(3)
- lasd_comp_grp(3)
- hpgvx(3)
- lapy2(3)
- lauu2(3)
- copy(3)
- getsqrhrt(3)
- stev_comp_grp(3)
- laev2(3)
- larfb_gett(3)
- trti2(3)
- laqz4(3)
- hegv_driver_grp(3)
- la_porfsx_extended(3)
- laruv(3)
- ggsvd_comp_grp(3)
- dot(3)
- gehd2(3)
- lanhf(3)
- hetri_rook(3)
- pfsv_comp(3)
- gbtrf(3)
- hpgst(3)
- getri(3)
- trevc(3)
- unmrz(3)
- hsein(3)
- lsamen(3)
- lasd6(3)
- trtri(3)
- ggglm(3)
- las2(3)
- latrs(3)
- lapll(3)
- gemlq(3)
- geqpf_comp_grp(3)
- stemr(3)
- rotm(3)
- disna(3)
- ggrqf(3)
- pptrf(3)
- lasd0(3)
- lals0(3)
- laqz2(3)
- hbev_driver2(3)
- geswlq_comp_grp(3)
- laqr0(3)
- trttp(3)
- stedc(3)
- lasq4(3)
- geev_comp_grp(3)
- ungbr(3)
- lanv2(3)
- hpsv(3)
- pprfs(3)
- gehrd(3)
- ppsv(3)
- lagtm(3)
- hpgv(3)
- trsv_comp(3)
- larfx(3)
- gesv_driver(3)
- gerfsx(3)
- la_geamv(3)
- laed9(3)
- tpqrt2(3)
- uncsd(3)
- gecs_comp_grp(3)
- bdsqr(3)
- hegv_comp_grp(3)
- labad(3)
- geqp3(3)
- gesvdq(3)
- tfttp(3)
- laln2(3)
- uncsd2by1(3)
- blas2_like_grp(3)
- latbs(3)
- hbgst(3)
- larrv(3)
- ilaenv2stage(3)
- bdsvdx(3)
- hegs2(3)
- lasq_comp_grp(3)
- hpr2(3)
- laqhe(3)
- larra(3)
- gemqrt(3)
- hbmv(3)
- hpsv_driver(3)
- lacp2(3)
- lapmt(3)
- gecon(3)
- unbdb5(3)
- la_gerpvgrw(3)
- tgex2(3)
- laqhp(3)
- tftri(3)
- getrf2(3)
- porfs(3)
- lartg(3)
- lagts(3)
- ggev_comp_grp(3)
- lasd3(3)
- geqr_comp2(3)
- laqz_group(3)
- pftri(3)
- hetri2x(3)
- lahef_aa(3)
- svd_driver_grp(3)
- gbsv_driver(3)
- hesv_comp_aasen2(3)
- laqtr(3)
- lag2(3)
- la_porcond(3)
- hbev(3)
- pbtrf(3)
- lascl(3)
- larr_comp_grp(3)
- hecon(3)
- pttrs(3)
- lasd8(3)
- lsame(3)
- unm2l(3)
- potrs(3)
- tptrs(3)
- lartv(3)
- trtrs(3)
- gsvj1(3)
- sum1(3)
- larrj(3)
- gbmv(3)
- posv(3)
- gghd3(3)
- geev_top(3)
- geqr_comp_grp(3)
- laset(3)
- hesvxx(3)
- posv_comp_grp(3)
- lahef_rk(3)
- lasd1(3)
- tprfb(3)
- potf2(3)
- laein(3)
- lamc4(3)
- stevd(3)
- gtsv_driver(3)
- gesvd_comp_grp(3)
- la_constants(3)
- gesvx(3)
- hseqr(3)
- launhr_col_getrfnp2(3)
- trcon(3)
- larre(3)
- gelsy(3)
- ptsv(3)
- lacon(3)
- laed_comp_grp(3)
- hpsvx(3)
- gemm(3)
- poequ(3)
- laesy(3)
- lagtf(3)
- trrfs(3)
- ggev3(3)
- pbstf(3)
- poequb(3)
- heevr(3)
- lanhp(3)
- unbdb3(3)
- tgsyl(3)
- lamc5(3)
- geqr2p(3)
- ungqr(3)
- laqz3(3)
- imax1(3)
- gels_top(3)
- hesv(3)
- gelqt(3)
- pfsv_driver(3)
- stegr(3)
- gerqf(3)
- laisnan(3)
- ilatrans(3)
- gbsv_comp(3)
- pbrfs(3)
- lascl2(3)
- larz(3)
- la_hercond(3)
- tgexc(3)
- ggesx(3)
- unbdb6(3)
- ungl2(3)
- laed_comp2(3)
- rscl(3)
- hegv(3)
- gelst(3)
- gbtrs(3)
- pftrf(3)
- langb(3)
- lantr(3)
- laqgb(3)
- ggsvp3(3)
- bdsdc(3)
- ladiv(3)
- laqge(3)
- iparmq(3)
- ggbal(3)
- hb2st_kernels(3)
- lartgs(3)
- lartgp(3)
- rot(3)
- ppequ(3)
- laed3(3)
- her(3)
- hptri(3)
- stevx(3)
- upgtr(3)
- lar2v(3)
- hbev_2stage(3)
- gejsv(3)
- ppsv_driver(3)
- unm22(3)
- gesvxx(3)
- laqz0(3)
- unmtr(3)
- laed5(3)
- tptri(3)
- laed0(3)
- heev_driver2(3)
- hpcon(3)
- lasd4(3)
- hetrf_aa(3)
- geqr_comp3(3)
- rot_aux_grp(3)
- aux_grp(3)
- laebz(3)
- trsyl3(3)
- gges(3)
- gesdd(3)
- trexc(3)
- ung2l(3)
- gesv(3)
- laed4(3)
- md__r_e_a_d_m_e(3)
- blas3_like_grp(3)
- laed1(3)
- larcm(3)
- hbevx(3)
- hesv_driver_grp(3)
- hetrs(3)
- hbevd_2stage(3)
- blas1_grp(3)
- laic1(3)
- geql_comp_grp(3)
- heev_2stage(3)
- hpmv(3)
- pbtf2(3)
- hetrf_aa_2stage(3)
- hbgv(3)
- pptrs(3)
- lapmr(3)
- tpqr_comp_grp(3)
- larfy(3)
- gedmd(3)
- lasr(3)
- hetrd_2stage(3)
- gerfs(3)
- ungtr(3)
- porfsx(3)
- tpmv(3)
- lasd_comp2(3)
- unmbr(3)
- tbtrs(3)
- hetd2(3)
- trsv_comp_grp(3)
- lapy3(3)
- ptts2(3)
- unmhr(3)
- hbev_driver(3)
- lalsa(3)
- tbsv_comp(3)
- hesv_comp_v1(3)
- geql2(3)
- sterf(3)
- larrd(3)
- larft(3)
- lagv2(3)
- gttrf(3)
- tpqrt(3)
- la_lin_berr(3)
- rotg(3)
- solve_top(3)
- lacgv(3)
- larrf(3)
- tbmv(3)
- trsyl(3)
- geequ(3)
- upmtr(3)
- hpgv_driver(3)
- tbsv(3)
- hesvx(3)
- latrz(3)
- tfttr(3)
- gesv_comp_grp(3)
- xerbla_grp(3)
- tpsv(3)
- blas3_grp(3)
- gesvd_driver(3)
- geqr_comp1(3)
- ggev_driver_grp(3)
- la_gbamv(3)
- tpmlqt(3)
- trttf(3)
- larzb(3)
- unmr3(3)
- hecon_rook(3)
- stebz(3)
- lantp(3)
- laqz1(3)
- hesv_rk(3)
- tbcon(3)
- xerbla(3)
- posv_mixed(3)
- latps(3)
- hesv_aa_driver(3)
- gemqr(3)
- larrr(3)
- gebrd(3)
- tgsna(3)
- la_gercond(3)
- gbsv(3)
- hesv_comp_grp(3)
- gesv_mixed(3)
- gghrd(3)
- gbrfs(3)
- tpmqrt(3)
- lasq3(3)
- tpsv_comp(3)
- largv(3)
- gelsd(3)
- pftrs(3)
- asum(3)
- launhr_col_getrfnp(3)
- hptrf(3)
- lacpy(3)
- gesc2(3)
- lasda(3)
- second(3)
- hprfs(3)
- hpsv_comp(3)
- lamrg(3)
- pbsv_comp(3)
- hegv_2stage(3)
- gerq2(3)
- lasdt(3)
- abs1(3)
- hbevd(3)
- hbev_comp(3)
- trsv(3)
- la_porpvgrw(3)
- la_gbrpvgrw(3)
- hbgv_driver(3)
- tgsja(3)
- gebd2(3)
- geqr2(3)
- unm2r(3)
- unmql(3)
- la_gbrfsx_extended(3)
- gelq_comp2(3)
- iparam2stage(3)
- ger(3)
- larf(3)
- ilaprec(3)
- labrd(3)
- unbdb1(3)
- unmlq(3)
- geequb(3)
- la_herfsx_extended(3)
- unbdb2(3)
- lapack_top(3)
- ptsv_driver(3)
- hetrs2(3)
- geqr_comp4(3)
- pbsv(3)
- posv_driver(3)
- steqr(3)
- gels(3)
- lar1v(3)
- hemv(3)
- la_transtype(3)
- hesv_aa(3)
- lacrm(3)
- stevr(3)
- hetf2_rk(3)
- blas2_banded(3)
- stein(3)
- unmrq(3)
- larrk(3)
- hetri2(3)
- hesv_aa_2stage(3)
- pttrf(3)
- gelss(3)
- pbsv_driver(3)
- lasq5(3)
- heevx_2stage(3)
- hetri(3)
- lasd2(3)
- laed2(3)
- pbcon(3)
- ptcon(3)
- laed7(3)
- gels_aux_grp(3)
- hpgvd(3)
- hetf2(3)
- tzrzf(3)
- hpr(3)
- unitary_top(3)
- latsqr(3)
- ungql(3)
- her2(3)
- hetri_3x(3)
- hetrd_hb2st(3)
- tgsen(3)
- ggsvd3(3)
- lasq6(3)
- set_grp(3)
- larfgp(3)
- gels_driver_grp(3)
- pbtrs(3)
- lamswlq(3)
- lanht(3)
- gbsvxx(3)
- tgevc(3)
- ilaenv(3)
- swap(3)
- lae2(3)
- iladiag(3)
- lasq2(3)
- la_heamv(3)
- blas_like_top(3)
- la_gerfsx_extended(3)
- hegst(3)
- tfsm(3)
- gesvd(3)
- ungr2(3)
- ggev(3)
- aux_top(3)
- blas2_packed(3)
- geqlf(3)
- hetrs_rook(3)
- gelq2(3)
- geqrfp(3)
- gbequb(3)
- stev(3)
- lauum(3)
- potrf2(3)
- lamc3(3)
- gbrfsx(3)
- gerq_comp_grp(3)
- pocon(3)
- tbrfs(3)
- heswapr(3)
- lamc2(3)
- hpevd(3)
- hesv_comp_aasen(3)
- scalar_grp(3)
- gemv(3)
- lasv2(3)
- lanhs(3)
- svd_top(3)
- gbsvx(3)
- gesvdx(3)
- tplq_comp_grp(3)
- hesv_driver(3)
- hesv_comp_v2(3)
- trsen(3)
- syconvf(3)
- lasd7(3)
- gbcon(3)
- unbdb(3)
- heev_driver_grp(3)
- ggqrf(3)
- heevx(3)
- gtsvx(3)
- lahef_rook(3)
- hetrf_rook(3)
- hetrf(3)
- trsna(3)
- gebak(3)
- larnv(3)
- ptsv_comp(3)
- laswlq(3)
- lags2(3)
- laed8(3)
- laswp(3)
- hptrs(3)
- unglq(3)
- la_wwaddw(3)
- getrf(3)
- gees(3)
- gbtf2(3)
- hegvx(3)
- latrs3(3)
- roundup_lwork(3)
- unghr(3)
- iamax(3)
- larzt(3)
- pteqr(3)
- ilaver(3)
- trmv(3)
- la_gbrcond(3)
- blas0_like_grp(3)
- nrm2(3)
- heev_top(3)
- gtcon(3)
- heevr_2stage(3)
- pstrf(3)
- rot_comp(3)
- laqr5(3)
- heevd_2stage(3)
- getsls(3)
- hetrd_he2hb(3)
- heequb(3)
- laqp2(3)
- axpy(3)
- blast_aux(3)
- rotmg(3)
- pbsvx(3)
- ilauplo(3)
- herfsx(3)
- laqr2(3)
- blas1_like_grp(3)
- lassq(3)
- larrb(3)
- stev_driver(3)
- geevx(3)
- tpttf(3)
- scal(3)
- laneg(3)
- posv_driver_grp(3)
- lasq1(3)
- hetrs_3(3)
- geqrt2(3)
- gbbrd(3)
- ilalr(3)
- hetri_3(3)
apt-get install liblapack-doc
Manual
gesvdq
NAMESYNOPSIS
Functions
Detailed Description
Function Documentation
subroutine cgesvdq (character joba, character jobp, character jobr,character jobu, character jobv, integer m, integer n, complex,dimension( lda, * ) a, integer lda, real, dimension( * ) s, complex,dimension( ldu, * ) u, integer ldu, complex, dimension( ldv, * ) v,integer ldv, integer numrank, integer, dimension( * ) iwork, integerliwork, complex, dimension( * ) cwork, integer lcwork, real, dimension(* ) rwork, integer lrwork, integer info)
subroutine dgesvdq (character joba, character jobp, character jobr,character jobu, character jobv, integer m, integer n, double precision,dimension( lda, * ) a, integer lda, double precision, dimension( * ) s,double precision, dimension( ldu, * ) u, integer ldu, double precision,dimension( ldv, * ) v, integer ldv, integer numrank, integer,dimension( * ) iwork, integer liwork, double precision, dimension( * )work, integer lwork, double precision, dimension( * ) rwork, integerlrwork, integer info)
subroutine sgesvdq (character joba, character jobp, character jobr,character jobu, character jobv, integer m, integer n, real, dimension(lda, * ) a, integer lda, real, dimension( * ) s, real, dimension( ldu,* ) u, integer ldu, real, dimension( ldv, * ) v, integer ldv, integernumrank, integer, dimension( * ) iwork, integer liwork, real,dimension( * ) work, integer lwork, real, dimension( * ) rwork, integerlrwork, integer info)
subroutine zgesvdq (character joba, character jobp, character jobr,character jobu, character jobv, integer m, integer n, complex*16,dimension( lda, * ) a, integer lda, double precision, dimension( * ) s,complex*16, dimension( ldu, * ) u, integer ldu, complex*16, dimension(ldv, * ) v, integer ldv, integer numrank, integer, dimension( * )iwork, integer liwork, complex*16, dimension( * ) cwork, integerlcwork, double precision, dimension( * ) rwork, integer lrwork, integerinfo)
Author
NAME
gesvdq - gesvdq: SVD, QR with pivoting
SYNOPSIS
Functions
subroutine
cgesvdq
(joba, jobp, jobr, jobu, jobv, m, n, a, lda,
s, u, ldu, v, ldv, numrank, iwork, liwork, cwork, lcwork,
rwork, lrwork, info)
CGESVDQ computes the singular value decomposition (SVD) with
a QR-Preconditioned QR SVD Method for GE matrices
subroutine
dgesvdq
(joba, jobp, jobr, jobu, jobv, m,
n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, work,
lwork, rwork, lrwork, info)
DGESVDQ computes the singular value decomposition (SVD) with
a QR-Preconditioned QR SVD Method for GE matrices
subroutine
sgesvdq
(joba, jobp, jobr, jobu, jobv, m,
n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, work,
lwork, rwork, lrwork, info)
SGESVDQ computes the singular value decomposition (SVD) with
a QR-Preconditioned QR SVD Method for GE matrices
subroutine
zgesvdq
(joba, jobp, jobr, jobu, jobv, m,
n, a, lda, s, u, ldu, v, ldv, numrank, iwork, liwork, cwork,
lcwork, rwork, lrwork, info)
ZGESVDQ computes the singular value decomposition (SVD) with
a QR-Preconditioned QR SVD Method for GE matrices
Detailed Description
Function Documentation
subroutine cgesvdq (character joba, character jobp, character jobr,character jobu, character jobv, integer m, integer n, complex,dimension( lda, * ) a, integer lda, real, dimension( * ) s, complex,dimension( ldu, * ) u, integer ldu, complex, dimension( ldv, * ) v,integer ldv, integer numrank, integer, dimension( * ) iwork, integerliwork, complex, dimension( * ) cwork, integer lcwork, real, dimension(* ) rwork, integer lrwork, integer info)
CGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices
Purpose:
CGESVDQ
computes the singular value decomposition (SVD) of a complex
M-by-N matrix A, where M >= N. The SVD of A is written as
[++] [xx] [x0] [xx]
A = U * SIGMA * VΛ*, [++] = [xx] * [ox] * [xx]
[++] [xx]
where SIGMA is an N-by-N diagonal matrix, U is an M-by-N
orthonormal
matrix, and V is an N-by-N unitary matrix. The diagonal
elements
of SIGMA are the singular values of A. The columns of U and
V are the
left and the right singular vectors of A, respectively.
Parameters
JOBA
JOBA is
CHARACTER*1
Specifies the level of accuracy in the computed SVD
= βAβ The requested accuracy corresponds to
having the backward
error bounded by || delta A ||_F <= f(m,n) * EPS * || A
||_F,
where EPS = SLAMCH(βEpsilonβ). This authorises
CGESVDQ to
truncate the computed triangular factor in a rank revealing
QR factorization whenever the truncated part is below the
threshold of the order of EPS * ||A||_F. This is aggressive
truncation level.
= βMβ Similarly as with βAβ, but the
truncation is more gentle: it
is allowed only when there is a drop on the diagonal of the
triangular factor in the QR factorization. This is medium
truncation level.
= βHβ High accuracy requested. No numerical rank
determination based
on the rank revealing QR factorization is attempted.
= βEβ Same as βHβ, and in addition
the condition number of column
scaled A is estimated and returned in RWORK(1).
NΛ(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <=
NΛ(1/4)*RWORK(1)
JOBP
JOBP is
CHARACTER*1
= βPβ The rows of A are ordered in decreasing
order with respect to
||A(i,:)||_\infty. This enhances numerical accuracy at the
cost
of extra data movement. Recommended for numerical
robustness.
= βNβ No row pivoting.
JOBR
JOBR is
CHARACTER*1
= βTβ After the initial pivoted QR
factorization, CGESVD is applied to
the adjoint R**H of the computed triangular factor R. This
involves
some extra data movement (matrix transpositions). Useful for
experiments, research and development.
= βNβ The triangular factor R is given as input
to CGESVD. This may be
preferred as it involves less data movement.
JOBU
JOBU is
CHARACTER*1
= βAβ All M left singular vectors are computed
and returned in the
matrix U. See the description of U.
= βSβ or βUβ N = min(M,N) left
singular vectors are computed and returned
in the matrix U. See the description of U.
= βRβ Numerical rank NUMRANK is determined and
only NUMRANK left singular
vectors are computed and returned in the matrix U.
= βFβ The N left singular vectors are returned
in factored form as the
product of the Q factor from the initial QR factorization
and the
N left singular vectors of (R**H , 0)**H. If row pivoting is
used,
then the necessary information on the row pivoting is stored
in
IWORK(N+1:N+M-1).
= βNβ The left singular vectors are not
computed.
JOBV
JOBV is
CHARACTER*1
= βAβ, βVβ All N right singular
vectors are computed and returned in
the matrix V.
= βRβ Numerical rank NUMRANK is determined and
only NUMRANK right singular
vectors are computed and returned in the matrix V. This
option is
allowed only if JOBU = βRβ or JOBU =
βNβ; otherwise it is illegal.
= βNβ The right singular vectors are not
computed.
M
M is INTEGER
The number of rows of the input matrix A. M >= 0.
N
N is INTEGER
The number of columns of the input matrix A. M >= N >=
0.
A
A is COMPLEX
array of dimensions LDA x N
On entry, the input matrix A.
On exit, if JOBU .NE. βNβ or JOBV .NE.
βNβ, the lower triangle of A contains
the Householder vectors as stored by CGEQP3. If JOBU =
βFβ, these Householder
vectors together with CWORK(1:N) can be used to restore the
Q factors from
the initial pivoted QR factorization of A. See the
description of U.
LDA
LDA is INTEGER.
The leading dimension of the array A. LDA >=
max(1,M).
S
S is REAL array
of dimension N.
The singular values of A, ordered so that S(i) >=
S(i+1).
U
U is COMPLEX
array, dimension
LDU x M if JOBU = βAβ; see the description of
LDU. In this case,
on exit, U contains the M left singular vectors.
LDU x N if JOBU = βSβ, βUβ,
βRβ ; see the description of LDU. In this
case, U contains the leading N or the leading NUMRANK left
singular vectors.
LDU x N if JOBU = βFβ ; see the description of
LDU. In this case U
contains N x N unitary matrix that can be used to form the
left
singular vectors.
If JOBU = βNβ, U is not referenced.
LDU
LDU is INTEGER.
The leading dimension of the array U.
If JOBU = βAβ, βSβ, βUβ,
βRβ, LDU >= max(1,M).
If JOBU = βFβ, LDU >= max(1,N).
Otherwise, LDU >= 1.
V
V is COMPLEX
array, dimension
LDV x N if JOBV = βAβ, βVβ,
βRβ or if JOBA = βEβ .
If JOBV = βAβ, or βVβ, V contains
the N-by-N unitary matrix V**H;
If JOBV = βRβ, V contains the first NUMRANK rows
of V**H (the right
singular vectors, stored rowwise, of the NUMRANK largest
singular values).
If JOBV = βNβ and JOBA = βEβ, V is
used as a workspace.
If JOBV = βNβ, and JOBA.NE.βEβ, V is
not referenced.
LDV
LDV is INTEGER
The leading dimension of the array V.
If JOBV = βAβ, βVβ, βRβ,
or JOBA = βEβ, LDV >= max(1,N).
Otherwise, LDV >= 1.
NUMRANK
NUMRANK is
INTEGER
NUMRANK is the numerical rank first determined after the
rank
revealing QR factorization, following the strategy specified
by the
value of JOBA. If JOBV = βRβ and JOBU =
βRβ, only NUMRANK
leading singular values and vectors are then requested in
the call
of CGESVD. The final value of NUMRANK might be further
reduced if
some singular values are computed as zeros.
IWORK
IWORK is
INTEGER array, dimension (max(1, LIWORK)).
On exit, IWORK(1:N) contains column pivoting permutation of
the
rank revealing QR factorization.
If JOBP = βPβ, IWORK(N+1:N+M-1) contains the
indices of the sequence
of row swaps used in row pivoting. These can be used to
restore the
left singular vectors in the case JOBU =
βFβ.
If LIWORK,
LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
IWORK(1) returns the minimal LIWORK.
LIWORK
LIWORK is
INTEGER
The dimension of the array IWORK.
LIWORK >= N + M - 1, if JOBP = βPβ;
LIWORK >= N if JOBP = βNβ.
If LIWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the CWORK, IWORK, and RWORK arrays, and no error
message related to LCWORK is issued by XERBLA.
CWORK
CWORK is
COMPLEX array, dimension (max(2, LCWORK)), used as a
workspace.
On exit, if, on entry, LCWORK.NE.-1, CWORK(1:N) contains
parameters
needed to recover the Q factor from the QR factorization
computed by
CGEQP3.
If LIWORK,
LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
CWORK(1) returns the optimal LCWORK, and
CWORK(2) returns the minimal LCWORK.
LCWORK
LCWORK is
INTEGER
The dimension of the array CWORK. It is determined as
follows:
Let LWQP3 = N+1, LWCON = 2*N, and let
LWUNQ = { MAX( N, 1 ), if JOBU = βRβ,
βSβ, or βUβ
{ MAX( M, 1 ), if JOBU = βAβ
LWSVD = MAX( 3*N, 1 )
LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 3*(N/2), 1 ), LWUNLQ =
MAX( N, 1 ),
LWQRF = MAX( N/2, 1 ), LWUNQ2 = MAX( N, 1 )
Then the minimal value of LCWORK is:
= MAX( N + LWQP3, LWSVD ) if only the singular values are
needed;
= MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values
are needed,
and a scaled condition estimate requested;
= N + MAX(
LWQP3, LWSVD, LWUNQ ) if the singular values and the left
singular vectors are requested;
= N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the singular
values and the left
singular vectors are requested, and also
a scaled condition estimate requested;
= N + MAX(
LWQP3, LWSVD ) if the singular values and the right
singular vectors are requested;
= N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and
the right
singular vectors are requested, and also
a scaled condition etimate requested;
= N + MAX(
LWQP3, LWSVD, LWUNQ ) if the full SVD is requested with JOBV
= βRβ;
independent of JOBR;
= N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the full SVD is
requested,
JOBV = βRβ and, also a scaled condition
estimate requested; independent of JOBR;
= MAX( N + MAX( LWQP3, LWSVD, LWUNQ ),
N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ) )
if the
full SVD is requested with JOBV = βAβ or
βVβ, and
JOBR =βNβ
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ,
LWUNQ ) )
if the full SVD is requested with JOBV = βAβ or
βVβ, and
JOBR =βNβ, and also a scaled condition number
estimate
requested.
= MAX( N + MAX( LWQP3, LWSVD, LWUNQ ),
N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) )
if the
full SVD is requested with JOBV = βAβ,
βVβ, and JOBR =βTβ
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2,
LWUNQ ) )
if the full SVD is requested with JOBV = βAβ,
βVβ and
JOBR =βTβ, and also a scaled condition number
estimate
requested.
Finally, LCWORK must be at least two: LCWORK = MAX( 2,
LCWORK ).
If LCWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the CWORK, IWORK, and RWORK arrays, and no error
message related to LCWORK is issued by XERBLA.
RWORK
RWORK is REAL
array, dimension (max(1, LRWORK)).
On exit,
1. If JOBA = βEβ, RWORK(1) contains an estimate
of the condition
number of column scaled A. If A = C * D where D is diagonal
and C
has unit columns in the Euclidean norm, then, assuming full
column rank,
NΛ(-1/4) * RWORK(1) <= ||pinv(C)||_2 <=
NΛ(1/4) * RWORK(1).
Otherwise, RWORK(1) = -1.
2. RWORK(2) contains the number of singular values computed
as
exact zeros in CGESVD applied to the upper triangular or
trapezoidal
R (from the initial QR factorization). In case of early exit
(no call to
CGESVD, such as in the case of zero matrix) RWORK(2) =
-1.
If LIWORK,
LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
RWORK(1) returns the minimal LRWORK.
LRWORK
LRWORK is
INTEGER.
The dimension of the array RWORK.
If JOBP =βPβ, then LRWORK >= MAX(2, M, 5*N);
Otherwise, LRWORK >= MAX(2, 5*N).
If LRWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the CWORK, IWORK, and RWORK arrays, and no error
message related to LCWORK is issued by XERBLA.
INFO
INFO is INTEGER
= 0: successful exit.
< 0: if INFO = -i, the i-th argument had an illegal
value.
> 0: if CBDSQR did not converge, INFO specifies how many
superdiagonals
of an intermediate bidiagonal form B (computed in CGESVD)
did not
converge to zero.
Further Details:
1. The data
movement (matrix transpose) is coded using simple nested
DO-loops because BLAS and LAPACK do not provide
corresponding subroutines.
Those DO-loops are easily identified in this source code -
by the CONTINUE
statements labeled with 11**. In an optimized version of
this code, the
nested DO loops should be replaced with calls to an
optimized subroutine.
2. This code scales A by 1/SQRT(M) if the largest
ABS(A(i,j)) could cause
column norm overflow. This is the minial precaution and it
is left to the
SVD routine (CGESVD) to do its own preemptive scaling if
potential over-
or underflows are detected. To avoid repeated scanning of
the array A,
an optimal implementation would do all necessary scaling
before calling
CGESVD and the scaling in CGESVD can be switched off.
3. Other comments related to code optimization are given in
comments in the
code, enclosed in [[double brackets]].
Bugs, examples and comments
Please report
all bugs and send interesting examples and/or comments to
drmac@math.hr. Thank you.
References
[1] Zlatko
Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
Computing the SVD with High Accuracy. ACM Trans. Math.
Softw.
44(1): 11:1-11:30 (2017)
SIGMA library,
xGESVDQ section updated February 2016.
Developed and coded by Zlatko Drmac, Department of
Mathematics
University of Zagreb, Croatia, drmac@math.hr
Contributors:
Developed and
coded by Zlatko Drmac, Department of Mathematics
University of Zagreb, Croatia, drmac@math.hr
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
subroutine dgesvdq (character joba, character jobp, character jobr,character jobu, character jobv, integer m, integer n, double precision,dimension( lda, * ) a, integer lda, double precision, dimension( * ) s,double precision, dimension( ldu, * ) u, integer ldu, double precision,dimension( ldv, * ) v, integer ldv, integer numrank, integer,dimension( * ) iwork, integer liwork, double precision, dimension( * )work, integer lwork, double precision, dimension( * ) rwork, integerlrwork, integer info)
DGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices
Purpose:
DGESVDQ
computes the singular value decomposition (SVD) of a real
M-by-N matrix A, where M >= N. The SVD of A is written as
[++] [xx] [x0] [xx]
A = U * SIGMA * VΛ*, [++] = [xx] * [ox] * [xx]
[++] [xx]
where SIGMA is an N-by-N diagonal matrix, U is an M-by-N
orthonormal
matrix, and V is an N-by-N orthogonal matrix. The diagonal
elements
of SIGMA are the singular values of A. The columns of U and
V are the
left and the right singular vectors of A, respectively.
Parameters
JOBA
JOBA is
CHARACTER*1
Specifies the level of accuracy in the computed SVD
= βAβ The requested accuracy corresponds to
having the backward
error bounded by || delta A ||_F <= f(m,n) * EPS * || A
||_F,
where EPS = DLAMCH(βEpsilonβ). This authorises
DGESVDQ to
truncate the computed triangular factor in a rank revealing
QR factorization whenever the truncated part is below the
threshold of the order of EPS * ||A||_F. This is aggressive
truncation level.
= βMβ Similarly as with βAβ, but the
truncation is more gentle: it
is allowed only when there is a drop on the diagonal of the
triangular factor in the QR factorization. This is medium
truncation level.
= βHβ High accuracy requested. No numerical rank
determination based
on the rank revealing QR factorization is attempted.
= βEβ Same as βHβ, and in addition
the condition number of column
scaled A is estimated and returned in RWORK(1).
NΛ(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <=
NΛ(1/4)*RWORK(1)
JOBP
JOBP is
CHARACTER*1
= βPβ The rows of A are ordered in decreasing
order with respect to
||A(i,:)||_\infty. This enhances numerical accuracy at the
cost
of extra data movement. Recommended for numerical
robustness.
= βNβ No row pivoting.
JOBR
JOBR is
CHARACTER*1
= βTβ After the initial pivoted QR
factorization, DGESVD is applied to
the transposed R**T of the computed triangular factor R.
This involves
some extra data movement (matrix transpositions). Useful for
experiments, research and development.
= βNβ The triangular factor R is given as input
to DGESVD. This may be
preferred as it involves less data movement.
JOBU
JOBU is
CHARACTER*1
= βAβ All M left singular vectors are computed
and returned in the
matrix U. See the description of U.
= βSβ or βUβ N = min(M,N) left
singular vectors are computed and returned
in the matrix U. See the description of U.
= βRβ Numerical rank NUMRANK is determined and
only NUMRANK left singular
vectors are computed and returned in the matrix U.
= βFβ The N left singular vectors are returned
in factored form as the
product of the Q factor from the initial QR factorization
and the
N left singular vectors of (R**T , 0)**T. If row pivoting is
used,
then the necessary information on the row pivoting is stored
in
IWORK(N+1:N+M-1).
= βNβ The left singular vectors are not
computed.
JOBV
JOBV is
CHARACTER*1
= βAβ, βVβ All N right singular
vectors are computed and returned in
the matrix V.
= βRβ Numerical rank NUMRANK is determined and
only NUMRANK right singular
vectors are computed and returned in the matrix V. This
option is
allowed only if JOBU = βRβ or JOBU =
βNβ; otherwise it is illegal.
= βNβ The right singular vectors are not
computed.
M
M is INTEGER
The number of rows of the input matrix A. M >= 0.
N
N is INTEGER
The number of columns of the input matrix A. M >= N >=
0.
A
A is DOUBLE
PRECISION array of dimensions LDA x N
On entry, the input matrix A.
On exit, if JOBU .NE. βNβ or JOBV .NE.
βNβ, the lower triangle of A contains
the Householder vectors as stored by DGEQP3. If JOBU =
βFβ, these Householder
vectors together with WORK(1:N) can be used to restore the Q
factors from
the initial pivoted QR factorization of A. See the
description of U.
LDA
LDA is INTEGER.
The leading dimension of the array A. LDA >=
max(1,M).
S
S is DOUBLE
PRECISION array of dimension N.
The singular values of A, ordered so that S(i) >=
S(i+1).
U
U is DOUBLE
PRECISION array, dimension
LDU x M if JOBU = βAβ; see the description of
LDU. In this case,
on exit, U contains the M left singular vectors.
LDU x N if JOBU = βSβ, βUβ,
βRβ ; see the description of LDU. In this
case, U contains the leading N or the leading NUMRANK left
singular vectors.
LDU x N if JOBU = βFβ ; see the description of
LDU. In this case U
contains N x N orthogonal matrix that can be used to form
the left
singular vectors.
If JOBU = βNβ, U is not referenced.
LDU
LDU is INTEGER.
The leading dimension of the array U.
If JOBU = βAβ, βSβ, βUβ,
βRβ, LDU >= max(1,M).
If JOBU = βFβ, LDU >= max(1,N).
Otherwise, LDU >= 1.
V
V is DOUBLE
PRECISION array, dimension
LDV x N if JOBV = βAβ, βVβ,
βRβ or if JOBA = βEβ .
If JOBV = βAβ, or βVβ, V contains
the N-by-N orthogonal matrix V**T;
If JOBV = βRβ, V contains the first NUMRANK rows
of V**T (the right
singular vectors, stored rowwise, of the NUMRANK largest
singular values).
If JOBV = βNβ and JOBA = βEβ, V is
used as a workspace.
If JOBV = βNβ, and JOBA.NE.βEβ, V is
not referenced.
LDV
LDV is INTEGER
The leading dimension of the array V.
If JOBV = βAβ, βVβ, βRβ,
or JOBA = βEβ, LDV >= max(1,N).
Otherwise, LDV >= 1.
NUMRANK
NUMRANK is
INTEGER
NUMRANK is the numerical rank first determined after the
rank
revealing QR factorization, following the strategy specified
by the
value of JOBA. If JOBV = βRβ and JOBU =
βRβ, only NUMRANK
leading singular values and vectors are then requested in
the call
of DGESVD. The final value of NUMRANK might be further
reduced if
some singular values are computed as zeros.
IWORK
IWORK is
INTEGER array, dimension (max(1, LIWORK)).
On exit, IWORK(1:N) contains column pivoting permutation of
the
rank revealing QR factorization.
If JOBP = βPβ, IWORK(N+1:N+M-1) contains the
indices of the sequence
of row swaps used in row pivoting. These can be used to
restore the
left singular vectors in the case JOBU =
βFβ.
If LIWORK,
LWORK, or LRWORK = -1, then on exit, if INFO = 0,
IWORK(1) returns the minimal LIWORK.
LIWORK
LIWORK is
INTEGER
The dimension of the array IWORK.
LIWORK >= N + M - 1, if JOBP = βPβ and JOBA
.NE. βEβ;
LIWORK >= N if JOBP = βNβ and JOBA .NE.
βEβ;
LIWORK >= N + M - 1 + N, if JOBP = βPβ and
JOBA = βEβ;
LIWORK >= N + N if JOBP = βNβ and JOBA =
βEβ.
If LIWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the WORK, IWORK, and RWORK arrays, and no error
message related to LWORK is issued by XERBLA.
WORK
WORK is DOUBLE
PRECISION array, dimension (max(2, LWORK)), used as a
workspace.
On exit, if, on entry, LWORK.NE.-1, WORK(1:N) contains
parameters
needed to recover the Q factor from the QR factorization
computed by
DGEQP3.
If LIWORK,
LWORK, or LRWORK = -1, then on exit, if INFO = 0,
WORK(1) returns the optimal LWORK, and
WORK(2) returns the minimal LWORK.
LWORK
LWORK is
INTEGER
The dimension of the array WORK. It is determined as
follows:
Let LWQP3 = 3*N+1, LWCON = 3*N, and let
LWORQ = { MAX( N, 1 ), if JOBU = βRβ,
βSβ, or βUβ
{ MAX( M, 1 ), if JOBU = βAβ
LWSVD = MAX( 5*N, 1 )
LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 5*(N/2), 1 ), LWORLQ =
MAX( N, 1 ),
LWQRF = MAX( N/2, 1 ), LWORQ2 = MAX( N, 1 )
Then the minimal value of LWORK is:
= MAX( N + LWQP3, LWSVD ) if only the singular values are
needed;
= MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values
are needed,
and a scaled condition estimate requested;
= N + MAX(
LWQP3, LWSVD, LWORQ ) if the singular values and the left
singular vectors are requested;
= N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the singular
values and the left
singular vectors are requested, and also
a scaled condition estimate requested;
= N + MAX(
LWQP3, LWSVD ) if the singular values and the right
singular vectors are requested;
= N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and
the right
singular vectors are requested, and also
a scaled condition etimate requested;
= N + MAX(
LWQP3, LWSVD, LWORQ ) if the full SVD is requested with JOBV
= βRβ;
independent of JOBR;
= N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the full SVD is
requested,
JOBV = βRβ and, also a scaled condition
estimate requested; independent of JOBR;
= MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ) )
if the
full SVD is requested with JOBV = βAβ or
βVβ, and
JOBR =βNβ
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ,
LWORQ ) )
if the full SVD is requested with JOBV = βAβ or
βVβ, and
JOBR =βNβ, and also a scaled condition number
estimate
requested.
= MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) )
if the
full SVD is requested with JOBV = βAβ,
βVβ, and JOBR =βTβ
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2,
LWORQ ) )
if the full SVD is requested with JOBV = βAβ or
βVβ, and
JOBR =βTβ, and also a scaled condition number
estimate
requested.
Finally, LWORK must be at least two: LWORK = MAX( 2, LWORK
).
If LWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the WORK, IWORK, and RWORK arrays, and no error
message related to LWORK is issued by XERBLA.
RWORK
RWORK is DOUBLE
PRECISION array, dimension (max(1, LRWORK)).
On exit,
1. If JOBA = βEβ, RWORK(1) contains an estimate
of the condition
number of column scaled A. If A = C * D where D is diagonal
and C
has unit columns in the Euclidean norm, then, assuming full
column rank,
NΛ(-1/4) * RWORK(1) <= ||pinv(C)||_2 <=
NΛ(1/4) * RWORK(1).
Otherwise, RWORK(1) = -1.
2. RWORK(2) contains the number of singular values computed
as
exact zeros in DGESVD applied to the upper triangular or
trapezoidal
R (from the initial QR factorization). In case of early exit
(no call to
DGESVD, such as in the case of zero matrix) RWORK(2) =
-1.
If LIWORK,
LWORK, or LRWORK = -1, then on exit, if INFO = 0,
RWORK(1) returns the minimal LRWORK.
LRWORK
LRWORK is
INTEGER.
The dimension of the array RWORK.
If JOBP =βPβ, then LRWORK >= MAX(2, M).
Otherwise, LRWORK >= 2
If LRWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the WORK, IWORK, and RWORK arrays, and no error
message related to LWORK is issued by XERBLA.
INFO
INFO is INTEGER
= 0: successful exit.
< 0: if INFO = -i, the i-th argument had an illegal
value.
> 0: if DBDSQR did not converge, INFO specifies how many
superdiagonals
of an intermediate bidiagonal form B (computed in DGESVD)
did not
converge to zero.
Further Details:
1. The data
movement (matrix transpose) is coded using simple nested
DO-loops because BLAS and LAPACK do not provide
corresponding subroutines.
Those DO-loops are easily identified in this source code -
by the CONTINUE
statements labeled with 11**. In an optimized version of
this code, the
nested DO loops should be replaced with calls to an
optimized subroutine.
2. This code scales A by 1/SQRT(M) if the largest
ABS(A(i,j)) could cause
column norm overflow. This is the minial precaution and it
is left to the
SVD routine (CGESVD) to do its own preemptive scaling if
potential over-
or underflows are detected. To avoid repeated scanning of
the array A,
an optimal implementation would do all necessary scaling
before calling
CGESVD and the scaling in CGESVD can be switched off.
3. Other comments related to code optimization are given in
comments in the
code, enclosed in [[double brackets]].
Bugs, examples and comments
Please report
all bugs and send interesting examples and/or comments to
drmac@math.hr. Thank you.
References
[1] Zlatko
Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
Computing the SVD with High Accuracy. ACM Trans. Math.
Softw.
44(1): 11:1-11:30 (2017)
SIGMA library,
xGESVDQ section updated February 2016.
Developed and coded by Zlatko Drmac, Department of
Mathematics
University of Zagreb, Croatia, drmac@math.hr
Contributors:
Developed and
coded by Zlatko Drmac, Department of Mathematics
University of Zagreb, Croatia, drmac@math.hr
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
subroutine sgesvdq (character joba, character jobp, character jobr,character jobu, character jobv, integer m, integer n, real, dimension(lda, * ) a, integer lda, real, dimension( * ) s, real, dimension( ldu,* ) u, integer ldu, real, dimension( ldv, * ) v, integer ldv, integernumrank, integer, dimension( * ) iwork, integer liwork, real,dimension( * ) work, integer lwork, real, dimension( * ) rwork, integerlrwork, integer info)
SGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices
Purpose:
SGESVDQ
computes the singular value decomposition (SVD) of a real
M-by-N matrix A, where M >= N. The SVD of A is written as
[++] [xx] [x0] [xx]
A = U * SIGMA * VΛ*, [++] = [xx] * [ox] * [xx]
[++] [xx]
where SIGMA is an N-by-N diagonal matrix, U is an M-by-N
orthonormal
matrix, and V is an N-by-N orthogonal matrix. The diagonal
elements
of SIGMA are the singular values of A. The columns of U and
V are the
left and the right singular vectors of A, respectively.
Parameters
JOBA
JOBA is
CHARACTER*1
Specifies the level of accuracy in the computed SVD
= βAβ The requested accuracy corresponds to
having the backward
error bounded by || delta A ||_F <= f(m,n) * EPS * || A
||_F,
where EPS = SLAMCH(βEpsilonβ). This authorises
CGESVDQ to
truncate the computed triangular factor in a rank revealing
QR factorization whenever the truncated part is below the
threshold of the order of EPS * ||A||_F. This is aggressive
truncation level.
= βMβ Similarly as with βAβ, but the
truncation is more gentle: it
is allowed only when there is a drop on the diagonal of the
triangular factor in the QR factorization. This is medium
truncation level.
= βHβ High accuracy requested. No numerical rank
determination based
on the rank revealing QR factorization is attempted.
= βEβ Same as βHβ, and in addition
the condition number of column
scaled A is estimated and returned in RWORK(1).
NΛ(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <=
NΛ(1/4)*RWORK(1)
JOBP
JOBP is
CHARACTER*1
= βPβ The rows of A are ordered in decreasing
order with respect to
||A(i,:)||_\infty. This enhances numerical accuracy at the
cost
of extra data movement. Recommended for numerical
robustness.
= βNβ No row pivoting.
JOBR
JOBR is
CHARACTER*1
= βTβ After the initial pivoted QR
factorization, SGESVD is applied to
the transposed R**T of the computed triangular factor R.
This involves
some extra data movement (matrix transpositions). Useful for
experiments, research and development.
= βNβ The triangular factor R is given as input
to SGESVD. This may be
preferred as it involves less data movement.
JOBU
JOBU is
CHARACTER*1
= βAβ All M left singular vectors are computed
and returned in the
matrix U. See the description of U.
= βSβ or βUβ N = min(M,N) left
singular vectors are computed and returned
in the matrix U. See the description of U.
= βRβ Numerical rank NUMRANK is determined and
only NUMRANK left singular
vectors are computed and returned in the matrix U.
= βFβ The N left singular vectors are returned
in factored form as the
product of the Q factor from the initial QR factorization
and the
N left singular vectors of (R**T , 0)**T. If row pivoting is
used,
then the necessary information on the row pivoting is stored
in
IWORK(N+1:N+M-1).
= βNβ The left singular vectors are not
computed.
JOBV
JOBV is
CHARACTER*1
= βAβ, βVβ All N right singular
vectors are computed and returned in
the matrix V.
= βRβ Numerical rank NUMRANK is determined and
only NUMRANK right singular
vectors are computed and returned in the matrix V. This
option is
allowed only if JOBU = βRβ or JOBU =
βNβ; otherwise it is illegal.
= βNβ The right singular vectors are not
computed.
M
M is INTEGER
The number of rows of the input matrix A. M >= 0.
N
N is INTEGER
The number of columns of the input matrix A. M >= N >=
0.
A
A is REAL array
of dimensions LDA x N
On entry, the input matrix A.
On exit, if JOBU .NE. βNβ or JOBV .NE.
βNβ, the lower triangle of A contains
the Householder vectors as stored by SGEQP3. If JOBU =
βFβ, these Householder
vectors together with WORK(1:N) can be used to restore the Q
factors from
the initial pivoted QR factorization of A. See the
description of U.
LDA
LDA is INTEGER.
The leading dimension of the array A. LDA >=
max(1,M).
S
S is REAL array
of dimension N.
The singular values of A, ordered so that S(i) >=
S(i+1).
U
U is REAL
array, dimension
LDU x M if JOBU = βAβ; see the description of
LDU. In this case,
on exit, U contains the M left singular vectors.
LDU x N if JOBU = βSβ, βUβ,
βRβ ; see the description of LDU. In this
case, U contains the leading N or the leading NUMRANK left
singular vectors.
LDU x N if JOBU = βFβ ; see the description of
LDU. In this case U
contains N x N orthogonal matrix that can be used to form
the left
singular vectors.
If JOBU = βNβ, U is not referenced.
LDU
LDU is INTEGER.
The leading dimension of the array U.
If JOBU = βAβ, βSβ, βUβ,
βRβ, LDU >= max(1,M).
If JOBU = βFβ, LDU >= max(1,N).
Otherwise, LDU >= 1.
V
V is REAL
array, dimension
LDV x N if JOBV = βAβ, βVβ,
βRβ or if JOBA = βEβ .
If JOBV = βAβ, or βVβ, V contains
the N-by-N orthogonal matrix V**T;
If JOBV = βRβ, V contains the first NUMRANK rows
of V**T (the right
singular vectors, stored rowwise, of the NUMRANK largest
singular values).
If JOBV = βNβ and JOBA = βEβ, V is
used as a workspace.
If JOBV = βNβ, and JOBA.NE.βEβ, V is
not referenced.
LDV
LDV is INTEGER
The leading dimension of the array V.
If JOBV = βAβ, βVβ, βRβ,
or JOBA = βEβ, LDV >= max(1,N).
Otherwise, LDV >= 1.
NUMRANK
NUMRANK is
INTEGER
NUMRANK is the numerical rank first determined after the
rank
revealing QR factorization, following the strategy specified
by the
value of JOBA. If JOBV = βRβ and JOBU =
βRβ, only NUMRANK
leading singular values and vectors are then requested in
the call
of SGESVD. The final value of NUMRANK might be further
reduced if
some singular values are computed as zeros.
IWORK
IWORK is
INTEGER array, dimension (max(1, LIWORK)).
On exit, IWORK(1:N) contains column pivoting permutation of
the
rank revealing QR factorization.
If JOBP = βPβ, IWORK(N+1:N+M-1) contains the
indices of the sequence
of row swaps used in row pivoting. These can be used to
restore the
left singular vectors in the case JOBU =
βFβ.
If LIWORK,
LWORK, or LRWORK = -1, then on exit, if INFO = 0,
IWORK(1) returns the minimal LIWORK.
LIWORK
LIWORK is
INTEGER
The dimension of the array IWORK.
LIWORK >= N + M - 1, if JOBP = βPβ and JOBA
.NE. βEβ;
LIWORK >= N if JOBP = βNβ and JOBA .NE.
βEβ;
LIWORK >= N + M - 1 + N, if JOBP = βPβ and
JOBA = βEβ;
LIWORK >= N + N if JOBP = βNβ and JOBA =
βEβ.
If LIWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the WORK, IWORK, and RWORK arrays, and no error
message related to LWORK is issued by XERBLA.
WORK
WORK is REAL
array, dimension (max(2, LWORK)), used as a workspace.
On exit, if, on entry, LWORK.NE.-1, WORK(1:N) contains
parameters
needed to recover the Q factor from the QR factorization
computed by
SGEQP3.
If LIWORK,
LWORK, or LRWORK = -1, then on exit, if INFO = 0,
WORK(1) returns the optimal LWORK, and
WORK(2) returns the minimal LWORK.
LWORK
LWORK is
INTEGER
The dimension of the array WORK. It is determined as
follows:
Let LWQP3 = 3*N+1, LWCON = 3*N, and let
LWORQ = { MAX( N, 1 ), if JOBU = βRβ,
βSβ, or βUβ
{ MAX( M, 1 ), if JOBU = βAβ
LWSVD = MAX( 5*N, 1 )
LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 5*(N/2), 1 ), LWORLQ =
MAX( N, 1 ),
LWQRF = MAX( N/2, 1 ), LWORQ2 = MAX( N, 1 )
Then the minimal value of LWORK is:
= MAX( N + LWQP3, LWSVD ) if only the singular values are
needed;
= MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values
are needed,
and a scaled condition estimate requested;
= N + MAX(
LWQP3, LWSVD, LWORQ ) if the singular values and the left
singular vectors are requested;
= N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the singular
values and the left
singular vectors are requested, and also
a scaled condition estimate requested;
= N + MAX(
LWQP3, LWSVD ) if the singular values and the right
singular vectors are requested;
= N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and
the right
singular vectors are requested, and also
a scaled condition etimate requested;
= N + MAX(
LWQP3, LWSVD, LWORQ ) if the full SVD is requested with JOBV
= βRβ;
independent of JOBR;
= N + MAX( LWQP3, LWCON, LWSVD, LWORQ ) if the full SVD is
requested,
JOBV = βRβ and, also a scaled condition
estimate requested; independent of JOBR;
= MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ, LWORQ) )
if the
full SVD is requested with JOBV = βAβ or
βVβ, and
JOBR =βNβ
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWORLQ,
LWORQ ) )
if the full SVD is requested with JOBV = βAβ or
βVβ, and
JOBR =βNβ, and also a scaled condition number
estimate
requested.
= MAX( N + MAX( LWQP3, LWSVD, LWORQ ),
N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2, LWORQ ) )
if the
full SVD is requested with JOBV = βAβ,
βVβ, and JOBR =βTβ
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWORQ ),
N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWORQ2,
LWORQ ) )
if the full SVD is requested with JOBV = βAβ or
βVβ, and
JOBR =βTβ, and also a scaled condition number
estimate
requested.
Finally, LWORK must be at least two: LWORK = MAX( 2, LWORK
).
If LWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the WORK, IWORK, and RWORK arrays, and no error
message related to LWORK is issued by XERBLA.
RWORK
RWORK is REAL
array, dimension (max(1, LRWORK)).
On exit,
1. If JOBA = βEβ, RWORK(1) contains an estimate
of the condition
number of column scaled A. If A = C * D where D is diagonal
and C
has unit columns in the Euclidean norm, then, assuming full
column rank,
NΛ(-1/4) * RWORK(1) <= ||pinv(C)||_2 <=
NΛ(1/4) * RWORK(1).
Otherwise, RWORK(1) = -1.
2. RWORK(2) contains the number of singular values computed
as
exact zeros in SGESVD applied to the upper triangular or
trapezoidal
R (from the initial QR factorization). In case of early exit
(no call to
SGESVD, such as in the case of zero matrix) RWORK(2) =
-1.
If LIWORK,
LWORK, or LRWORK = -1, then on exit, if INFO = 0,
RWORK(1) returns the minimal LRWORK.
LRWORK
LRWORK is
INTEGER.
The dimension of the array RWORK.
If JOBP =βPβ, then LRWORK >= MAX(2, M).
Otherwise, LRWORK >= 2
If LRWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the WORK, IWORK, and RWORK arrays, and no error
message related to LWORK is issued by XERBLA.
INFO
INFO is INTEGER
= 0: successful exit.
< 0: if INFO = -i, the i-th argument had an illegal
value.
> 0: if SBDSQR did not converge, INFO specifies how many
superdiagonals
of an intermediate bidiagonal form B (computed in SGESVD)
did not
converge to zero.
Further Details:
1. The data
movement (matrix transpose) is coded using simple nested
DO-loops because BLAS and LAPACK do not provide
corresponding subroutines.
Those DO-loops are easily identified in this source code -
by the CONTINUE
statements labeled with 11**. In an optimized version of
this code, the
nested DO loops should be replaced with calls to an
optimized subroutine.
2. This code scales A by 1/SQRT(M) if the largest
ABS(A(i,j)) could cause
column norm overflow. This is the minial precaution and it
is left to the
SVD routine (CGESVD) to do its own preemptive scaling if
potential over-
or underflows are detected. To avoid repeated scanning of
the array A,
an optimal implementation would do all necessary scaling
before calling
CGESVD and the scaling in CGESVD can be switched off.
3. Other comments related to code optimization are given in
comments in the
code, enclosed in [[double brackets]].
Bugs, examples and comments
Please report
all bugs and send interesting examples and/or comments to
drmac@math.hr. Thank you.
References
[1] Zlatko
Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
Computing the SVD with High Accuracy. ACM Trans. Math.
Softw.
44(1): 11:1-11:30 (2017)
SIGMA library,
xGESVDQ section updated February 2016.
Developed and coded by Zlatko Drmac, Department of
Mathematics
University of Zagreb, Croatia, drmac@math.hr
Contributors:
Developed and
coded by Zlatko Drmac, Department of Mathematics
University of Zagreb, Croatia, drmac@math.hr
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
subroutine zgesvdq (character joba, character jobp, character jobr,character jobu, character jobv, integer m, integer n, complex*16,dimension( lda, * ) a, integer lda, double precision, dimension( * ) s,complex*16, dimension( ldu, * ) u, integer ldu, complex*16, dimension(ldv, * ) v, integer ldv, integer numrank, integer, dimension( * )iwork, integer liwork, complex*16, dimension( * ) cwork, integerlcwork, double precision, dimension( * ) rwork, integer lrwork, integerinfo)
ZGESVDQ computes the singular value decomposition (SVD) with a QR-Preconditioned QR SVD Method for GE matrices
Purpose:
ZCGESVDQ
computes the singular value decomposition (SVD) of a complex
M-by-N matrix A, where M >= N. The SVD of A is written as
[++] [xx] [x0] [xx]
A = U * SIGMA * VΛ*, [++] = [xx] * [ox] * [xx]
[++] [xx]
where SIGMA is an N-by-N diagonal matrix, U is an M-by-N
orthonormal
matrix, and V is an N-by-N unitary matrix. The diagonal
elements
of SIGMA are the singular values of A. The columns of U and
V are the
left and the right singular vectors of A, respectively.
Parameters
JOBA
JOBA is
CHARACTER*1
Specifies the level of accuracy in the computed SVD
= βAβ The requested accuracy corresponds to
having the backward
error bounded by || delta A ||_F <= f(m,n) * EPS * || A
||_F,
where EPS = DLAMCH(βEpsilonβ). This authorises
ZGESVDQ to
truncate the computed triangular factor in a rank revealing
QR factorization whenever the truncated part is below the
threshold of the order of EPS * ||A||_F. This is aggressive
truncation level.
= βMβ Similarly as with βAβ, but the
truncation is more gentle: it
is allowed only when there is a drop on the diagonal of the
triangular factor in the QR factorization. This is medium
truncation level.
= βHβ High accuracy requested. No numerical rank
determination based
on the rank revealing QR factorization is attempted.
= βEβ Same as βHβ, and in addition
the condition number of column
scaled A is estimated and returned in RWORK(1).
NΛ(-1/4)*RWORK(1) <= ||pinv(A_scaled)||_2 <=
NΛ(1/4)*RWORK(1)
JOBP
JOBP is
CHARACTER*1
= βPβ The rows of A are ordered in decreasing
order with respect to
||A(i,:)||_\infty. This enhances numerical accuracy at the
cost
of extra data movement. Recommended for numerical
robustness.
= βNβ No row pivoting.
JOBR
JOBR is
CHARACTER*1
= βTβ After the initial pivoted QR
factorization, ZGESVD is applied to
the adjoint R**H of the computed triangular factor R. This
involves
some extra data movement (matrix transpositions). Useful for
experiments, research and development.
= βNβ The triangular factor R is given as input
to CGESVD. This may be
preferred as it involves less data movement.
JOBU
JOBU is
CHARACTER*1
= βAβ All M left singular vectors are computed
and returned in the
matrix U. See the description of U.
= βSβ or βUβ N = min(M,N) left
singular vectors are computed and returned
in the matrix U. See the description of U.
= βRβ Numerical rank NUMRANK is determined and
only NUMRANK left singular
vectors are computed and returned in the matrix U.
= βFβ The N left singular vectors are returned
in factored form as the
product of the Q factor from the initial QR factorization
and the
N left singular vectors of (R**H , 0)**H. If row pivoting is
used,
then the necessary information on the row pivoting is stored
in
IWORK(N+1:N+M-1).
= βNβ The left singular vectors are not
computed.
JOBV
JOBV is
CHARACTER*1
= βAβ, βVβ All N right singular
vectors are computed and returned in
the matrix V.
= βRβ Numerical rank NUMRANK is determined and
only NUMRANK right singular
vectors are computed and returned in the matrix V. This
option is
allowed only if JOBU = βRβ or JOBU =
βNβ; otherwise it is illegal.
= βNβ The right singular vectors are not
computed.
M
M is INTEGER
The number of rows of the input matrix A. M >= 0.
N
N is INTEGER
The number of columns of the input matrix A. M >= N >=
0.
A
A is COMPLEX*16
array of dimensions LDA x N
On entry, the input matrix A.
On exit, if JOBU .NE. βNβ or JOBV .NE.
βNβ, the lower triangle of A contains
the Householder vectors as stored by ZGEQP3. If JOBU =
βFβ, these Householder
vectors together with CWORK(1:N) can be used to restore the
Q factors from
the initial pivoted QR factorization of A. See the
description of U.
LDA
LDA is INTEGER.
The leading dimension of the array A. LDA >=
max(1,M).
S
S is DOUBLE
PRECISION array of dimension N.
The singular values of A, ordered so that S(i) >=
S(i+1).
U
U is COMPLEX*16
array, dimension
LDU x M if JOBU = βAβ; see the description of
LDU. In this case,
on exit, U contains the M left singular vectors.
LDU x N if JOBU = βSβ, βUβ,
βRβ ; see the description of LDU. In this
case, U contains the leading N or the leading NUMRANK left
singular vectors.
LDU x N if JOBU = βFβ ; see the description of
LDU. In this case U
contains N x N unitary matrix that can be used to form the
left
singular vectors.
If JOBU = βNβ, U is not referenced.
LDU
LDU is INTEGER.
The leading dimension of the array U.
If JOBU = βAβ, βSβ, βUβ,
βRβ, LDU >= max(1,M).
If JOBU = βFβ, LDU >= max(1,N).
Otherwise, LDU >= 1.
V
V is COMPLEX*16
array, dimension
LDV x N if JOBV = βAβ, βVβ,
βRβ or if JOBA = βEβ .
If JOBV = βAβ, or βVβ, V contains
the N-by-N unitary matrix V**H;
If JOBV = βRβ, V contains the first NUMRANK rows
of V**H (the right
singular vectors, stored rowwise, of the NUMRANK largest
singular values).
If JOBV = βNβ and JOBA = βEβ, V is
used as a workspace.
If JOBV = βNβ, and JOBA.NE.βEβ, V is
not referenced.
LDV
LDV is INTEGER
The leading dimension of the array V.
If JOBV = βAβ, βVβ, βRβ,
or JOBA = βEβ, LDV >= max(1,N).
Otherwise, LDV >= 1.
NUMRANK
NUMRANK is
INTEGER
NUMRANK is the numerical rank first determined after the
rank
revealing QR factorization, following the strategy specified
by the
value of JOBA. If JOBV = βRβ and JOBU =
βRβ, only NUMRANK
leading singular values and vectors are then requested in
the call
of CGESVD. The final value of NUMRANK might be further
reduced if
some singular values are computed as zeros.
IWORK
IWORK is
INTEGER array, dimension (max(1, LIWORK)).
On exit, IWORK(1:N) contains column pivoting permutation of
the
rank revealing QR factorization.
If JOBP = βPβ, IWORK(N+1:N+M-1) contains the
indices of the sequence
of row swaps used in row pivoting. These can be used to
restore the
left singular vectors in the case JOBU =
βFβ.
If LIWORK,
LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
IWORK(1) returns the minimal LIWORK.
LIWORK
LIWORK is
INTEGER
The dimension of the array IWORK.
LIWORK >= N + M - 1, if JOBP = βPβ;
LIWORK >= N if JOBP = βNβ.
If LIWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the CWORK, IWORK, and RWORK arrays, and no error
message related to LCWORK is issued by XERBLA.
CWORK
CWORK is
COMPLEX*12 array, dimension (max(2, LCWORK)), used as a
workspace.
On exit, if, on entry, LCWORK.NE.-1, CWORK(1:N) contains
parameters
needed to recover the Q factor from the QR factorization
computed by
ZGEQP3.
If LIWORK,
LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
CWORK(1) returns the optimal LCWORK, and
CWORK(2) returns the minimal LCWORK.
LCWORK
LCWORK is
INTEGER
The dimension of the array CWORK. It is determined as
follows:
Let LWQP3 = N+1, LWCON = 2*N, and let
LWUNQ = { MAX( N, 1 ), if JOBU = βRβ,
βSβ, or βUβ
{ MAX( M, 1 ), if JOBU = βAβ
LWSVD = MAX( 3*N, 1 )
LWLQF = MAX( N/2, 1 ), LWSVD2 = MAX( 3*(N/2), 1 ), LWUNLQ =
MAX( N, 1 ),
LWQRF = MAX( N/2, 1 ), LWUNQ2 = MAX( N, 1 )
Then the minimal value of LCWORK is:
= MAX( N + LWQP3, LWSVD ) if only the singular values are
needed;
= MAX( N + LWQP3, LWCON, LWSVD ) if only the singular values
are needed,
and a scaled condition estimate requested;
= N + MAX(
LWQP3, LWSVD, LWUNQ ) if the singular values and the left
singular vectors are requested;
= N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the singular
values and the left
singular vectors are requested, and also
a scaled condition estimate requested;
= N + MAX(
LWQP3, LWSVD ) if the singular values and the right
singular vectors are requested;
= N + MAX( LWQP3, LWCON, LWSVD ) if the singular values and
the right
singular vectors are requested, and also
a scaled condition etimate requested;
= N + MAX(
LWQP3, LWSVD, LWUNQ ) if the full SVD is requested with JOBV
= βRβ;
independent of JOBR;
= N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ) if the full SVD is
requested,
JOBV = βRβ and, also a scaled condition
estimate requested; independent of JOBR;
= MAX( N + MAX( LWQP3, LWSVD, LWUNQ ),
N + MAX( LWQP3, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ, LWUNQ) )
if the
full SVD is requested with JOBV = βAβ or
βVβ, and
JOBR =βNβ
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
N + MAX( LWQP3, LWCON, N/2+LWLQF, N/2+LWSVD2, N/2+LWUNLQ,
LWUNQ ) )
if the full SVD is requested with JOBV = βAβ or
βVβ, and
JOBR =βNβ, and also a scaled condition number
estimate
requested.
= MAX( N + MAX( LWQP3, LWSVD, LWUNQ ),
N + MAX( LWQP3, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2, LWUNQ ) )
if the
full SVD is requested with JOBV = βAβ,
βVβ, and JOBR =βTβ
= MAX( N + MAX( LWQP3, LWCON, LWSVD, LWUNQ ),
N + MAX( LWQP3, LWCON, N/2+LWQRF, N/2+LWSVD2, N/2+LWUNQ2,
LWUNQ ) )
if the full SVD is requested with JOBV = βAβ,
βVβ and
JOBR =βTβ, and also a scaled condition number
estimate
requested.
Finally, LCWORK must be at least two: LCWORK = MAX( 2,
LCWORK ).
If LCWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the CWORK, IWORK, and RWORK arrays, and no error
message related to LCWORK is issued by XERBLA.
RWORK
RWORK is DOUBLE
PRECISION array, dimension (max(1, LRWORK)).
On exit,
1. If JOBA = βEβ, RWORK(1) contains an estimate
of the condition
number of column scaled A. If A = C * D where D is diagonal
and C
has unit columns in the Euclidean norm, then, assuming full
column rank,
NΛ(-1/4) * RWORK(1) <= ||pinv(C)||_2 <=
NΛ(1/4) * RWORK(1).
Otherwise, RWORK(1) = -1.
2. RWORK(2) contains the number of singular values computed
as
exact zeros in ZGESVD applied to the upper triangular or
trapezoidal
R (from the initial QR factorization). In case of early exit
(no call to
ZGESVD, such as in the case of zero matrix) RWORK(2) =
-1.
If LIWORK,
LCWORK, or LRWORK = -1, then on exit, if INFO = 0,
RWORK(1) returns the minimal LRWORK.
LRWORK
LRWORK is
INTEGER.
The dimension of the array RWORK.
If JOBP =βPβ, then LRWORK >= MAX(2, M, 5*N);
Otherwise, LRWORK >= MAX(2, 5*N).
If LRWORK = -1,
then a workspace query is assumed; the routine
only calculates and returns the optimal and minimal sizes
for the CWORK, IWORK, and RWORK arrays, and no error
message related to LCWORK is issued by XERBLA.
INFO
INFO is INTEGER
= 0: successful exit.
< 0: if INFO = -i, the i-th argument had an illegal
value.
> 0: if ZBDSQR did not converge, INFO specifies how many
superdiagonals
of an intermediate bidiagonal form B (computed in ZGESVD)
did not
converge to zero.
Further Details:
1. The data
movement (matrix transpose) is coded using simple nested
DO-loops because BLAS and LAPACK do not provide
corresponding subroutines.
Those DO-loops are easily identified in this source code -
by the CONTINUE
statements labeled with 11**. In an optimized version of
this code, the
nested DO loops should be replaced with calls to an
optimized subroutine.
2. This code scales A by 1/SQRT(M) if the largest
ABS(A(i,j)) could cause
column norm overflow. This is the minial precaution and it
is left to the
SVD routine (CGESVD) to do its own preemptive scaling if
potential over-
or underflows are detected. To avoid repeated scanning of
the array A,
an optimal implementation would do all necessary scaling
before calling
CGESVD and the scaling in CGESVD can be switched off.
3. Other comments related to code optimization are given in
comments in the
code, enclosed in [[double brackets]].
Bugs, examples and comments
Please report
all bugs and send interesting examples and/or comments to
drmac@math.hr. Thank you.
References
[1] Zlatko
Drmac, Algorithm 977: A QR-Preconditioned QR SVD Method for
Computing the SVD with High Accuracy. ACM Trans. Math.
Softw.
44(1): 11:1-11:30 (2017)
SIGMA library,
xGESVDQ section updated February 2016.
Developed and coded by Zlatko Drmac, Department of
Mathematics
University of Zagreb, Croatia, drmac@math.hr
Contributors:
Developed and
coded by Zlatko Drmac, Department of Mathematics
University of Zagreb, Croatia, drmac@math.hr
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Author
Generated automatically by Doxygen for LAPACK from the source code.