/* comcof.f -- translated by f2c (version of 3 February 1990  3:36:42).
   You must link the resulting object file with the libraries:
	-lF77 -lI77 -lm -lc   (in that order)
*/

#include "f2c.h"

/* Common Block Declarations */

struct {
    integer ielmnt, isbckt, nsbckt, iunsat, nunsat, itemps, numtem, isens, 
	    nsens, ifour, nfour, ifield, icode, idelim, icolum, insize, 
	    junode, lsbkpt, numbkp, iorder, jmnode, iur, iuc, ilc, ilr, 
	    numoff, isr, nmoffc, iseq, iseq1, neqn, nodevs, ndiag, iswap, 
	    iequa, macins, lvnim1, lx0, lvn, lynl, lyu, lyl, lx1, lx2, lx3, 
	    lx4, lx5, lx6, lx7, ld0, ld1, ltd, imynl, imvn, lcvn, nsnod, 
	    nsmat, nsval, icnod, icmat, icval, loutpt, lpol, lzer, irswpf, 
	    irswpr, icswpf, icswpr, irpt, jcpt, irowno, jcolno, nttbr, nttar, 
	    lvntmp;
} tabinf_;

#define tabinf_1 tabinf_

struct {
    doublereal omega, time, delta, delold[7], ag[7], vt, xni, egfet, xmu, 
	    sfactr;
    integer mode, modedc, icalc, initf, method, iord, maxord, noncon, iterno, 
	    itemno, nosolv, modac, ipiv, ivmflg, ipostp, iscrch, iofile;
} status_;

#define status_1 status_

struct {
    integer iprnta, iprntl, iprntm, iprntn, iprnto, limtim, limpts, lvlcod, 
	    lvltim, itl1, itl2, itl3, itl4, itl5, itl6, igoof, nogo, keof;
} flags_;

#define flags_1 flags_

struct {
    doublereal value[200000];
} blank_;

#define blank_1 blank_

/*<       subroutine comcof >*/
/* Subroutine */ int comcof_()
{
    /* System generated locals */
    integer i_1, i_2, i_3;

    /* Local variables */
    static doublereal gmat[49]	/* was [7][7] */;
    extern /* Subroutine */ int zero8_();
    static integer i, j, k, istop, ir;
#define nodplc ((integer *)&blank_1)
#define cvalue ((complex *)&blank_1)
    static integer jstart;
    static doublereal arg, arg1;

/*<       implicit double precision (a-h,o-z) >*/

/*     this routine calculates the timestep-dependent terms used in the */

/* numerical integration. */

/* spice version 2g.6  sccsid=tabinf 3/15/83 */
/*<       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem, >*/
/*<      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize, >*/
/*<      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr, >*/
/*<      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1, >*/
/*<      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd, >*/
/*<      5   imynl,imvn,lcvn,nsnod,nsmat,nsval,icnod,icmat,icval, >*/
/*<      6   loutpt,lpol,lzer,irswpf,irswpr,icswpf,icswpr,irpt,jcpt, >*/
/*<      7   irowno,jcolno,nttbr,nttar,lvntmp >*/
/* spice version 2g.6  sccsid=status 3/15/83 */
/*<       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet, >*/
/*<      1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon, >*/
/*<      2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile >*/
/* spice version 2g.6  sccsid=flags 3/15/83 */
/*<       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts, >*/
/*<      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,itl6,igoof,nogo,keof >*/
/* spice version 2g.6  sccsid=blank 3/15/83 */
/*<       common /blank/ value(200000) >*/
/*<       integer nodplc(64) >*/
/*<       complex cvalue(32) >*/
/*<       equivalence (value(1),nodplc(1),cvalue(1)) >*/
/*<       dimension gmat(7,7) >*/

/*  compute coefficients for particular integration method */

/*<       if (method.ne.1) go to 5 >*/
    if (status_1.method != 1) {
	goto L5;
    }
/*<       if (iord.eq.1) go to 5 >*/
    if (status_1.iord == 1) {
	goto L5;
    }
/* ...  trapezoidal method */
/*<       ag(1)=1.0d0/delta/(1.0d0-xmu) >*/
    status_1.ag[0] = 1. / status_1.delta / (1. - status_1.xmu);
/*<       ag(2)=xmu/(1.0d0-xmu) >*/
    status_1.ag[1] = status_1.xmu / (1. - status_1.xmu);
/*<       go to 200 >*/
    goto L200;

/*  construct gear coefficient matrix */

/*<     5 istop=iord+1 >*/
L5:
    istop = status_1.iord + 1;
/*<       call zero8(ag,istop) >*/
    zero8_(status_1.ag, &istop);
/*<       ag(2)=-1.0d0 >*/
    status_1.ag[1] = -1.;
/*<       do 10 i=1,istop >*/
    i_1 = istop;
    for (i = 1; i <= i_1; ++i) {
/*<       gmat(1,i)=1.0d0 >*/
	gmat[i * 7 - 7] = 1.;
/*<    10 continue >*/
/* L10: */
    }
/*<       do 20 i=2,istop >*/
    i_1 = istop;
    for (i = 2; i <= i_1; ++i) {
/*<       gmat(i,1)=0.0d0 >*/
	gmat[i - 1] = 0.;
/*<    20 continue >*/
/* L20: */
    }
/*<       arg=0.0d0 >*/
    arg = 0.;
/*<       do 40 i=2,istop >*/
    i_1 = istop;
    for (i = 2; i <= i_1; ++i) {
/*<       arg=arg+delold(i-1) >*/
	arg += status_1.delold[i - 2];
/*<       arg1=1.0d0 >*/
	arg1 = 1.;
/*<       do 30 j=2,istop >*/
	i_2 = istop;
	for (j = 2; j <= i_2; ++j) {
/*<       arg1=arg1*arg >*/
	    arg1 *= arg;
/*<       gmat(j,i)=arg1 >*/
	    gmat[j + i * 7 - 8] = arg1;
/*<    30 continue >*/
/* L30: */
	}
/*<    40 continue >*/
/* L40: */
    }

/*  solve for gear coefficients ag(*) */


/*  lu decomposition */

/*<       do 70 i=2,istop >*/
    i_1 = istop;
    for (i = 2; i <= i_1; ++i) {
/*<       jstart=i+1 >*/
	jstart = i + 1;
/*<       if (jstart.gt.istop) go to 70 >*/
	if (jstart > istop) {
	    goto L70;
	}
/*<       do 60 j=jstart,istop >*/
	i_2 = istop;
	for (j = jstart; j <= i_2; ++j) {
/*<       gmat(j,i)=gmat(j,i)/gmat(i,i) >*/
	    gmat[j + i * 7 - 8] /= gmat[i + i * 7 - 8];
/*<       do 50 k=jstart,istop >*/
	    i_3 = istop;
	    for (k = jstart; k <= i_3; ++k) {
/*<       gmat(j,k)=gmat(j,k)-gmat(j,i)*gmat(i,k) >*/
		gmat[j + k * 7 - 8] -= gmat[j + i * 7 - 8] * gmat[i + k * 7 - 
			8];
/*<    50 continue >*/
/* L50: */
	    }
/*<    60 continue >*/
/* L60: */
	}
/*<    70 continue >*/
L70:
    ;}

/*  forward substitution */

/*<       do 90 i=2,istop >*/
    i_1 = istop;
    for (i = 2; i <= i_1; ++i) {
/*<       jstart=i+1 >*/
	jstart = i + 1;
/*<       if (jstart.gt.istop) go to 90 >*/
	if (jstart > istop) {
	    goto L90;
	}
/*<       do 80 j=jstart,istop >*/
	i_2 = istop;
	for (j = jstart; j <= i_2; ++j) {
/*<       ag(j)=ag(j)-gmat(j,i)*ag(i) >*/
	    status_1.ag[j - 1] -= gmat[j + i * 7 - 8] * status_1.ag[i - 1];
/*<    80 continue >*/
/* L80: */
	}
/*<    90 continue >*/
L90:
    ;}

/*  backward substitution */

/*<       ag(istop)=ag(istop)/gmat(istop,istop) >*/
    status_1.ag[istop - 1] /= gmat[istop + istop * 7 - 8];
/*<       ir=istop >*/
    ir = istop;
/*<       do 110 i=2,istop >*/
    i_1 = istop;
    for (i = 2; i <= i_1; ++i) {
/*<       jstart=ir >*/
	jstart = ir;
/*<       ir=ir-1 >*/
	--ir;
/*<       do 100 j=jstart,istop >*/
	i_2 = istop;
	for (j = jstart; j <= i_2; ++j) {
/*<       ag(ir)=ag(ir)-gmat(ir,j)*ag(j) >*/
	    status_1.ag[ir - 1] -= gmat[ir + j * 7 - 8] * status_1.ag[j - 1];
/*<   100 continue >*/
/* L100: */
	}
/*<       ag(ir)=ag(ir)/gmat(ir,ir) >*/
	status_1.ag[ir - 1] /= gmat[ir + ir * 7 - 8];
/*<   110 continue >*/
/* L110: */
    }

/*  finished */

/*<   200 return >*/
L200:
    return 0;
/*<       end >*/
} /* comcof_ */

#undef cvalue
#undef nodplc


