slasd2.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__1 = 1;
static real c_b30 = 0.f;

/* Subroutine */ int slasd2_(integer *nl, integer *nr, integer *sqre, integer 
	*k, real *d__, real *z__, real *alpha, real *beta, real *u, integer *
	ldu, real *vt, integer *ldvt, real *dsigma, real *u2, integer *ldu2, 
	real *vt2, integer *ldvt2, integer *idxp, integer *idx, integer *idxc,
	 integer *idxq, integer *coltyp, integer *info)
{
    /* System generated locals */
    integer u_dim1, u_offset, u2_dim1, u2_offset, vt_dim1, vt_offset, 
	    vt2_dim1, vt2_offset, i__1;
    real r__1, r__2;

    /* Local variables */
    static integer idxi, idxj, ctot[4];
    extern /* Subroutine */ int srot_(integer *, real *, integer *, real *, 
	    integer *, real *, real *);
    static real c__;
    static integer i__, j, m, n;
    static real s;
    static integer idxjp, jprev, k2;
    extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, 
	    integer *);
    static real z1;
    extern doublereal slapy2_(real *, real *);
    static integer ct, jp;
    extern doublereal slamch_(char *);
    extern /* Subroutine */ int xerbla_(char *, integer *), slamrg_(
	    integer *, integer *, real *, integer *, integer *, integer *);
    static real hlftol;
    extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *, 
	    integer *, real *, integer *), slaset_(char *, integer *, 
	    integer *, real *, real *, real *, integer *);
    static real eps, tau, tol;
    static integer psm[4], nlp1, nlp2;


#define u_ref(a_1,a_2) u[(a_2)*u_dim1 + a_1]
#define u2_ref(a_1,a_2) u2[(a_2)*u2_dim1 + a_1]
#define vt_ref(a_1,a_2) vt[(a_2)*vt_dim1 + a_1]
#define vt2_ref(a_1,a_2) vt2[(a_2)*vt2_dim1 + a_1]


/*  -- LAPACK auxiliary routine (instrumented to count ops, version 3.0) --   
       Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,   
       Courant Institute, NAG Ltd., and Rice University   
       October 31, 1999   


    Purpose   
    =======   

    SLASD2 merges the two sets of singular values together into a single   
    sorted set.  Then it tries to deflate the size of the problem.   
    There are two ways in which deflation can occur:  when two or more   
    singular values are close together or if there is a tiny entry in the   
    Z vector.  For each such occurrence the order of the related secular   
    equation problem is reduced by one.   

    SLASD2 is called from SLASD1.   

    Arguments   
    =========   

    NL     (input) INTEGER   
           The row dimension of the upper block.  NL >= 1.   

    NR     (input) INTEGER   
           The row dimension of the lower block.  NR >= 1.   

    SQRE   (input) INTEGER   
           = 0: the lower block is an NR-by-NR square matrix.   
           = 1: the lower block is an NR-by-(NR+1) rectangular matrix.   

           The bidiagonal matrix has N = NL + NR + 1 rows and   
           M = N + SQRE >= N columns.   

    K      (output) INTEGER   
           Contains the dimension of the non-deflated matrix,   
           This is the order of the related secular equation. 1 <= K <=N.   

    D      (input/output) REAL array, dimension(N)   
           On entry D contains the singular values of the two submatrices   
           to be combined.  On exit D contains the trailing (N-K) updated   
           singular values (those which were deflated) sorted into   
           increasing order.   

    ALPHA  (input) REAL   
           Contains the diagonal element associated with the added row.   

    BETA   (input) REAL   
           Contains the off-diagonal element associated with the added   
           row.   

    U      (input/output) REAL array, dimension(LDU,N)   
           On entry U contains the left singular vectors of two   
           submatrices in the two square blocks with corners at (1,1),   
           (NL, NL), and (NL+2, NL+2), (N,N).   
           On exit U contains the trailing (N-K) updated left singular   
           vectors (those which were deflated) in its last N-K columns.   

    LDU    (input) INTEGER   
           The leading dimension of the array U.  LDU >= N.   

    Z      (output) REAL array, dimension(N)   
           On exit Z contains the updating row vector in the secular   
           equation.   

    DSIGMA (output) REAL array, dimension (N)   
           Contains a copy of the diagonal elements (K-1 singular values   
           and one zero) in the secular equation.   

    U2     (output) REAL array, dimension(LDU2,N)   
           Contains a copy of the first K-1 left singular vectors which   
           will be used by SLASD3 in a matrix multiply (SGEMM) to solve   
           for the new left singular vectors. U2 is arranged into four   
           blocks. The first block contains a column with 1 at NL+1 and   
           zero everywhere else; the second block contains non-zero   
           entries only at and above NL; the third contains non-zero   
           entries only below NL+1; and the fourth is dense.   

    LDU2   (input) INTEGER   
           The leading dimension of the array U2.  LDU2 >= N.   

    VT     (input/output) REAL array, dimension(LDVT,M)   
           On entry VT' contains the right singular vectors of two   
           submatrices in the two square blocks with corners at (1,1),   
           (NL+1, NL+1), and (NL+2, NL+2), (M,M).   
           On exit VT' contains the trailing (N-K) updated right singular   
           vectors (those which were deflated) in its last N-K columns.   
           In case SQRE =1, the last row of VT spans the right null   
           space.   

    LDVT   (input) INTEGER   
           The leading dimension of the array VT.  LDVT >= M.   

    VT2    (output) REAL array, dimension(LDVT2,N)   
           VT2' contains a copy of the first K right singular vectors   
           which will be used by SLASD3 in a matrix multiply (SGEMM) to   
           solve for the new right singular vectors. VT2 is arranged into   
           three blocks. The first block contains a row that corresponds   
           to the special 0 diagonal element in SIGMA; the second block   
           contains non-zeros only at and before NL +1; the third block   
           contains non-zeros only at and after  NL +2.   

    LDVT2  (input) INTEGER   
           The leading dimension of the array VT2.  LDVT2 >= M.   

    IDXP   (workspace) INTEGER array, dimension(N)   
           This will contain the permutation used to place deflated   
           values of D at the end of the array. On output IDXP(2:K)   
           points to the nondeflated D-values and IDXP(K+1:N)   
           points to the deflated singular values.   

    IDX    (workspace) INTEGER array, dimension(N)   
           This will contain the permutation used to sort the contents of   
           D into ascending order.   

    IDXC   (output) INTEGER array, dimension(N)   
           This will contain the permutation used to arrange the columns   
           of the deflated U matrix into three groups:  the first group   
           contains non-zero entries only at and above NL, the second   
           contains non-zero entries only below NL+2, and the third is   
           dense.   

    COLTYP (workspace/output) INTEGER array, dimension(N)   
           As workspace, this will contain a label which will indicate   
           which of the following types a column in the U2 matrix or a   
           row in the VT2 matrix is:   
           1 : non-zero in the upper half only   
           2 : non-zero in the lower half only   
           3 : dense   
           4 : deflated   

           On exit, it is an array of dimension 4, with COLTYP(I) being   
           the dimension of the I-th type columns.   

    IDXQ   (input) INTEGER array, dimension(N)   
           This contains the permutation which separately sorts the two   
           sub-problems in D into ascending order.  Note that entries in   
           the first hlaf of this permutation must first be moved one   
           position backward; and entries in the second half   
           must first have NL+1 added to their values.   

    INFO   (output) INTEGER   
            = 0:  successful exit.   
            < 0:  if INFO = -i, the i-th argument had an illegal value.   

    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__;
    --z__;
    u_dim1 = *ldu;
    u_offset = 1 + u_dim1 * 1;
    u -= u_offset;
    vt_dim1 = *ldvt;
    vt_offset = 1 + vt_dim1 * 1;
    vt -= vt_offset;
    --dsigma;
    u2_dim1 = *ldu2;
    u2_offset = 1 + u2_dim1 * 1;
    u2 -= u2_offset;
    vt2_dim1 = *ldvt2;
    vt2_offset = 1 + vt2_dim1 * 1;
    vt2 -= vt2_offset;
    --idxp;
    --idx;
    --idxc;
    --idxq;
    --coltyp;

    /* Function Body */
    *info = 0;

    if (*nl < 1) {
	*info = -1;
    } else if (*nr < 1) {
	*info = -2;
    } else if (*sqre != 1 && *sqre != 0) {
	*info = -3;
    }

    n = *nl + *nr + 1;
    m = n + *sqre;

    if (*ldu < n) {
	*info = -10;
    } else if (*ldvt < m) {
	*info = -12;
    } else if (*ldu2 < n) {
	*info = -15;
    } else if (*ldvt2 < m) {
	*info = -17;
    }
    if (*info != 0) {
	i__1 = -(*info);
	xerbla_("SLASD2", &i__1);
	return 0;
    }

    nlp1 = *nl + 1;
    nlp2 = *nl + 2;

/*     Generate the first part of the vector Z; and move the singular   
       values in the first part of D one position backward. */

    latime_1.ops += (real) (*nl + 1);
    z1 = *alpha * vt_ref(nlp1, nlp1);
    z__[1] = z1;
    for (i__ = *nl; i__ >= 1; --i__) {
	z__[i__ + 1] = *alpha * vt_ref(i__, nlp1);
	d__[i__ + 1] = d__[i__];
	idxq[i__ + 1] = idxq[i__] + 1;
/* L10: */
    }

/*     Generate the second part of the vector Z. */

    latime_1.ops += (real) (m - nlp2 + 1);
    i__1 = m;
    for (i__ = nlp2; i__ <= i__1; ++i__) {
	z__[i__] = *beta * vt_ref(i__, nlp2);
/* L20: */
    }

/*     Initialize some reference arrays. */

    i__1 = nlp1;
    for (i__ = 2; i__ <= i__1; ++i__) {
	coltyp[i__] = 1;
/* L30: */
    }
    i__1 = n;
    for (i__ = nlp2; i__ <= i__1; ++i__) {
	coltyp[i__] = 2;
/* L40: */
    }

/*     Sort the singular values into increasing order */