/* Subroutine */ int dlasd8_(integer *icompq, integer *k, doublereal *d__,
doublereal *z__, doublereal *vf, doublereal *vl, doublereal *difl,
doublereal *difr, integer *lddifr, doublereal *dsigma, doublereal *
work, integer *info)
{
/* System generated locals */
integer difr_dim1, difr_offset, i__1, i__2;
doublereal d__1, d__2;
/* Builtin functions */
double sqrt(doublereal), d_sign(doublereal *, doublereal *);
/* Local variables */
integer i__, j;
doublereal dj, rho;
integer iwk1, iwk2, iwk3;
extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *,
integer *);
doublereal temp;
extern doublereal dnrm2_(integer *, doublereal *, integer *);
integer iwk2i, iwk3i;
doublereal diflj, difrj, dsigj;
extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
doublereal *, integer *);
extern doublereal dlamc3_(doublereal *, doublereal *);
extern /* Subroutine */ int dlasd4_(integer *, integer *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, integer *), dlascl_(char *, integer *, integer *,
doublereal *, doublereal *, integer *, integer *, doublereal *,
integer *, integer *), dlaset_(char *, integer *, integer
*, doublereal *, doublereal *, doublereal *, integer *),
xerbla_(char *, integer *);
doublereal dsigjp;
/* -- LAPACK auxiliary routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DLASD8 finds the square roots of the roots of the secular equation, */
/* as defined by the values in DSIGMA and Z. It makes the appropriate */
/* calls to DLASD4, and stores, for each element in D, the distance */
/* to its two nearest poles (elements in DSIGMA). It also updates */
/* the arrays VF and VL, the first and last components of all the */
/* right singular vectors of the original bidiagonal matrix. */
/* DLASD8 is called from DLASD6. */
/* Arguments */
/* ========= */
/* ICOMPQ (input) INTEGER */
/* Specifies whether singular vectors are to be computed in */
/* factored form in the calling routine: */
/* = 0: Compute singular values only. */
/* = 1: Compute singular vectors in factored form as well. */
/* K (input) INTEGER */
/* The number of terms in the rational function to be solved */
/* by DLASD4. K >= 1. */
/* D (output) DOUBLE PRECISION array, dimension ( K ) */
/* On output, D contains the updated singular values. */
/* Z (input) DOUBLE PRECISION array, dimension ( K ) */
/* The first K elements of this array contain the components */
/* of the deflation-adjusted updating row vector. */
/* VF (input/output) DOUBLE PRECISION array, dimension ( K ) */
/* On entry, VF contains information passed through DBEDE8. */
/* On exit, VF contains the first K components of the first */
/* components of all right singular vectors of the bidiagonal */
/* matrix. */
/* VL (input/output) DOUBLE PRECISION array, dimension ( K ) */
/* On entry, VL contains information passed through DBEDE8. */
/* On exit, VL contains the first K components of the last */
/* components of all right singular vectors of the bidiagonal */
/* matrix. */
/* DIFL (output) DOUBLE PRECISION array, dimension ( K ) */
/* On exit, DIFL(I) = D(I) - DSIGMA(I). */
/* DIFR (output) DOUBLE PRECISION array, */
/* dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and */
/* dimension ( K ) if ICOMPQ = 0. */
/* On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not */
/* defined and will not be referenced. */
/* If ICOMPQ = 1, DIFR(1:K,2) is an array containing the */
/* normalizing factors for the right singular vector matrix. */
//.........这里部分代码省略.........
doublereal dqrt14_(char *trans, integer *m, integer *n, integer *nrhs,
doublereal *a, integer *lda, doublereal *x, integer *ldx, doublereal *
work, integer *lwork)
{
/* System generated locals */
integer a_dim1, a_offset, x_dim1, x_offset, i__1, i__2, i__3;
doublereal ret_val, d__1, d__2, d__3;
/* Local variables */
integer i__, j;
doublereal err;
integer info;
doublereal anrm;
logical tpsd;
doublereal xnrm;
extern logical lsame_(char *, char *);
doublereal rwork[1];
extern /* Subroutine */ int dgelq2_(integer *, integer *, doublereal *,
integer *, doublereal *, doublereal *, integer *), dgeqr2_(
integer *, integer *, doublereal *, integer *, doublereal *,
doublereal *, integer *);
extern doublereal dlamch_(char *), dlange_(char *, integer *,
integer *, doublereal *, integer *, doublereal *);
extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
doublereal *, doublereal *, integer *, integer *, doublereal *,
integer *, integer *), dlacpy_(char *, integer *, integer
*, doublereal *, integer *, doublereal *, integer *),
xerbla_(char *, integer *);
integer ldwork;
/* -- LAPACK test routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DQRT14 checks whether X is in the row space of A or A'. It does so */
/* by scaling both X and A such that their norms are in the range */
/* [sqrt(eps), 1/sqrt(eps)], then computing a QR factorization of [A,X] */
/* (if TRANS = 'T') or an LQ factorization of [A',X]' (if TRANS = 'N'), */
/* and returning the norm of the trailing triangle, scaled by */
/* MAX(M,N,NRHS)*eps. */
/* Arguments */
/* ========= */
/* TRANS (input) CHARACTER*1 */
/* = 'N': No transpose, check for X in the row space of A */
/* = 'T': Transpose, check for X in the row space of A'. */
/* M (input) INTEGER */
/* The number of rows of the matrix A. */
/* N (input) INTEGER */
/* The number of columns of the matrix A. */
/* NRHS (input) INTEGER */
/* The number of right hand sides, i.e., the number of columns */
/* of X. */
/* A (input) DOUBLE PRECISION array, dimension (LDA,N) */
/* The M-by-N matrix A. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. */
/* X (input) DOUBLE PRECISION array, dimension (LDX,NRHS) */
/* If TRANS = 'N', the N-by-NRHS matrix X. */
/* IF TRANS = 'T', the M-by-NRHS matrix X. */
/* LDX (input) INTEGER */
/* The leading dimension of the array X. */
/* WORK (workspace) DOUBLE PRECISION array dimension (LWORK) */
/* LWORK (input) INTEGER */
/* length of workspace array required */
/* If TRANS = 'N', LWORK >= (M+NRHS)*(N+2); */
/* if TRANS = 'T', LWORK >= (N+NRHS)*(M+2). */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. Local Arrays .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
//.........这里部分代码省略.........
/* Subroutine */ int dlasd6_(integer *icompq, integer *nl, integer *nr,
integer *sqre, doublereal *d__, doublereal *vf, doublereal *vl,
doublereal *alpha, doublereal *beta, integer *idxq, integer *perm,
integer *givptr, integer *givcol, integer *ldgcol, doublereal *givnum,
integer *ldgnum, doublereal *poles, doublereal *difl, doublereal *
difr, doublereal *z__, integer *k, doublereal *c__, doublereal *s,
doublereal *work, integer *iwork, integer *info)
{
/* System generated locals */
integer givcol_dim1, givcol_offset, givnum_dim1, givnum_offset,
poles_dim1, poles_offset, i__1;
doublereal d__1, d__2;
/* Local variables */
static integer idxc, idxp, ivfw, ivlw, i__, m, n;
extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
doublereal *, integer *);
static integer n1, n2;
extern /* Subroutine */ int dlasd7_(integer *, integer *, integer *,
integer *, integer *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, integer *, integer *,
integer *, integer *, integer *, integer *, integer *, doublereal
*, integer *, doublereal *, doublereal *, integer *), dlasd8_(
integer *, integer *, doublereal *, doublereal *, doublereal *,
doublereal *, doublereal *, doublereal *, integer *, doublereal *,
doublereal *, integer *);
static integer iw;
extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
doublereal *, doublereal *, integer *, integer *, doublereal *,
integer *, integer *), dlamrg_(integer *, integer *,
doublereal *, integer *, integer *, integer *);
static integer isigma;
extern /* Subroutine */ int xerbla_(char *, integer *);
static doublereal orgnrm;
static integer idx;
#define poles_ref(a_1,a_2) poles[(a_2)*poles_dim1 + a_1]
/* -- LAPACK auxiliary routine (instrumented to count ops, version 3.0) --
Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
Courant Institute, Argonne National Lab, and Rice University
June 30, 1999
Purpose
=======
DLASD6 computes the SVD of an updated upper bidiagonal matrix B
obtained by merging two smaller ones by appending a row. This
routine is used only for the problem which requires all singular
values and optionally singular vector matrices in factored form.
B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE.
A related subroutine, DLASD1, handles the case in which all singular
values and singular vectors of the bidiagonal matrix are desired.
DLASD6 computes the SVD as follows:
( D1(in) 0 0 0 )
B = U(in) * ( Z1' a Z2' b ) * VT(in)
( 0 0 D2(in) 0 )
= U(out) * ( D(out) 0) * VT(out)
where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M
with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
elsewhere; and the entry b is empty if SQRE = 0.
The singular values of B can be computed using D1, D2, the first
components of all the right singular vectors of the lower block, and
the last components of all the right singular vectors of the upper
block. These components are stored and updated in VF and VL,
respectively, in DLASD6. Hence U and VT are not explicitly
referenced.
The singular values are stored in D. The algorithm consists of two
stages:
The first stage consists of deflating the size of the problem
when there are multiple singular values or if there is a zero
in the Z vector. For each such occurence the dimension of the
secular equation problem is reduced by one. This stage is
performed by the routine DLASD7.
The second stage consists of calculating the updated
singular values. This is done by finding the roots of the
secular equation via the routine DLASD4 (as called by DLASD8).
This routine also updates VF and VL and computes the distances
between the updated singular values and the old singular
values.
DLASD6 is called from DLASDA.
Arguments
=========
ICOMPQ (input) INTEGER
Specifies whether singular vectors are to be computed in
//.........这里部分代码省略.........
/* Subroutine */ int dlasq1_(integer *n, doublereal *d__, doublereal *e,
doublereal *work, integer *info)
{
/* System generated locals */
integer i__1, i__2;
doublereal d__1, d__2, d__3;
/* Builtin functions */
double sqrt(doublereal);
/* Local variables */
integer i__;
doublereal eps;
extern /* Subroutine */ int dlas2_(doublereal *, doublereal *, doublereal
*, doublereal *, doublereal *);
doublereal scale;
integer iinfo;
doublereal sigmn;
extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
doublereal *, integer *);
doublereal sigmx;
extern /* Subroutine */ int dlasq2_(integer *, doublereal *, integer *);
extern doublereal dlamch_(char *);
extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
doublereal *, doublereal *, integer *, integer *, doublereal *,
integer *, integer *);
doublereal safmin;
extern /* Subroutine */ int xerbla_(char *, integer *), dlasrt_(
char *, integer *, doublereal *, integer *);
/* -- LAPACK routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DLASQ1 computes the singular values of a real N-by-N bidiagonal */
/* matrix with diagonal D and off-diagonal E. The singular values */
/* are computed to high relative accuracy, in the absence of */
/* denormalization, underflow and overflow. The algorithm was first */
/* presented in */
/* "Accurate singular values and differential qd algorithms" by K. V. */
/* Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230, */
/* 1994, */
/* and the present implementation is described in "An implementation of */
/* the dqds Algorithm (Positive Case)", LAPACK Working Note. */
/* Arguments */
/* ========= */
/* N (input) INTEGER */
/* The number of rows and columns in the matrix. N >= 0. */
/* D (input/output) DOUBLE PRECISION array, dimension (N) */
/* On entry, D contains the diagonal elements of the */
/* bidiagonal matrix whose SVD is desired. On normal exit, */
/* D contains the singular values in decreasing order. */
/* E (input/output) DOUBLE PRECISION array, dimension (N) */
/* On entry, elements E(1:N-1) contain the off-diagonal elements */
/* of the bidiagonal matrix whose SVD is desired. */
/* On exit, E is overwritten. */
/* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* > 0: the algorithm failed */
/* = 1, a split was marked by a positive value in E */
/* = 2, current block of Z not diagonalized after 30*N */
/* iterations (in inner while loop) */
/* = 3, termination criterion of outer while loop not met */
/* (program created more than N unreduced blocks) */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Parameter adjustments */
--work;
//.........这里部分代码省略.........
客服电话
APP下载
官方微信
请发表评论