/* Subroutine */ int PASTEF77(s,rotg)(real *sa, real *sb, real *c__, real *s)
{
    /* System generated locals */
    real r__1, r__2;

    /* Builtin functions */
    double sqrt(doublereal), r_sign(real *, real *);

    /* Local variables */
    real r__, scale, z__, roe;


/*     construct givens plane rotation. */
/*     jack dongarra, linpack, 3/11/78. */


    roe = *sb;
    if (abs(*sa) > abs(*sb)) {
	roe = *sa;
    }
    scale = abs(*sa) + abs(*sb);
    if (scale != 0.f) {
	goto L10;
    }
    *c__ = 1.f;
    *s = 0.f;
    r__ = 0.f;
    z__ = 0.f;
    goto L20;
L10:
/* Computing 2nd power */
    r__1 = *sa / scale;
/* Computing 2nd power */
    r__2 = *sb / scale;
    r__ = scale * sqrt(r__1 * r__1 + r__2 * r__2);
    r__ = r_sign(&sc_b4, &roe) * r__;
    *c__ = *sa / r__;
    *s = *sb / r__;
    z__ = 1.f;
    if (abs(*sa) > abs(*sb)) {
	z__ = *s;
    }
    if (abs(*sb) >= abs(*sa) && *c__ != 0.f) {
	z__ = 1.f / *c__;
    }
L20:
    *sa = r__;
    *sb = z__;
    return 0;
} /* srotg_ */

/*< 	S >*/
/* Subroutine */ int onset_0_(int n__, real *pebuf, integer *osbuf, integer *
	osptr, integer *oslen, integer *sbufl, integer *sbufh, integer *
	lframe)
{
    /* Initialized data */

    static real n = 0.f;
    static real d__ = 1.f;
    static real l2buf[16] = { 0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f };
    static real l2sum1 = 0.f;
    static integer l2ptr1 = 1;
    static integer l2ptr2 = 9;
    static logical hyst = FALSE_;

    /* System generated locals */
    integer i__1;

    /* Builtin functions */
    double r_sign(real *, real *);

    /* Local variables */
    integer i__;
    static integer lasti;
    real l2sum2;
    static real fpc;

/*< 	INCLUDE 'config.fh' >*/
/*< 	INTEGER OSLEN, SBUFL, SBUFH, LFRAME >*/
/*       Arguments */
/* $Log: onset.c,v $
/* Revision 1.2  2001-01-25 23:45:49  jpoehlmann
/* Version 1.7c. Identical with files on the ftp Server ftp.franken.de.
/* (+ 1 patch in cli.c, wich is on the server too)
/* Not compiled now
/* */
/* Revision 1.3  1996/03/29  22:03:47  jaf */
/* Removed definitions for any constants that were no longer used. */

/* Revision 1.2  1996/03/26  19:34:33  jaf */
/* Added comments indicating which constants are not needed in an */
/* application that uses the LPC-10 coder. */

/* Revision 1.1  1996/02/07  14:43:51  jaf */
/* Initial revision */

/*   LPC Configuration parameters: */
/* Frame size, Prediction order, Pitch period */
/*< 	parameter (MAXFRM = 180, MAXORD = 10, MAXPIT = 156) >*/
/*< 	REAL PEBUF(SBUFL:SBUFH) >*/
/*< 	INTEGER OSBUF(OSLEN), OSPTR >*/
/*       Parameters/constants */
/*   Parameters for onset detection algorithm: */
/*    L2		Threshold for filtered slope of FPC (function of L2WID!) */
/*    L2LAG	Lag due to both filters which compute filtered slope of FPC */
/*    L2WID	Width of the filter which computes the slope of FPC */
/*    OSHYST	The number of samples of slope(FPC) which must be below */
/* 	        the threshold before a new onset may be declared. */
/*< 	INTEGER L2LAG, L2WID, OSHYST, TEMP >*/
/*< 	REAL L2 >*/
/*< 	PARAMETER (L2=1.7, L2LAG=9, L2WID=16, OSHYST=10) >*/
/*< 	PARAMETER (TEMP=1+L2WID/2) >*/
/*       Local variables that need not be saved */
/*< 	INTEGER I >*/
/*< 	REAL L2SUM2 >*/
/*       Local state */
/*   Variables */
/*    N, D       Numerator and denominator of prediction filters */
/*    FPC        Current prediction coefs */
/*    L2BUF, L2SUM1, L2SUM2    State of slope filter */
/*       The only "significant" change I've made is to change L2SUM2 out 
*/
/*       of the list of local variables that need to be saved, since it */
/*       didn't need to be. */
/*       L2SUM1 need not be, but avoiding saving it would require a small 
*/
/*       change to the body of the code.  See comments below for an */
/*       example of how the code could be changed to avoid saving L2SUM1. 
*/
/*       FPC and LASTI are saved from one invocation to the next, but */
/*       they are not given initial values.  This is acceptable, because 
*/
/*       FPC will be assigned a value the first time that this function */
/*       is called after D is initialized to 1, since the formula to */
/*       change D will not change it to 0 in one step, and the IF (D */
/*       .NE. 0) statement will execute its THEN part, initializing FPC. 
*/

/*       LASTI's value will not be used until HYST is .TRUE., and */
/*       whenever HYST is changed from its initial value of .FALSE., */
/*       LASTI is assigned a value. */
/*       In a C version of this coder, it would be nice if all of these */
/*       saved things, in this and all other subroutines, could be stored 
*/
/*       in a single struct lpc10_coder_state_t, initialized with a call 
*/
/*       to a function like lpc10_init(&lpc10_coder_state).  In this way, 
*/
/*       a program that used these functions could conveniently alternate 
//.........这里部分代码省略.........

//.........这里部分代码省略.........

/*        If remaining matrix is 2-by-2, use SLAE2 or SLAEV2 */
/*        to compute its eigensystem. */

	if (m == l + 1) {
	    if (icompz > 0) {
		slaev2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
		work[l] = c__;
		work[*n - 1 + l] = s;
		slasr_("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], &
			z__[l * z_dim1 + 1], ldz);
	    } else {
		slae2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
	    }
	    d__[l] = rt1;
	    d__[l + 1] = rt2;
	    e[l] = 0.f;
	    l += 2;
	    if (l <= lend) {
		goto L40;
	    }
	    goto L140;
	}

	if (jtot == nmaxit) {
	    goto L140;
	}
	++jtot;

/*        Form shift. */

	g = (d__[l + 1] - p) / (e[l] * 2.f);
	r__ = slapy2_(&g, &c_b10);
	g = d__[m] - p + e[l] / (g + r_sign(&r__, &g));

	s = 1.f;
	c__ = 1.f;
	p = 0.f;

/*        Inner loop */

	mm1 = m - 1;
	i__1 = l;
	for (i__ = mm1; i__ >= i__1; --i__) {
	    f = s * e[i__];
	    b = c__ * e[i__];
	    slartg_(&g, &f, &c__, &s, &r__);
	    if (i__ != m - 1) {
		e[i__ + 1] = r__;
	    }
	    g = d__[i__ + 1] - p;
	    r__ = (d__[i__] - g) * s + c__ * 2.f * b;
	    p = s * r__;
	    d__[i__ + 1] = g + p;
	    g = c__ * r__ - b;

/*           If eigenvectors are desired, then save rotations. */

	    if (icompz > 0) {
		work[i__] = c__;
		work[*n - 1 + i__] = -s;
	    }

	}

/*        If eigenvectors are desired, then apply saved rotations. */

/* Subroutine */
int cla_heamv_(integer *uplo, integer *n, real *alpha, complex *a, integer *lda, complex *x, integer *incx, real *beta, real *y, integer *incy)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2, i__3;
    real r__1, r__2;
    /* Builtin functions */
    double r_imag(complex *), r_sign(real *, real *);
    /* Local variables */
    integer i__, j;
    logical symb_zero__;
    integer iy, jx, kx, ky, info;
    real temp, safe1;
    extern real slamch_(char *);
    extern /* Subroutine */
    int xerbla_(char *, integer *);
    extern integer ilauplo_(char *);
    /* -- LAPACK computational routine (version 3.4.2) -- */
    /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
    /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
    /* September 2012 */
    /* .. Scalar Arguments .. */
    /* .. */
    /* .. Array Arguments .. */
    /* .. */
    /* ===================================================================== */
    /* .. Parameters .. */
    /* .. */
    /* .. Local Scalars .. */
    /* .. */
    /* .. External Subroutines .. */
    /* .. */
    /* .. External Functions .. */
    /* .. */
    /* .. Intrinsic Functions .. */
    /* .. */
    /* .. Statement Functions .. */
    /* .. */
    /* .. Statement Function Definitions .. */
    /* .. */
    /* .. Executable Statements .. */
    /* Test the input parameters. */
    /* Parameter adjustments */
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    --x;
    --y;
    /* Function Body */
    info = 0;
    if (*uplo != ilauplo_("U") && *uplo != ilauplo_("L") )
    {
        info = 1;
    }
    else if (*n < 0)
    {
        info = 2;
    }
    else if (*lda < max(1,*n))
    {
        info = 5;
    }
    else if (*incx == 0)
    {
        info = 7;
    }
    else if (*incy == 0)
    {
        info = 10;
    }
    if (info != 0)
    {
        xerbla_("CHEMV ", &info);
        return 0;
    }
    /* Quick return if possible. */
    if (*n == 0 || *alpha == 0.f && *beta == 1.f)
    {
        return 0;
    }
    /* Set up the start points in X and Y. */
    if (*incx > 0)
    {
        kx = 1;
    }
    else
    {
        kx = 1 - (*n - 1) * *incx;
    }
    if (*incy > 0)
    {
        ky = 1;
    }
    else
    {
        ky = 1 - (*n - 1) * *incy;
    }
    /* Set SAFE1 essentially to be the underflow threshold times the */
    /* number of additions in each row. */
    safe1 = slamch_("Safe minimum");
//.........这里部分代码省略.........

/* Subroutine */ int sdrvgg_(integer *nsizes, integer *nn, integer *ntypes, 
	logical *dotype, integer *iseed, real *thresh, real *thrshn, integer *
	nounit, real *a, integer *lda, real *b, real *s, real *t, real *s2, 
	real *t2, real *q, integer *ldq, real *z__, real *alphr1, real *
	alphi1, real *beta1, real *alphr2, real *alphi2, real *beta2, real *
	vl, real *vr, real *work, integer *lwork, real *result, integer *info)
{
    /* Initialized data */

    static integer kclass[26] = { 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,
	    2,2,2,3 };
    static integer kbmagn[26] = { 1,1,1,1,1,1,1,1,3,2,3,2,2,3,1,1,1,1,1,1,1,3,
	    2,3,2,1 };
    static integer ktrian[26] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,
	    1,1,1,1 };
    static integer iasign[26] = { 0,0,0,0,0,0,2,0,2,2,0,0,2,2,2,0,2,0,0,0,2,2,
	    2,2,2,0 };
    static integer ibsign[26] = { 0,0,0,0,0,0,0,2,0,0,2,2,0,0,2,0,2,0,0,0,0,0,
	    0,0,0,0 };
    static integer kz1[6] = { 0,1,2,1,3,3 };
    static integer kz2[6] = { 0,0,1,2,1,1 };
    static integer kadd[6] = { 0,0,0,0,3,2 };
    static integer katype[26] = { 0,1,0,1,2,3,4,1,4,4,1,1,4,4,4,2,4,5,8,7,9,4,
	    4,4,4,0 };
    static integer kbtype[26] = { 0,0,1,1,2,-3,1,4,1,1,4,4,1,1,-4,2,-4,8,8,8,
	    8,8,8,8,8,0 };
    static integer kazero[26] = { 1,1,1,1,1,1,2,1,2,2,1,1,2,2,3,1,3,5,5,5,5,3,
	    3,3,3,1 };
    static integer kbzero[26] = { 1,1,1,1,1,1,1,2,1,1,2,2,1,1,4,1,4,6,6,6,6,4,
	    4,4,4,1 };
    static integer kamagn[26] = { 1,1,1,1,1,1,1,1,2,3,2,3,2,3,1,1,1,1,1,1,1,2,
	    3,3,2,1 };

    /* Format strings */
    static char fmt_9999[] = "(\002 SDRVGG: \002,a,\002 returned INFO=\002,i"
	    "6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002, ISEED="
	    "(\002,3(i5,\002,\002),i5,\002)\002)";
    static char fmt_9997[] = "(\002 SDRVGG: SGET53 returned INFO=\002,i1,"
	    "\002 for eigenvalue \002,i6,\002.\002,/9x,\002N=\002,i6,\002, JT"
	    "YPE=\002,i6,\002, ISEED=(\002,3(i5,\002,\002),i5,\002)\002)";
    static char fmt_9996[] = "(\002 SDRVGG: S not in Schur form at eigenvalu"
	    "e \002,i6,\002.\002,/9x,\002N=\002,i6,\002, JTYPE=\002,i6,\002, "
	    "ISEED=(\002,3(i5,\002,\002),i5,\002)\002)";
    static char fmt_9998[] = "(\002 SDRVGG: \002,a,\002 Eigenvectors from"
	    " \002,a,\002 incorrectly \002,\002normalized.\002,/\002 Bits of "
	    "error=\002,0p,g10.3,\002,\002,9x,\002N=\002,i6,\002, JTYPE=\002,"
	    "i6,\002, ISEED=(\002,3(i5,\002,\002),i5,\002)\002)";
    static char fmt_9995[] = "(/1x,a3,\002 -- Real Generalized eigenvalue pr"
	    "oblem driver\002)";
    static char fmt_9994[] = "(\002 Matrix types (see SDRVGG for details):"
	    " \002)";
    static char fmt_9993[] = "(\002 Special Matrices:\002,23x,\002(J'=transp"
	    "osed Jordan block)\002,/\002   1=(0,0)  2=(I,0)  3=(0,I)  4=(I,I"
	    ")  5=(J',J')  \002,\0026=(diag(J',I), diag(I,J'))\002,/\002 Diag"
	    "onal Matrices:  ( \002,\002D=diag(0,1,2,...) )\002,/\002   7=(D,"
	    "I)   9=(large*D, small*I\002,\002)  11=(large*I, small*D)  13=(l"
	    "arge*D, large*I)\002,/\002   8=(I,D)  10=(small*D, large*I)  12="
	    "(small*I, large*D) \002,\002 14=(small*D, small*I)\002,/\002  15"
	    "=(D, reversed D)\002)";
    static char fmt_9992[] = "(\002 Matrices Rotated by Random \002,a,\002 M"
	    "atrices U, V:\002,/\002  16=Transposed Jordan Blocks            "
	    " 19=geometric \002,\002alpha, beta=0,1\002,/\002  17=arithm. alp"
	    "ha&beta             \002,\002      20=arithmetic alpha, beta=0,"
	    "1\002,/\002  18=clustered \002,\002alpha, beta=0,1            21"
	    "=random alpha, beta=0,1\002,/\002 Large & Small Matrices:\002,"
	    "/\002  22=(large, small)   \002,\00223=(small,large)    24=(smal"
	    "l,small)    25=(large,large)\002,/\002  26=random O(1) matrices"
	    ".\002)";
    static char fmt_9991[] = "(/\002 Tests performed:  (S is Schur, T is tri"
	    "angular, \002,\002Q and Z are \002,a,\002,\002,/20x,\002l and r "
	    "are the appropriate left and right\002,/19x,\002eigenvectors, re"
	    "sp., a is alpha, b is beta, and\002,/19x,a,\002 means \002,a,"
	    "\002.)\002,/\002 1 = | A - Q S Z\002,a,\002 | / ( |A| n ulp )   "
	    "   2 = | B - Q T Z\002,a,\002 | / ( |B| n ulp )\002,/\002 3 = | "
	    "I - QQ\002,a,\002 | / ( n ulp )             4 = | I - ZZ\002,a"
	    ",\002 | / ( n ulp )\002,/\002 5 = difference between (alpha,beta"
	    ") and diagonals of\002,\002 (S,T)\002,/\002 6 = max | ( b A - a "
	    "B )\002,a,\002 l | / const.   7 = max | ( b A - a B ) r | / cons"
	    "t.\002,/1x)";
    static char fmt_9990[] = "(\002 Matrix order=\002,i5,\002, type=\002,i2"
	    ",\002, seed=\002,4(i4,\002,\002),\002 result \002,i3,\002 is\002"
	    ",0p,f8.2)";
    static char fmt_9989[] = "(\002 Matrix order=\002,i5,\002, type=\002,i2"
	    ",\002, seed=\002,4(i4,\002,\002),\002 result \002,i3,\002 is\002"
	    ",1p,e10.3)";

    /* System generated locals */
    integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, s_dim1, 
	    s_offset, s2_dim1, s2_offset, t_dim1, t_offset, t2_dim1, 
	    t2_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, z_dim1, 
	    z_offset, i__1, i__2, i__3, i__4;
    real r__1, r__2, r__3, r__4, r__5, r__6, r__7, r__8, r__9, r__10;

    /* Builtin functions */
    double r_sign(real *, real *);
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    integer j, n, i1, n1, jc, nb, in, jr, ns, nbz;
    real ulp;
//.........这里部分代码省略.........

/* Subroutine */ int slaed3_(integer *k, integer *n, integer *n1, real *d__, 
	real *q, integer *ldq, real *rho, real *dlamda, real *q2, integer *
	indx, integer *ctot, real *w, real *s, integer *info)
{
    /* System generated locals */
    integer q_dim1, q_offset, i__1, i__2;
    real r__1;

    /* Builtin functions */
    double sqrt(doublereal), r_sign(real *, real *);

    /* Local variables */
    integer i__, j, n2, n12, ii, n23, iq2;
    real temp;
    extern doublereal snrm2_(integer *, real *, integer *);
    extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 
	    integer *, real *, real *, integer *, real *, integer *, real *, 
	    real *, integer *), scopy_(integer *, real *, 
	    integer *, real *, integer *), slaed4_(integer *, integer *, real 
	    *, real *, real *, real *, real *, integer *);
    extern doublereal slamc3_(real *, real *);
    extern /* Subroutine */ int xerbla_(char *, integer *), slacpy_(
	    char *, integer *, integer *, real *, integer *, real *, integer *
), slaset_(char *, integer *, integer *, real *, real *, 
	    real *, integer *);


/*  -- LAPACK routine (version 3.2) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  SLAED3 finds the roots of the secular equation, as defined by the */
/*  values in D, W, and RHO, between 1 and K.  It makes the */
/*  appropriate calls to SLAED4 and then updates the eigenvectors by */
/*  multiplying the matrix of eigenvectors of the pair of eigensystems */
/*  being combined by the matrix of eigenvectors of the K-by-K system */
/*  which is solved here. */

/*  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 X-MP, Cray Y-MP, Cray C-90, or Cray-2. */
/*  It could conceivably fail on hexadecimal or decimal machines */
/*  without guard digits, but we know of none. */

/*  Arguments */
/*  ========= */

/*  K       (input) INTEGER */
/*          The number of terms in the rational function to be solved by */
/*          SLAED4.  K >= 0. */

/*  N       (input) INTEGER */
/*          The number of rows and columns in the Q matrix. */
/*          N >= K (deflation may result in N>K). */

/*  N1      (input) INTEGER */
/*          The location of the last eigenvalue in the leading submatrix. */
/*          min(1,N) <= N1 <= N/2. */

/*  D       (output) REAL array, dimension (N) */
/*          D(I) contains the updated eigenvalues for */
/*          1 <= I <= K. */

/*  Q       (output) REAL array, dimension (LDQ,N) */
/*          Initially the first K columns are used as workspace. */
/*          On output the columns 1 to K contain */
/*          the updated eigenvectors. */

/*  LDQ     (input) INTEGER */
/*          The leading dimension of the array Q.  LDQ >= max(1,N). */

/*  RHO     (input) REAL */
/*          The value of the parameter in the rank one update equation. */
/*          RHO >= 0 required. */

/*  DLAMDA  (input/output) REAL 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. May be changed on output by */
/*          having lowest order bit set to zero on Cray X-MP, Cray Y-MP, */
/*          Cray-2, or Cray C-90, as described above. */

/*  Q2      (input) REAL array, dimension (LDQ2, N) */
/*          The first K columns of this matrix contain the non-deflated */
/*          eigenvectors for the split problem. */

/*  INDX    (input) INTEGER array, dimension (N) */
/*          The permutation used to arrange the columns of the deflated */
/*          Q matrix into three groups (see SLAED2). */
/*          The rows of the eigenvectors found by SLAED4 must be likewise */
/*          permuted before the matrix multiply can take place. */
//.........这里部分代码省略.........

//.........这里部分代码省略.........
/*                 linear interpolation, Constructive Aspects of the */
/*                 Fundamental Theorem of Algebra, edited by B. Dejon */
/*                 and P. Henrici, Wiley-Interscience, 1969. */
/* ***ROUTINES CALLED  R1MACH */
/* ***REVISION HISTORY  (YYMMDD) */
/*   700901  DATE WRITTEN */
/*   890531  Changed all specific intrinsics to generic.  (WRB) */
/*   890531  REVISION DATE from Version 3.2 */
/*   891214  Prologue converted to Version 4.0 format.  (BAB) */
/*   920501  Reformatted the REFERENCES section.  (WRB) */
/* ***END PROLOGUE  FZERO */
/* ***FIRST EXECUTABLE STATEMENT  FZERO */

/*   ER is two times the computer unit roundoff value which is defined */
/*   here by the function R1MACH. */

    er = r1mach_(&c__4) * 2.f;

/*   Initialize. */

    z__ = *r__;
    if (*r__ <= dmin(*b,*c__) || *r__ >= dmax(*b,*c__)) {
	z__ = *c__;
    }
    rw = dmax(*re,er);
    aw = dmax(*ae,0.f);
    ic = 0;
    t = z__;
    fz = (*f)(&t);
    fc = fz;
    t = *b;
    fb = (*f)(&t);
    kount = 2;
    if (r_sign(&c_b3, &fz) == r_sign(&c_b3, &fb)) {
	goto L1;
    }
    *c__ = z__;
    goto L2;
L1:
    if (z__ == *c__) {
	goto L2;
    }
    t = *c__;
    fc = (*f)(&t);
    kount = 3;
    if (r_sign(&c_b3, &fz) == r_sign(&c_b3, &fc)) {
	goto L2;
    }
    *b = z__;
    fb = fz;
L2:
    a = *c__;
    fa = fc;
    acbs = (r__1 = *b - *c__, dabs(r__1));
/* Computing MAX */
    r__1 = dabs(fb), r__2 = dabs(fc);
    fx = dmax(r__1,r__2);

L3:
    if (dabs(fc) >= dabs(fb)) {
	goto L4;
    }

/*   Perform interchange. */

    a = *b;

/* Subroutine */ int slaic1_(integer *job, integer *j, real *x, real *sest, 
	real *w, real *gamma, real *sestpr, real *s, real *c__)
{
    /* System generated locals */
    real r__1, r__2, r__3, r__4;

    /* Builtin functions */
    double sqrt(doublereal), r_sign(real *, real *);

    /* Local variables */
    static real sine;
    extern doublereal sdot_(integer *, real *, integer *, real *, integer *);
    static real test, zeta1, zeta2, b, t, alpha, norma, s1, s2, absgam, 
	    absalp;
    extern doublereal slamch_(char *);
    static real cosine, absest, eps, tmp;


/*  -- 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   
    =======   

    SLAIC1 applies one step of incremental condition estimation in   
    its simplest version:   

    Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j   
    lower triangular matrix L, such that   
             twonorm(L*x) = sest   
    Then SLAIC1 computes sestpr, s, c such that   
    the vector   
                    [ s*x ]   
             xhat = [  c  ]   
    is an approximate singular vector of   
                    [ L     0  ]   
             Lhat = [ w' gamma ]   
    in the sense that   
             twonorm(Lhat*xhat) = sestpr.   

    Depending on JOB, an estimate for the largest or smallest singular   
    value is computed.   

    Note that [s c]' and sestpr**2 is an eigenpair of the system   

        diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]   
                                              [ gamma ]   

    where  alpha =  x'*w.   

    Arguments   
    =========   

    JOB     (input) INTEGER   
            = 1: an estimate for the largest singular value is computed.   
            = 2: an estimate for the smallest singular value is computed.   

    J       (input) INTEGER   
            Length of X and W   

    X       (input) REAL array, dimension (J)   
            The j-vector x.   

    SEST    (input) REAL   
            Estimated singular value of j by j matrix L   

    W       (input) REAL array, dimension (J)   
            The j-vector w.   

    GAMMA   (input) REAL   
            The diagonal element gamma.   

    SESTPR  (output) REAL   
            Estimated singular value of (j+1) by (j+1) matrix Lhat.   

    S       (output) REAL   
            Sine needed in forming xhat.   

    C       (output) REAL   
            Cosine needed in forming xhat.   

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


       Parameter adjustments */
    --w;
    --x;

    /* Function Body */
    eps = slamch_("Epsilon");
    alpha = sdot_(j, &x[1], &c__1, &w[1], &c__1);

    absalp = dabs(alpha);
    absgam = dabs(*gamma);
    absest = dabs(*sest);

    if (*job == 1) {
//.........这里部分代码省略.........

/* Subroutine */ int slanv2_(real *a, real *b, real *c__, real *d__, real *
	rt1r, real *rt1i, real *rt2r, real *rt2i, real *cs, real *sn)
{
    /* System generated locals */
    real r__1;

    /* Builtin functions */
    double r_sign(real *, real *), sqrt(doublereal);

    /* Local variables */
    static real p, aa, bb, cc, dd, cs1, sn1, sab, sac, tau, temp, sigma;
    extern doublereal slapy2_(real *, real *);


/*  -- LAPACK auxiliary routine (version 2.0) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/*     Courant Institute, Argonne National Lab, and Rice University */
/*     September 30, 1994 */

/*     .. Scalar Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric */
/*  matrix in standard form: */

/*       [ A  B ] = [ CS -SN ] [ AA  BB ] [ CS  SN ] */
/*       [ C  D ]   [ SN  CS ] [ CC  DD ] [-SN  CS ] */

/*  where either */
/*  1) CC = 0 so that AA and DD are real eigenvalues of the matrix, or */
/*  2) AA = DD and BB*CC < 0, so that AA + or - sqrt(BB*CC) are complex */
/*  conjugate eigenvalues. */

/*  Arguments */
/*  ========= */

/*  A       (input/output) REAL */
/*  B       (input/output) REAL */
/*  C       (input/output) REAL */
/*  D       (input/output) REAL */
/*          On entry, the elements of the input matrix. */
/*          On exit, they are overwritten by the elements of the */
/*          standardised Schur form. */

/*  RT1R    (output) REAL */
/*  RT1I    (output) REAL */
/*  RT2R    (output) REAL */
/*  RT2I    (output) REAL */
/*          The real and imaginary parts of the eigenvalues. If the */
/*          eigenvalues are both real, abs(RT1R) >= abs(RT2R); if the */
/*          eigenvalues are a complex conjugate pair, RT1I > 0. */

/*  CS      (output) REAL */
/*  SN      (output) REAL */
/*          Parameters of the rotation matrix. */

/*  ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

/*     Initialize CS and SN */

    *cs = 1.f;
    *sn = 0.f;

    if (*c__ == 0.f) {
	goto L10;

    } else if (*b == 0.f) {

/*        Swap rows and columns */

	*cs = 0.f;
	*sn = 1.f;
	temp = *d__;
	*d__ = *a;
	*a = temp;
	*b = -(*c__);
	*c__ = 0.f;
	goto L10;
    } else if (*a - *d__ == 0.f && r_sign(&c_b3, b) != r_sign(&c_b3, c__)) {
	goto L10;
    } else {

/*        Make diagonal elements equal */

	temp = *a - *d__;
	p = temp * .5f;
	sigma = *b + *c__;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	c2i = dabs(c2i);
	c2m = dmax(c2r,c2i);
	if (c2m <= ascle) {
	    goto L120;
	}
	++kflag;
	ascle = bry[kflag - 1];
	q__1.r = s1.r * c1.r - s1.i * c1.i, q__1.i = s1.r * c1.i + s1.i * 
		c1.r;
	s1.r = q__1.r, s1.i = q__1.i;
	s2.r = c2.r, s2.i = c2.i;
	i__2 = kflag - 1;
	q__1.r = s1.r * css[i__2].r - s1.i * css[i__2].i, q__1.i = s1.r * css[
		i__2].i + s1.i * css[i__2].r;
	s1.r = q__1.r, s1.i = q__1.i;
	i__2 = kflag - 1;
	q__1.r = s2.r * css[i__2].r - s2.i * css[i__2].i, q__1.i = s2.r * css[
		i__2].i + s2.i * css[i__2].r;
	s2.r = q__1.r, s2.i = q__1.i;
	i__2 = kflag - 1;
	c1.r = csr[i__2].r, c1.i = csr[i__2].i;
L120:
	;
    }
L160:
    if (*mr == 0) {
	return 0;
    }
/* ----------------------------------------------------------------------- */
/*     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0E0 */
/* ----------------------------------------------------------------------- */
    *nz = 0;
    fmr = (real) (*mr);
    sgn = -r_sign(&pi, &fmr);
/* ----------------------------------------------------------------------- */
/*     CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP. */
/* ----------------------------------------------------------------------- */
    q__1.r = 0.f, q__1.i = sgn;
    csgn.r = q__1.r, csgn.i = q__1.i;
    inu = *fnu;
    fnf = *fnu - inu;
    ifn = inu + *n - 1;
    ang = fnf * sgn;
    cpn = cos(ang);
    spn = sin(ang);
    q__1.r = cpn, q__1.i = spn;
    cspn.r = q__1.r, cspn.i = q__1.i;
    if (ifn % 2 == 1) {
	q__1.r = -cspn.r, q__1.i = -cspn.i;
	cspn.r = q__1.r, cspn.i = q__1.i;
    }
    asc = bry[0];
    kk = *n;
    iuf = 0;
    kdflg = 1;
    --ib;
    ic = ib - 1;
    i__1 = *n;
    for (k = 1; k <= i__1; ++k) {
	fn = *fnu + (kk - 1);
/* ----------------------------------------------------------------------- */
/*     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K */
/*     FUNCTION ABOVE */
/* ----------------------------------------------------------------------- */
	m = 3;
	if (*n > 2) {

static void vparms(int32_t vwin[],
                   float *inbuf,
                   float *lpbuf,
                   const int32_t buflim[],
                   int32_t half,
                   float *dither,
                   int32_t *mintau,
                   int32_t *zc, 
                   int32_t *lbe,
                   int32_t *fbe,
                   float *qs,
                   float *rc1,
                   float *ar_b,
                   float *ar_f)
{
    int32_t inbuf_offset;
    int32_t lpbuf_offset;
    int32_t vlen;
    int32_t stop;
    int32_t i;
    int32_t start;
    float r1;
    float r2;
    float e_pre;
    float ap_rms;
    float e_0;
    float oldsgn;
    float lp_rms;
    float e_b;
    float e_f;
    float r_b;
    float r_f;
    float e0ap;

    /* Calculate zero crossings (ZC) and several energy and correlation */
    /* measures on low band and full band speech.  Each measure is taken */
    /* over either the first or the second half of the voicing window, */
    /* depending on the variable HALF. */
    lpbuf_offset = buflim[2];
    lpbuf -= lpbuf_offset;
    inbuf_offset = buflim[0];
    inbuf -= inbuf_offset;

    lp_rms = 0.0f;
    ap_rms = 0.0f;
    e_pre = 0.0f;
    e0ap = 0.0f;
    *rc1 = 0.0f;
    e_0 = 0.0f;
    e_b = 0.0f;
    e_f = 0.0f;
    r_f = 0.0f;
    r_b = 0.0f;
    *zc = 0;
    vlen = vwin[1] - vwin[0] + 1;
    start = vwin[0] + half*vlen/2 + 1;
    stop = start + vlen/2 - 1;

    /* I'll use the symbol HVL in the table below to represent the value */
    /* VLEN/2.  Note that if VLEN is odd, then HVL should be rounded down, */
    /* i.e., HVL = (VLEN-1)/2. */

    /* HALF  START          STOP */

    /* 1     VWIN(1)+1      VWIN(1)+HVL */
    /* 2     VWIN(1)+HVL+1  VWIN(1)+2*HVL */
    oldsgn = r_sign(1.0f, inbuf[start - 1] - *dither);
    for (i = start;  i <= stop;  i++)
    {
        lp_rms += fabsf(lpbuf[i]);
        ap_rms += fabsf(inbuf[i]);
        e_pre += fabsf(inbuf[i] - inbuf[i - 1]);
        r1 = inbuf[i];
        e0ap += r1*r1;
        *rc1 += inbuf[i]*inbuf[i - 1];
        r1 = lpbuf[i];
        e_0 += r1*r1;
        r1 = lpbuf[i - *mintau];
        e_b += r1*r1;
        r1 = lpbuf[i + *mintau];
        e_f += r1*r1;
        r_f += lpbuf[i]*lpbuf[i + *mintau];
        r_b += lpbuf[i]*lpbuf[i - *mintau];
        r1 = inbuf[i] + *dither;
        if (r_sign(1.0f, r1) != oldsgn)
        {
            ++(*zc);
            oldsgn = -oldsgn;
        }
        *dither = -(*dither);
    }
    /* Normalized short-term autocovariance coefficient at unit sample delay */
    *rc1 /= max(e0ap, 1.0f);
    /* Ratio of the energy of the first difference signal (6 dB/oct preemphasis)*/
    /* to the energy of the full band signal */
    /* Computing MAX */
    r1 = ap_rms*2.0f;
    *qs = e_pre/max(r1, 1.0f);
    /* aR_b is the product of the forward and reverse prediction gains, */
    /* looking backward in time (the causal case). */
//.........这里部分代码省略.........

/* Subroutine */ int ssterf_(integer *n, real *d__, real *e, integer *info)
{
    /* System generated locals */
    integer i__1;
    real r__1, r__2, r__3;

    /* Builtin functions */
    double sqrt(doublereal), r_sign(real *, real *);

    /* Local variables */
    real c__;
    integer i__, l, m;
    real p, r__, s;
    integer l1;
    real bb, rt1, rt2, eps, rte;
    integer lsv;
    real eps2, oldc;
    integer lend, jtot;
    extern /* Subroutine */ int slae2_(real *, real *, real *, real *, real *)
	    ;
    real gamma, alpha, sigma, anorm;
    extern doublereal slapy2_(real *, real *);
    integer iscale;
    real oldgam;
    extern doublereal slamch_(char *);
    real safmin;
    extern /* Subroutine */ int xerbla_(char *, integer *);
    real safmax;
    extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *, 
	    real *, integer *, integer *, real *, integer *, integer *);
    integer lendsv;
    real ssfmin;
    integer nmaxit;
    real ssfmax;
    extern doublereal slanst_(char *, integer *, real *, real *);
    extern /* Subroutine */ int slasrt_(char *, integer *, real *, integer *);


/*  -- LAPACK routine (version 3.1) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  SSTERF computes all eigenvalues of a symmetric tridiagonal matrix */
/*  using the Pal-Walker-Kahan variant of the QL or QR algorithm. */

/*  Arguments */
/*  ========= */

/*  N       (input) INTEGER */
/*          The order of the matrix.  N >= 0. */

/*  D       (input/output) REAL array, dimension (N) */
/*          On entry, the n diagonal elements of the tridiagonal matrix. */
/*          On exit, if INFO = 0, the eigenvalues in ascending order. */

/*  E       (input/output) REAL array, dimension (N-1) */
/*          On entry, the (n-1) subdiagonal elements of the tridiagonal */
/*          matrix. */
/*          On exit, E has been destroyed. */

/*  INFO    (output) INTEGER */
/*          = 0:  successful exit */
/*          < 0:  if INFO = -i, the i-th argument had an illegal value */
/*          > 0:  the algorithm failed to find all of the eigenvalues in */
/*                a total of 30*N iterations; if INFO = i, then i */
/*                elements of E have not converged to zero. */

/*  ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Executable Statements .. */

/*     Test the input parameters. */

    /* Parameter adjustments */
    --e;
    --d__;

    /* Function Body */
    *info = 0;

/*     Quick return if possible */

//.........这里部分代码省略.........

/* Subroutine */ int sla_syamv__(integer *uplo, integer *n, real *alpha, real
                                 *a, integer *lda, real *x, integer *incx, real *beta, real *y,
                                 integer *incy)
{
    /* System generated locals */
    integer a_dim1, a_offset, i__1, i__2;
    real r__1;

    /* Builtin functions */
    double r_sign(real *, real *);

    /* Local variables */
    integer i__, j;
    logical symb_zero__;
    integer iy, jx, kx, ky, info;
    real temp, safe1;
    extern doublereal slamch_(char *);
    extern /* Subroutine */ int xerbla_(char *, integer *);
    extern integer ilauplo_(char *);


    /*     -- LAPACK routine (version 3.2)                                 -- */
    /*     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- */
    /*     -- Jason Riedy of Univ. of California Berkeley.                 -- */
    /*     -- November 2008                                                -- */

    /*     -- LAPACK is a software package provided by Univ. of Tennessee, -- */
    /*     -- Univ. of California Berkeley and NAG Ltd.                    -- */

    /*     .. */
    /*     .. Scalar Arguments .. */
    /*     .. */
    /*     .. Array Arguments .. */
    /*     .. */

    /*  Purpose */
    /*  ======= */

    /*  SLA_SYAMV  performs the matrix-vector operation */

    /*          y := alpha*abs(A)*abs(x) + beta*abs(y), */

    /*  where alpha and beta are scalars, x and y are vectors and A is an */
    /*  n by n symmetric matrix. */

    /*  This function is primarily used in calculating error bounds. */
    /*  To protect against underflow during evaluation, components in */
    /*  the resulting vector are perturbed away from zero by (N+1) */
    /*  times the underflow threshold.  To prevent unnecessarily large */
    /*  errors for block-structure embedded in general matrices, */
    /*  "symbolically" zero components are not perturbed.  A zero */
    /*  entry is considered "symbolic" if all multiplications involved */
    /*  in computing that entry have at least one zero multiplicand. */

    /*  Parameters */
    /*  ========== */

    /*  UPLO   - INTEGER */
    /*           On entry, UPLO specifies whether the upper or lower */
    /*           triangular part of the array A is to be referenced as */
    /*           follows: */

    /*              UPLO = BLAS_UPPER   Only the upper triangular part of A */
    /*                                  is to be referenced. */

    /*              UPLO = BLAS_LOWER   Only the lower triangular part of A */
    /*                                  is to be referenced. */

    /*           Unchanged on exit. */

    /*  N      - INTEGER. */
    /*           On entry, N specifies the number of columns of the matrix A. */
    /*           N must be at least zero. */
    /*           Unchanged on exit. */

    /*  ALPHA  - REAL            . */
    /*           On entry, ALPHA specifies the scalar alpha. */
    /*           Unchanged on exit. */

    /*  A      - REAL             array of DIMENSION ( LDA, n ). */
    /*           Before entry, the leading m by n part of the array A must */
    /*           contain the matrix of coefficients. */
    /*           Unchanged on exit. */

    /*  LDA    - INTEGER. */
    /*           On entry, LDA specifies the first dimension of A as declared */
    /*           in the calling (sub) program. LDA must be at least */
    /*           max( 1, n ). */
    /*           Unchanged on exit. */

    /*  X      - REAL             array of DIMENSION at least */
    /*           ( 1 + ( n - 1 )*abs( INCX ) ) */
    /*           Before entry, the incremented array X must contain the */
    /*           vector x. */
    /*           Unchanged on exit. */

    /*  INCX   - INTEGER. */
    /*           On entry, INCX specifies the increment for the elements of */
    /*           X. INCX must not be zero. */
    /*           Unchanged on exit. */
//.........这里部分代码省略.........

/* Subroutine */ int slaed3_(integer *k, integer *n, integer *n1, real *d__, 
	real *q, integer *ldq, real *rho, real *dlamda, real *q2, integer *
	indx, integer *ctot, real *w, real *s, integer *info)
{
    /* System generated locals */
    integer q_dim1, q_offset, i__1, i__2;
    real r__1;

    /* Builtin functions */
    double sqrt(doublereal), r_sign(real *, real *);

    /* Local variables */
    static real temp;
    extern doublereal snrm2_(integer *, real *, integer *);
    static integer i__, j;
    extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 
	    integer *, real *, real *, integer *, real *, integer *, real *, 
	    real *, integer *), scopy_(integer *, real *, 
	    integer *, real *, integer *);
    static integer n2;
    extern /* Subroutine */ int slaed4_(integer *, integer *, real *, real *, 
	    real *, real *, real *, integer *);
    extern doublereal slamc3_(real *, real *);
    static integer n12, ii, n23;
    extern /* Subroutine */ int xerbla_(char *, integer *), slacpy_(
	    char *, integer *, integer *, real *, integer *, real *, integer *
	    ), slaset_(char *, integer *, integer *, real *, real *, 
	    real *, integer *);
    static integer iq2;


#define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1]


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

       Common block to return operation count and iteration count   
       ITCNT is unchanged, OPS is only incremented   

    Purpose   
    =======   

    SLAED3 finds the roots of the secular equation, as defined by the   
    values in D, W, and RHO, between 1 and K.  It makes the   
    appropriate calls to SLAED4 and then updates the eigenvectors by   
    multiplying the matrix of eigenvectors of the pair of eigensystems   
    being combined by the matrix of eigenvectors of the K-by-K system   
    which is solved here.   

    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 X-MP, Cray Y-MP, Cray C-90, or Cray-2.   
    It could conceivably fail on hexadecimal or decimal machines   
    without guard digits, but we know of none.   

    Arguments   
    =========   

    K       (input) INTEGER   
            The number of terms in the rational function to be solved by   
            SLAED4.  K >= 0.   

    N       (input) INTEGER   
            The number of rows and columns in the Q matrix.   
            N >= K (deflation may result in N>K).   

    N1      (input) INTEGER   
            The location of the last eigenvalue in the leading submatrix.   
            min(1,N) <= N1 <= N/2.   

    D       (output) REAL array, dimension (N)   
            D(I) contains the updated eigenvalues for   
            1 <= I <= K.   

    Q       (output) REAL array, dimension (LDQ,N)   
            Initially the first K columns are used as workspace.   
            On output the columns 1 to K contain   
            the updated eigenvectors.   

    LDQ     (input) INTEGER   
            The leading dimension of the array Q.  LDQ >= max(1,N).   

    RHO     (input) REAL   
            The value of the parameter in the rank one update equation.   
            RHO >= 0 required.   

    DLAMDA  (input/output) REAL 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. May be changed on output by   
            having lowest order bit set to zero on Cray X-MP, Cray Y-MP,   
            Cray-2, or Cray C-90, as described above.   

    Q2      (input) REAL array, dimension (LDQ2, N)   
            The first K columns of this matrix contain the non-deflated   
            eigenvectors for the split problem.   
//.........这里部分代码省略.........

 int sbdsdc_(char *uplo, char *compq, int *n, float *d__, 
	float *e, float *u, int *ldu, float *vt, int *ldvt, float *q, 
	int *iq, float *work, int *iwork, int *info)
{
    /* System generated locals */
    int u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
    float r__1;

    /* Builtin functions */
    double r_sign(float *, float *), log(double);

    /* Local variables */
    int i__, j, k;
    float p, r__;
    int z__, ic, ii, kk;
    float cs;
    int is, iu;
    float sn;
    int nm1;
    float eps;
    int ivt, difl, difr, ierr, perm, mlvl,