292 lines
		
	
	
		
			8.0 KiB
		
	
	
	
		
			C
		
	
	
	
	
	
			
		
		
	
	
			292 lines
		
	
	
		
			8.0 KiB
		
	
	
	
		
			C
		
	
	
	
	
	
| /* dlasd0.f -- translated by f2c (version 20061008).
 | |
|    You must link the resulting object file with libf2c:
 | |
| 	on Microsoft Windows system, link with libf2c.lib;
 | |
| 	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
 | |
| 	or, if you install libf2c.a in a standard place, with -lf2c -lm
 | |
| 	-- in that order, at the end of the command line, as in
 | |
| 		cc *.o -lf2c -lm
 | |
| 	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
 | |
| 
 | |
| 		http://www.netlib.org/f2c/libf2c.zip
 | |
| */
 | |
| 
 | |
| #include "clapack.h"
 | |
| 
 | |
| 
 | |
| /* Table of constant values */
 | |
| 
 | |
| static integer c__0 = 0;
 | |
| static integer c__2 = 2;
 | |
| 
 | |
| /* Subroutine */ int dlasd0_(integer *n, integer *sqre, doublereal *d__, 
 | |
| 	doublereal *e, doublereal *u, integer *ldu, doublereal *vt, integer *
 | |
| 	ldvt, integer *smlsiz, integer *iwork, doublereal *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 */
 | |
|     integer i__, j, m, i1, ic, lf, nd, ll, nl, nr, im1, ncc, nlf, nrf, iwk, 
 | |
| 	    lvl, ndb1, nlp1, nrp1;
 | |
|     doublereal beta;
 | |
|     integer idxq, nlvl;
 | |
|     doublereal alpha;
 | |
|     integer inode, ndiml, idxqc, ndimr, itemp, sqrei;
 | |
|     extern /* Subroutine */ int dlasd1_(integer *, integer *, integer *, 
 | |
| 	    doublereal *, doublereal *, doublereal *, doublereal *, integer *, 
 | |
| 	     doublereal *, integer *, integer *, integer *, doublereal *, 
 | |
| 	    integer *), dlasdq_(char *, integer *, integer *, integer *, 
 | |
| 	    integer *, integer *, doublereal *, doublereal *, doublereal *, 
 | |
| 	    integer *, doublereal *, integer *, doublereal *, integer *, 
 | |
| 	    doublereal *, integer *), dlasdt_(integer *, integer *, 
 | |
| 	    integer *, integer *, integer *, integer *, integer *), xerbla_(
 | |
| 	    char *, integer *);
 | |
| 
 | |
| 
 | |
| /*  -- LAPACK auxiliary routine (version 3.2) -- */
 | |
| /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
 | |
| /*     November 2006 */
 | |
| 
 | |
| /*     .. Scalar Arguments .. */
 | |
| /*     .. */
 | |
| /*     .. Array Arguments .. */
 | |
| /*     .. */
 | |
| 
 | |
| /*  Purpose */
 | |
| /*  ======= */
 | |
| 
 | |
| /*  Using a divide and conquer approach, DLASD0 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, DLASDA, 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) DOUBLE PRECISION 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) DOUBLE PRECISION array, dimension (M-1) */
 | |
| /*         Contains the subdiagonal entries of the bidiagonal matrix. */
 | |
| /*         On exit, E has been destroyed. */
 | |
| 
 | |
| /*  U      (output) DOUBLE PRECISION 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) DOUBLE PRECISION 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  (workspace) INTEGER work array. */
 | |
| /*         Dimension must be at least (8 * N) */
 | |
| 
 | |
| /*  WORK   (workspace) DOUBLE PRECISION 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 */
 | |
| 
 | |
| /*  ===================================================================== */
 | |
| 
 | |
| /*     .. Local Scalars .. */
 | |
| /*     .. */
 | |
| /*     .. External Subroutines .. */
 | |
| /*     .. */
 | |
| /*     .. Executable Statements .. */
 | |
| 
 | |
| /*     Test the input parameters. */
 | |
| 
 | |
|     /* Parameter adjustments */
 | |
|     --d__;
 | |
|     --e;
 | |
|     u_dim1 = *ldu;
 | |
|     u_offset = 1 + u_dim1;
 | |
|     u -= u_offset;
 | |
|     vt_dim1 = *ldvt;
 | |
|     vt_offset = 1 + vt_dim1;
 | |
|     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_("DLASD0", &i__1);
 | |
| 	return 0;
 | |
|     }
 | |
| 
 | |
| /*     If the input matrix is too small, call DLASDQ to find the SVD. */
 | |
| 
 | |
|     if (*n <= *smlsiz) {
 | |
| 	dlasdq_("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;
 | |
|     dlasdt_(n, &nlvl, &nd, &iwork[inode], &iwork[ndiml], &iwork[ndimr], 
 | |
| 	    smlsiz);
 | |
| 
 | |
| /*     For the nodes on bottom level of the tree, solve */
 | |
| /*     their subproblems by DLASDQ. */
 | |
| 
 | |
|     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;
 | |
| 	dlasdq_("U", &sqrei, &nl, &nlp1, &nl, &ncc, &d__[nlf], &e[nlf], &vt[
 | |
| 		nlf + nlf * vt_dim1], ldvt, &u[nlf + nlf * u_dim1], ldu, &u[
 | |
| 		nlf + nlf * u_dim1], 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;
 | |
| 	dlasdq_("U", &sqrei, &nr, &nrp1, &nr, &ncc, &d__[nrf], &e[nrf], &vt[
 | |
| 		nrf + nrf * vt_dim1], ldvt, &u[nrf + nrf * u_dim1], ldu, &u[
 | |
| 		nrf + nrf * u_dim1], 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];
 | |
| 	    dlasd1_(&nl, &nr, &sqrei, &d__[nlf], &alpha, &beta, &u[nlf + nlf *
 | |
| 		     u_dim1], ldu, &vt[nlf + nlf * vt_dim1], ldvt, &iwork[
 | |
| 		    idxqc], &iwork[iwk], &work[1], info);
 | |
| 	    if (*info != 0) {
 | |
| 		return 0;
 | |
| 	    }
 | |
| /* L40: */
 | |
| 	}
 | |
| /* L50: */
 | |
|     }
 | |
| 
 | |
|     return 0;
 | |
| 
 | |
| /*     End of DLASD0 */
 | |
| 
 | |
| } /* dlasd0_ */
 | 
