/* DECK INDXC */
/* Subroutine */ int indxc_(integer *i__, integer *ir, integer *idxc, integer 
	*nc)
{
/* ***BEGIN PROLOGUE  INDXC */
/* ***SUBSIDIARY */
/* ***PURPOSE  Subsidiary to BLKTRI */
/* ***LIBRARY   SLATEC */
/* ***TYPE      INTEGER (INDXC-I) */
/* ***AUTHOR  (UNKNOWN) */
/* ***SEE ALSO  BLKTRI */
/* ***ROUTINES CALLED  (NONE) */
/* ***COMMON BLOCKS    CBLKT */
/* ***REVISION HISTORY  (YYMMDD) */
/*   801001  DATE WRITTEN */
/*   891214  Prologue converted to Version 4.0 format.  (BAB) */
/*   900402  Added TYPE section.  (WRB) */
/* ***END PROLOGUE  INDXC */
/* ***FIRST EXECUTABLE STATEMENT  INDXC */
    *nc = pow_ii(&c__2, ir);
    *idxc = *i__;
    if (*idxc + *nc - 1 - cblkt_1.nm <= 0) {
	goto L102;
    } else {
	goto L101;
    }
L101:
    *nc = 0;
L102:
    return 0;
} /* indxc_ */

/*     ---------- ------ */
/* Subroutine */ void fcon_cntdif(integer n, double *d)
{
  /* Local variables */
  integer i, k, k1;
  double sc;
  integer np1;

/* Generates the coefficients of the central difference formula for */
/* Nth derivative at uniformly spaced points */

  d[0] = 1.;
  if (n == 0) {
    return;
  }

  for (i = 0; i < n; ++i) {
    d[i + 1] = 0.;
    for (k = 1; k <= i + 1; ++k) {
      k1 = i + 2 - k;
      d[k1] = d[k1 - 1] - d[k1];
    }
    d[0] = -d[0];
  }

  /* Scale to [0,1]  : */

  sc = (double) pow_ii(n, n);
  np1 = n + 1;
  for (i = 0; i < np1; ++i) {
    d[i] = sc * d[i];
  }

  return;
} /* fcon_cntdif_ */

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

    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 Ren-Cang Li, Computer Science Division, University of   
         California at Berkeley, USA   
       Osni Marques, LBNL/NERSC, USA   

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


       Test the input parameters.   

       Parameter adjustments */
    /* Table of constant values */
    static real c_b9 = 1.f;
    static real c_b10 = 0.f;
    static integer c__2 = 2;
    
    /* System generated locals */
    integer givcol_dim1, givcol_offset, perm_dim1, perm_offset, difl_dim1, 
	    difl_offset, difr_dim1, difr_offset, givnum_dim1, givnum_offset, 
	    poles_dim1, poles_offset, u_dim1, u_offset, vt_dim1, vt_offset, 
	    z_dim1, z_offset, b_dim1, b_offset, bx_dim1, bx_offset, i__1, 
	    i__2, i__3, i__4, i__5, i__6;
    complex q__1;
    /* Builtin functions */
    double r_imag(complex *);
    integer pow_ii(integer *, integer *);
    /* Local variables */
    static integer jcol, nlvl, sqre, jrow, i__, j, jimag, jreal, inode, ndiml;
    extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *, 
	    integer *, real *, real *, integer *, real *, integer *, real *, 
	    real *, integer *);
    static integer ndimr;
    extern /* Subroutine */ int ccopy_(integer *, complex *, integer *, 
	    complex *, integer *);
    static integer i1;
    extern /* Subroutine */ int clals0_(integer *, integer *, integer *, 
	    integer *, integer *, complex *, integer *, complex *, integer *, 
	    integer *, integer *, integer *, integer *, real *, integer *, 
	    real *, real *, real *, real *, integer *, real *, real *, real *,
	     integer *);
    static integer ic, lf, nd, ll, nl, nr;
    extern /* Subroutine */ int xerbla_(char *, integer *), slasdt_(
	    integer *, integer *, integer *, integer *, integer *, integer *, 
	    integer *);
    static integer im1, nlf, nrf, lvl, ndb1, nlp1, lvl2, nrp1;
#define difl_ref(a_1,a_2) difl[(a_2)*difl_dim1 + a_1]
#define difr_ref(a_1,a_2) difr[(a_2)*difr_dim1 + a_1]
#define perm_ref(a_1,a_2) perm[(a_2)*perm_dim1 + a_1]
#define b_subscr(a_1,a_2) (a_2)*b_dim1 + a_1
#define b_ref(a_1,a_2) b[b_subscr(a_1,a_2)]
#define u_ref(a_1,a_2) u[(a_2)*u_dim1 + a_1]
#define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
#define poles_ref(a_1,a_2) poles[(a_2)*poles_dim1 + a_1]
#define bx_subscr(a_1,a_2) (a_2)*bx_dim1 + a_1
#define bx_ref(a_1,a_2) bx[bx_subscr(a_1,a_2)]
#define vt_ref(a_1,a_2) vt[(a_2)*vt_dim1 + a_1]
#define givcol_ref(a_1,a_2) givcol[(a_2)*givcol_dim1 + a_1]
#define givnum_ref(a_1,a_2) givnum[(a_2)*givnum_dim1 + a_1]

/* Subroutine */ int dlasda_(integer *icompq, integer *smlsiz, integer *n, 
	integer *sqre, doublereal *d__, doublereal *e, doublereal *u, integer 
	*ldu, doublereal *vt, integer *k, doublereal *difl, doublereal *difr, 
	doublereal *z__, doublereal *poles, integer *givptr, integer *givcol, 
	integer *ldgcol, integer *perm, doublereal *givnum, doublereal *c__, 
	doublereal *s, doublereal *work, integer *iwork, integer *info)
{
    /* System generated locals */
    integer givcol_dim1, givcol_offset, perm_dim1, perm_offset, difl_dim1, 
	    difl_offset, difr_dim1, difr_offset, givnum_dim1, givnum_offset, 
	    poles_dim1, poles_offset, u_dim1, u_offset, vt_dim1, vt_offset, 
	    z_dim1, z_offset, i__1, i__2;

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

    /* Local variables */
    static doublereal beta;
    static integer idxq, nlvl, i__, j, m;
    static doublereal alpha;
    static integer inode, ndiml, ndimr, idxqi, itemp;
    extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
	    doublereal *, integer *);
    static integer sqrei, i1;
    extern /* Subroutine */ int dlasd6_(integer *, integer *, integer *, 
	    integer *, doublereal *, doublereal *, doublereal *, doublereal *,
	     doublereal *, integer *, integer *, integer *, integer *, 
	    integer *, doublereal *, integer *, doublereal *, doublereal *, 
	    doublereal *, doublereal *, integer *, doublereal *, doublereal *,
	     doublereal *, integer *, integer *);
    static integer ic, nwork1, lf, nd, nwork2, ll, nl, vf, nr, vl;
    extern /* Subroutine */ int 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 *), dlaset_(
	    char *, integer *, integer *, doublereal *, doublereal *, 
	    doublereal *, integer *), xerbla_(char *, integer *);
    static integer im1, smlszp, ncc, nlf, nrf, vfi, iwk, vli, lvl, nru, ndb1, 
	    nlp1, lvl2, nrp1;


#define difl_ref(a_1,a_2) difl[(a_2)*difl_dim1 + a_1]
#define difr_ref(a_1,a_2) difr[(a_2)*difr_dim1 + a_1]
#define perm_ref(a_1,a_2) perm[(a_2)*perm_dim1 + a_1]
#define u_ref(a_1,a_2) u[(a_2)*u_dim1 + a_1]
#define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
#define poles_ref(a_1,a_2) poles[(a_2)*poles_dim1 + a_1]
#define vt_ref(a_1,a_2) vt[(a_2)*vt_dim1 + a_1]
#define givcol_ref(a_1,a_2) givcol[(a_2)*givcol_dim1 + a_1]
#define givnum_ref(a_1,a_2) givnum[(a_2)*givnum_dim1 + a_1]


/*  -- LAPACK auxiliary routine (version 3.0) --   
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
       Courant Institute, Argonne National Lab, and Rice University   
       October 31, 1999   


    Purpose   
    =======   

    Using a divide and conquer approach, DLASDA 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 the singular values in the SVD B = U * S * VT.   
    The orthogonal matrices U and VT are optionally computed in   
    compact form.   

    A related subroutine, DLASD0, computes the singular values and   
    the singular vectors in explicit form.   

    Arguments   
    =========   

    ICOMPQ (input) INTEGER   
           Specifies whether singular vectors are to be computed   
           in compact form, as follows   
           = 0: Compute singular values only.   
           = 1: Compute singular vectors of upper bidiagonal   
                matrix in compact form.   

    SMLSIZ (input) INTEGER   
           The maximum size of the subproblems at the bottom of the   
           computation tree.   

    N      (input) INTEGER   
           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 )   
//.........这里部分代码省略.........

 int claed0_(int *qsiz, int *n, float *d__, float *e, 
	complex *q, int *ldq, complex *qstore, int *ldqs, float *rwork, 
	 int *iwork, int *info)
{
    /* System generated locals */
    int q_dim1, q_offset, qstore_dim1, qstore_offset, i__1, i__2;
    float r__1;

    /* Builtin functions */
    double log(double);
    int pow_ii(int *, int *);

    /* Local variables */
    int i__, j, k, ll, iq, lgn, msd2, smm1, spm1, spm2;
    float temp;
    int curr, iperm;
    extern  int ccopy_(int *, complex *, int *, 
	    complex *, int *);
    int indxq, iwrem;
    extern  int scopy_(int *, float *, int *, float *, 
	    int *);
    int iqptr;
    extern  int claed7_(int *, int *, int *, 
	    int *, int *, int *, float *, complex *, int *, 
	    float *, int *, float *, int *, int *, int *, 
	    int *, int *, float *, complex *, float *, int *, 
	    int *);
    int tlvls;
    extern  int clacrm_(int *, int *, complex *, 
	    int *, float *, int *, complex *, int *, float *);
    int igivcl;
    extern  int xerbla_(char *, int *);
    extern int ilaenv_(int *, char *, char *, int *, int *, 
	    int *, int *);
    int igivnm, submat, curprb, subpbs, igivpt, curlvl, matsiz, iprmpt, 
	    smlsiz;
    extern  int ssteqr_(char *, int *, float *, float *, 
	    float *, int *, float *, int *);


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

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

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

/*  Using the divide and conquer method, CLAED0 computes all eigenvalues */
/*  of a symmetric tridiagonal matrix which is one diagonal block of */
/*  those from reducing a dense or band Hermitian matrix and */
/*  corresponding eigenvectors of the dense or band matrix. */

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

/*  QSIZ   (input) INTEGER */
/*         The dimension of the unitary matrix used to reduce */
/*         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1. */

/*  N      (input) INTEGER */
/*         The dimension of the symmetric tridiagonal matrix.  N >= 0. */

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

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

/*  Q      (input/output) COMPLEX array, dimension (LDQ,N) */
/*         On entry, Q must contain an QSIZ x N matrix whose columns */
/*         unitarily orthonormal. It is a part of the unitary matrix */
/*         that reduces the full dense Hermitian matrix to a */
/*         (reducible) symmetric tridiagonal matrix. */

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

/*  IWORK  (workspace) INTEGER array, */
/*         the dimension of IWORK must be at least */
/*                      6 + 6*N + 5*N*lg N */
/*                      ( lg( N ) = smallest int k */
/*                                  such that 2^k >= N ) */

/*  RWORK  (workspace) REAL array, */
/*                               dimension (1 + 3*N + 2*N*lg N + 3*N**2) */
/*                        ( lg( N ) = smallest int k */
/*                                    such that 2^k >= N ) */

/*  QSTORE (workspace) COMPLEX array, dimension (LDQS, N) */
/*         Used to store parts of */
/*         the eigenvector matrix when the updating matrix multiplies */
/*         take place. */

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

/* DECK MPBLAS */
/* Subroutine */ int mpblas_(integer *i1)
{
    /* System generated locals */
    integer i__1, i__2;

    /* Local variables */
    extern integer i1mach_(integer *);
    static integer mpbexp;
    extern /* Subroutine */ int xermsg_(char *, char *, char *, integer *, 
	    integer *, ftnlen, ftnlen, ftnlen);

/* ***BEGIN PROLOGUE  MPBLAS */
/* ***SUBSIDIARY */
/* ***PURPOSE  Subsidiary to DQDOTA and DQDOTI */
/* ***LIBRARY   SLATEC */
/* ***TYPE      ALL (MPBLAS-A) */
/* ***AUTHOR  (UNKNOWN) */
/* ***DESCRIPTION */

/*     This subroutine is called to set up Brent's 'mp' package */
/*     for use by the extended precision inner products from the BLAS. */

/*     In the SLATEC library we require the Extended Precision MP number */
/*     to have a mantissa twice as long as Double Precision numbers. */
/*     The calculation of MPT (and MPMXR which is the actual array size) */
/*     in this routine will give 2x (or slightly more) on the machine */
/*     that we are running on.  The INTEGER array size of 30 was chosen */
/*     to be slightly longer than the longest INTEGER array needed on */
/*     any machine that we are currently aware of. */

/* ***SEE ALSO  DQDOTA, DQDOTI */
/* ***REFERENCES  R. P. Brent, A Fortran multiple-precision arithmetic */
/*                 package, ACM Transactions on Mathematical Software 4, */
/*                 1 (March 1978), pp. 57-70. */
/*               R. P. Brent, MP, a Fortran multiple-precision arithmetic */
/*                 package, Algorithm 524, ACM Transactions on Mathema- */
/*                 tical Software 4, 1 (March 1978), pp. 71-81. */
/* ***ROUTINES CALLED  I1MACH, XERMSG */
/* ***COMMON BLOCKS    MPCOM */
/* ***REVISION HISTORY  (YYMMDD) */
/*   791001  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) */
/*   900402  Added TYPE section.  (WRB) */
/*   920501  Reformatted the REFERENCES section.  (WRB) */
/*   930124  Increased Array size in MPCON for SUN -r8, and calculate */
/*               size for Quad Precision for 2x DP.  (RWC) */
/* ***END PROLOGUE  MPBLAS */
/* ***FIRST EXECUTABLE STATEMENT  MPBLAS */
    *i1 = 1;

/*     For full extended precision accuracy, MPB should be as large as */
/*     possible, subject to the restrictions in Brent's paper. */

/*     Statements below are for an integer wordlength of  48, 36, 32, */
/*     24, 18, and 16.  Pick one, or generate a new one. */
/*       48     MPB = 4194304 */
/*       36     MPB =   65536 */
/*       32     MPB =   16384 */
/*       24     MPB =    1024 */
/*       18     MPB =     128 */
/*       16     MPB =      64 */

    mpbexp = i1mach_(&c__8) / 2 - 2;
    mpcom_1.mpb = pow_ii(&c__2, &mpbexp);

/*     Set up remaining parameters */
/*                  UNIT FOR ERROR MESSAGES */
    mpcom_1.mplun = i1mach_(&c__4);
/*                  NUMBER OF MP DIGITS */
    mpcom_1.mpt = ((i1mach_(&c__14) << 1) + mpbexp - 1) / mpbexp;
/*                  DIMENSION OF R */
    mpcom_1.mpmxr = mpcom_1.mpt + 4;

    if (mpcom_1.mpmxr > 30) {
	xermsg_("SLATEC", "MPBLAS", "Array space not sufficient for Quad Pre"
		"cision 2x Double Precision, Proceeding.", &c__1, &c__1, (
		ftnlen)6, (ftnlen)6, (ftnlen)78);
	mpcom_1.mpt = 26;
	mpcom_1.mpmxr = 30;
    }
/*                  EXPONENT RANGE */
/* Computing MIN */
    i__1 = 32767, i__2 = i1mach_(&c__9) / 4 - 1;
    mpcom_1.mpm = min(i__1,i__2);
    return 0;
} /* mpblas_ */

 int claed7_(int *n, int *cutpnt, int *qsiz, 
	int *tlvls, int *curlvl, int *curpbm, float *d__, complex *
	q, int *ldq, float *rho, int *indxq, float *qstore, int *
	qptr, int *prmptr, int *perm, int *givptr, int *
	givcol, float *givnum, complex *work, float *rwork, int *iwork, 
	int *info)
{
    /* System generated locals */
    int q_dim1, q_offset, i__1, i__2;

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

    /* Local variables */
    int i__, k, n1, n2, iq, iw, iz, ptr, indx, curr, indxc, indxp;
    extern  int claed8_(int *, int *, int *, 
	    complex *, int *, float *, float *, int *, float *, float *, 
	    complex *, int *, float *, int *, int *, int *, 
	    int *, int *, int *, float *, int *), slaed9_(
	    int *, int *, int *, int *, float *, float *, 
	    int *, float *, float *, float *, float *, int *, int *), 
	    slaeda_(int *, int *, int *, int *, int *, 
	    int *, int *, int *, float *, float *, int *, float *
, float *, int *);
    int idlmda;
    extern  int clacrm_(int *, int *, complex *, 
	    int *, float *, int *, complex *, int *, float *), 
	    xerbla_(char *, int *), slamrg_(int *, int *, 
	    float *, int *, int *, int *);
    int coltyp;


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

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

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

/*  CLAED7 computes the updated eigensystem of a diagonal */
/*  matrix after modification by a rank-one symmetric matrix. This */
/*  routine is used only for the eigenproblem which requires all */
/*  eigenvalues and optionally eigenvectors of a dense or banded */
/*  Hermitian matrix that has been reduced to tridiagonal form. */

/*    T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out) */

/*    where Z = Q'u, u is a vector of length N with ones in the */
/*    CUTPNT and CUTPNT + 1 th elements and zeros elsewhere. */

/*     The eigenvectors of the original matrix are stored in Q, and the */
/*     eigenvalues are in D.  The algorithm consists of three stages: */

/*        The first stage consists of deflating the size of the problem */
/*        when there are multiple eigenvalues or if there is a zero in */
/*        the Z vector.  For each such occurence the dimension of the */
/*        secular equation problem is reduced by one.  This stage is */
/*        performed by the routine SLAED2. */

/*        The second stage consists of calculating the updated */
/*        eigenvalues. This is done by finding the roots of the secular */
/*        equation via the routine SLAED4 (as called by SLAED3). */
/*        This routine also calculates the eigenvectors of the current */
/*        problem. */

/*        The final stage consists of computing the updated eigenvectors */
/*        directly using the updated eigenvalues.  The eigenvectors for */
/*        the current problem are multiplied with the eigenvectors from */
/*        the overall problem. */

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

/*  N      (input) INTEGER */
/*         The dimension of the symmetric tridiagonal matrix.  N >= 0. */

/*  CUTPNT (input) INTEGER */
/*         Contains the location of the last eigenvalue in the leading */
/*         sub-matrix.  MIN(1,N) <= CUTPNT <= N. */

/*  QSIZ   (input) INTEGER */
/*         The dimension of the unitary matrix used to reduce */
/*         the full matrix to tridiagonal form.  QSIZ >= N. */

/*  TLVLS  (input) INTEGER */
/*         The total number of merging levels in the overall divide and */
/*         conquer tree. */

/*  CURLVL (input) INTEGER */
/*         The current level in the overall merge routine, */
/*         0 <= curlvl <= tlvls. */

/*  CURPBM (input) INTEGER */
/*         The current problem in the current level in the overall */
/*         merge routine (counting from upper left to lower right). */
//.........这里部分代码省略.........

/* Subroutine */ int initsv_(integer *indeps)
{
    /* Initialized data */

    static doublereal rvdw[53] = { 1.08,1.,1.8,999.,999.,1.53,1.48,1.36,1.3,
	    999.,2.3,999.,2.05,2.1,1.75,1.7,1.65,999.,2.8,2.75,999.,999.,999.,
	    999.,999.,999.,999.,999.,999.,999.,999.,999.,999.,999.,1.8,999.,
	    999.,999.,999.,999.,999.,999.,999.,999.,999.,999.,999.,999.,999.,
	    999.,999.,999.,2.05 };

    /* System generated locals */
    integer i__1, i__2;
    doublereal d__1;

    /* Builtin functions */
    double sqrt(doublereal);
    integer s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
	    e_wsle(void);
    /* Subroutine */ int s_stop(char *, ftnlen);
    integer i_indx(char *, char *, ftnlen, ftnlen), i_dnnt(doublereal *);
    double log(doublereal);
    integer pow_ii(integer *, integer *);

    /* Local variables */
    static integer i__, n;
    static doublereal x;
    static integer i4;
    static doublereal x0, z3, z4;
#define iw ((integer *)&chanel_1 + 5)
    static integer iat;
    static doublereal epsi, avdw;
    extern doublereal reada_(char *, integer *, ftnlen);
    static doublereal delsc, disex;
#define dirsm ((doublereal *)&solv_1 + 1325)
    static doublereal rsolv;
    static integer indels, indise;
    extern /* Subroutine */ int dvfill_(integer *, doublereal *);
#define dirsmh ((doublereal *)&solv_1 + 4571)
    static integer maxnps, inrsol;
    static doublereal usevdw[53];

    /* Fortran I/O blocks */
    static cilist io___10 = { 0, 0, 0, 0, 0 };
    static cilist io___15 = { 0, 0, 0, 0, 0 };


/* COMDECK SIZES */
/* *********************************************************************** */
/*   THIS FILE CONTAINS ALL THE ARRAY SIZES FOR USE IN MOPAC. */

/*     THERE ARE ONLY 5 PARAMETERS THAT THE PROGRAMMER NEED SET: */
/*     MAXHEV = MAXIMUM NUMBER OF HEAVY ATOMS (HEAVY: NON-HYDROGEN ATOMS) */
/*     MAXLIT = MAXIMUM NUMBER OF HYDROGEN ATOMS. */
/*     MAXTIM = DEFAULT TIME FOR A JOB. (SECONDS) */
/*     MAXDMP = DEFAULT TIME FOR AUTOMATIC RESTART FILE GENERATION (SECS) */
/*     ISYBYL = 1 IF MOPAC IS TO BE USED IN THE SYBYL PACKAGE, =0 OTHERWISE */
/*     SEE ALSO NMECI, NPULAY AND MESP AT THE END OF THIS FILE */


/* *********************************************************************** */

/*   THE FOLLOWING CODE DOES NOT NEED TO BE ALTERED BY THE PROGRAMMER */

/* *********************************************************************** */

/*    ALL OTHER PARAMETERS ARE DERIVED FUNCTIONS OF THESE TWO PARAMETERS */

/*      NAME                   DEFINITION */
/*     NUMATM         MAXIMUM NUMBER OF ATOMS ALLOWED. */
/*     MAXORB         MAXIMUM NUMBER OF ORBITALS ALLOWED. */
/*     MAXPAR         MAXIMUM NUMBER OF PARAMETERS FOR OPTIMISATION. */
/*     N2ELEC         MAXIMUM NUMBER OF TWO ELECTRON INTEGRALS ALLOWED. */
/*     MPACK          AREA OF LOWER HALF TRIANGLE OF DENSITY MATRIX. */
/*     MORB2          SQUARE OF THE MAXIMUM NUMBER OF ORBITALS ALLOWED. */
/*     MAXHES         AREA OF HESSIAN MATRIX */
/*     MAXALL         LARGER THAN MAXORB OR MAXPAR. */
/* *********************************************************************** */

/* *********************************************************************** */
/* DECK MOPAC */
    for (i__ = 1; i__ <= 53; ++i__) {
/* L10: */
	usevdw[i__ - 1] = rvdw[i__ - 1];
    }
    epsi = reada_(keywrd_1.keywrd, indeps, (ftnlen)241);
    solv_1.fepsi = (epsi - 1.) / (epsi + .5);
    solvps_1.nps = 0;
    *iw = 6;
    solv_1.nden = molkst_1.norbs * 3 - (molkst_1.numat << 1);
    maxnps = sqrt(324000.25099999999f) - solv_1.nden - .5f;
    maxnps = min(maxnps,400);
/*     WRITE(IW,*) 'MAXIMUM NUMBER OF SEGMENTS ALLOWED:',MAXNPS */
    if (solv_1.nden * (solv_1.nden + 1) / 2 > 162000) {
	io___10.ciunit = *iw;
	s_wsle(&io___10);
	do_lio(&c__9, &c__1, "PARAMETER LENABC IS TOO SMALL FOR THIS SYSTEM", 
		(ftnlen)45);
	e_wsle();
	s_stop("PARAMETER LENABC IS TOO SMALL FOR THIS SYSTEM", (ftnlen)45);
    }
//.........这里部分代码省略.........

/* ----------------------------------------------------------------------| */
/* Subroutine */ int zgpadm_(integer *ideg, integer *m, doublereal *t, 
	doublecomplex *h__, integer *ldh, doublecomplex *wsp, integer *lwsp, 
	integer *ipiv, integer *iexph, integer *ns, integer *iflag)
{
    /* System generated locals */
    integer h_dim1, h_offset, i__1, i__2, i__3, i__4;
    doublereal d__1, d__2;
    doublecomplex z__1, z__2;

    /* Builtin functions */
    /* Subroutine */ int s_stop(char *, ftnlen);
    double z_abs(doublecomplex *), log(doublereal);
    integer pow_ii(integer *, integer *);

    /* Local variables */
    static integer i__, j, k;
    static doublecomplex cp, cq;
    static integer ip, mm, iq, ih2, iodd, iget, iput, icoef;
    static doublecomplex scale;
    static integer ifree, iused;
    extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, 
	    integer *, doublecomplex *, doublecomplex *, integer *, 
	    doublecomplex *, integer *, doublecomplex *, doublecomplex *, 
	    integer *, ftnlen, ftnlen);
    static doublereal hnorm;
    extern /* Subroutine */ int zgesv_(integer *, integer *, doublecomplex *, 
	    integer *, integer *, doublecomplex *, integer *, integer *);
    static doublecomplex scale2;
    extern /* Subroutine */ int zaxpy_(integer *, doublecomplex *, 
	    doublecomplex *, integer *, doublecomplex *, integer *), zdscal_(
	    integer *, doublereal *, doublecomplex *, integer *);

/* -----Purpose----------------------------------------------------------| */

/*     Computes exp(t*H), the matrix exponential of a general complex */
/*     matrix in full, using the irreducible rational Pade approximation */
/*     to the exponential exp(z) = r(z) = (+/-)( I + 2*(q(z)/p(z)) ), */
/*     combined with scaling-and-squaring. */

/* -----Arguments--------------------------------------------------------| */

/*     ideg      : (input) the degre of the diagonal Pade to be used. */
/*                 a value of 6 is generally satisfactory. */

/*     m         : (input) order of H. */

/*     H(ldh,m)  : (input) argument matrix. */

/*     t         : (input) time-scale (can be < 0). */

/*     wsp(lwsp) : (workspace/output) lwsp .ge. 4*m*m+ideg+1. */

/*     ipiv(m)   : (workspace) */

/* >>>> iexph     : (output) number such that wsp(iexph) points to exp(tH) */
/*                 i.e., exp(tH) is located at wsp(iexph ... iexph+m*m-1) */
/*                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ */
/*                 NOTE: if the routine was called with wsp(iptr), */
/*                       then exp(tH) will start at wsp(iptr+iexph-1). */

/*     ns        : (output) number of scaling-squaring used. */

/*     iflag     : (output) exit flag. */
/*                       0 - no problem */
/*                      <0 - problem */

/* ----------------------------------------------------------------------| */
/*     Roger B. Sidje ([email protected]) */
/*     EXPOKIT: Software Package for Computing Matrix Exponentials. */
/*     ACM - Transactions On Mathematical Software, 24(1):130-156, 1998 */
/* ----------------------------------------------------------------------| */

/* ---  check restrictions on input parameters ... */
    /* Parameter adjustments */
    --ipiv;
    h_dim1 = *ldh;
    h_offset = 1 + h_dim1;
    h__ -= h_offset;
    --wsp;

    /* Function Body */
    mm = *m * *m;
    *iflag = 0;
    if (*ldh < *m) {
	*iflag = -1;
    }
    if (*lwsp < (mm << 2) + *ideg + 1) {
	*iflag = -2;
    }
    if (*iflag != 0) {
	s_stop("bad sizes (in input of ZGPADM)", (ftnlen)30);
    }

/* ---  initialise pointers ... */

    icoef = 1;
    ih2 = icoef + (*ideg + 1);
    ip = ih2 + mm;
    iq = ip + mm;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                          6 + 6*N + 5*N*lg N.   
                          ( lg( N ) = smallest integer k   
                                      such that 2^k >= N )   
           If ICOMPQ = 2, the dimension of IWORK must be at least   
                          2 + 5*N.   

    INFO   (output) INTEGER   
            = 0:  successful exit.   
            < 0:  if INFO = -i, the i-th argument had an illegal value.   
            > 0:  The algorithm failed to compute an eigenvalue while   
                  working on the submatrix lying in rows and columns   
                  INFO/(N+1) through mod(INFO,N+1).   

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


       Test the input parameters.   

    
   Parameter adjustments   
       Function Body */
    /* Table of constant values */
    static integer c__2 = 2;
    static doublereal c_b16 = 1.;
    static doublereal c_b17 = 0.;
    static integer c__1 = 1;
    
    /* System generated locals */
    integer q_dim1, q_offset, qstore_dim1, qstore_offset, i__1, i__2;
    doublereal d__1;
    /* Builtin functions */
    double log(doublereal);
    integer pow_ii(integer *, integer *);
    /* Local variables */
    static doublereal temp;
    static integer curr, i, j, k;
    extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *, 
	    integer *, doublereal *, doublereal *, integer *, doublereal *, 
	    integer *, doublereal *, doublereal *, integer *);
    static integer iperm;
    extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
	    doublereal *, integer *);
    static integer indxq, iwrem;
    extern /* Subroutine */ int dlaed1_(integer *, doublereal *, doublereal *,
	     integer *, integer *, doublereal *, integer *, doublereal *, 
	    integer *, integer *);
    static integer iqptr;
    extern /* Subroutine */ int dlaed7_(integer *, integer *, integer *, 
	    integer *, integer *, integer *, doublereal *, doublereal *, 
	    integer *, integer *, doublereal *, integer *, doublereal *, 
	    integer *, integer *, integer *, integer *, integer *, doublereal 
	    *, doublereal *, integer *, integer *);
    static integer tlvls, iq;
    extern /* Subroutine */ int dlacpy_(char *, integer *, integer *, 
	    doublereal *, integer *, doublereal *, integer *);
    static integer igivcl;
    extern /* Subroutine */ int xerbla_(char *, integer *);
    static integer igivnm, submat, curprb, subpbs, igivpt;
    extern /* Subroutine */ int dsteqr_(char *, integer *, doublereal *, 
	    doublereal *, doublereal *, integer *, doublereal *, integer *);
    static integer curlvl, matsiz, iprmpt, lgn, msd2, smm1, spm1, spm2;



#define D(I) d[(I)-1]

/* Subroutine */ int dvfill_(integer *nppa, doublereal *dirvec)
{
    /* Initialized data */

    static integer kset[60]	/* was [2][30] */ = { 1,2,1,3,1,4,1,5,1,6,12,
	    11,12,10,12,9,12,8,12,7,2,3,3,4,4,5,5,6,6,2,7,8,8,9,9,10,10,11,11,
	    7,2,7,7,3,3,8,8,4,4,9,9,5,5,10,10,6,6,11,11,2 };
    static integer fset[60]	/* was [3][20] */ = { 1,2,3,1,3,4,1,4,5,1,5,6,
	    1,6,2,12,11,10,12,10,9,12,9,8,12,8,7,12,7,11,2,3,7,3,4,8,4,5,9,5,
	    6,10,6,2,11,7,8,3,8,9,4,9,10,5,10,11,6,11,7,2 };

    /* System generated locals */
    integer i__1, i__2;
    doublereal d__1;

    /* Builtin functions */
    double sqrt(doublereal), cos(doublereal), sin(doublereal);
    integer pow_ii(integer *, integer *);
    /* Subroutine */ int s_stop(char *, ftnlen);

    /* Local variables */
    static doublereal h__;
    static integer i__, j, k, l, m;
    static doublereal r__, t;
    static integer j1, j2, na, nb, nc, nd, kh, ix;
    static doublereal beta, dist;

/* 	FUELLEN DES FELDES DIRVEC */
    /* Parameter adjustments */
    dirvec -= 4;

    /* Function Body */
    dirvec[4] = -1.;
    dirvec[5] = 0.;
    dirvec[6] = 0.;
    nd = 1;
    r__ = sqrt(.8);
    h__ = sqrt(.2);
    for (i__ = -1; i__ <= 1; i__ += 2) {
	for (j = 1; j <= 5; ++j) {
	    ++nd;
	    beta = j * 1.25663706 + 1. + (i__ + 1) * .3141593;
	    dirvec[nd * 3 + 2] = r__ * cos(beta);
	    dirvec[nd * 3 + 3] = r__ * sin(beta);
	    dirvec[nd * 3 + 1] = i__ * h__;
/* L10: */
	}
    }
    dirvec[38] = 0.;
    dirvec[39] = 0.;
    dirvec[37] = 1.;
    nd = 12;
/*  NPPA=10*3**K*4**L+2 */
    m = (*nppa - 2) / 10;
    for (k = 0; k <= 10; ++k) {
	if (m / 3 * 3 != m) {
	    goto L30;
	}
/* L20: */
	m /= 3;
    }
L30:
    for (l = 0; l <= 10; ++l) {
	if (m / 4 << 2 != m) {
	    goto L50;
	}
/* L40: */
	m /= 4;
    }
L50:
    if (pow_ii(&c__3, &k) * 10 * pow_ii(&c__4, &l) + 2 != *nppa) {
	s_stop("VALUE OF NPPA NOT ALLOWED:      IT MUST BE 10*3**K*4**L+2", (
		ftnlen)57);
    }
    kh = k / 2;
    m = pow_ii(&c__2, &l) * pow_ii(&c__3, &kh);
/* CREATE ON EACH EDGE 2**L*3**KH-1 NEW POINTS */
    for (i__ = 1; i__ <= 30; ++i__) {
	na = kset[(i__ << 1) - 2];
	nb = kset[(i__ << 1) - 1];
	i__1 = m - 1;
	for (j = 1; j <= i__1; ++j) {
	    ++nd;
	    for (ix = 1; ix <= 3; ++ix) {
/* L60: */
		dirvec[ix + nd * 3] = dirvec[ix + na * 3] * (m - j) + dirvec[
			ix + nb * 3] * j;
	    }
/* L70: */
	}
    }
/* CREATE POINTS WITHIN EACH TRIANGLE */
    for (i__ = 1; i__ <= 20; ++i__) {
	na = fset[i__ * 3 - 3];
	nb = fset[i__ * 3 - 2];
	nc = fset[i__ * 3 - 1];
	i__1 = m - 1;
	for (j1 = 1; j1 <= i__1; ++j1) {
	    i__2 = m - j1 - 1;
	    for (j2 = 1; j2 <= i__2; ++j2) {
//.........这里部分代码省略.........

/* Subroutine */ int zlaed0_(integer *qsiz, integer *n, doublereal *d__, 
	doublereal *e, doublecomplex *q, integer *ldq, doublecomplex *qstore, 
	integer *ldqs, doublereal *rwork, integer *iwork, integer *info)
{
/*  -- LAPACK routine (version 3.0) --   
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
       Courant Institute, Argonne National Lab, and Rice University   
       September 30, 1994   


    Purpose   
    =======   

    Using the divide and conquer method, ZLAED0 computes all eigenvalues   
    of a symmetric tridiagonal matrix which is one diagonal block of   
    those from reducing a dense or band Hermitian matrix and   
    corresponding eigenvectors of the dense or band matrix.   

    Arguments   
    =========   

    QSIZ   (input) INTEGER   
           The dimension of the unitary matrix used to reduce   
           the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.   

    N      (input) INTEGER   
           The dimension of the symmetric tridiagonal matrix.  N >= 0.   

    D      (input/output) DOUBLE PRECISION array, dimension (N)   
           On entry, the diagonal elements of the tridiagonal matrix.   
           On exit, the eigenvalues in ascending order.   

    E      (input/output) DOUBLE PRECISION array, dimension (N-1)   
           On entry, the off-diagonal elements of the tridiagonal matrix.   
           On exit, E has been destroyed.   

    Q      (input/output) COMPLEX*16 array, dimension (LDQ,N)   
           On entry, Q must contain an QSIZ x N matrix whose columns   
           unitarily orthonormal. It is a part of the unitary matrix   
           that reduces the full dense Hermitian matrix to a   
           (reducible) symmetric tridiagonal matrix.   

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

    IWORK  (workspace) INTEGER array,   
           the dimension of IWORK must be at least   
                        6 + 6*N + 5*N*lg N   
                        ( lg( N ) = smallest integer k   
                                    such that 2^k >= N )   

    RWORK  (workspace) DOUBLE PRECISION array,   
                                 dimension (1 + 3*N + 2*N*lg N + 3*N**2)   
                          ( lg( N ) = smallest integer k   
                                      such that 2^k >= N )   

    QSTORE (workspace) COMPLEX*16 array, dimension (LDQS, N)   
           Used to store parts of   
           the eigenvector matrix when the updating matrix multiplies   
           take place.   

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

    INFO   (output) INTEGER   
            = 0:  successful exit.   
            < 0:  if INFO = -i, the i-th argument had an illegal value.   
            > 0:  The algorithm failed to compute an eigenvalue while   
                  working on the submatrix lying in rows and columns   
                  INFO/(N+1) through mod(INFO,N+1).   

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

    Warning:      N could be as big as QSIZ!   


       Test the input parameters.   

       Parameter adjustments */
    /* Table of constant values */
    static integer c__9 = 9;
    static integer c__0 = 0;
    static integer c__2 = 2;
    static integer c__1 = 1;
    
    /* System generated locals */
    integer q_dim1, q_offset, qstore_dim1, qstore_offset, i__1, i__2;
    doublereal d__1;
    /* Builtin functions */
    double log(doublereal);
    integer pow_ii(integer *, integer *);
    /* Local variables */
    static doublereal temp;
    static integer curr, i__, j, k, iperm;
    extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
	    doublereal *, integer *);
    static integer indxq, iwrem, iqptr, tlvls;
    extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, 
	    doublecomplex *, integer *), zlaed7_(integer *, integer *, 
//.........这里部分代码省略.........

/*<       subroutine gpfa3f(a,b,trigs,inc,jump,n,mm,lot,isign) >*/
/* Subroutine */ int gpfa3f_(real *a, real *b, real *trigs, integer *inc,
        integer *jump, integer *n, integer *mm, integer *lot, integer *isign)
{
    /* Initialized data */

    static real sin60 = (float).866025403784437; /* constant */
    static integer lvr = 128; /* constant */

    /* System generated locals */
    integer i__1, i__2, i__3, i__4, i__5, i__6, i__7, i__8, i__9, i__10;

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

    /* Local variables */
    integer j, k, l, m;
    real s, c1;
    integer n3;
    real t1, t2, t3, u1, u2, u3;
    integer ja, jb, la, jc, jd, nb, je, jf, jg, jh, mh, kk, ji, ll, mu, nu;
    real co1, co2, si1, si2, aja, ajb, ajc, bjb, bjc, bja, ajd, bjd, aje, ajf,
             ajh, bje, bjf, bjh, aji, ajg, bji, bjg;
    integer jjj, ink, inq, ninc, left, nvex, ipass, nblox, jstep, laincl,
            jstepl, istart, jstepx;

/*<       real a(*), b(*), trigs(*) >*/
/*<       integer inc, jump, n, mm, lot, isign >*/
/*<       real s, c1, t1, t2, t3, u1, u2, u3, co1, co2 >*/
/*<       real si1, si2, aja, ajb, ajc, bjb, bjc, bja, ajd, bjd >*/
/*<       real aje, ajf, ajh, bje, bjf, bjh, aji, ajg, bji, bjg >*/
/*<       data sin60/0.866025403784437/ >*/
    /* Parameter adjustments */
    --trigs;
    --b;
    --a;

    /* Function Body */
/*<       data lvr/128/ >*/

/*     *************************************************************** */
/*     *                                                             * */
/*     *  N.B. LVR = LENGTH OF VECTOR REGISTERS, SET TO 128 FOR C90. * */
/*     *  RESET TO 64 FOR OTHER CRAY MACHINES, OR TO ANY LARGE VALUE * */
/*     *  (GREATER THAN OR EQUAL TO LOT) FOR A SCALAR COMPUTER.      * */
/*     *                                                             * */
/*     *************************************************************** */

/*<       n3 = 3**mm >*/
    n3 = pow_ii(&c__3, mm);
/*<       inq = n/n3 >*/
    inq = *n / n3;
/*<       jstepx = (n3-n) * inc >*/
    jstepx = (n3 - *n) * *inc;
/*<       ninc = n * inc >*/
    ninc = *n * *inc;
/*<       ink = inc * inq >*/
    ink = *inc * inq;
/*<       mu = mod(inq,3) >*/
    mu = inq % 3;
/*<       if (isign.eq.-1) mu = 3-mu >*/
    if (*isign == -1) {
        mu = 3 - mu;
    }
/*<       m = mm >*/
    m = *mm;
/*<       mh = (m+1)/2 >*/
    mh = (m + 1) / 2;
/*<       s = float(isign) >*/
    s = (real) (*isign);
/*<       c1 = sin60 >*/
    c1 = sin60;
/*<       if (mu.eq.2) c1 = -c1 >*/
    if (mu == 2) {
        c1 = -c1;
    }

/*<       nblox = 1 + (lot-1)/lvr >*/
    nblox = (*lot - 1) / lvr + 1;
/*<       left = lot >*/
    left = *lot;
/*<       s = float(isign) >*/
    s = (real) (*isign);
/*<       istart = 1 >*/
    istart = 1;

/*  loop on blocks of lvr transforms */
/*  -------------------------------- */
/*<       do 500 nb = 1 , nblox >*/
    i__1 = nblox;
    for (nb = 1; nb <= i__1; ++nb) {

/*<       if (left.le.lvr) then >*/
        if (left <= lvr) {
/*<          nvex = left >*/
            nvex = left;
/*<       else if (left.lt.(2*lvr)) then >*/
        } else if (left < lvr << 1) {
/*<          nvex = left/2 >*/
            nvex = left / 2;
//.........这里部分代码省略.........

/* Subroutine */ int decode_(integer *ipitv, integer *irms, 
	integer *irc, integer *voice, integer *pitch, real *rms, real *rc,
			     struct lpc10_decoder_state *st)
{
    /* Initialized data */

    logical *first;
    static integer ethrs = 2048;
    static integer ethrs1 = 128;
    static integer ethrs2 = 1024;
    static integer ethrs3 = 2048;
    static integer ivtab[32] = { 24960,24960,24960,24960,25480,25480,25483,
	    25480,16640,1560,1560,1560,16640,1816,1563,1560,24960,24960,24859,
	    24856,26001,25881,25915,25913,1560,1560,7800,3640,1561,1561,3643,
	    3641 };
    static real corth[32]	/* was [4][8] */ = { 32767.f,10.f,5.f,0.f,
	    32767.f,8.f,4.f,0.f,32.f,6.4f,3.2f,0.f,32.f,6.4f,3.2f,0.f,32.f,
	    11.2f,6.4f,0.f,32.f,11.2f,6.4f,0.f,16.f,5.6f,3.2f,0.f,16.f,5.6f,
	    3.2f,0.f };
    static integer detau[128] = { 0,0,0,3,0,3,3,31,0,3,3,21,3,3,29,30,0,3,3,
	    20,3,25,27,26,3,23,58,22,3,24,28,3,0,3,3,3,3,39,33,32,3,37,35,36,
	    3,38,34,3,3,42,46,44,50,40,48,3,54,3,56,3,52,3,3,1,0,3,3,108,3,78,
	    100,104,3,84,92,88,156,80,96,3,3,74,70,72,66,76,68,3,62,3,60,3,64,
	    3,3,1,3,116,132,112,148,152,3,3,140,3,136,3,144,3,3,1,124,120,128,
	    3,3,3,3,1,3,3,3,1,3,1,1,1 };
    static integer rmst[64] = { 1024,936,856,784,718,656,600,550,502,460,420,
	    384,352,328,294,270,246,226,206,188,172,158,144,132,120,110,102,
	    92,84,78,70,64,60,54,50,46,42,38,34,32,30,26,24,22,20,18,17,16,15,
	    14,13,12,11,10,9,8,7,6,5,4,3,2,1,0 };
    static integer detab7[32] = { 4,11,18,25,32,39,46,53,60,66,72,77,82,87,92,
	    96,101,104,108,111,114,115,117,119,121,122,123,124,125,126,127,
	    127 };
    static real descl[8] = { .6953f,.625f,.5781f,.5469f,.5312f,.5391f,.4688f,
	    .3828f };
    integer *ivp2h;
    static integer deadd[8] = { 1152,-2816,-1536,-3584,-1280,-2432,768,-1920 }
	    ;
    static integer qb[8] = { 511,511,1023,1023,1023,1023,2047,4095 };
    static integer nbit[10] = { 8,8,5,5,4,4,4,4,3,2 };
    static integer zrc[10] = { 0,0,0,0,0,3,0,2,0,0 };
    static integer bit[5] = { 2,4,8,16,32 };
    integer *iovoic;
    integer *iavgp;
    integer *iptold;
    integer *erate;
    integer *drc;
    integer *dpit;
    integer *drms;

    /* System generated locals */
    integer i__1, i__2;

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

    /* Local variables */
    extern /* Subroutine */ int ham84_(integer *, integer *, integer *);
    integer ipit, iout, i__, icorf, index, ivoic, ixcor, i1, i2, i4;
    extern integer median_(integer *, integer *, integer *);
    integer ishift, errcnt, lsb;

/* $Log: decode_.c,v $
/* Revision 1.1  2004/05/04 11:16:42  csoutheren
/* Initial version
/*
/* Revision 1.2  2002/02/15 03:57:55  yurik
/* Warnings removed during compilation, patch courtesy of Jehan Bing, [email protected]
/*
/* Revision 1.1  2000/06/05 04:45:12  robertj
/* Added LPC-10 2400bps codec
/*
 * Revision 1.2  1996/08/20  20:22:39  jaf
 * Removed all static local variables that were SAVE'd in the Fortran
 * code, and put them in struct lpc10_decoder_state that is passed as an
 * argument.
 *
 * Removed init function, since all initialization is now done in
 * init_lpc10_decoder_state().
 *
 * Revision 1.1  1996/08/19  22:32:38  jaf
 * Initial revision
 * */
/* 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 */
/*       Arguments */
/* $Log: decod