dlasd3.c

著名的LAPACK矩阵计算软件包, 是比较新的版本, 一般用到矩阵分解的朋友也许会用到

#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 {
    doublereal ops, itcnt;
} latime_;

#define latime_1 latime_

/* Table of constant values */

static integer c__1 = 1;
static integer c__0 = 0;
static doublereal c_b13 = 1.;
static doublereal c_b27 = 0.;

/* Subroutine */ int dlasd3_(integer *nl, integer *nr, integer *sqre, integer 
	*k, doublereal *d__, doublereal *q, integer *ldq, doublereal *dsigma, 
	doublereal *u, integer *ldu, doublereal *u2, integer *ldu2, 
	doublereal *vt, integer *ldvt, doublereal *vt2, integer *ldvt2, 
	integer *idxc, integer *ctot, doublereal *z__, integer *info)
{
    /* System generated locals */
    integer q_dim1, q_offset, u_dim1, u_offset, u2_dim1, u2_offset, vt_dim1, 
	    vt_offset, vt2_dim1, vt2_offset, i__1, i__2;
    doublereal d__1, d__2;

    /* Builtin functions */
    double sqrt(doublereal), d_sign(doublereal *, doublereal *);

    /* Local variables */
    static doublereal temp;
    extern doublereal dnrm2_(integer *, doublereal *, integer *);
    static integer i__, j, m, n;
    extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 
	    integer *, doublereal *, doublereal *, integer *, doublereal *, 
	    integer *, doublereal *, doublereal *, integer *);
    static integer ctemp;
    extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
	    doublereal *, integer *);
    static integer ktemp;
    extern doublereal dlamc3_(doublereal *, doublereal *);
    extern /* Subroutine */ int dlasd4_(integer *, integer *, doublereal *, 
	    doublereal *, doublereal *, doublereal *, doublereal *, 
	    doublereal *, integer *);
    extern doublereal dopbl3_(char *, integer *, integer *, integer *)
	    ;
    static integer jc;
    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 *);
    static doublereal rho;
    static integer nlp1, nlp2, nrp1;


#define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1]
#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   
    =======   

    DLASD3 finds all the square roots of the roots of the secular   
    equation, as defined by the values in D and Z.  It makes the   
    appropriate calls to DLASD4 and then updates the singular   
    vectors by matrix multiplication.   

    This code makes very mild assumptions about floating point   
    arithmetic. It will work on machines with a guard digit in   
    add/subtract, or on those binary machines without guard digits   
    which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.   
    It could conceivably fail on hexadecimal or decimal machines   
    without guard digits, but we know of none.   

    DLASD3 is called from DLASD1.   

    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      (input) INTEGER   
           The size of the secular equation, 1 =< K = < N.   

    D      (output) DOUBLE PRECISION array, dimension(K)   
           On exit the square roots of the roots of the secular equation,   
           in ascending order.   

    Q      (workspace) DOUBLE PRECISION array,   
                       dimension at least (LDQ,K).   

    LDQ    (input) INTEGER   
           The leading dimension of the array Q.  LDQ >= K.   

    DSIGMA (input) DOUBLE PRECISION array, dimension(K)   
           The first K elements of this array contain the old roots   
           of the deflated updating problem.  These are the poles   
           of the secular equation.   

    U      (input) DOUBLE PRECISION array, dimension (LDU, N)   
           The last N - K columns of this matrix contain the deflated   
           left singular vectors.   

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

    U2     (input) DOUBLE PRECISION array, dimension (LDU2, N)   
           The first K columns of this matrix contain the non-deflated   
           left singular vectors for the split problem.   

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

    VT     (input) DOUBLE PRECISION array, dimension (LDVT, M)   
           The last M - K columns of VT' contain the deflated   
           right singular vectors.   

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

    VT2    (input) DOUBLE PRECISION array, dimension (LDVT2, N)   
           The first K columns of VT2' contain the non-deflated   
           right singular vectors for the split problem.   

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

    IDXC   (input) INTEGER array, dimension ( N )   
           The permutation used to arrange the columns of U (and rows of   
           VT) into three groups:  the first group contains non-zero   
           entries only at and above (or before) NL +1; the second   
           contains non-zero entries only at and below (or after) NL+2;   
           and the third is dense. The first column of U and the row of   
           VT are treated separately, however.   

           The rows of the singular vectors found by DLASD4   
           must be likewise permuted before the matrix multiplies can   
           take place.   

    CTOT   (input) INTEGER array, dimension ( 4 )   
           A count of the total number of the various types of columns   
           in U (or rows in VT), as described in IDXC. The fourth column   
           type is any column which has been deflated.   

    Z      (input) DOUBLE PRECISION array, dimension (K)   
           The first K elements of this array contain the components   
           of the deflation-adjusted updating row vector.   

    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__;
    q_dim1 = *ldq;
    q_offset = 1 + q_dim1 * 1;
    q -= q_offset;
    --dsigma;
    u_dim1 = *ldu;
    u_offset = 1 + u_dim1 * 1;
    u -= u_offset;
    u2_dim1 = *ldu2;
    u2_offset = 1 + u2_dim1 * 1;
    u2 -= u2_offset;
    vt_dim1 = *ldvt;
    vt_offset = 1 + vt_dim1 * 1;
    vt -= vt_offset;
    vt2_dim1 = *ldvt2;
    vt2_offset = 1 + vt2_dim1 * 1;
    vt2 -= vt2_offset;
    --idxc;
    --ctot;
    --z__;

    /* 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;
    nlp1 = *nl + 1;
    nlp2 = *nl + 2;

    if (*k < 1 || *k > n) {
	*info = -4;
    } else if (*ldq < *k) {
	*info = -7;
    } else if (*ldu < n) {