slasd0.c

提供矩阵类的函数库

#include "blaswrap.h"
/*  -- translated by f2c (version 19990503).
   You must link the resulting object file with the libraries:
	-lf2c -lm   (in that order)
*/

#include "f2c.h"

/* Common Block Declarations */

struct {
    real ops, itcnt;
} latime_;

#define latime_1 latime_

/* Table of constant values */

static integer c__0 = 0;
static integer c__2 = 2;

/* Subroutine */ int slasd0_(integer *n, integer *sqre, real *d__, real *e, 
	real *u, integer *ldu, real *vt, integer *ldvt, integer *smlsiz, 
	integer *iwork, real *work, integer *info)
{
    /* System generated locals */
    integer u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;

    /* Builtin functions */
    integer pow_ii(integer *, integer *);

    /* Local variables */
    static real beta;
    static integer idxq, nlvl, i__, j, m;
    static real alpha;
    static integer inode, ndiml, idxqc, ndimr, itemp, sqrei, i1;
    extern /* Subroutine */ int slasd1_(integer *, integer *, integer *, real 
	    *, real *, real *, real *, integer *, real *, integer *, integer *
	    , integer *, real *, integer *);
    static integer ic, lf, nd, ll, nl, nr;
    extern /* Subroutine */ int xerbla_(char *, integer *), slasdq_(
	    char *, integer *, integer *, integer *, integer *, integer *, 
	    real *, real *, real *, integer *, real *, integer *, real *, 
	    integer *, real *, integer *), slasdt_(integer *, integer 
	    *, integer *, integer *, integer *, integer *, integer *);
    static integer im1, ncc, nlf, nrf, iwk, lvl, ndb1, nlp1, nrp1;


#define u_ref(a_1,a_2) u[(a_2)*u_dim1 + a_1]
#define vt_ref(a_1,a_2) vt[(a_2)*vt_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   
    =======   

    Using a divide and conquer approach, SLASD0 computes the singular   
    value decomposition (SVD) of a real upper bidiagonal N-by-M   
    matrix B with diagonal D and offdiagonal E, where M = N + SQRE.   
    The algorithm computes orthogonal matrices U and VT such that   
    B = U * S * VT. The singular values S are overwritten on D.   

    A related subroutine, SLASDA, computes only the singular values,   
    and optionally, the singular vectors in compact form.   

    Arguments   
    =========   

    N      (input) INTEGER   
           On entry, the row dimension of the upper bidiagonal matrix.   
           This is also the dimension of the main diagonal array D.   

    SQRE   (input) INTEGER   
           Specifies the column dimension of the bidiagonal matrix.   
           = 0: The bidiagonal matrix has column dimension M = N;   
           = 1: The bidiagonal matrix has column dimension M = N+1;   

    D      (input/output) REAL array, dimension (N)   
           On entry D contains the main diagonal of the bidiagonal   
           matrix.   
           On exit D, if INFO = 0, contains its singular values.   

    E      (input) REAL array, dimension (M-1)   
           Contains the subdiagonal entries of the bidiagonal matrix.   
           On exit, E has been destroyed.   

    U      (output) REAL array, dimension at least (LDQ, N)   
           On exit, U contains the left singular vectors.   

    LDU    (input) INTEGER   
           On entry, leading dimension of U.   

    VT     (output) REAL array, dimension at least (LDVT, M)   
           On exit, VT' contains the right singular vectors.   

    LDVT   (input) INTEGER   
           On entry, leading dimension of VT.   

    SMLSIZ (input) INTEGER   
           On entry, maximum size of the subproblems at the   
           bottom of the computation tree.   

    IWORK  INTEGER work array.   
           Dimension must be at least (8 * N)   

    WORK   REAL work array.   
           Dimension must be at least (3 * M**2 + 2 * M)   

    INFO   (output) INTEGER   
            = 0:  successful exit.   
            < 0:  if INFO = -i, the i-th argument had an illegal value.   
            > 0:  if INFO = 1, an singular value did not converge   

    Further Details   
    ===============   

    Based on contributions by   
       Ming Gu and Huan Ren, Computer Science Division, University of   
       California at Berkeley, USA   

    =====================================================================   


       Test the input parameters.   

       Parameter adjustments */
    --d__;
    --e;
    u_dim1 = *ldu;
    u_offset = 1 + u_dim1 * 1;
    u -= u_offset;
    vt_dim1 = *ldvt;
    vt_offset = 1 + vt_dim1 * 1;
    vt -= vt_offset;
    --iwork;
    --work;

    /* Function Body */
    *info = 0;

    if (*n < 0) {
	*info = -1;
    } else if (*sqre < 0 || *sqre > 1) {
	*info = -2;
    }

    m = *n + *sqre;

    if (*ldu < *n) {
	*info = -6;
    } else if (*ldvt < m) {
	*info = -8;
    } else if (*smlsiz < 3) {
	*info = -9;
    }
    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("SLASD0", &i__1);
	return 0;
    }

/*     If the input matrix is too small, call SLASDQ to find the SVD. */

    if (*n <= *smlsiz) {
	slasdq_("U", sqre, n, &m, n, &c__0, &d__[1], &e[1], &vt[vt_offset], 
		ldvt, &u[u_offset], ldu, &u[u_offset], ldu, &work[1], info);
	return 0;
    }

/*     Set up the computation tree. */

    inode = 1;
    ndiml = inode + *n;
    ndimr = ndiml + *n;
    idxq = ndimr + *n;
    iwk = idxq + *n;
    slasdt_(n, &nlvl, &nd, &iwork[inode], &iwork[ndiml], &iwork[ndimr], 
	    smlsiz);

/*     For the nodes on bottom level of the tree, solve   
       their subproblems by SLASDQ. */

    ndb1 = (nd + 1) / 2;
    ncc = 0;
    i__1 = nd;
    for (i__ = ndb1; i__ <= i__1; ++i__) {

/*     IC : center row of each node   
       NL : number of rows of left  subproblem   
       NR : number of rows of right subproblem   
       NLF: starting row of the left   subproblem   
       NRF: starting row of the right  subproblem */

	i1 = i__ - 1;
	ic = iwork[inode + i1];
	nl = iwork[ndiml + i1];
	nlp1 = nl + 1;
	nr = iwork[ndimr + i1];
	nrp1 = nr + 1;
	nlf = ic - nl;
	nrf = ic + 1;
	sqrei = 1;
	slasdq_("U", &sqrei, &nl, &nlp1, &nl, &ncc, &d__[nlf], &e[nlf], &
		vt_ref(nlf, nlf), ldvt, &u_ref(nlf, nlf), ldu, &u_ref(nlf, 
		nlf), ldu, &work[1], info);
	if (*info != 0) {
	    return 0;
	}
	itemp = idxq + nlf - 2;
	i__2 = nl;
	for (j = 1; j <= i__2; ++j) {
	    iwork[itemp + j] = j;
/* L10: */
	}
	if (i__ == nd) {
	    sqrei = *sqre;
	} else {
	    sqrei = 1;
	}
	nrp1 = nr + sqrei;
	slasdq_("U", &sqrei, &nr, &nrp1, &nr, &ncc, &d__[nrf], &e[nrf], &
		vt_ref(nrf, nrf), ldvt, &u_ref(nrf, nrf), ldu, &u_ref(nrf, 
		nrf), ldu, &work[1], info);
	if (*info != 0) {
	    return 0;
	}
	itemp = idxq + ic;
	i__2 = nr;
	for (j = 1; j <= i__2; ++j) {
	    iwork[itemp + j - 1] = j;
/* L20: */
	}
/* L30: */
    }

/*     Now conquer each subproblem bottom-up. */

    for (lvl = nlvl; lvl >= 1; --lvl) {

/*        Find the first node LF and last node LL on the   
          current level LVL. */

	if (lvl == 1) {
	    lf = 1;
	    ll = 1;
	} else {
	    i__1 = lvl - 1;
	    lf = pow_ii(&c__2, &i__1);
	    ll = (lf << 1) - 1;
	}
	i__1 = ll;
	for (i__ = lf; i__ <= i__1; ++i__) {
	    im1 = i__ - 1;
	    ic = iwork[inode + im1];
	    nl = iwork[ndiml + im1];
	    nr = iwork[ndimr + im1];
	    nlf = ic - nl;
	    if (*sqre == 0 && i__ == ll) {
		sqrei = *sqre;
	    } else {
		sqrei = 1;
	    }
	    idxqc = idxq + nlf - 1;
	    alpha = d__[ic];
	    beta = e[ic];
	    slasd1_(&nl, &nr, &sqrei, &d__[nlf], &alpha, &beta, &u_ref(nlf, 
		    nlf), ldu, &vt_ref(nlf, nlf), ldvt, &iwork[idxqc], &iwork[
		    iwk], &work[1], info);
	    if (*info != 0) {
		return 0;
	    }
/* L40: */
	}
/* L50: */
    }

    return 0;

/*     End of SLASD0 */

} /* slasd0_ */

#undef vt_ref
#undef u_ref