/* autlib2.f -- translated by f2c (version of 28 December 1990  16:16:33).
   You must link the resulting object file with the libraries:
	-lF77 -lI77 -lm -lc   (in that order)
*/

#include "f2c.h"
#include "autlim.h"
int FLOWK;
/* Common Block Declarations */

typedef struct {
  int irot;
  int nrot[1000];
  double torper;
} ROTCHK;

extern ROTCHK blrtn;
struct {
    integer ndim, ips, irs, ilp, icp[20];
    doublereal par[20];
} blbcn_;

#define blbcn_1 blbcn_

struct {
    doublereal rdsold, a, rl[20], rlold[20], rldot[20];
} blcrl_;

#define blcrl_1 blcrl_

struct {
    integer ntst, ncol, iad, isp, isw, iplt, nbc, nint;
} blcde_;

#define blcde_1 blcde_

struct {
    integer ndm, ndmp1, nrow, nclm, nrc, ncc, npar, nfpar, nbc0, nint0;
} blicn_;

#define blicn_1 blicn_

struct {
    integer itpst, itpsp, ibrsp;
} blitp_;

#define blitp_1 blitp_

struct {
    doublereal ds, dsmin, dsmax;
    integer iads;
} bldls_;

#define bldls_1 bldls_

struct {
    doublereal detge;
    integer nins;
} bldet_;

#define bldet_1 bldet_

struct {
    integer npr, mxbf, iid, itmx, itnw, nwtn, jac;
} blmax_;

#define blmax_1 blmax_

struct {
    doublereal half, zero, one, two, hmach, rsmall, rlarge;
} blrcn_;

#define blrcn_1 blrcn_

struct {
    integer nmx, nuzr;
    doublereal rl0, rl1, a0, a1;
} bllim_;

#define bllim_1 bllim_

struct {
    integer iuzr;
} blusz_;

#define blusz_1 blusz_

struct {
    integer ndimp1, ndirc, ntstp1, ndcc, ndrhs, ndbc, nuicd, ndicd, nwbr, 
	    niwbr;
} bldim_;

#define bldim_1 bldim_

struct {
    doublereal u0xx[N3AUTO], u1xx[N3AUTO], u2xx[N3AUTO], f1xx[N3AUTO], f2xx[N3AUTO], dfuxx[NFAUTO]	/* was [50][120] */, dfpxx[NPAUTO]	/* was [50][
	    20] */, dduxx[NAUTO], ddpxx[20];
} blwif_;


/*
 *The COMMON blocks BLWIF and BLDIF are used for Interface Workspace.
 *Array dimensions in BLWIF and BLDIF are:
 *
 * BLWIF:
 *       U0XX(3n+2),U1XX(3n+2),U2XX(3n+2),F1XX(3n+2),F2XX(3n+2),
 *       DFUXX(n,2n+20),DFPXX(n,20),DDUXX(n),DDPXX(20)
 * BLDIF:
 *       U1ZZ(n),U2ZZ(n),F1ZZ(n),F2ZZ(n)
 *
 *  where n is the largest value of NDIM or NINT or NBC.
 */

#define blwif_1 blwif_

struct {
    doublereal u1zz[NAUTO], u2zz[NAUTO], f1zz[NAUTO], f2zz[NAUTO];
} bldif_;

#define bldif_1 bldif_

struct {
    doublereal epsl[20], epsu, epss;
} bleps_;

#define bleps_1 bleps_

struct {
    doublereal delref;
} blref_;

#define blref_1 blref_

struct {
    doublereal thetal[20], thetau;
} bltht_;

#define bltht_1 bltht_

struct {
    doublereal w[56]	/* was [8][7] */, wp[56]	/* was [8][7] */, wh[
	    8], wi[8];
} blwts_;

#define blwts_1 blwts_

struct {
    doublereal tsetub, tconpa, tconrh, tinfpa, treduc, twr8;
} bltim_;

#define bltim_1 bltim_

struct {
    integer nam1, nap1, nxe;
} blbnd_;

#define blbnd_1 blbnd_

struct {
    integer ndecom, nbcksb;
} blcnt_;

#define blcnt_1 blcnt_

/* Table of constant values */

static integer c__0 = 0;
static integer c__1 = 1;
static integer c__3 = 3;
static integer c_n1 = -1;
static integer c__2 = 2;

/* ------------------------------------------------------------------- */
/* ------------------------------------------------------------------- */
/*                    General Boundary Value Problems */
/* ------------------------------------------------------------------- */
/* ------------------------------------------------------------------- */

/*     ---------- ------ */
/* Subroutine */ int cnrlbv_(funi, bcni, icni, stpnt, fnbpbv, ibr, m1aa, m2aa,
	 aa, m1bb, m2bb, bb, m1cc, cc, m1dd, dd, wbrbd, m1u, ups, uoldps, 
	upoldp, udotps, rhsa, rhsd, tint, uint, dups, eqf, uneq, tm, dtm, tm2,
	 u, f, m1df, dfdu, dfdp, itm, ial, ubc0, ubc1, m1bc, dbc, uicd, ficd, 
	m1ic, dicd, ir, ic, iwbrbd, p0, p1, poin, ev, wkev, ndim2, smat, 
	rnllv)
/* Subroutine */ int (*funi) (), (*bcni) (), (*icni) (), (*stpnt) (), (*
	fnbpbv) ();
integer *ibr, *m1aa, *m2aa;
doublereal *aa;
integer *m1bb, *m2bb;
doublereal *bb;
integer *m1cc;
doublereal *cc;
integer *m1dd;
doublereal *dd, *wbrbd;
integer *m1u;
doublereal *ups, *uoldps, *upoldp, *udotps, *rhsa, *rhsd, *tint, *uint, *dups,
	 *eqf, *uneq, *tm, *dtm, *tm2, *u, *f;
integer *m1df;
doublereal *dfdu, *dfdp;
integer *itm, *ial;
doublereal *ubc0, *ubc1;
integer *m1bc;
doublereal *dbc, *uicd, *ficd;
integer *m1ic;
doublereal *dicd;
integer *ir, *ic, *iwbrbd;
doublereal *p0, *p1, *poin;
doublecomplex *ev;
doublereal *wkev;
integer *ndim2;
doublereal *smat, *rnllv;
{
    /* System generated locals */
    integer aa_dim1, aa_dim2, aa_offset, bb_dim1, bb_dim2, bb_offset, cc_dim1,
	     cc_offset, dd_dim1, dd_offset, ups_dim1, ups_offset, uoldps_dim1,
	     uoldps_offset, upoldp_dim1, upoldp_offset, udotps_dim1, 
	    udotps_offset, rhsa_dim1, rhsa_offset, uint_dim1, uint_offset, 
	    dups_dim1, dups_offset, dfdu_dim1, dfdu_offset, dfdp_dim1, 
	    dfdp_offset, dbc_dim1, dbc_offset, dicd_dim1, dicd_offset, 
	    p0_dim1, p0_offset, p1_dim1, p1_offset, poin_dim1, poin_offset, 
	    smat_dim1, smat_offset, i__1, i__2;

    /* Local variables */
    extern /* Subroutine */ int sthd_();
    static integer ntot, i, k;
    extern /* Subroutine */ int adapt_();
    static integer nodir, nitps, istop;
    extern /* Subroutine */ int adptds_();
    extern /* Subroutine */ doublereal fnlpbv_();
    extern /* Subroutine */ int lcspbv_(), contbv_();
    static logical limpnt;
    extern /* Subroutine */ int stdrbv_();
    extern /* Subroutine */ doublereal fnuzbv_();
    extern /* Subroutine */ int stplbv_(), extrbv_(), solvbv_(), rsptbv_(), 
	    tpspbv_();
    static doublereal sp1;
    static integer lab;
    static doublereal rds, rlp;
    static integer itp;
    static doublereal uzr[20];


/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Controls the computation of solution branches. */




/* SGLE COMPLEX  EV(NDIM) */


/* INITIALIZE COMPUTATION OF BRANCH */

    /* Parameter adjustments */
    --rnllv;
    smat_dim1 = *ndim2;
    smat_offset = smat_dim1 + 1;
    smat -= smat_offset;
    --wkev;
    --ev;
    poin_dim1 = blbcn_1.ndim;
    poin_offset = poin_dim1 + 1;
    poin -= poin_offset;
    p1_dim1 = blbcn_1.ndim;
    p1_offset = p1_dim1 + 1;
    p1 -= p1_offset;
    p0_dim1 = blbcn_1.ndim;
    p0_offset = p0_dim1 + 1;
    p0 -= p0_offset;
    --iwbrbd;
    --ic;
    --ir;
    dicd_dim1 = *m1ic;
    dicd_offset = dicd_dim1 + 1;
    dicd -= dicd_offset;
    --ficd;
    --uicd;
    dbc_dim1 = *m1bc;
    dbc_offset = dbc_dim1 + 1;
    dbc -= dbc_offset;
    --ubc1;
    --ubc0;
    --ial;
    --itm;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    --tm2;
    --dtm;
    --tm;
    --uneq;
    --eqf;
    dups_dim1 = *m1u;
    dups_offset = dups_dim1 + 1;
    dups -= dups_offset;
    uint_dim1 = *m1u;
    uint_offset = uint_dim1 + 1;
    uint -= uint_offset;
    --tint;
    --rhsd;
    rhsa_dim1 = *m1u;
    rhsa_offset = rhsa_dim1 + 1;
    rhsa -= rhsa_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;
    --wbrbd;
    dd_dim1 = *m1dd;
    dd_offset = dd_dim1 + 1;
    dd -= dd_offset;
    cc_dim1 = *m1cc;
    cc_offset = cc_dim1 + 1;
    cc -= cc_offset;
    bb_dim1 = *m1bb;
    bb_dim2 = *m2bb;
    bb_offset = bb_dim1 * (bb_dim2 + 1) + 1;
    bb -= bb_offset;
    aa_dim1 = *m1aa;
    aa_dim2 = *m2aa;
    aa_offset = aa_dim1 * (aa_dim2 + 1) + 1;
    aa -= aa_offset;

    /* Function Body */
    rds = bldls_1.ds;
    blcrl_1.rdsold = rds;
    if (blcde_1.isp < 0) {
	blcde_1.isp = -blcde_1.isp;
    }
    sp1 = blrcn_1.zero;
    rlp = blrcn_1.zero;
    if (bllim_1.nuzr > 0) {
	i__1 = bllim_1.nuzr;
	for (i = 1; i <= i__1; ++i) {
	    uzr[i - 1] = blrcn_1.zero;
/* L15: */
	}
    }
    nitps = 0;
    ntot = 0;
    istop = 0;
    limpnt = FALSE_;

    rsptbv_(funi, stpnt, &rds, &istop, &ntot, &lab, ibr, m1u, &ups[ups_offset]
	    , &uoldps[uoldps_offset], &udotps[udotps_offset], &upoldp[
	    upoldp_offset], &tint[1], &uint[uint_offset], &eqf[1], &uneq[1], &
	    dups[dups_offset], &tm[1], &dtm[1], &tm2[1], &itm[1], &ial[1], &u[
	    1], &f[1], m1df, &dfdu[dfdu_offset], &dfdp[dfdp_offset], &ev[1], 
	    ndim2, &smat[smat_offset], &rnllv[1], &ir[1], &ic[1], &nodir);

     setrtn_( blbcn_1.ndim,blcde_1.ntst,*m1u,ups);
    if (nodir == 1 && blcde_1.isw > 0) {
	stdrbv_(funi, bcni, icni, &rds, m1aa, m2aa, &aa[aa_offset], m1bb, 
		m2bb, &bb[bb_offset], m1cc, &cc[cc_offset], m1dd, &dd[
		dd_offset], &wbrbd[1], m1u, &ups[ups_offset], &uoldps[
		uoldps_offset], &udotps[udotps_offset], &upoldp[upoldp_offset]
		, &rhsa[rhsa_offset], &rhsd[1], &dups[dups_offset], &dtm[1], &
		u[1], &f[1], m1df, &dfdu[dfdu_offset], &dfdp[dfdp_offset], &
		ubc0[1], &ubc1[1], m1bc, &dbc[dbc_offset], &uicd[1], &ficd[1],
		 m1ic, &dicd[dicd_offset], &ir[1], &ic[1], &iwbrbd[1], &c__0);

    } else if (blbcn_1.irs != 0 && blcde_1.isw < 0) {
	stdrbv_(funi, bcni, icni, &rds, m1aa, m2aa, &aa[aa_offset], m1bb, 
		m2bb, &bb[bb_offset], m1cc, &cc[cc_offset], m1dd, &dd[
		dd_offset], &wbrbd[1], m1u, &ups[ups_offset], &uoldps[
		uoldps_offset], &udotps[udotps_offset], &upoldp[upoldp_offset]
		, &rhsa[rhsa_offset], &rhsd[1], &dups[dups_offset], &dtm[1], &
		u[1], &f[1], m1df, &dfdu[dfdu_offset], &dfdp[dfdp_offset], &
		ubc0[1], &ubc1[1], m1bc, &dbc[dbc_offset], &uicd[1], &ficd[1],
		 m1ic, &dicd[dicd_offset], &ir[1], &ic[1], &iwbrbd[1], &c__1);

    }

/* Store plotting data for restart point : */

    sthd_();
    if (blbcn_1.irs == 0) {
	itp = blitp_1.itpst * 10 + 9;
    } else {
	itp = 0;
    }
    istop = 0;
    stplbv_(&istop, &itp, &c__0, &ntot, &lab, ibr, m1u, &ups[ups_offset], &
	    udotps[udotps_offset], &tm[1], &dtm[1], m1df);
    if (istop == 1) {
	goto L7;
    }

    extrbv_(funi, &rds, m1u, &ups[ups_offset], &uoldps[uoldps_offset], &
	    udotps[udotps_offset], &upoldp[upoldp_offset], &u[1], &f[1], m1df,
	     &dfdu[dfdu_offset], &dfdp[dfdp_offset], &dtm[1]);

    itp = 0;
    goto L4;

L1:
    itp = 0;

/* Adapt the mesh to the solution. */

    if (blcde_1.iad == 0) {
	goto L2;
    }
    if (ntot % blcde_1.iad == 0) {
	adapt_(&blcde_1.ntst, &blcde_1.ncol, &blcde_1.ntst, &blcde_1.ncol, &
		tm[1], &dtm[1], m1u, &ups[ups_offset], &uoldps[uoldps_offset],
		 &tint[1], &uint[uint_offset], &eqf[1], &uneq[1], &dups[
		dups_offset], &tm2[1], &itm[1], &ial[1]);
    }
L2:

/* Adapt the stepsize along the branch. */

    if (bldls_1.iads == 0) {
	goto L3;
    }
    if (ntot % bldls_1.iads == 0) {
	adptds_(&rds, &nitps);
    }

/* Provide initial approximation and determine next point. */

L3:
    contbv_(funi, &rds, m1u, &ups[ups_offset], &uoldps[uoldps_offset], &
	    udotps[udotps_offset], &upoldp[upoldp_offset], &u[1], &f[1], m1df,
	     &dfdu[dfdu_offset], &dfdp[dfdp_offset], &dtm[1]);
L4:
    solvbv_(funi, bcni, icni, &istop, &rds, &nitps, ibr, &ntot, m1aa, m2aa, &
	    aa[aa_offset], m1bb, m2bb, &bb[bb_offset], m1cc, &cc[cc_offset], 
	    m1dd, &dd[dd_offset], &wbrbd[1], m1u, &ups[ups_offset], &uoldps[
	    uoldps_offset], &udotps[udotps_offset], &upoldp[upoldp_offset], &
	    rhsa[rhsa_offset], &rhsd[1], &dups[dups_offset], &tm[1], &dtm[1], 
	    &u[1], &f[1], m1df, &dfdu[dfdu_offset], &dfdp[dfdp_offset], &ubc0[
	    1], &ubc1[1], m1bc, &dbc[dbc_offset], &uicd[1], &ficd[1], m1ic, &
	    dicd[dicd_offset], &ir[1], &ic[1], &iwbrbd[1]);
    if (istop == 1) {
	goto L6;
    }

/* Check for limit point (fold). */

    if (blbcn_1.ilp == 0) {
	goto L5;
    }
    lcspbv_(fnlpbv_, funi, bcni, icni, &istop, &itp, &rlp, &nitps, ibr, &ntot,
	     m1aa, m2aa, &aa[aa_offset], m1bb, m2bb, &bb[bb_offset], m1cc, &
	    cc[cc_offset], m1dd, &dd[dd_offset], &wbrbd[1], m1u, &ups[
	    ups_offset], &uoldps[uoldps_offset], &udotps[udotps_offset], &
	    upoldp[upoldp_offset], &rhsa[rhsa_offset], &rhsd[1], &dups[
	    dups_offset], &tm[1], &dtm[1], &u[1], &f[1], m1df, &dfdu[
	    dfdu_offset], &dfdp[dfdp_offset], &ubc0[1], &ubc1[1], m1bc, &dbc[
	    dbc_offset], &uicd[1], &ficd[1], m1ic, &dicd[dicd_offset], &ir[1],
	     &ic[1], &iwbrbd[1], &p0[p0_offset], &p1[p1_offset], &poin[
	    poin_offset], &ev[1], &wkev[1]);
    if (istop == 1) {
	goto L6;
    }

    if (itp == -1) {
	itp = blitp_1.itpst * 10 + 5;
	rlp = blrcn_1.zero;
	sp1 = blrcn_1.zero;
	limpnt = TRUE_;
    } else {
	limpnt = FALSE_;
    }

/* Check for bifurcation. */

L5:
    if (blcde_1.isp == 0 || blcde_1.isp == 1 && blbcn_1.ips == 4) {
	goto L55;
    }

    lcspbv_(fnbpbv, funi, bcni, icni, &istop, &itp, &sp1, &nitps, ibr, &ntot, 
	    m1aa, m2aa, &aa[aa_offset], m1bb, m2bb, &bb[bb_offset], m1cc, &cc[
	    cc_offset], m1dd, &dd[dd_offset], &wbrbd[1], m1u, &ups[ups_offset]
	    , &uoldps[uoldps_offset], &udotps[udotps_offset], &upoldp[
	    upoldp_offset], &rhsa[rhsa_offset], &rhsd[1], &dups[dups_offset], 
	    &tm[1], &dtm[1], &u[1], &f[1], m1df, &dfdu[dfdu_offset], &dfdp[
	    dfdp_offset], &ubc0[1], &ubc1[1], m1bc, &dbc[dbc_offset], &uicd[1]
	    , &ficd[1], m1ic, &dicd[dicd_offset], &ir[1], &ic[1], &iwbrbd[1], 
	    &p0[p0_offset], &p1[p1_offset], &poin[poin_offset], &ev[1], &wkev[
	    1]);
    if (istop == 1) {
	goto L6;
    }

    if (limpnt) {
	if (blbcn_1.ips == 2 || blbcn_1.ips == 3) {
	    sp1 = 0.;
	}
	limpnt = FALSE_;
    }

    if (itp == -1) {
	if (blbcn_1.ips != 2 && blbcn_1.ips != 3 && blbcn_1.ips != 12 && 
		blbcn_1.ips != 13 && blbcn_1.ips != 6) {
/*          **BIFURCATION POINT */
	    itp = blitp_1.itpst * 10 + 6;
	} else {
/*          **SECONDARY PERIODIC BIFURCATION : DETERMINE TYPE */
	    tpspbv_(&ev[1], &itp);
	}
	rlp = blrcn_1.zero;
	sp1 = blrcn_1.zero;
    }

/* Check for zero(es) of user supplied function(s) USZR. */

L55:
    if (bllim_1.nuzr <= 0) {
	goto L6;
    }

    i__1 = bllim_1.nuzr;
    for (i = 1; i <= i__1; ++i) {
	blusz_1.iuzr = i;
	lcspbv_(fnuzbv_, funi, bcni, icni, &istop, &itp, &uzr[i - 1], &nitps, 
		ibr, &ntot, m1aa, m2aa, &aa[aa_offset], m1bb, m2bb, &bb[
		bb_offset], m1cc, &cc[cc_offset], m1dd, &dd[dd_offset], &
		wbrbd[1], m1u, &ups[ups_offset], &uoldps[uoldps_offset], &
		udotps[udotps_offset], &upoldp[upoldp_offset], &rhsa[
		rhsa_offset], &rhsd[1], &dups[dups_offset], &tm[1], &dtm[1], &
		u[1], &f[1], m1df, &dfdu[dfdu_offset], &dfdp[dfdp_offset], &
		ubc0[1], &ubc1[1], m1bc, &dbc[dbc_offset], &uicd[1], &ficd[1],
		 m1ic, &dicd[dicd_offset], &ir[1], &ic[1], &iwbrbd[1], &p0[
		p0_offset], &p1[p1_offset], &poin[poin_offset], &ev[1], &wkev[
		1]);
	if (istop == 1) {
	    goto L6;
	}

	if (itp == -1) {
	    itp = -4 - blitp_1.itpst * 10;
	    i__2 = bllim_1.nuzr;
	    for (k = 1; k <= i__2; ++k) {
		uzr[k - 1] = blrcn_1.zero;
/* L56: */
	    }
	    goto L6;
	}
/* L57: */
    }

/* Store plotting data. */

L6:
    stplbv_(&istop, &itp, &nitps, &ntot, &lab, ibr, m1u, &ups[ups_offset], &
	    udotps[udotps_offset], &tm[1], &dtm[1], m1df);
    if (istop == 1) {
	goto L7;
    }

/* COMPUTE THE NEXT POINT ON THE BRANCH */

    if (istop == 0) {
	goto L1;
    }

L7:
    return 0;
} /* cnrlbv_ */


/* Subroutine */ int setrtn_(n, ntst, ndx, ups)
integer n, ntst, ndx;
doublereal *ups;
{
    /* System generated locals */
    integer ups_dim1, ups_offset, i__1,j;
    doublereal d__1;

    /* Builtin functions */
    integer i_dnnt();

    /* Local variables */
    static integer i__;




    ups_dim1 = ndx;
    ups_offset = ups_dim1 + 1;


    /* Function Body */
    blrtn.irot = 0;
    /* for(j=0;j<=ntst;j++)
      printf("%g %g %g \n",ups[j+ups_offset],ups[j+ups_dim1+ups_offset],
      ups[j+ups_offset+2*ups_dim1] ); */
    for (i__ = 0; i__ < n; ++i__) {

	d__1 = (ups[ntst+ i__ * ups_dim1+ups_offset] - ups[i__ * ups_dim1+ups_offset]) / 
		blrtn.torper;
	blrtn.nrot[i__ ] = i_dnnt(&d__1);
	/* printf(" i=%d nrot=%d  tp=%g \n",i__,blrtn.nrot[i__],blrtn.torper); 
         */
	if (blrtn.nrot[i__] != 0) {
	    blrtn.irot = 1;
	}
    }

    return 0;
} /* setrtn_ */



/*     ---------- ------ */
/* Subroutine */ int contbv_(funi, rds, m1u, ups, uoldps, udotps, upoldp, u, 
	f, m1df, dfdu, dfdp, dtm)
/* Subroutine */ int (*funi) ();
doublereal *rds;
integer *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *u, *f;
integer *m1df;
doublereal *dfdu, *dfdp, *dtm;
{
    /* System generated locals */
    integer ups_dim1, ups_offset, udotps_dim1, udotps_offset, upoldp_dim1, 
	    upoldp_offset, uoldps_dim1, uoldps_offset, dfdu_dim1, dfdu_offset,
	     dfdp_dim1, dfdp_offset, i__1, i__2;

    /* Local variables */
    static integer i, j;
    extern /* Subroutine */ int scalebb_(), extrbv_(), stupbv_();
    static doublereal dds;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Determines an initial approximation to the next solution point, */
/* by extrapolating from the two preceding points. */
/* The stepsize used in the preceding step has been stored in RDSOLD. */




/* Compute rate of change (along branch) of PAR(ICP(1)) and U : */

    /* Parameter adjustments */
    --dtm;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;

    /* Function Body */
    dds = blrcn_1.one / blcrl_1.rdsold;
    i__1 = blcde_1.ntst + 1;
    for (j = 1; j <= i__1; ++j) {
	i__2 = blicn_1.nrow;
	for (i = 1; i <= i__2; ++i) {
	    udotps[j + i * udotps_dim1] = (ups[j + i * ups_dim1] - uoldps[j + 
		    i * uoldps_dim1]) * dds;
/* L1: */
	}
/* L2: */
    }
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	blcrl_1.rldot[i - 1] = (blcrl_1.rl[i - 1] - blcrl_1.rlold[i - 1]) * 
		dds;
/* L3: */
    }
/*        Rescale, to set the norm of (UDOTPS,RLDOT) equal to 1. */
    scalebb_(m1u, &udotps[udotps_offset], blcrl_1.rldot, &dtm[1]);

/* Extrapolate to get initial approximation to next solution point. */

    extrbv_(funi, rds, m1u, &ups[ups_offset], &uoldps[uoldps_offset], &udotps[
	    udotps_offset], &upoldp[upoldp_offset], &u[1], &f[1], m1df, &dfdu[
	    dfdu_offset], &dfdp[dfdp_offset], &dtm[1]);

/* Store time-derivative. */

    stupbv_(funi, m1u, &ups[ups_offset], &uoldps[uoldps_offset], &upoldp[
	    upoldp_offset], &u[1], &f[1], m1df, &dfdu[dfdu_offset], &dfdp[
	    dfdp_offset]);

    return 0;
} /* contbv_ */


/*     ---------- ------ */
/* Subroutine */ int extrbv_(funi, rds, m1u, ups, uoldps, udotps, upoldp, u, 
	f, m1df, dfdu, dfdp, dtm)
/* Subroutine */ int (*funi) ();
doublereal *rds;
integer *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *u, *f;
integer *m1df;
doublereal *dfdu, *dfdp, *dtm;
{
    /* System generated locals */
    integer ups_dim1, ups_offset, udotps_dim1, udotps_offset, upoldp_dim1, 
	    upoldp_offset, uoldps_dim1, uoldps_offset, dfdu_dim1, dfdu_offset,
	     dfdp_dim1, dfdp_offset, i__1, i__2;

    /* Local variables */
    static integer i, j;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Determines an initial approximation to the next solution by */
/* extrapolating from the two preceding points. */
/* The stepsize used in the preceding step has been stored in RDSOLD. */




    /* Parameter adjustments */
    --dtm;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;

    /* Function Body */
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	blcrl_1.rlold[i - 1] = blcrl_1.rl[i - 1];
	blcrl_1.rl[i - 1] += *rds * blcrl_1.rldot[i - 1];
/* L1: */
    }
    i__1 = blcde_1.ntst + 1;
    for (j = 1; j <= i__1; ++j) {
	i__2 = blicn_1.nrow;
	for (i = 1; i <= i__2; ++i) {
	    uoldps[j + i * uoldps_dim1] = ups[j + i * ups_dim1];
	    ups[j + i * ups_dim1] += *rds * udotps[j + i * udotps_dim1];
/* L2: */
	}
/* L3: */
    }

    return 0;
} /* extrbv_ */


/*     ---------- ------ */
/* Subroutine */ int stupbv_(funi, m1u, ups, uoldps, upoldp, u, f, m1df, dfdu,
	 dfdp)
/* Subroutine */ int (*funi) ();
integer *m1u;
doublereal *ups, *uoldps, *upoldp, *u, *f;
integer *m1df;
doublereal *dfdu, *dfdp;
{
    /* System generated locals */
    integer ups_dim1, ups_offset, uoldps_dim1, uoldps_offset, upoldp_dim1, 
	    upoldp_offset, dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, 
	    i__1, i__2, i__3;

    /* Local variables */
    static integer i, j, k, n1, nc1;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Stores U-prime (derivative with respect to T) in UPOLDP. */



/* -----------------------------------------------------------------------
 */
/* The COMMON blocks BLWIF and BLDIF are used for Interface Workspace. */
/* Array dimensions in BLWIF and BLDIF are: */

/* BLWIF: */
/*       U0XX(3n+2),U1XX(3n+2),U2XX(3n+2),F1XX(3n+2),F2XX(3n+2), */
/*       DFUXX(n,2n+20),DFPXX(n,20),DDUXX(n),DDPXX(20) */
/* BLDIF: */
/*       U1ZZ(n),U2ZZ(n),F1ZZ(n),F2ZZ(n) */

/*  where n is the largest value of NDIM or NINT or NBC. */

/* -----------------------------------------------------------------------
 */


    /* Parameter adjustments */
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;

    /* Function Body */
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	blbcn_1.par[blbcn_1.icp[i - 1] - 1] = blcrl_1.rlold[i - 1];
/* L1: */
    }

    i__1 = blcde_1.ntst + 1;
    for (j = 1; j <= i__1; ++j) {
	i__2 = blbcn_1.ndim;
	for (i = 1; i <= i__2; ++i) {
	    u[i] = uoldps[j + i * uoldps_dim1];
	    if (blbcn_1.ips == 14) {
		blwif_1.u0xx[i - 1] = uoldps[j + i * uoldps_dim1] * 2 - ups[j 
			+ i * ups_dim1];
	    }
/* L2: */
	}
	(*funi)(&blbcn_1.ndim, &u[1], blwif_1.u0xx, blbcn_1.icp, blbcn_1.par, 
		&c__0, &f[1], &dfdu[dfdu_offset], &dfdp[dfdp_offset]);
	i__2 = blbcn_1.ndim;
	for (i = 1; i <= i__2; ++i) {
	    upoldp[j + i * upoldp_dim1] = f[i];
/* L3: */
	}
/* L4: */
    }

    nc1 = blcde_1.ncol - 1;
    i__1 = nc1;
    for (k = 1; k <= i__1; ++k) {
	n1 = k * blbcn_1.ndim;
	i__2 = blcde_1.ntst;
	for (j = 1; j <= i__2; ++j) {
	    i__3 = blbcn_1.ndim;
	    for (i = 1; i <= i__3; ++i) {
		u[i] = uoldps[j + (n1 + i) * uoldps_dim1];
		if (blbcn_1.ips == 14) {
		    blwif_1.u0xx[i - 1] = uoldps[j + (n1 + i) * uoldps_dim1] *
			     2 - ups[j + (n1 + i) * ups_dim1];
		}
/* L5: */
	    }
	    (*funi)(&blbcn_1.ndim, &u[1], blwif_1.u0xx, blbcn_1.icp, 
		    blbcn_1.par, &c__0, &f[1], &dfdu[dfdu_offset], &dfdp[
		    dfdp_offset]);
	    i__3 = blbcn_1.ndim;
	    for (i = 1; i <= i__3; ++i) {
		upoldp[j + (n1 + i) * upoldp_dim1] = f[i];
/* L6: */
	    }
/* L7: */
	}
/* L8: */
    }

    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	blbcn_1.par[blbcn_1.icp[i - 1] - 1] = blcrl_1.rl[i - 1];
/* L9: */
    }

    return 0;
} /* stupbv_ */


/*     ---------- ------ */
/* Subroutine */ int solvbv_(funi, bcni, icni, istop, rds, nitps, ibr, ntot, 
	m1aa, m2aa, aa, m1bb, m2bb, bb, m1cc, cc, m1dd, dd, wbrbd, m1u, ups, 
	uoldps, udotps, upoldp, rhsa, rhsd, dups, tm, dtm, u, f, m1df, dfdu, 
	dfdp, ubc0, ubc1, m1bc, dbc, uicd, ficd, m1ic, dicd, ir, ic, iwbrbd)
/* Subroutine */ int (*funi) (), (*bcni) (), (*icni) ();
integer *istop;
doublereal *rds;
integer *nitps, *ibr, *ntot, *m1aa, *m2aa;
doublereal *aa;
integer *m1bb, *m2bb;
doublereal *bb;
integer *m1cc;
doublereal *cc;
integer *m1dd;
doublereal *dd, *wbrbd;
integer *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *rhsa, *rhsd, *dups, *tm, *dtm, *
	u, *f;
integer *m1df;
doublereal *dfdu, *dfdp, *ubc0, *ubc1;
integer *m1bc;
doublereal *dbc, *uicd, *ficd;
integer *m1ic;
doublereal *dicd;
integer *ir, *ic, *iwbrbd;
{
    /* Format strings */
    static char fmt_101[] = "(\002 *** NO CONVERGENCE USING FIXED STEPSIZ\
E\002)";
    static char fmt_102[] = "(\002 STEP FAILED : TRY AGAIN WITH HALF THE STE\
PSIZE\002)";
    static char fmt_103[] = "(\002 *** NO CONVERGENCE USING MINIMUM STEPSIZ\
E\002)";

    /* System generated locals */
    integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, aa_dim1, aa_dim2, 
	    aa_offset, bb_dim1, bb_dim2, bb_offset, cc_dim1, cc_offset, 
	    dd_dim1, dd_offset, ups_dim1, ups_offset, uoldps_dim1, 
	    uoldps_offset, udotps_dim1, udotps_offset, upoldp_dim1, 
	    upoldp_offset, rhsa_dim1, rhsa_offset, dups_dim1, dups_offset, 
	    dbc_dim1, dbc_offset, dicd_dim1, dicd_offset, i__1, i__2, i__3;
    doublereal d__1, d__2;

    /* Builtin functions */
    integer s_wsfe(), e_wsfe();

    /* Local variables */
    extern /* Subroutine */ int brbd_();
    static logical done;
    static integer iibr;
    static doublereal adrl, rdrl;
    static integer ifst;
    static doublereal dumx;
    static integer i, j;
    static doublereal rdumx, au;
    extern /* Subroutine */ int wrtbv9_();
    static doublereal delmax;
    extern /* Subroutine */ int adptds_(), setrbv_(), setubv_();
    static doublereal adu;
    static integer mxt;
    static doublereal umx;
    static integer nit1;

    /* Fortran I/O blocks */
    static cilist io___38 = { 0, 9, 0, fmt_101, 0 };
    static cilist io___40 = { 0, 9, 0, fmt_102, 0 };
    static cilist io___41 = { 0, 9, 0, fmt_103, 0 };



/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Controls the solution of the nonlinear equations (by Newton's method) 
*/
/* for the next solution (PAR(ICP(*)) , U) on a branch of solutions. */




    /* Parameter adjustments */
    --iwbrbd;
    --ic;
    --ir;
    dicd_dim1 = *m1ic;
    dicd_offset = dicd_dim1 + 1;
    dicd -= dicd_offset;
    --ficd;
    --uicd;
    dbc_dim1 = *m1bc;
    dbc_offset = dbc_dim1 + 1;
    dbc -= dbc_offset;
    --ubc1;
    --ubc0;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    --dtm;
    --tm;
    dups_dim1 = *m1u;
    dups_offset = dups_dim1 + 1;
    dups -= dups_offset;
    --rhsd;
    rhsa_dim1 = *m1u;
    rhsa_offset = rhsa_dim1 + 1;
    rhsa -= rhsa_offset;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;
    --wbrbd;
    dd_dim1 = *m1dd;
    dd_offset = dd_dim1 + 1;
    dd -= dd_offset;
    cc_dim1 = *m1cc;
    cc_offset = cc_dim1 + 1;
    cc -= cc_offset;
    bb_dim1 = *m1bb;
    bb_dim2 = *m2bb;
    bb_offset = bb_dim1 * (bb_dim2 + 1) + 1;
    bb -= bb_offset;
    aa_dim1 = *m1aa;
    aa_dim2 = *m2aa;
    aa_offset = aa_dim1 * (aa_dim2 + 1) + 1;
    aa -= aa_offset;

    /* Function Body */
L1:
    blcrl_1.rdsold = *rds;
    *nitps = 0;

/* Write additional output on unit 9 if requested. */

    wrtbv9_(nitps, ibr, ntot, m1u, &ups[ups_offset], &tm[1], &dtm[1]);

/* Generate the Jacobian matrix and the right hand side. */

    i__1 = blmax_1.itnw;
    for (nit1 = 1; nit1 <= i__1; ++nit1) {

	*nitps = nit1;
	ifst = 0;
	if (*nitps <= blmax_1.nwtn) {
	    ifst = 1;
	}

	if (ifst == 1) {
	    setubv_(funi, bcni, icni, rds, m1aa, m2aa, &aa[aa_offset], m1bb, 
		    m2bb, &bb[bb_offset], m1cc, &cc[cc_offset], m1dd, &dd[
		    dd_offset], m1u, &ups[ups_offset], &uoldps[uoldps_offset],
		     &udotps[udotps_offset], &upoldp[upoldp_offset], &rhsa[
		    rhsa_offset], &rhsd[1], &dups[dups_offset], &dtm[1], &u[1]
		    , &f[1], m1df, &dfdu[dfdu_offset], &dfdp[dfdp_offset], &
		    ubc0[1], &ubc1[1], m1bc, &dbc[dbc_offset], &uicd[1], &
		    ficd[1], m1ic, &dicd[dicd_offset]);

/* Generate the right hand side only. */

	} else {

	    setrbv_(funi, bcni, icni, rds, m1u, &ups[ups_offset], &uoldps[
		    uoldps_offset], &udotps[udotps_offset], &upoldp[
		    upoldp_offset], &rhsa[rhsa_offset], &rhsd[1], &dups[
		    dups_offset], &dtm[1], &u[1], &f[1], m1df, &dfdu[
		    dfdu_offset], &dfdp[dfdp_offset], &ubc0[1], &ubc1[1], 
		    m1bc, &dbc[dbc_offset], &uicd[1], &ficd[1], m1ic, &dicd[
		    dicd_offset]);
	}

/* Solve the linearized system. */

	if (blmax_1.iid < 4) {
	    iibr = 0;
	} else if (blmax_1.iid == 4) {
	    iibr = 1;
	} else {
	    iibr = 2;
	}
	brbd_(&blcde_1.ntst, &blicn_1.nrow, &blicn_1.nclm, m1aa, m2aa, &aa[
		aa_offset], &blicn_1.nfpar, m1bb, m2bb, &bb[bb_offset], &
		blicn_1.nrc, m1cc, &cc[cc_offset], m1dd, &dd[dd_offset], m1u, 
		&rhsa[rhsa_offset], &rhsd[1], &wbrbd[1], &iibr, &ifst, &ir[1],
		 &ic[1], &iwbrbd[1], &c__0);

/* Add Newton increments. */

	i__2 = blbcn_1.ndim;
	for (i = 1; i <= i__2; ++i) {
	    ups[blcde_1.ntst + 1 + i * ups_dim1] += rhsd[i];
/* L2: */
	}
	i__2 = blicn_1.nfpar;
	for (i = 1; i <= i__2; ++i) {
	    blcrl_1.rl[i - 1] += rhsd[blbcn_1.ndim + i];
	    blbcn_1.par[blbcn_1.icp[i - 1] - 1] = blcrl_1.rl[i - 1];
/* L3: */
	}

	dumx = blrcn_1.zero;
	umx = blrcn_1.zero;
	i__2 = blcde_1.ntst;
	for (j = 1; j <= i__2; ++j) {
	    i__3 = blicn_1.nrow;
	    for (i = 1; i <= i__3; ++i) {
		adu = (d__1 = rhsa[j + i * rhsa_dim1], abs(d__1));
/* SGLE        ADU= ABS(RHSA(J,I)) */
		if (adu > dumx) {
		    dumx = adu;
		}
		au = (d__1 = ups[j + i * ups_dim1], abs(d__1));
/* SGLE        AU= ABS(UPS(J,I)) */
		if (au > umx) {
		    umx = au;
		}
		ups[j + i * ups_dim1] += rhsa[j + i * rhsa_dim1];
/* L4: */
	    }
/* L5: */
	}

	wrtbv9_(nitps, ibr, ntot, m1u, &ups[ups_offset], &tm[1], &dtm[1]);

/* Check whether user-supplied error tolerances have been met : */

	done = TRUE_;
	rdrl = blrcn_1.zero;
	i__2 = blicn_1.nfpar;
	for (i = 1; i <= i__2; ++i) {
	    adrl = (d__1 = rhsd[blbcn_1.ndim + i], abs(d__1)) / (blrcn_1.one 
		    + (d__2 = blcrl_1.rl[i - 1], abs(d__2)));
/* SGLE      ADRL= ABS(RHSD(NDIM+I))/(ONE+ ABS(RL(I))) */
	    if (adrl > bleps_1.epsl[i - 1]) {
		done = FALSE_;
	    }
	    if (adrl > rdrl) {
		rdrl = adrl;
	    }
/* L6: */
	}
	rdumx = dumx / (blrcn_1.one + umx);
	if (done && rdumx < bleps_1.epsu) {
	    return 0;
	}

	if (*nitps == 1) {
	    blref_1.delref = max(rdrl,rdumx) * 20;
/* SGLE      DELREF=20*AMAX1(RDRL,RDUMX) */
	} else {
	    delmax = max(rdrl,rdumx);
/* SGLE      DELMAX=AMAX1(RDRL,RDUMX) */
	    if (delmax > blref_1.delref) {
		goto L8;
	    }
	}

/* L7: */
    }

/* Maximum number of iterations reached. */

L8:

    if (bldls_1.iads == 0) {
	s_wsfe(&io___38);
	e_wsfe();
    }
    if (bldls_1.iads == 0) {
	goto L13;
    }

/* Reduce stepsize and try again. */

    mxt = blmax_1.itnw;
    adptds_(rds, &mxt);
    if (abs(*rds) < bldls_1.dsmin) {
	goto L12;
    }
/* SGLE  IF( ABS(RDS).LT.DSMIN)GOTO 12 */
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	blcrl_1.rl[i - 1] = blcrl_1.rlold[i - 1] + *rds * blcrl_1.rldot[i - 1]
		;
/* L9: */
    }
    i__1 = blcde_1.ntst + 1;
    for (j = 1; j <= i__1; ++j) {
	i__2 = blicn_1.nrow;
	for (i = 1; i <= i__2; ++i) {
	    ups[j + i * ups_dim1] = uoldps[j + i * uoldps_dim1] + *rds * 
		    udotps[j + i * udotps_dim1];
/* L10: */
	}
/* L11: */
    }
    if (blmax_1.iid >= 2) {
	s_wsfe(&io___40);
	e_wsfe();
    }
    goto L1;

/* Minimum stepsize reached. */

L12:
    s_wsfe(&io___41);
    e_wsfe();
L13:
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	blcrl_1.rl[i - 1] = blcrl_1.rlold[i - 1];
	blbcn_1.par[blbcn_1.icp[i - 1] - 1] = blcrl_1.rl[i - 1];
/* L14: */
    }
    i__1 = blcde_1.ntst + 1;
    for (j = 1; j <= i__1; ++j) {
	i__2 = blicn_1.nrow;
	for (i = 1; i <= i__2; ++i) {
	    ups[j + i * ups_dim1] = uoldps[j + i * uoldps_dim1];
/* L15: */
	}
/* L16: */
    }
    *istop = 1;


    return 0;
} /* solvbv_ */


/*     ---------- ------ */
/* Subroutine */ int setubv_(funi, bcni, icni, rds, m1aa, m2aa, aa, m1bb, 
	m2bb, bb, m1cc, cc, m1dd, dd, m1u, ups, uoldps, udotps, upoldp, rhsa, 
	rhsd, dups, dtm, u, f, m1df, dfdu, dfdp, ubc0, ubc1, m1bc, dbc, uicd, 
	ficd, m1ic, dicd)
/* Subroutine */ int (*funi) (), (*bcni) (), (*icni) ();
doublereal *rds;
integer *m1aa, *m2aa;
doublereal *aa;
integer *m1bb, *m2bb;
doublereal *bb;
integer *m1cc;
doublereal *cc;
integer *m1dd;
doublereal *dd;
integer *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *rhsa, *rhsd, *dups, *dtm, *u, *f;

integer *m1df;
doublereal *dfdu, *dfdp, *ubc0, *ubc1;
integer *m1bc;
doublereal *dbc, *uicd, *ficd;
integer *m1ic;
doublereal *dicd;
{
    /* System generated locals */
    integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, aa_dim1, aa_dim2, 
	    aa_offset, bb_dim1, bb_dim2, bb_offset, cc_dim1, cc_offset, 
	    dd_dim1, dd_offset, ups_dim1, ups_offset, uoldps_dim1, 
	    uoldps_offset, udotps_dim1, udotps_offset, upoldp_dim1, 
	    upoldp_offset, rhsa_dim1, rhsa_offset, dups_dim1, dups_offset, 
	    dbc_dim1, dbc_offset, dicd_dim1, dicd_offset, i__1, i__2, i__3, 
	    i__4, i__5;
    doublereal d__1;

    /* Local variables */
    static doublereal time0, time1;
    static integer i, j, k, l, m;
    static doublereal wploc[56]	/* was [8][7] */;
    static integer i1, j1, l1, k1;
    extern doublereal rinpr_();
    static doublereal rlsum;
    extern /* Subroutine */ int autim0_(), autim1_();
    static integer ib, ic;
    static doublereal sc, dt;
    static integer ir, ib1, ic1, jp1;
    static doublereal ddt;
    static integer ncp1;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/*                               ! AA  BB !             ! H1 ! */
/* Generate the Jacobian matrix  !        ! and the rhs !    ! . */
/*                               ! CC  DD !             ! H2 ! */




/* Set initial time (for timing of this subroutine). */

    /* Parameter adjustments */
    dicd_dim1 = *m1ic;
    dicd_offset = dicd_dim1 + 1;
    dicd -= dicd_offset;
    --ficd;
    --uicd;
    dbc_dim1 = *m1bc;
    dbc_offset = dbc_dim1 + 1;
    dbc -= dbc_offset;
    --ubc1;
    --ubc0;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    --dtm;
    dups_dim1 = *m1u;
    dups_offset = dups_dim1 + 1;
    dups -= dups_offset;
    --rhsd;
    rhsa_dim1 = *m1u;
    rhsa_offset = rhsa_dim1 + 1;
    rhsa -= rhsa_offset;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;
    dd_dim1 = *m1dd;
    dd_offset = dd_dim1 + 1;
    dd -= dd_offset;
    cc_dim1 = *m1cc;
    cc_offset = cc_dim1 + 1;
    cc -= cc_offset;
    bb_dim1 = *m1bb;
    bb_dim2 = *m2bb;
    bb_offset = bb_dim1 * (bb_dim2 + 1) + 1;
    bb -= bb_offset;
    aa_dim1 = *m1aa;
    aa_dim2 = *m2aa;
    aa_offset = aa_dim1 * (aa_dim2 + 1) + 1;
    aa -= aa_offset;

    /* Function Body */
    autim0_(&time0);

/* Initialize to zero. */

    i__1 = blicn_1.nrc;
    for (i = 1; i <= i__1; ++i) {
	rhsd[i] = blrcn_1.zero;
	i__2 = blicn_1.nfpar;
	for (k = 1; k <= i__2; ++k) {
	    dd[i + k * dd_dim1] = blrcn_1.zero;
/* L1: */
	}
/* L2: */
    }

    i__1 = blcde_1.ntst;
    for (j = 1; j <= i__1; ++j) {
	i__2 = blicn_1.nrow;
	for (i = 1; i <= i__2; ++i) {
	    i__3 = blicn_1.nfpar;
	    for (k = 1; k <= i__3; ++k) {
		bb[j + (i + k * bb_dim2) * bb_dim1] = blrcn_1.zero;
/* L3: */
	    }
/* L4: */
	}
/* L5: */
    }

    i__1 = blicn_1.nrc;
    for (i = 1; i <= i__1; ++i) {
	i__2 = blicn_1.ncc;
	for (j = 1; j <= i__2; ++j) {
	    cc[j + i * cc_dim1] = blrcn_1.zero;
/* L6: */
	}
/* L7: */
    }

    i__1 = blicn_1.nclm;
    for (ic = 1; ic <= i__1; ++ic) {
	i__2 = blicn_1.nrow;
	for (ir = 1; ir <= i__2; ++ir) {
	    i__3 = blcde_1.ntst;
	    for (j = 1; j <= i__3; ++j) {
		aa[j + (ir + ic * aa_dim2) * aa_dim1] = blrcn_1.zero;
/* L8: */
	    }
/* L9: */
	}
/* L10: */
    }

/* Set constants. */

    ncp1 = blcde_1.ncol + 1;

/* Generate AA , BB and H1 : */

    i__1 = blcde_1.ntst;
    for (j = 1; j <= i__1; ++j) {
	jp1 = j + 1;
	dt = dtm[j];
	ddt = blrcn_1.one / dt;
	i__2 = blcde_1.ncol;
	for (ic = 1; ic <= i__2; ++ic) {
	    i__3 = ncp1;
	    for (ib = 1; ib <= i__3; ++ib) {
		wploc[ib + (ic << 3) - 9] = ddt * blwts_1.wp[ib + (ic << 3) - 
			9];
/* L11: */
	    }
/* L12: */
	}
	i__2 = blcde_1.ncol;
	for (ic = 1; ic <= i__2; ++ic) {
	    i__3 = blbcn_1.ndim;
	    for (k = 1; k <= i__3; ++k) {
		u[k] = blwts_1.w[ncp1 + (ic << 3) - 9] * ups[jp1 + k * 
			ups_dim1];
		i__4 = blcde_1.ncol;
		for (l = 1; l <= i__4; ++l) {
		    l1 = (l - 1) * blbcn_1.ndim + k;
		    u[k] += blwts_1.w[l + (ic << 3) - 9] * ups[j + l1 * 
			    ups_dim1];
/* L13: */
		}
/* L14: */
	    }
	    if (blbcn_1.ips == 14) {
/*          ** Time evolution computations (parabolic systems)
 */
		i__3 = blbcn_1.ndim;
		for (k = 1; k <= i__3; ++k) {
		    blwif_1.u0xx[k - 1] = blwts_1.w[ncp1 + (ic << 3) - 9] * 
			    uoldps[jp1 + k * uoldps_dim1];
		    i__4 = blcde_1.ncol;
		    for (l = 1; l <= i__4; ++l) {
			l1 = (l - 1) * blbcn_1.ndim + k;
			blwif_1.u0xx[k - 1] += blwts_1.w[l + (ic << 3) - 9] * 
				uoldps[j + l1 * uoldps_dim1];
/* L131: */
		    }
/* L141: */
		}
	    }
	    i__3 = blicn_1.nfpar;
	    for (i = 1; i <= i__3; ++i) {
		blbcn_1.par[blbcn_1.icp[i - 1] - 1] = blcrl_1.rl[i - 1];
/* L15: */
	    }
	    (*funi)(&blbcn_1.ndim, &u[1], blwif_1.u0xx, blbcn_1.icp, 
		    blbcn_1.par, &c__1, &f[1], &dfdu[dfdu_offset], &dfdp[
		    dfdp_offset]);
	    ic1 = (ic - 1) * blbcn_1.ndim;
	    i__3 = ncp1;
	    for (ib = 1; ib <= i__3; ++ib) {
		ib1 = (ib - 1) * blbcn_1.ndim;
		i__4 = blbcn_1.ndim;
		for (i = 1; i <= i__4; ++i) {
		    aa[j + (ic1 + i + (ib1 + i) * aa_dim2) * aa_dim1] = wploc[
			    ib + (ic << 3) - 9];
		    i__5 = blbcn_1.ndim;
		    for (k = 1; k <= i__5; ++k) {
			aa[j + (ic1 + i + (ib1 + k) * aa_dim2) * aa_dim1] -= 
				blwts_1.w[ib + (ic << 3) - 9] * dfdu[i + k * 
				dfdu_dim1];
/* L16: */
		    }
/* L17: */
		}
/* L18: */
	    }
	    i__3 = blbcn_1.ndim;
	    for (i = 1; i <= i__3; ++i) {
		i__4 = blicn_1.nfpar;
		for (k = 1; k <= i__4; ++k) {
		    bb[j + (ic1 + i + k * bb_dim2) * bb_dim1] = -dfdp[i + 
			    blbcn_1.icp[k - 1] * dfdp_dim1];
/* L19: */
		}
		rhsa[j + (ic1 + i) * rhsa_dim1] = f[i] - wploc[ncp1 + (ic << 
			3) - 9] * ups[jp1 + i * ups_dim1];
		i__4 = blcde_1.ncol;
		for (k = 1; k <= i__4; ++k) {
		    k1 = (k - 1) * blbcn_1.ndim + i;
		    rhsa[j + (ic1 + i) * rhsa_dim1] -= wploc[k + (ic << 3) - 
			    9] * ups[j + k1 * ups_dim1];
/* L20: */
		}
/* L21: */
	    }
/* L22: */
	}
/* L23: */
    }

/* Generate CC, DD and H2 : */

/*    Boundary conditions. */

    if (blcde_1.nbc <= 0) {
	goto L30;
    }
    i__1 = blbcn_1.ndim;
    for (i = 1; i <= i__1; ++i) {
	ubc0[i] = ups[i * ups_dim1 + 1];
	ubc1[i] = ups[blcde_1.ntst + 1 + i * ups_dim1];
/* L24: */
    }

    (*bcni)(&blbcn_1.ndim, blbcn_1.par, blbcn_1.icp, &blcde_1.nbc, &ubc0[1], &
	    ubc1[1], &rhsd[1], &c__1, &dbc[dbc_offset]);

    i__1 = blcde_1.nbc;
    for (i = 1; i <= i__1; ++i) {
	rhsd[i] = -rhsd[i];
	i__2 = blbcn_1.ndim;
	for (k = 1; k <= i__2; ++k) {
	    cc[k + i * cc_dim1] = dbc[i + k * dbc_dim1];
	    cc[blicn_1.ncc - blbcn_1.ndim + k + i * cc_dim1] = dbc[i + (
		    blbcn_1.ndim + k) * dbc_dim1];
/* L25: */
	}
	i__2 = blicn_1.nfpar;
	for (k = 1; k <= i__2; ++k) {
	    dd[i + k * dd_dim1] = dbc[i + ((blbcn_1.ndim << 1) + blbcn_1.icp[
		    k - 1]) * dbc_dim1];
/* L26: */
	}
/* L27: */
    }

/*   Save difference : */

    i__1 = blcde_1.ntst + 1;
    for (j = 1; j <= i__1; ++j) {
	i__2 = blicn_1.nrow;
	for (i = 1; i <= i__2; ++i) {
	    dups[j + i * dups_dim1] = ups[j + i * ups_dim1] - uoldps[j + i * 
		    uoldps_dim1];
/* L28: */
	}
/* L29: */
    }

/*   Integral constraints : */

L30:
    if (blcde_1.nint <= 0) {
	goto L37;
    }

    i__1 = blcde_1.ntst;
    for (j = 1; j <= i__1; ++j) {
	jp1 = j + 1;
	i__2 = ncp1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = blbcn_1.ndim;
	    for (i = 1; i <= i__3; ++i) {
		i1 = (k - 1) * blbcn_1.ndim + i;
		j1 = j;
		if (k == ncp1) {
		    i1 = i;
		}
		if (k == ncp1) {
		    j1 = jp1;
		}
		uicd[i] = ups[j1 + i1 * ups_dim1];
		uicd[blbcn_1.ndim + i] = uoldps[j1 + i1 * uoldps_dim1];
		uicd[(blbcn_1.ndim << 1) + i] = udotps[j1 + i1 * udotps_dim1];

		uicd[blbcn_1.ndim * 3 + i] = upoldp[j1 + i1 * upoldp_dim1];
/* L31: */
	    }
	    (*icni)(&blbcn_1.ndim, blbcn_1.par, blbcn_1.icp, &blcde_1.nint, &
		    uicd[1], &uicd[blbcn_1.ndim + 1], &uicd[(blbcn_1.ndim << 
		    1) + 1], &uicd[blbcn_1.ndim * 3 + 1], &ficd[1], &c__1, &
		    dicd[dicd_offset]);
	    i__3 = blcde_1.nint;
	    for (m = 1; m <= i__3; ++m) {
		i__4 = blbcn_1.ndim;
		for (i = 1; i <= i__4; ++i) {
		    k1 = (k - 1) * blbcn_1.ndim + i;
		    l1 = (j - 1) * blicn_1.nrow + k1;
		    cc[l1 + (blcde_1.nbc + m) * cc_dim1] += dtm[j] * 
			    blwts_1.wi[k - 1] * dicd[m + i * dicd_dim1];
/* L32: */
		}
		i__4 = blicn_1.nfpar;
		for (i = 1; i <= i__4; ++i) {
		    dd[blcde_1.nbc + m + i * dd_dim1] += dtm[j] * blwts_1.wi[
			    k - 1] * dicd[m + (blbcn_1.ndim + blbcn_1.icp[i - 
			    1]) * dicd_dim1];
/* L33: */
		}
		rhsd[blcde_1.nbc + m] -= dtm[j] * blwts_1.wi[k - 1] * ficd[m];

/* L34: */
	    }
/* L35: */
	}
/* L36: */
    }

L37:

/*   Pseudo-arclength equation : */

    i__1 = blcde_1.ntst;
    for (j = 1; j <= i__1; ++j) {
	jp1 = j + 1;
/* Computing 2nd power */
	d__1 = bltht_1.thetau;
	sc = dtm[j] * (d__1 * d__1);
	i__2 = blbcn_1.ndim;
	for (i = 1; i <= i__2; ++i) {
	    i__3 = blcde_1.ncol;
	    for (k = 1; k <= i__3; ++k) {
		k1 = (k - 1) * blbcn_1.ndim + i;
		l1 = (j - 1) * blicn_1.nrow + k1;
		cc[l1 + blicn_1.nrc * cc_dim1] += sc * blwts_1.wi[k - 1] * 
			udotps[j + k1 * udotps_dim1];
/* L38: */
	    }
	    cc[j * blicn_1.nrow + i + blicn_1.nrc * cc_dim1] += sc * 
		    blwts_1.wi[ncp1 - 1] * udotps[jp1 + i * udotps_dim1];
/* L39: */
	}
/* L40: */
    }

    rlsum = blrcn_1.zero;
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
/* Computing 2nd power */
	d__1 = bltht_1.thetal[i - 1];
	dd[blicn_1.nrc + i * dd_dim1] = d__1 * d__1 * blcrl_1.rldot[i - 1];
/* Computing 2nd power */
	d__1 = bltht_1.thetal[i - 1];
	rlsum += d__1 * d__1 * (blcrl_1.rl[i - 1] - blcrl_1.rlold[i - 1]) * 
		blcrl_1.rldot[i - 1];
/* L41: */
    }

/* Computing 2nd power */
    d__1 = bltht_1.thetau;
    rhsd[blicn_1.nrc] = *rds - d__1 * d__1 * rinpr_(&blbcn_1.ndim, m1u, &
	    udotps[udotps_offset], &dups[dups_offset], &dtm[1]) - rlsum;

/* Determine the time spent in this subroutine. */

    autim1_(&time1);
    bltim_1.tsetub = bltim_1.tsetub + time1 - time0;

    return 0;
} /* setubv_ */


/*     ---------- ------ */
/* Subroutine */ int setrbv_(funi, bcni, icni, rds, m1u, ups, uoldps, udotps, 
	upoldp, rhsa, rhsd, dups, dtm, u, f, m1df, dfdu, dfdp, ubc0, ubc1, 
	m1bc, dbc, uicd, ficd, m1ic, dicd)
/* Subroutine */ int (*funi) (), (*bcni) (), (*icni) ();
doublereal *rds;
integer *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *rhsa, *rhsd, *dups, *dtm, *u, *f;

integer *m1df;
doublereal *dfdu, *dfdp, *ubc0, *ubc1;
integer *m1bc;
doublereal *dbc, *uicd, *ficd;
integer *m1ic;
doublereal *dicd;
{
    /* System generated locals */
    integer ups_dim1, ups_offset, uoldps_dim1, uoldps_offset, udotps_dim1, 
	    udotps_offset, upoldp_dim1, upoldp_offset, dfdu_dim1, dfdu_offset,
	     dfdp_dim1, dfdp_offset, rhsa_dim1, rhsa_offset, dups_dim1, 
	    dups_offset, dbc_dim1, dbc_offset, dicd_dim1, dicd_offset, i__1, 
	    i__2, i__3, i__4;
    doublereal d__1;

    /* Local variables */
    static integer i, j, k, l, m;
    static doublereal wploc[56]	/* was [8][7] */;
    static integer i1, j1, l1, k1;
    extern doublereal rinpr_();
    static doublereal rlsum;
    static integer ib, ic;
    static doublereal dt;
    static integer ic1, jp1;
    static doublereal ddt;
    static integer ncp1;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/*  Generate the right hand side only of the Newton equations. */




    /* Parameter adjustments */
    dicd_dim1 = *m1ic;
    dicd_offset = dicd_dim1 + 1;
    dicd -= dicd_offset;
    --ficd;
    --uicd;
    dbc_dim1 = *m1bc;
    dbc_offset = dbc_dim1 + 1;
    dbc -= dbc_offset;
    --ubc1;
    --ubc0;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    --dtm;
    dups_dim1 = *m1u;
    dups_offset = dups_dim1 + 1;
    dups -= dups_offset;
    --rhsd;
    rhsa_dim1 = *m1u;
    rhsa_offset = rhsa_dim1 + 1;
    rhsa -= rhsa_offset;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;

    /* Function Body */
    i__1 = blicn_1.nrc;
    for (i = 1; i <= i__1; ++i) {
	rhsd[i] = blrcn_1.zero;
/* L1: */
    }

/* SET CONSTANTS */

    ncp1 = blcde_1.ncol + 1;

    i__1 = blcde_1.ntst;
    for (j = 1; j <= i__1; ++j) {
	jp1 = j + 1;
	dt = dtm[j];
	ddt = blrcn_1.one / dt;
	i__2 = blcde_1.ncol;
	for (ic = 1; ic <= i__2; ++ic) {
	    i__3 = ncp1;
	    for (ib = 1; ib <= i__3; ++ib) {
		wploc[ib + (ic << 3) - 9] = ddt * blwts_1.wp[ib + (ic << 3) - 
			9];
/* L2: */
	    }
/* L3: */
	}
	i__2 = blcde_1.ncol;
	for (ic = 1; ic <= i__2; ++ic) {
	    i__3 = blbcn_1.ndim;
	    for (k = 1; k <= i__3; ++k) {
		u[k] = blwts_1.w[ncp1 + (ic << 3) - 9] * ups[jp1 + k * 
			ups_dim1];
		i__4 = blcde_1.ncol;
		for (l = 1; l <= i__4; ++l) {
		    l1 = (l - 1) * blbcn_1.ndim + k;
		    u[k] += blwts_1.w[l + (ic << 3) - 9] * ups[j + l1 * 
			    ups_dim1];
/* L4: */
		}
/* L5: */
	    }
	    if (blbcn_1.ips == 14) {
/*          ** TIME EVOLUTION COMPUTATIONS (PARABOLIC SYSTEMS)
 */
		i__3 = blbcn_1.ndim;
		for (k = 1; k <= i__3; ++k) {
		    blwif_1.u0xx[k - 1] = blwts_1.w[ncp1 + (ic << 3) - 9] * 
			    uoldps[jp1 + k * uoldps_dim1];
		    i__4 = blcde_1.ncol;
		    for (l = 1; l <= i__4; ++l) {
			l1 = (l - 1) * blbcn_1.ndim + k;
			blwif_1.u0xx[k - 1] += blwts_1.w[l + (ic << 3) - 9] * 
				uoldps[j + l1 * uoldps_dim1];
/* L41: */
		    }
/* L51: */
		}
	    }
	    i__3 = blicn_1.nfpar;
	    for (i = 1; i <= i__3; ++i) {
		blbcn_1.par[blbcn_1.icp[i - 1] - 1] = blcrl_1.rl[i - 1];
/* L6: */
	    }
	    (*funi)(&blbcn_1.ndim, &u[1], blwif_1.u0xx, blbcn_1.icp, 
		    blbcn_1.par, &c__0, &f[1], &dfdu[dfdu_offset], &dfdp[
		    dfdp_offset]);
	    ic1 = (ic - 1) * blbcn_1.ndim;
	    i__3 = blbcn_1.ndim;
	    for (i = 1; i <= i__3; ++i) {
		rhsa[j + (ic1 + i) * rhsa_dim1] = f[i] - wploc[ncp1 + (ic << 
			3) - 9] * ups[jp1 + i * ups_dim1];
		i__4 = blcde_1.ncol;
		for (k = 1; k <= i__4; ++k) {
		    k1 = (k - 1) * blbcn_1.ndim + i;
		    rhsa[j + (ic1 + i) * rhsa_dim1] -= wploc[k + (ic << 3) - 
			    9] * ups[j + k1 * ups_dim1];
/* L7: */
		}
/* L8: */
	    }
/* L9: */
	}
/* L10: */
    }

/*    Boundary conditions. */

    if (blcde_1.nbc == 0) {
	goto L15;
    }
    i__1 = blbcn_1.ndim;
    for (i = 1; i <= i__1; ++i) {
	ubc0[i] = ups[i * ups_dim1 + 1];
	ubc1[i] = ups[blcde_1.ntst + 1 + i * ups_dim1];
/* L11: */
    }

    (*bcni)(&blbcn_1.ndim, blbcn_1.par, blbcn_1.icp, &blcde_1.nbc, &ubc0[1], &
	    ubc1[1], &rhsd[1], &c__0, &dbc[dbc_offset]);

    i__1 = blcde_1.nbc;
    for (i = 1; i <= i__1; ++i) {
	rhsd[i] = -rhsd[i];
/* L12: */
    }

/*   Save difference : */

    i__1 = blcde_1.ntst + 1;
    for (j = 1; j <= i__1; ++j) {
	i__2 = blicn_1.nrow;
	for (i = 1; i <= i__2; ++i) {
	    dups[j + i * dups_dim1] = ups[j + i * ups_dim1] - uoldps[j + i * 
		    uoldps_dim1];
/* L13: */
	}
/* L14: */
    }

/*   Integral constraints : */

L15:
    if (blcde_1.nint <= 0) {
	goto L20;
    }

    i__1 = blcde_1.ntst;
    for (j = 1; j <= i__1; ++j) {
	jp1 = j + 1;
	i__2 = ncp1;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = blbcn_1.ndim;
	    for (i = 1; i <= i__3; ++i) {
		i1 = (k - 1) * blbcn_1.ndim + i;
		j1 = j;
		if (k == ncp1) {
		    i1 = i;
		}
		if (k == ncp1) {
		    j1 = jp1;
		}
		uicd[i] = ups[j1 + i1 * ups_dim1];
		uicd[blbcn_1.ndim + i] = uoldps[j1 + i1 * uoldps_dim1];
		uicd[(blbcn_1.ndim << 1) + i] = udotps[j1 + i1 * udotps_dim1];

		uicd[blbcn_1.ndim * 3 + i] = upoldp[j1 + i1 * upoldp_dim1];
/* L16: */
	    }
	    (*icni)(&blbcn_1.ndim, blbcn_1.par, blbcn_1.icp, &blcde_1.nint, &
		    uicd[1], &uicd[blbcn_1.ndim + 1], &uicd[(blbcn_1.ndim << 
		    1) + 1], &uicd[blbcn_1.ndim * 3 + 1], &ficd[1], &c__0, &
		    dicd[dicd_offset]);
	    i__3 = blcde_1.nint;
	    for (m = 1; m <= i__3; ++m) {
		rhsd[blcde_1.nbc + m] -= dtm[j] * blwts_1.wi[k - 1] * ficd[m];

/* L17: */
	    }
/* L18: */
	}
/* L19: */
    }

/*   Pseudo-arclength equation : */

L20:
    rlsum = blrcn_1.zero;
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
/* Computing 2nd power */
	d__1 = bltht_1.thetal[i - 1];
	rlsum += d__1 * d__1 * (blcrl_1.rl[i - 1] - blcrl_1.rlold[i - 1]) * 
		blcrl_1.rldot[i - 1];
/* L21: */
    }

/* Computing 2nd power */
    d__1 = bltht_1.thetau;
    rhsd[blicn_1.nrc] = *rds - d__1 * d__1 * rinpr_(&blbcn_1.ndim, m1u, &
	    udotps[udotps_offset], &dups[dups_offset], &dtm[1]) - rlsum;

    return 0;
} /* setrbv_ */


/* ------------------------------------------------------------------- */
/* ------------------------------------------------------------------- */
/*      Restart of Solution Branches ( Differential Equations ) */
/* ------------------------------------------------------------------- */
/* ------------------------------------------------------------------- */

/*     ---------- ------ */
/* Subroutine */ int rsptbv_(funi, stpnt, rds, istop, ntot, lab, ibr, m1u, 
	ups, uoldps, udotps, upoldp, tint, uint, eqf, uneq, dups, tm, dtm, 
	tm2, itm, ial, u, f, m1df, dfdu, dfdp, ev, ndim2, smat, rnllv, ir, ic,
	 nodir)
/* Subroutine */ int (*funi) (), (*stpnt) ();
doublereal *rds;
integer *istop, *ntot, *lab, *ibr, *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *tint, *uint, *eqf, *uneq, *dups, 
	*tm, *dtm, *tm2;
integer *itm, *ial;
doublereal *u, *f;
integer *m1df;
doublereal *dfdu, *dfdp;
doublecomplex *ev;
integer *ndim2;
doublereal *smat, *rnllv;
integer *ir, *ic, *nodir;
{
    /* System generated locals */
    integer ups_dim1, ups_offset, uoldps_dim1, uoldps_offset, udotps_dim1, 
	    udotps_offset, upoldp_dim1, upoldp_offset, uint_dim1, uint_offset,
	     dups_dim1, dups_offset, dfdu_dim1, dfdu_offset, dfdp_dim1, 
	    dfdp_offset, smat_dim1, smat_offset, i__1, i__2;

    /* Local variables */
    static integer i, j;
    extern /* Subroutine */ int adapt_(), newlab_();
    static integer ncolrs;
    extern /* Subroutine */ int stupbv_();
    static integer ntstrs;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Restarts computation of a branch of solutions at point labelled IRS. */

/* The output written on unit 8 by a previous run is now expected as */
/* input on unit 3. The label IRS, where computation is to resume, must */

/* be specified in the user-supplied subroutine INIT. */
/* If IRS=0 then the starting point must be provided analytically in the 
*/
/* user-supplied subroutine STPNT. */



/* SGLE COMPLEX  EV(NDIM) */


/* Get restart data : */

    /* Parameter adjustments */
    --ic;
    --ir;
    --rnllv;
    smat_dim1 = *ndim2;
    smat_offset = smat_dim1 + 1;
    smat -= smat_offset;
    --ev;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    --ial;
    --itm;
    --tm2;
    --dtm;
    --tm;
    dups_dim1 = *m1u;
    dups_offset = dups_dim1 + 1;
    dups -= dups_offset;
    --uneq;
    --eqf;
    uint_dim1 = *m1u;
    uint_offset = uint_dim1 + 1;
    uint -= uint_offset;
    --tint;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;

    /* Function Body */
    (*stpnt)(&ntstrs, &ncolrs, lab, ibr, m1u, &u[1], &ups[ups_offset], &
	    udotps[udotps_offset], &upoldp[upoldp_offset], &tm[1], &dtm[1], 
	    ndim2, &smat[smat_offset], &rnllv[1], &ir[1], &ic[1], &f[1], &
	    dfdu[dfdu_offset], &dfdp[dfdp_offset], nodir);

/* Determine a suitable starting label and branch number. */

    if (blbcn_1.irs > 0) {
	newlab_(&blcde_1.isw, ibr, lab);
    }

    i__1 = ntstrs;
    for (j = 1; j <= i__1; ++j) {
	dtm[j] = tm[j + 1] - tm[j];
/* L1: */
    }

/* Adapt mesh if necessary : */

    if (blcde_1.ntst != ntstrs || blcde_1.ncol != ncolrs) {
	adapt_(&ntstrs, &ncolrs, &blcde_1.ntst, &blcde_1.ncol, &tm[1], &dtm[1]
		, m1u, &ups[ups_offset], &udotps[udotps_offset], &tint[1], &
		uint[uint_offset], &eqf[1], &uneq[1], &dups[dups_offset], &
		tm2[1], &itm[1], &ial[1]);
    }

/* Set UOLDPS, RLOLD. */

    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	blcrl_1.rl[i - 1] = blbcn_1.par[blbcn_1.icp[i - 1] - 1];
	blcrl_1.rlold[i - 1] = blcrl_1.rl[i - 1];
/* L2: */
    }

    i__1 = blicn_1.nrow;
    for (i = 1; i <= i__1; ++i) {
	i__2 = blcde_1.ntst + 1;
	for (j = 1; j <= i__2; ++j) {
	    uoldps[j + i * uoldps_dim1] = ups[j + i * ups_dim1];
/* L3: */
	}
/* L4: */
    }

/* Store U-prime (derivative with respect to time or space variable). */

    if (*nodir == -1) {
/*        ** Restart from a Hopf bifurcation. */
	*nodir = 0;
	blcde_1.isw = 1;
    } else {
/*        ** Restart from orbit. */
	stupbv_(funi, m1u, &ups[ups_offset], &uoldps[uoldps_offset], &upoldp[
		upoldp_offset], &u[1], &f[1], m1df, &dfdu[dfdu_offset], &dfdp[
		dfdp_offset]);
    }

    return 0;
} /* rsptbv_ */


/*     ---------- ------ */
/* Subroutine */ int stpnbv_(ntstrs, ncolrs, lab, ibr, m1u, u, ups, udotps, 
	upoldp, tm, dtm, ndim2, smat, rnllv, ir, ic, f, dfdu, dfdp, nodir)
integer *ntstrs, *ncolrs, *lab, *ibr, *m1u;
doublereal *u, *ups, *udotps, *upoldp, *tm, *dtm;
integer *ndim2;
doublereal *smat, *rnllv;
integer *ir, *ic;
doublereal *f, *dfdu, *dfdp;
integer *nodir;
{
    /* Format strings */
    static char fmt_101[] = "(4x,1p7e18.10)";

    /* System generated locals */
    integer ups_dim1, ups_offset, udotps_dim1, udotps_offset, upoldp_dim1, 
	    upoldp_offset, smat_dim1, smat_offset, dfdu_dim1, dfdu_offset, 
	    dfdp_dim1, dfdp_offset, i__1, i__2, i__3;

    /* Builtin functions */
    integer s_rsle(), do_lio(), e_rsle(), s_rsfe(), do_fio(), e_rsfe();

    /* Local variables */
    static doublereal temp[7];
    extern /* Subroutine */ int skip3_();
    static integer ntpl1, nrsp1, ntot1, i, j, k;
    extern /* Subroutine */ int pdble_();
    static logical found;
    static integer icprs[20], k1, k2;
    extern /* Subroutine */ int findl3_();
    static integer nfpar1, nskip1, nskipr;
    static logical eof3;
    static integer nar1, itp1, num1, num2, isw1;

    /* Fortran I/O blocks */
    static cilist io___90 = { 0, 3, 0, 0, 0 };
    static cilist io___105 = { 0, 3, 0, fmt_101, 0 };
    static cilist io___109 = { 0, 3, 0, fmt_101, 0 };
    static cilist io___110 = { 0, 3, 0, fmt_101, 0 };
    static cilist io___111 = { 0, 3, 0, fmt_101, 0 };
    static cilist io___112 = { 0, 3, 0, fmt_101, 0 };
    static cilist io___113 = { 0, 3, 0, fmt_101, 0 };



/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* This subroutine locates and retrieves the information required to */
/* restart computation at the point with label IRS. */
/* This information is expected on unit 3. */




    /* Parameter adjustments */
    dfdp_dim1 = blbcn_1.ndim;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = blbcn_1.ndim;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --ic;
    --ir;
    --rnllv;
    smat_dim1 = *ndim2;
    smat_offset = smat_dim1 + 1;
    smat -= smat_offset;
    --dtm;
    --tm;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;
    --u;

    /* Function Body */
    findl3_(&blbcn_1.irs, &itp1, &nfpar1, &found);
    s_rsle(&io___90);
    do_lio(&c__3, &c__1, (char *)&(*ibr), (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&ntot1, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&itp1, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&(*lab), (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&nfpar1, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&isw1, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&ntpl1, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&nar1, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&nskip1, (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&(*ntstrs), (ftnlen)sizeof(integer));
    do_lio(&c__3, &c__1, (char *)&(*ncolrs), (ftnlen)sizeof(integer));
    i__1 = nfpar1;
    for (i = 1; i <= i__1; ++i) {
	do_lio(&c__3, &c__1, (char *)&icprs[i - 1], (ftnlen)sizeof(integer));
    }
    e_rsle();
    nrsp1 = *ntstrs + 1;

    num1 = blbcn_1.ndim / 7 + 1;
    num2 = nar1 / 8 + 1;
    nskipr = num2 - num1;

    i__1 = *ntstrs;
    for (j = 1; j <= i__1; ++j) {
	i__2 = *ncolrs;
	for (i = 1; i <= i__2; ++i) {
	    k1 = (i - 1) * blbcn_1.ndim + 1;
	    k2 = i * blbcn_1.ndim;
	    s_rsfe(&io___105);
	    do_fio(&c__1, (char *)&temp[i - 1], (ftnlen)sizeof(doublereal));
	    i__3 = k2;
	    for (k = k1; k <= i__3; ++k) {
		do_fio(&c__1, (char *)&ups[j + k * ups_dim1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_rsfe();
	    if (nskipr > 0) {
		skip3_(&nskipr, &eof3);
	    }
/* L1: */
	}
	tm[j] = temp[0];
/* L2: */
    }
    s_rsfe(&io___109);
    do_fio(&c__1, (char *)&tm[nrsp1], (ftnlen)sizeof(doublereal));
    i__1 = blbcn_1.ndim;
    for (k = 1; k <= i__1; ++k) {
	do_fio(&c__1, (char *)&ups[nrsp1 + k * ups_dim1], (ftnlen)sizeof(
		doublereal));
    }
    e_rsfe();
    if (nskipr > 0) {
	skip3_(&nskipr, &eof3);
    }

    s_rsfe(&io___110);
    i__1 = nfpar1;
    for (i = 1; i <= i__1; ++i) {
	do_fio(&c__1, (char *)&blcrl_1.rldot[i - 1], (ftnlen)sizeof(
		doublereal));
    }
    e_rsfe();

/* Read U-dot (deriv. with respect to arclength along solution branch). */


    num1 = blbcn_1.ndim / 8 + 1;
    num2 = nar1 / 9 + 1;
    nskipr = num2 - num1;

    i__1 = *ntstrs;
    for (j = 1; j <= i__1; ++j) {
	i__2 = *ncolrs;
	for (i = 1; i <= i__2; ++i) {
	    k1 = (i - 1) * blbcn_1.ndim + 1;
	    k2 = i * blbcn_1.ndim;
	    s_rsfe(&io___111);
	    i__3 = k2;
	    for (k = k1; k <= i__3; ++k) {
		do_fio(&c__1, (char *)&udotps[j + k * udotps_dim1], (ftnlen)
			sizeof(doublereal));
	    }
	    e_rsfe();
	    if (nskipr > 0) {
		skip3_(&nskipr, &eof3);
	    }
/* L3: */
	}
/* L4: */
    }
    s_rsfe(&io___112);
    i__1 = blbcn_1.ndim;
    for (k = 1; k <= i__1; ++k) {
	do_fio(&c__1, (char *)&udotps[nrsp1 + k * udotps_dim1], (ftnlen)
		sizeof(doublereal));
    }
    e_rsfe();
    if (nskipr > 0) {
	skip3_(&nskipr, &eof3);
    }

/* Read the parameter values. */

    s_rsfe(&io___113);
    i__1 = blicn_1.npar;
    for (i = 1; i <= i__1; ++i) {
	do_fio(&c__1, (char *)&blbcn_1.par[i - 1], (ftnlen)sizeof(doublereal))
		;
    }
    e_rsfe();
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	blcrl_1.rl[i - 1] = blbcn_1.par[blbcn_1.icp[i - 1] - 1];
/* L5: */
    }

/* Special case : Preprocess restart data in case of branch switching */
/* at a period doubling bifurcation. */

    if ((blbcn_1.ips == 2 || blbcn_1.ips == 3 || (real) blbcn_1.ips == (float)
	    6. || blbcn_1.ips == 12 || blbcn_1.ips == 13) && blcde_1.isw == 
	    -1 && itp1 == 7) {
	pdble_(&blbcn_1.ndim, ntstrs, ncolrs, m1u, &ups[ups_offset], &udotps[
		udotps_offset], &tm[1], &blbcn_1.par[10]);
	return 0;
    }

/* Take care of the case where the free parameters have been changed at */

/* the restart point. */

    if (nfpar1 != blicn_1.nfpar) {
	*nodir = 1;
	return 0;
    }
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	if (icprs[i - 1] != blbcn_1.icp[i - 1]) {
	    *nodir = 1;
	    return 0;
	}
/* L6: */
    }


    return 0;
} /* stpnbv_ */


/*     ---------- ------ */
/* Subroutine */ int stpnub_(ntstrs, ncolrs, lab, ibr, m1u, u, ups, udotps, 
	upoldp, tm, dtm, ndim2, smat, rnllv, ir, ic, f, dfdu, dfdp, nodir)
integer *ntstrs, *ncolrs, *lab, *ibr, *m1u;
doublereal *u, *ups, *udotps, *upoldp, *tm, *dtm;
integer *ndim2;
doublereal *smat, *rnllv;
integer *ir, *ic;
doublereal *f, *dfdu, *dfdp;
integer *nodir;
{
    /* System generated locals */
    integer ups_dim1, ups_offset, udotps_dim1, udotps_offset, upoldp_dim1, 
	    upoldp_offset, smat_dim1, smat_offset, dfdu_dim1, dfdu_offset, 
	    dfdp_dim1, dfdp_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer ncol1, i, j, k;
    static doublereal t;
    static integer k1, k2;
    extern /* Subroutine */ int stpnt_();
    static doublereal dt;
    extern /* Subroutine */ int msh_();


/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Generates a starting point for the continuation of a branch of */
/* of solutions to general boundary value problems by calling the user */
/* supplied subroutine STPNT where an analytical solution is given. */



/* Generate the (initially uniform) mesh. */

    /* Parameter adjustments */
    dfdp_dim1 = blbcn_1.ndim;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = blbcn_1.ndim;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --ic;
    --ir;
    --rnllv;
    smat_dim1 = *ndim2;
    smat_offset = smat_dim1 + 1;
    smat -= smat_offset;
    --dtm;
    --tm;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;
    --u;

    /* Function Body */
    msh_(&tm[1]);
    dt = blrcn_1.one / (blcde_1.ntst * blcde_1.ncol);

    i__1 = blcde_1.ntst + 1;
    for (j = 1; j <= i__1; ++j) {
	if (j == blcde_1.ntst + 1) {
	    ncol1 = 1;
	} else {
	    ncol1 = blcde_1.ncol;
	}
	i__2 = ncol1;
	for (i = 1; i <= i__2; ++i) {
	    t = tm[j] + (i - 1) * dt;
	    k1 = (i - 1) * blbcn_1.ndim + 1;
	    k2 = i * blbcn_1.ndim;
	    stpnt_(&blbcn_1.ndim, &u[1], blbcn_1.par, &t);
	    i__3 = k2;
	    for (k = k1; k <= i__3; ++k) {
		ups[j + k * ups_dim1] = u[k - k1 + 1];
/* L1: */
	    }
/* L2: */
	}
/* L3: */
    }

    *ntstrs = blcde_1.ntst;
    *ncolrs = blcde_1.ncol;
    *ibr = 1;
    *lab = 0;

    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	blcrl_1.rl[i - 1] = blbcn_1.par[blbcn_1.icp[i - 1] - 1];
/* L4: */
    }

    *nodir = 1;

    return 0;
} /* stpnub_ */


/*     ---------- ------ */
/* Subroutine */ int stdrbv_(funi, bcni, icni, rds, m1aa, m2aa, aa, m1bb, 
	m2bb, bb, m1cc, cc, m1dd, dd, wbrbd, m1u, ups, uoldps, udotps, upoldp,
	 rhsa, rhsd, dups, dtm, u, f, m1df, dfdu, dfdp, ubc0, ubc1, m1bc, dbc,
	 uicd, ficd, m1ic, dicd, ir, ic, iwbrbd, iperp)
/* Subroutine */ int (*funi) (), (*bcni) (), (*icni) ();
doublereal *rds;
integer *m1aa, *m2aa;
doublereal *aa;
integer *m1bb, *m2bb;
doublereal *bb;
integer *m1cc;
doublereal *cc;
integer *m1dd;
doublereal *dd, *wbrbd;
integer *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *rhsa, *rhsd, *dups, *dtm, *u, *f;

integer *m1df;
doublereal *dfdu, *dfdp, *ubc0, *ubc1;
integer *m1bc;
doublereal *dbc, *uicd, *ficd;
integer *m1ic;
doublereal *dicd;
integer *ir, *ic, *iwbrbd, *iperp;
{
    /* Format strings */
    static char fmt_101[] = "(\002 STARTING DIRECTION OF THE FREE PARAMETER(\
S) : \002,/,10e11.3,/,10e11.3)";

    /* System generated locals */
    integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, aa_dim1, aa_dim2, 
	    aa_offset, bb_dim1, bb_dim2, bb_offset, cc_dim1, cc_offset, 
	    dd_dim1, dd_offset, ups_dim1, ups_offset, uoldps_dim1, 
	    uoldps_offset, udotps_dim1, udotps_offset, upoldp_dim1, 
	    upoldp_offset, rhsa_dim1, rhsa_offset, dups_dim1, dups_offset, 
	    dbc_dim1, dbc_offset, dicd_dim1, dicd_offset, i__1, i__2;

    /* Builtin functions */
    integer s_wsfe(), do_fio(), e_wsfe();

    /* Local variables */
    extern /* Subroutine */ int brbd_();
    static integer iibr, ifst;
    static doublereal rdsw;
    static integer i, j;
    extern /* Subroutine */ int scalebb_(), setubv_();

    /* Fortran I/O blocks */
    static cilist io___127 = { 0, 9, 0, fmt_101, 0 };



/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Generates a direction vector (UDOTPS,RLDOT) that is needed to start */
/* the computation of a branch when no direction vector is given. */




/* Generate the Jacobian matrix with zero direction vector. */
/* (Then the last row of the Jacobian will be zero) */
/* in case the starting direction is to be determined. */

    /* Parameter adjustments */
    --iwbrbd;
    --ic;
    --ir;
    dicd_dim1 = *m1ic;
    dicd_offset = dicd_dim1 + 1;
    dicd -= dicd_offset;
    --ficd;
    --uicd;
    dbc_dim1 = *m1bc;
    dbc_offset = dbc_dim1 + 1;
    dbc -= dbc_offset;
    --ubc1;
    --ubc0;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    --dtm;
    dups_dim1 = *m1u;
    dups_offset = dups_dim1 + 1;
    dups -= dups_offset;
    --rhsd;
    rhsa_dim1 = *m1u;
    rhsa_offset = rhsa_dim1 + 1;
    rhsa -= rhsa_offset;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;
    --wbrbd;
    dd_dim1 = *m1dd;
    dd_offset = dd_dim1 + 1;
    dd -= dd_offset;
    cc_dim1 = *m1cc;
    cc_offset = cc_dim1 + 1;
    cc -= cc_offset;
    bb_dim1 = *m1bb;
    bb_dim2 = *m2bb;
    bb_offset = bb_dim1 * (bb_dim2 + 1) + 1;
    bb -= bb_offset;
    aa_dim1 = *m1aa;
    aa_dim2 = *m2aa;
    aa_offset = aa_dim1 * (aa_dim2 + 1) + 1;
    aa -= aa_offset;

    /* Function Body */
    if (*iperp == 0) {
	i__1 = blcde_1.ntst + 1;
	for (j = 1; j <= i__1; ++j) {
	    i__2 = blicn_1.nrow;
	    for (i = 1; i <= i__2; ++i) {
		udotps[j + i * udotps_dim1] = blrcn_1.zero;
/* L1: */
	    }
/* L2: */
	}
	i__1 = blicn_1.nfpar;
	for (i = 1; i <= i__1; ++i) {
	    blcrl_1.rldot[i - 1] = blrcn_1.zero;
/* L3: */
	}
    }

    rdsw = blrcn_1.zero;
    setubv_(funi, bcni, icni, &rdsw, m1aa, m2aa, &aa[aa_offset], m1bb, m2bb, &
	    bb[bb_offset], m1cc, &cc[cc_offset], m1dd, &dd[dd_offset], m1u, &
	    ups[ups_offset], &uoldps[uoldps_offset], &udotps[udotps_offset], &
	    upoldp[upoldp_offset], &rhsa[rhsa_offset], &rhsd[1], &dups[
	    dups_offset], &dtm[1], &u[1], &f[1], m1df, &dfdu[dfdu_offset], &
	    dfdp[dfdp_offset], &ubc0[1], &ubc1[1], m1bc, &dbc[dbc_offset], &
	    uicd[1], &ficd[1], m1ic, &dicd[dicd_offset]);

/* Find the null vector. */

    ifst = 1;
    if (blmax_1.iid < 4) {
	iibr = 0;
    } else if (blmax_1.iid == 4) {
	iibr = 1;
    } else {
	iibr = 2;
    }
    brbd_(&blcde_1.ntst, &blicn_1.nrow, &blicn_1.nclm, m1aa, m2aa, &aa[
	    aa_offset], &blicn_1.nfpar, m1bb, m2bb, &bb[bb_offset], &
	    blicn_1.nrc, m1cc, &cc[cc_offset], m1dd, &dd[dd_offset], m1u, &
	    rhsa[rhsa_offset], &rhsd[1], &wbrbd[1], &iibr, &ifst, &ir[1], &ic[
	    1], &iwbrbd[1], &c__1);

/* Compute the starting direction. */

    i__1 = blbcn_1.ndim;
    for (i = 1; i <= i__1; ++i) {
	udotps[blcde_1.ntst + 1 + i * udotps_dim1] = rhsd[i];
/* L4: */
    }
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	blcrl_1.rldot[i - 1] = rhsd[blbcn_1.ndim + i];
	blbcn_1.par[blbcn_1.icp[i - 1] - 1] = blcrl_1.rl[i - 1];
/* L5: */
    }

    i__1 = blcde_1.ntst;
    for (j = 1; j <= i__1; ++j) {
	i__2 = blicn_1.nrow;
	for (i = 1; i <= i__2; ++i) {
	    udotps[j + i * udotps_dim1] = rhsa[j + i * rhsa_dim1];
/* L6: */
	}
/* L7: */
    }

/* Scale the starting direction. */

    scalebb_(m1u, &udotps[udotps_offset], blcrl_1.rldot, &dtm[1]);

/* Make sure that RLDOT(1) is positive (unless zero). */

    if (blcrl_1.rldot[0] < blrcn_1.zero) {
	i__1 = blicn_1.nfpar;
	for (i = 1; i <= i__1; ++i) {
	    blcrl_1.rldot[i - 1] = -blcrl_1.rldot[i - 1];
/* L8: */
	}
	i__1 = blcde_1.ntst;
	for (j = 1; j <= i__1; ++j) {
	    i__2 = blicn_1.nrow;
	    for (i = 1; i <= i__2; ++i) {
		udotps[j + i * udotps_dim1] = -udotps[j + i * udotps_dim1];
/* L9: */
	    }
/* L10: */
	}
    }

    if (blmax_1.iid >= 2) {
	s_wsfe(&io___127);
	i__1 = blicn_1.nfpar;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&blcrl_1.rldot[i - 1], (ftnlen)sizeof(
		    doublereal));
	}
	e_wsfe();
    }

    return 0;
} /* stdrbv_ */


/* ------------------------------------------------------------------- */
/* ------------------------------------------------------------------- */
/*  Detection and Location of Bifurcations in Boundary Value Problems */
/* ------------------------------------------------------------------- */
/* ------------------------------------------------------------------- */

/*     ---------- ------ */
/* Subroutine */ int lcspbv_(fncs, funi, bcni, icni, istop, itp, sp1, nitps, 
	ibr, ntot, m1aa, m2aa, aa, m1bb, m2bb, bb, m1cc, cc, m1dd, dd, wbrbd, 
	m1u, ups, uoldps, udotps, upoldp, rhsa, rhsd, dups, tm, dtm, u, f, 
	m1df, dfdu, dfdp, ubc0, ubc1, m1bc, dbc, uicd, ficd, m1ic, dicd, ir, 
	ic, iwbrbd, p0, p1, poin, ev, wkev)
doublereal (*fncs) ();
/* Subroutine */ int (*funi) (), (*bcni) (), (*icni) ();
integer *istop, *itp;
doublereal *sp1;
integer *nitps, *ibr, *ntot, *m1aa, *m2aa;
doublereal *aa;
integer *m1bb, *m2bb;
doublereal *bb;
integer *m1cc;
doublereal *cc;
integer *m1dd;
doublereal *dd, *wbrbd;
integer *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *rhsa, *rhsd, *dups, *tm, *dtm, *
	u, *f;
integer *m1df;
doublereal *dfdu, *dfdp, *ubc0, *ubc1;
integer *m1bc;
doublereal *dbc, *uicd, *ficd;
integer *m1ic;
doublereal *dicd;
integer *ir, *ic, *iwbrbd;
doublereal *p0, *p1, *poin;
doublecomplex *ev;
doublereal *wkev;
{
    /* Format strings */
    static char fmt_102[] = "(\002 * DETECTION OF SINGULAR POINT : ITERATI\
ON \002,i3,\002 STEPSIZE =\002,e11.3)";
    static char fmt_101[] = "(\002 *** POSSIBLE SINGULAR POINT (BRANCH \002,\
i3,\002  POINT \002,i4,\002)\002)";

    /* System generated locals */
    integer aa_dim1, aa_dim2, aa_offset, bb_dim1, bb_dim2, bb_offset, cc_dim1,
	     cc_offset, dd_dim1, dd_offset, ups_dim1, ups_offset, udotps_dim1,
	     udotps_offset, upoldp_dim1, upoldp_offset, uoldps_dim1, 
	    uoldps_offset, rhsa_dim1, rhsa_offset, dups_dim1, dups_offset, 
	    dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, dbc_dim1, 
	    dbc_offset, dicd_dim1, dicd_offset, p0_dim1, p0_offset, p1_dim1, 
	    p1_offset, poin_dim1, poin_offset;

    /* Builtin functions */
    integer s_wsfe(), do_fio(), e_wsfe();

    /* Local variables */
    static logical chng;
    static doublereal rrds;
    static integer nitsp1, ntotp1;
    extern /* Subroutine */ int contbv_(), solvbv_();
    static doublereal sp10, rds, dsp1, psp1;

    /* Fortran I/O blocks */
    static cilist io___135 = { 0, 9, 0, fmt_102, 0 };
    static cilist io___137 = { 0, 9, 0, fmt_101, 0 };



/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* This subroutine uses the Secant method to accurately locate limit */
/* points, bifurcation points, and zero(es) of function(s) from USZR. */
/* Such points are located as points on a solution branch where the */
/* EXTERNAL function FNCS changes sign. */
/* It involves calling the basic solution subroutines CONTBV and SOLVBV */

/* with decreasing values of RDS (stepsize along branch). */
/* The point is assumed to have been found with sufficient accuracy if */
/* the ratio between RDS and the user supplied value of DS is less than */

/* the user-supplied tolerance EPSS. */
/* This subroutine is called from CNRLB, which controls the computation */

/* of branches of solutions to general boundary value problems. */



/* SGLE COMPLEX  EV(NDIM) */



    /* Parameter adjustments */
    --wkev;
    --ev;
    poin_dim1 = blbcn_1.ndim;
    poin_offset = poin_dim1 + 1;
    poin -= poin_offset;
    p1_dim1 = blbcn_1.ndim;
    p1_offset = p1_dim1 + 1;
    p1 -= p1_offset;
    p0_dim1 = blbcn_1.ndim;
    p0_offset = p0_dim1 + 1;
    p0 -= p0_offset;
    --iwbrbd;
    --ic;
    --ir;
    dicd_dim1 = *m1ic;
    dicd_offset = dicd_dim1 + 1;
    dicd -= dicd_offset;
    --ficd;
    --uicd;
    dbc_dim1 = *m1bc;
    dbc_offset = dbc_dim1 + 1;
    dbc -= dbc_offset;
    --ubc1;
    --ubc0;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    --dtm;
    --tm;
    dups_dim1 = *m1u;
    dups_offset = dups_dim1 + 1;
    dups -= dups_offset;
    --rhsd;
    rhsa_dim1 = *m1u;
    rhsa_offset = rhsa_dim1 + 1;
    rhsa -= rhsa_offset;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;
    --wbrbd;
    dd_dim1 = *m1dd;
    dd_offset = dd_dim1 + 1;
    dd -= dd_offset;
    cc_dim1 = *m1cc;
    cc_offset = cc_dim1 + 1;
    cc -= cc_offset;
    bb_dim1 = *m1bb;
    bb_dim2 = *m2bb;
    bb_offset = bb_dim1 * (bb_dim2 + 1) + 1;
    bb -= bb_offset;
    aa_dim1 = *m1aa;
    aa_dim2 = *m2aa;
    aa_offset = aa_dim1 * (aa_dim2 + 1) + 1;
    aa -= aa_offset;

    /* Function Body */
    sp10 = *sp1;

/* Check for zero. */

    *sp1 = (*fncs)(&chng, funi, bcni, icni, istop, itp, nitps, &p0[p0_offset],
	     &p1[p1_offset], &poin[poin_offset], &ev[1], &wkev[1], ibr, ntot, 
	    m1aa, m2aa, &aa[aa_offset], m1bb, m2bb, &bb[bb_offset], m1cc, &cc[
	    cc_offset], m1dd, &dd[dd_offset], &wbrbd[1], m1u, &ups[ups_offset]
	    , &uoldps[uoldps_offset], &udotps[udotps_offset], &upoldp[
	    upoldp_offset], &rhsa[rhsa_offset], &rhsd[1], &dups[dups_offset], 
	    &tm[1], &dtm[1], &u[1], &f[1], m1df, &dfdu[dfdu_offset], &dfdp[
	    dfdp_offset], &ubc0[1], &ubc1[1], m1bc, &dbc[dbc_offset], &uicd[1]
	    , &ficd[1], m1ic, &dicd[dicd_offset], &ir[1], &ic[1], &iwbrbd[1]);


    psp1 = sp10 * *sp1;
    ntotp1 = *ntot + 1;
    if (psp1 >= blrcn_1.zero || ! chng) {
	return 0;
    }

/* Compute next RDS by a perturbed Secant method : */

    rds = blcrl_1.rdsold;
    nitsp1 = 0;
L1:
    dsp1 = sp10 - *sp1;
    if (dsp1 == blrcn_1.zero) {
	rds = blrcn_1.zero;
    }
    if (dsp1 != blrcn_1.zero) {
	rds = *sp1 / dsp1 * rds;
    }
    rds = (blrcn_1.one + blrcn_1.hmach) * rds;

/* If requested write additional output on unit 9 : */

    if (blmax_1.iid >= 2) {
	s_wsfe(&io___135);
	do_fio(&c__1, (char *)&nitsp1, (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&rds, (ftnlen)sizeof(doublereal));
	e_wsfe();
    }

/* Return if tolerance has been met : */

    rrds = abs(rds) / (blrcn_1.one + abs(bldls_1.ds));
/* SGLE  RRDS= ABS(RDS)/(ONE+ ABS(DS)) */
    if (rrds < bleps_1.epss) {
	*itp = -1;
	return 0;
    }

    contbv_(funi, &rds, m1u, &ups[ups_offset], &uoldps[uoldps_offset], &
	    udotps[udotps_offset], &upoldp[upoldp_offset], &u[1], &f[1], m1df,
	     &dfdu[dfdu_offset], &dfdp[dfdp_offset], &dtm[1]);
    solvbv_(funi, bcni, icni, istop, &rds, nitps, ibr, ntot, m1aa, m2aa, &aa[
	    aa_offset], m1bb, m2bb, &bb[bb_offset], m1cc, &cc[cc_offset], 
	    m1dd, &dd[dd_offset], &wbrbd[1], m1u, &ups[ups_offset], &uoldps[
	    uoldps_offset], &udotps[udotps_offset], &upoldp[upoldp_offset], &
	    rhsa[rhsa_offset], &rhsd[1], &dups[dups_offset], &tm[1], &dtm[1], 
	    &u[1], &f[1], m1df, &dfdu[dfdu_offset], &dfdp[dfdp_offset], &ubc0[
	    1], &ubc1[1], m1bc, &dbc[dbc_offset], &uicd[1], &ficd[1], m1ic, &
	    dicd[dicd_offset], &ir[1], &ic[1], &iwbrbd[1]);
    if (*istop != 0) {
	*sp1 = blrcn_1.zero;
	return 0;
    }

/* Check for zero. */

    sp10 = *sp1;
    *sp1 = (*fncs)(&chng, funi, bcni, icni, istop, itp, nitps, &p0[p0_offset],
	     &p1[p1_offset], &poin[poin_offset], &ev[1], &wkev[1], ibr, ntot, 
	    m1aa, m2aa, &aa[aa_offset], m1bb, m2bb, &bb[bb_offset], m1cc, &cc[
	    cc_offset], m1dd, &dd[dd_offset], &wbrbd[1], m1u, &ups[ups_offset]
	    , &uoldps[uoldps_offset], &udotps[udotps_offset], &upoldp[
	    upoldp_offset], &rhsa[rhsa_offset], &rhsd[1], &dups[dups_offset], 
	    &tm[1], &dtm[1], &u[1], &f[1], m1df, &dfdu[dfdu_offset], &dfdp[
	    dfdp_offset], &ubc0[1], &ubc1[1], m1bc, &dbc[dbc_offset], &uicd[1]
	    , &ficd[1], m1ic, &dicd[dicd_offset], &ir[1], &ic[1], &iwbrbd[1]);


    ++nitsp1;
    if (nitsp1 <= blmax_1.itmx) {
	goto L1;
    }

    s_wsfe(&io___137);
    do_fio(&c__1, (char *)&(*ibr), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&ntotp1, (ftnlen)sizeof(integer));
    e_wsfe();
    *sp1 = blrcn_1.zero;

    return 0;
} /* lcspbv_ */


/*     ------ --------- */
doublereal fnlpbv_(chng, funi, bcni, icni, istop, itp, nitps, p0, p1, poin, 
	ev, wkev, ibr, ntot, m1aa, m2aa, aa, m1bb, m2bb, bb, m1cc, cc, m1dd, 
	dd, wbrbd, m1u, ups, uoldps, udotps, upoldp, rhsa, rhsd, dups, tm, 
	dtm, u, f, m1df, dfdu, dfdp, ubc0, ubc1, m1bc, dbc, uicd, ficd, m1ic, 
	dicd, ir, ic, iwbrbd)
logical *chng;
/* Subroutine */ int (*funi) (), (*bcni) (), (*icni) ();
integer *istop, *itp, *nitps;
doublereal *p0, *p1, *poin;
doublecomplex *ev;
doublereal *wkev;
integer *ibr, *ntot, *m1aa, *m2aa;
doublereal *aa;
integer *m1bb, *m2bb;
doublereal *bb;
integer *m1cc;
doublereal *cc;
integer *m1dd;
doublereal *dd, *wbrbd;
integer *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *rhsa, *rhsd, *dups, *tm, *dtm, *
	u, *f;
integer *m1df;
doublereal *dfdu, *dfdp, *ubc0, *ubc1;
integer *m1bc;
doublereal *dbc, *uicd, *ficd;
integer *m1ic;
doublereal *dicd;
integer *ir, *ic, *iwbrbd;
{
    /* Format strings */
    static char fmt_101[] = "(\002 LIMIT POINT FUNCTION = \002,e11.3)";

    /* System generated locals */
    integer aa_dim1, aa_dim2, aa_offset, bb_dim1, bb_dim2, bb_offset, cc_dim1,
	     cc_offset, dd_dim1, dd_offset, ups_dim1, ups_offset, udotps_dim1,
	     udotps_offset, upoldp_dim1, upoldp_offset, uoldps_dim1, 
	    uoldps_offset, rhsa_dim1, rhsa_offset, dups_dim1, dups_offset, 
	    dbc_dim1, dbc_offset, dicd_dim1, dicd_offset, dfdu_dim1, 
	    dfdu_offset, dfdp_dim1, dfdp_offset, p0_dim1, p0_offset, p1_dim1, 
	    p1_offset, poin_dim1, poin_offset, i__1, i__2;
    doublereal ret_val;

    /* Builtin functions */
    integer s_wsfe(), do_fio(), e_wsfe();

    /* Local variables */
    extern /* Subroutine */ int brbd_();
    static integer iibr, ifst, i, j;
    extern /* Subroutine */ int scalebb_();

    /* Fortran I/O blocks */
    static cilist io___142 = { 0, 9, 0, fmt_101, 0 };


/* SGLE REAL */

/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* RETURNS A QUANTITY THAT CHANGES SIGN AT A LIMIT POINT (BVP) */


/* SGLE COMPLEX  EV(NDIM) */




/* Find the direction vector. */

    /* Parameter adjustments */
    --iwbrbd;
    --ic;
    --ir;
    dicd_dim1 = *m1ic;
    dicd_offset = dicd_dim1 + 1;
    dicd -= dicd_offset;
    --ficd;
    --uicd;
    dbc_dim1 = *m1bc;
    dbc_offset = dbc_dim1 + 1;
    dbc -= dbc_offset;
    --ubc1;
    --ubc0;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    --dtm;
    --tm;
    dups_dim1 = *m1u;
    dups_offset = dups_dim1 + 1;
    dups -= dups_offset;
    --rhsd;
    rhsa_dim1 = *m1u;
    rhsa_offset = rhsa_dim1 + 1;
    rhsa -= rhsa_offset;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;
    --wbrbd;
    dd_dim1 = *m1dd;
    dd_offset = dd_dim1 + 1;
    dd -= dd_offset;
    cc_dim1 = *m1cc;
    cc_offset = cc_dim1 + 1;
    cc -= cc_offset;
    bb_dim1 = *m1bb;
    bb_dim2 = *m2bb;
    bb_offset = bb_dim1 * (bb_dim2 + 1) + 1;
    bb -= bb_offset;
    aa_dim1 = *m1aa;
    aa_dim2 = *m2aa;
    aa_offset = aa_dim1 * (aa_dim2 + 1) + 1;
    aa -= aa_offset;
    --wkev;
    --ev;
    poin_dim1 = blbcn_1.ndim;
    poin_offset = poin_dim1 + 1;
    poin -= poin_offset;
    p1_dim1 = blbcn_1.ndim;
    p1_offset = p1_dim1 + 1;
    p1 -= p1_offset;
    p0_dim1 = blbcn_1.ndim;
    p0_offset = p0_dim1 + 1;
    p0 -= p0_offset;

    /* Function Body */
    ifst = 0;
    if (blmax_1.iid < 4) {
	iibr = 0;
    } else if (blmax_1.iid == 4) {
	iibr = 1;
    } else {
	iibr = 2;
    }
    brbd_(&blcde_1.ntst, &blicn_1.nrow, &blicn_1.nclm, m1aa, m2aa, &aa[
	    aa_offset], &blicn_1.nfpar, m1bb, m2bb, &bb[bb_offset], &
	    blicn_1.nrc, m1cc, &cc[cc_offset], m1dd, &dd[dd_offset], m1u, &
	    rhsa[rhsa_offset], &rhsd[1], &wbrbd[1], &iibr, &ifst, &ir[1], &ic[
	    1], &iwbrbd[1], &c_n1);

    i__1 = blbcn_1.ndim;
    for (i = 1; i <= i__1; ++i) {
	udotps[blcde_1.ntst + 1 + i * udotps_dim1] = rhsd[i];
/* L1: */
    }
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	blcrl_1.rldot[i - 1] = rhsd[blbcn_1.ndim + i];
/* L2: */
    }

    i__1 = blcde_1.ntst;
    for (j = 1; j <= i__1; ++j) {
	i__2 = blicn_1.nrow;
	for (i = 1; i <= i__2; ++i) {
	    udotps[j + i * udotps_dim1] = rhsa[j + i * rhsa_dim1];
/* L3: */
	}
/* L4: */
    }

/* Scale the direction vector. */

    scalebb_(m1u, &udotps[udotps_offset], blcrl_1.rldot, &dtm[1]);
    if (blmax_1.iid >= 2) {
	s_wsfe(&io___142);
	do_fio(&c__1, (char *)&blcrl_1.rldot[0], (ftnlen)sizeof(doublereal));
	e_wsfe();
    }

/* Set the quantity to be returned. */

    ret_val = blcrl_1.rldot[0];
    if (abs(blcde_1.isp) == 3) {
	ret_val = blcrl_1.rldot[1];
    }
    *chng = TRUE_;


    return ret_val;
} /* fnlpbv_ */


/*     ------ --------- */
doublereal fnbpbv_(chng, funi, bcni, icni, istop, itp, nitps, p0, p1, poin, 
	ev, wkev, ibr, ntot, m1aa, m2aa, aa, m1bb, m2bb, bb, m1cc, cc, m1dd, 
	dd, wbrbd, m1u, ups, uoldps, udotps, upoldp, rhsa, rhsd, dups, tm, 
	dtm, u, f, m1df, dfdu, dfdp, ubc0, ubc1, m1bc, dbc, uicd, ficd, m1ic, 
	dicd, ir, ic, iwbrbd)
logical *chng;
/* Subroutine */ int (*funi) (), (*bcni) (), (*icni) ();
integer *istop, *itp, *nitps;
doublereal *p0, *p1, *poin;
doublecomplex *ev;
doublereal *wkev;
integer *ibr, *ntot, *m1aa, *m2aa;
doublereal *aa;
integer *m1bb, *m2bb;
doublereal *bb;
integer *m1cc;
doublereal *cc;
integer *m1dd;
doublereal *dd, *wbrbd;
integer *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *rhsa, *rhsd, *dups, *tm, *dtm, *
	u, *f;
integer *m1df;
doublereal *dfdu, *dfdp, *ubc0, *ubc1;
integer *m1bc;
doublereal *dbc, *uicd, *ficd;
integer *m1ic;
doublereal *dicd;
integer *ir, *ic, *iwbrbd;
{
    /* Format strings */
    static char fmt_101[] = "(\002 BIFURCATION FUNCTION = \002,e11.3)";

    /* System generated locals */
    integer aa_dim1, aa_dim2, aa_offset, bb_dim1, bb_dim2, bb_offset, cc_dim1,
	     cc_offset, dd_dim1, dd_offset, ups_dim1, ups_offset, udotps_dim1,
	     udotps_offset, upoldp_dim1, upoldp_offset, uoldps_dim1, 
	    uoldps_offset, rhsa_dim1, rhsa_offset, dups_dim1, dups_offset, 
	    dbc_dim1, dbc_offset, dicd_dim1, dicd_offset, dfdu_dim1, 
	    dfdu_offset, dfdp_dim1, dfdp_offset, p0_dim1, p0_offset, p1_dim1, 
	    p1_offset, poin_dim1, poin_offset;
    doublereal ret_val;

    /* Builtin functions */
    integer s_wsfe(), do_fio(), e_wsfe();

    /* Local variables */
    extern /* Subroutine */ int ge_();
    static doublereal det;

    /* Fortran I/O blocks */
    static cilist io___144 = { 0, 9, 0, fmt_101, 0 };


/* SGLE REAL */

/* SGLE IMPLICIT REAL             (A-H,O-Z) */


/* SGLE COMPLEX  EV(NDIM) */




/* Save the determinant of the reduced system. */

    /* Parameter adjustments */
    --iwbrbd;
    --ic;
    --ir;
    dicd_dim1 = *m1ic;
    dicd_offset = dicd_dim1 + 1;
    dicd -= dicd_offset;
    --ficd;
    --uicd;
    dbc_dim1 = *m1bc;
    dbc_offset = dbc_dim1 + 1;
    dbc -= dbc_offset;
    --ubc1;
    --ubc0;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    --dtm;
    --tm;
    dups_dim1 = *m1u;
    dups_offset = dups_dim1 + 1;
    dups -= dups_offset;
    --rhsd;
    rhsa_dim1 = *m1u;
    rhsa_offset = rhsa_dim1 + 1;
    rhsa -= rhsa_offset;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;
    --wbrbd;
    dd_dim1 = *m1dd;
    dd_offset = dd_dim1 + 1;
    dd -= dd_offset;
    cc_dim1 = *m1cc;
    cc_offset = cc_dim1 + 1;
    cc -= cc_offset;
    bb_dim1 = *m1bb;
    bb_dim2 = *m2bb;
    bb_offset = bb_dim1 * (bb_dim2 + 1) + 1;
    bb -= bb_offset;
    aa_dim1 = *m1aa;
    aa_dim2 = *m2aa;
    aa_offset = aa_dim1 * (aa_dim2 + 1) + 1;
    aa -= aa_offset;
    --wkev;
    --ev;
    poin_dim1 = blbcn_1.ndim;
    poin_offset = poin_dim1 + 1;
    poin -= poin_offset;
    p1_dim1 = blbcn_1.ndim;
    p1_offset = p1_dim1 + 1;
    p1 -= p1_offset;
    p0_dim1 = blbcn_1.ndim;
    p0_offset = p0_dim1 + 1;
    p0 -= p0_offset;

    /* Function Body */
    det = bldet_1.detge;

/* Compute the determinant of P1. */
    ge_(&blbcn_1.ndim, &blbcn_1.ndim, &p1[p1_offset], &c__0, &c__1, &u[1], &
	    c__1, &f[1], &ir[1], &ic[1]);

/* Set the determinant of the normalized reduced system. */

    if (bldet_1.detge != blrcn_1.zero) {
	ret_val = det / bldet_1.detge;
	*chng = TRUE_;
    } else {
	ret_val = blrcn_1.zero;
	*chng = FALSE_;
    }

    if (blmax_1.iid >= 2) {
	s_wsfe(&io___144);
	do_fio(&c__1, (char *)&ret_val, (ftnlen)sizeof(doublereal));
	e_wsfe();
    }

    return ret_val;
} /* fnbpbv_ */


/*     ------ --------- */
doublereal fnspbv_(chng, funi, bcni, icni, istop, itp, nitps, p0, p1, poin, 
	ev, wkev, ibr, ntot, m1aa, m2aa, aa, m1bb, m2bb, bb, m1cc, cc, m1dd, 
	dd, wbrbd, m1u, ups, uoldps, udotps, upoldp, rhsa, rhsd, dups, tm, 
	dtm, u, f, m1df, dfdu, dfdp, ubc0, ubc1, m1bc, dbc, uicd, ficd, m1ic, 
	dicd, ir, ic, iwbrbd)
logical *chng;
/* Subroutine */ int (*funi) (), (*bcni) (), (*icni) ();
integer *istop, *itp, *nitps;
doublereal *p0, *p1, *poin;
doublecomplex *ev;
doublereal *wkev;
integer *ibr, *ntot, *m1aa, *m2aa;
doublereal *aa;
integer *m1bb, *m2bb;
doublereal *bb;
integer *m1cc;
doublereal *cc;
integer *m1dd;
doublereal *dd, *wbrbd;
integer *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *rhsa, *rhsd, *dups, *tm, *dtm, *
	u, *f;
integer *m1df;
doublereal *dfdu, *dfdp, *ubc0, *ubc1;
integer *m1bc;
doublereal *dbc, *uicd, *ficd;
integer *m1ic;
doublereal *dicd;
integer *ir, *ic, *iwbrbd;
{
    /* Format strings */
    static char fmt_104[] = "(\002 FLOQUET MULTIPLIERS :\002)";
    static char fmt_105[] = "(\002 BRANCH \002,i3,\002  POINT \002,i4,2(2x,2\
e12.5),50(/,23x,2(2x,2e12.5)))";
    static char fmt_101[] = "(\002 *** FLOQUET MULTIPLIER AT 1 INACCURATE\
\002)";
    static char fmt_102[] = "(\002 *** FLOQUET MULTIPLIER AT 1 ACCURATE AG\
AIN\002)";
    static char fmt_103[] = "(\002 BIFURCATION FUNCTION = \002,e11.3,\002 # \
OF MULTIPLIERS IN UNIT CIRCLE =\002,i3)";

    /* System generated locals */
    integer aa_dim1, aa_dim2, aa_offset, bb_dim1, bb_dim2, bb_offset, cc_dim1,
	     cc_offset, dd_dim1, dd_offset, ups_dim1, ups_offset, udotps_dim1,
	     udotps_offset, upoldp_dim1, upoldp_offset, uoldps_dim1, 
	    uoldps_offset, rhsa_dim1, rhsa_offset, dups_dim1, dups_offset, 
	    dbc_dim1, dbc_offset, dicd_dim1, dicd_offset, dfdu_dim1, 
	    dfdu_offset, dfdp_dim1, dfdp_offset, p0_dim1, p0_offset, p1_dim1, 
	    p1_offset, poin_dim1, poin_offset, i__1, i__2, i__3,i_1,i_2;
    doublereal ret_val;
    doublecomplex z__1;

    /* Builtin functions */
    double z_abs();
    integer s_wsfe(), e_wsfe(), do_fio();

    /* Local variables */
    static doublereal amin;
    static doublecomplex ztmp;
    static integer nins1, nins2;
    static doublereal d;
    static integer ntot1, i, j;
    extern /* Subroutine */ int poinc_();
    static doublereal ad;
    extern /* Subroutine */ int eig_();
    extern  int flowkm_() ;
    static integer loc, ier;
    static doublereal azm1;
    static integer isp1;

    /* Fortran I/O blocks */
    static cilist io___157 = { 0, 9, 0, fmt_104, 0 };
    static cilist io___158 = { 0, 9, 0, fmt_105, 0 };
    static cilist io___159 = { 0, 9, 0, fmt_101, 0 };
    static cilist io___160 = { 0, 9, 0, fmt_102, 0 };
    static cilist io___161 = { 0, 9, 0, fmt_104, 0 };
    static cilist io___162 = { 0, 9, 0, fmt_105, 0 };
    static cilist io___164 = { 0, 9, 0, fmt_103, 0 };
    static cilist io___165 = { 0, 9, 0, fmt_104, 0 };
    static cilist io___166 = { 0, 9, 0, fmt_105, 0 };


/* SGLE REAL */

/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* This function returns a quantity that changes sign when a complex */
/* pair of eigenvalues of the linearized Poincare map moves in or out */
/* of the unit circle. */


/* SGLE COMPLEX  EV(NDIM),ZTMP */




/* Initialize. */

    /* Parameter adjustments */
    --iwbrbd;
    --ic;
    --ir;
    dicd_dim1 = *m1ic;
    dicd_offset = dicd_dim1 + 1;
    dicd -= dicd_offset;
    --ficd;
    --uicd;
    dbc_dim1 = *m1bc;
    dbc_offset = dbc_dim1 + 1;
    dbc -= dbc_offset;
    --ubc1;
    --ubc0;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    --dtm;
    --tm;
    dups_dim1 = *m1u;
    dups_offset = dups_dim1 + 1;
    dups -= dups_offset;
    --rhsd;
    rhsa_dim1 = *m1u;
    rhsa_offset = rhsa_dim1 + 1;
    rhsa -= rhsa_offset;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;
    --wbrbd;
    dd_dim1 = *m1dd;
    dd_offset = dd_dim1 + 1;
    dd -= dd_offset;
    cc_dim1 = *m1cc;
    cc_offset = cc_dim1 + 1;
    cc -= cc_offset;
    bb_dim1 = *m1bb;
    bb_dim2 = *m2bb;
    bb_offset = bb_dim1 * (bb_dim2 + 1) + 1;
    bb -= bb_offset;
    aa_dim1 = *m1aa;
    aa_dim2 = *m2aa;
    aa_offset = aa_dim1 * (aa_dim2 + 1) + 1;
    aa -= aa_offset;
    --wkev;
    --ev;
    poin_dim1 = blbcn_1.ndim;
    poin_offset = poin_dim1 + 1;
    poin -= poin_offset;
    p1_dim1 = blbcn_1.ndim;
    p1_offset = p1_dim1 + 1;
    p1 -= p1_offset;
    p0_dim1 = blbcn_1.ndim;
    p0_offset = p0_dim1 + 1;
    p0 -= p0_offset;

    /* Function Body */
    ret_val = blrcn_1.zero;
    d = blrcn_1.zero;
    *chng = FALSE_;
   if(FLOWK==0)  {
/* Compute the linearization of the Poincare map. */

    poinc_(&blbcn_1.ndim, &p0[p0_offset], &p1[p1_offset], &poin[poin_offset], 
	    &blmax_1.iid, &ir[1], &ic[1]);

/* Compute the Floquet multipliers. */

    eig_(&blbcn_1.ndim, &blbcn_1.ndim, &poin[poin_offset], &ev[1], &wkev[1], &
	    ier);
}   else   {


 	i_1 = blbcn_1.ndim;
	for (j = 1; j <= i_1; ++j) {
	    i_2 = blbcn_1.ndim;
	    for (i = 1; i <= i_2; ++i) {
		p1[i + j * p1_dim1] = -p1[i + j * p1_dim1];

		}
		}

	flowkm_(&blbcn_1.ndim, &p0[p0_offset], &p1[p1_offset], &blmax_1.iid, &
		ir[1], &ic[1], &poin[poin_offset], &ev[1]);
		
      }


/* Order the Floquet multipliers by distance from z=1. */
    send_mult(*ibr,*ntot+1,blbcn_1.ndim,&ev[1]);
    i__1 = blbcn_1.ndim - 1;
    for (i = 1; i <= i__1; ++i) {
	amin = blrcn_1.rlarge;
	i__2 = blbcn_1.ndim;
	for (j = i; j <= i__2; ++j) {
	    i__3 = j;
	    z__1.r = ev[i__3].r - blrcn_1.one, z__1.i = ev[i__3].i;
	    azm1 = z_abs(&z__1);
/* SGLE      AZM1= CABS( EV(J) - ONE ) */
	    if (azm1 > amin) {
		goto L1;
	    }
	    amin = azm1;
	    loc = j;
L1:
	    ;
	}
	if (loc != i) {
	    i__2 = loc;
	    ztmp.r = ev[i__2].r, ztmp.i = ev[i__2].i;
	    i__2 = loc;
	    i__3 = i;
	    ev[i__2].r = ev[i__3].r, ev[i__2].i = ev[i__3].i;
	    i__2 = i;
	    ev[i__2].r = ztmp.r, ev[i__2].i = ztmp.i;
	}
/* L2: */
    }

/* Find eigenvalue closest to the unit circle */
/* (excluding the eigenvalue at z=1). */

    nins1= abs(blcde_1.isp) - 1;
    if (nins1 < 1) {
	nins1 = 1;
    }
    if (nins1 > blbcn_1.ndim) {
	nins1 = blbcn_1.ndim;
    }

    if (nins1 < blbcn_1.ndim) {
	amin = blrcn_1.rlarge;
	i__1 = blbcn_1.ndim;
	for (i = nins1 + 1; i <= i__1; ++i) {
	    d = z_abs(&ev[i]) - blrcn_1.one;
/* SGLE      D= CABS(EV(I)) - ONE */
	    ad = abs(d);
/* SGLE      AD= ABS(D) */
	    if (ad > amin) {
		goto L3;
	    }
	    amin = ad;
	    loc = i;
L3:
	    ;
	}
/*        Interchange, to put eigenvalue in ISP'th position. */
	isp1 = blcde_1.isp;
	if (blcde_1.isp == 1) {
	    isp1 = 2;
	}
	if (loc != isp1) {
	    i__1 = loc;
	    ztmp.r = ev[i__1].r, ztmp.i = ev[i__1].i;
	    i__1 = loc;
	    i__2 = isp1;
	    ev[i__1].r = ev[i__2].r, ev[i__1].i = ev[i__2].i;
	    i__1 = isp1;
	    ev[i__1].r = ztmp.r, ev[i__1].i = ztmp.i;
	}
    }

/* Print error message if the Floquet multiplier at z=1 is inaccurate. */
/* (ISP is set to negative and detection of bifurations is discontinued) 
*/

    z__1.r = ev[1].r - blrcn_1.one, z__1.i = ev[1].i;
    amin = z_abs(&z__1);
/* SGLE  AMIN= CABS( EV(1) - ONE ) */

    if (amin > .05 && blcde_1.isp > 1) {
	ntot1 = *ntot + 1;
	if (blmax_1.iid >= 2) {
	    s_wsfe(&io___157);
	    e_wsfe();
	}
	s_wsfe(&io___158);
	do_fio(&c__1, (char *)&(*ibr), (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&ntot1, (ftnlen)sizeof(integer));
	i__1 = blbcn_1.ndim;

	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__2, (char *)&ev[i], (ftnlen)sizeof(doublereal));
	}
	e_wsfe();
	bldet_1.nins = 0;
	s_wsfe(&io___159);
	e_wsfe();
	blcde_1.isp = -blcde_1.isp;
	return ret_val;
    }

/* Restart automatic detection if the Floquet multiplier at z=1 is */
/* sufficiently accurate again. */

    if (blcde_1.isp < 0) {
	if (amin < .02) {
	    s_wsfe(&io___160);
	    e_wsfe();
	    blcde_1.isp = -blcde_1.isp;
	} else {
	    ntot1 = *ntot + 1;
	    if (blmax_1.iid >= 2) {
		s_wsfe(&io___161);
		e_wsfe();
	    }
	    s_wsfe(&io___162);
	    do_fio(&c__1, (char *)&(*ibr), (ftnlen)sizeof(integer));
	    do_fio(&c__1, (char *)&ntot1, (ftnlen)sizeof(integer));
	    i__1 = blbcn_1.ndim;

	    for (i = 1; i <= i__1; ++i) {
		do_fio(&c__2, (char *)&ev[i], (ftnlen)sizeof(doublereal));
	    }
	    e_wsfe();
	    return ret_val;
	}
    }

/* Count the number of Floquet multipliers inside the unit circle. */

    if (nins1 == blbcn_1.ndim) {
	d = blrcn_1.zero;
	ret_val = d;
	goto L5;
    }

    nins2 = nins1;
    i__1 = blbcn_1.ndim;
    for (i = nins2 + 1; i <= i__1; ++i) {
	if (z_abs(&ev[i]) <= blrcn_1.one) {
	    ++nins1;
	}
/* SGLE    IF( CABS(EV(I)).LE.ONE)NINS1=NINS1+1 */
/* L4: */
    }

    if (blcde_1.isp >= 2) {
	d = z_abs(&ev[blcde_1.isp]) - blrcn_1.one;
/* SGLE    D= CABS(EV(ISP)) - ONE */
	ret_val = d;
	if (nins1 != bldet_1.nins) {
	    *chng = TRUE_;
	}
    }
L5:
    bldet_1.nins = nins1;
    if (blmax_1.iid >= 2 && blcde_1.isp >= 1) {
	s_wsfe(&io___164);
	do_fio(&c__1, (char *)&d, (ftnlen)sizeof(doublereal));
	do_fio(&c__1, (char *)&bldet_1.nins, (ftnlen)sizeof(integer));
	e_wsfe();
    }

/* Print the Floquet multipliers. */

    ntot1 = *ntot + 1;
    if (bldet_1.nins == blbcn_1.ndim) {
	ntot1 = -ntot1;
    }
    if (blmax_1.iid >= 2) {
	s_wsfe(&io___165);
	e_wsfe();
    }
    s_wsfe(&io___166);
    do_fio(&c__1, (char *)&(*ibr), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&ntot1, (ftnlen)sizeof(integer));
    i__1 = blbcn_1.ndim;
    for (i = 1; i <= i__1; ++i) {
	do_fio(&c__2, (char *)&ev[i], (ftnlen)sizeof(doublereal));
    }
    e_wsfe();


    return ret_val;
} /* fnspbv_ */


/*     ------ --------- */
doublereal fnuzbv_(chng, funi, bcni, icni, istop, itp, nitps, p0, p1, poin, 
	ev, wkev, ibr, ntot, m1aa, m2aa, aa, m1bb, m2bb, bb, m1cc, cc, m1dd, 
	dd, wbrbd, m1u, ups, uoldps, udotps, upoldp, rhsa, rhsd, dups, tm, 
	dtm, u, f, m1df, dfdu, dfdp, ubc0, ubc1, m1bc, dbc, uicd, ficd, m1ic, 
	dicd, ir, ic, iwbrbd)
logical *chng;
/* Subroutine */ int (*funi) (), (*bcni) (), (*icni) ();
integer *istop, *itp, *nitps;
doublereal *p0, *p1, *poin;
doublecomplex *ev;
doublereal *wkev;
integer *ibr, *ntot, *m1aa, *m2aa;
doublereal *aa;
integer *m1bb, *m2bb;
doublereal *bb;
integer *m1cc;
doublereal *cc;
integer *m1dd;
doublereal *dd, *wbrbd;
integer *m1u;
doublereal *ups, *uoldps, *udotps, *upoldp, *rhsa, *rhsd, *dups, *tm, *dtm, *
	u, *f;
integer *m1df;
doublereal *dfdu, *dfdp, *ubc0, *ubc1;
integer *m1bc;
doublereal *dbc, *uicd, *ficd;
integer *m1ic;
doublereal *dicd;
integer *ir, *ic, *iwbrbd;
{
    /* Format strings */
    static char fmt_101[] = "(\002 USZR        FUNCTION = \002,e11.3)";

    /* System generated locals */
    integer aa_dim1, aa_dim2, aa_offset, bb_dim1, bb_dim2, bb_offset, cc_dim1,
	     cc_offset, dd_dim1, dd_offset, ups_dim1, ups_offset, udotps_dim1,
	     udotps_offset, upoldp_dim1, upoldp_offset, uoldps_dim1, 
	    uoldps_offset, rhsa_dim1, rhsa_offset, dups_dim1, dups_offset, 
	    dbc_dim1, dbc_offset, dicd_dim1, dicd_offset, dfdu_dim1, 
	    dfdu_offset, dfdp_dim1, dfdp_offset, p0_dim1, p0_offset, p1_dim1, 
	    p1_offset, poin_dim1, poin_offset;
    doublereal ret_val;

    /* Builtin functions */
    integer s_wsfe(), do_fio(), e_wsfe();

    /* Local variables */
    extern doublereal uszr_();

    /* Fortran I/O blocks */
    static cilist io___167 = { 0, 9, 0, fmt_101, 0 };


/* SGLE REAL */

/* SGLE IMPLICIT REAL             (A-H,O-Z) */


/* SGLE COMPLEX  EV(NDIM) */




    /* Parameter adjustments */
    --iwbrbd;
    --ic;
    --ir;
    dicd_dim1 = *m1ic;
    dicd_offset = dicd_dim1 + 1;
    dicd -= dicd_offset;
    --ficd;
    --uicd;
    dbc_dim1 = *m1bc;
    dbc_offset = dbc_dim1 + 1;
    dbc -= dbc_offset;
    --ubc1;
    --ubc0;
    dfdp_dim1 = *m1df;
    dfdp_offset = dfdp_dim1 + 1;
    dfdp -= dfdp_offset;
    dfdu_dim1 = *m1df;
    dfdu_offset = dfdu_dim1 + 1;
    dfdu -= dfdu_offset;
    --f;
    --u;
    --dtm;
    --tm;
    dups_dim1 = *m1u;
    dups_offset = dups_dim1 + 1;
    dups -= dups_offset;
    --rhsd;
    rhsa_dim1 = *m1u;
    rhsa_offset = rhsa_dim1 + 1;
    rhsa -= rhsa_offset;
    upoldp_dim1 = *m1u;
    upoldp_offset = upoldp_dim1 + 1;
    upoldp -= upoldp_offset;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    uoldps_dim1 = *m1u;
    uoldps_offset = uoldps_dim1 + 1;
    uoldps -= uoldps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;
    --wbrbd;
    dd_dim1 = *m1dd;
    dd_offset = dd_dim1 + 1;
    dd -= dd_offset;
    cc_dim1 = *m1cc;
    cc_offset = cc_dim1 + 1;
    cc -= cc_offset;
    bb_dim1 = *m1bb;
    bb_dim2 = *m2bb;
    bb_offset = bb_dim1 * (bb_dim2 + 1) + 1;
    bb -= bb_offset;
    aa_dim1 = *m1aa;
    aa_dim2 = *m2aa;
    aa_offset = aa_dim1 * (aa_dim2 + 1) + 1;
    aa -= aa_offset;
    --wkev;
    --ev;
    poin_dim1 = blbcn_1.ndim;
    poin_offset = poin_dim1 + 1;
    poin -= poin_offset;
    p1_dim1 = blbcn_1.ndim;
    p1_offset = p1_dim1 + 1;
    p1 -= p1_offset;
    p0_dim1 = blbcn_1.ndim;
    p0_offset = p0_dim1 + 1;
    p0 -= p0_offset;

    /* Function Body */
    ret_val = uszr_(&blusz_1.iuzr, &bllim_1.nuzr, blbcn_1.par);
    *chng = TRUE_;

    if (blmax_1.iid >= 2) {
	s_wsfe(&io___167);
	do_fio(&c__1, (char *)&ret_val, (ftnlen)sizeof(doublereal));
	e_wsfe();
    }

    return ret_val;
} /* fnuzbv_ */


/*     ---------- ----- */
/* Subroutine */ int poinc_(ndim, p0, p1, poin, iid, ir, ic)
integer *ndim;
doublereal *p0, *p1, *poin;
integer *iid, *ir, *ic;
{
    /* Format strings */
    static char fmt_101[] = "(\002 LINEARIZED POINCARE MAP\002)";
    static char fmt_102[] = "(1x,6e21.14)";

    /* System generated locals */
    integer p0_dim1, p0_offset, p1_dim1, p1_offset, poin_dim1, poin_offset, 
	    i__1, i__2;

    /* Builtin functions */
    integer s_wsfe(), e_wsfe(), do_fio();

    /* Local variables */
    static integer i, j;
    extern /* Subroutine */ int ge_();

    /* Fortran I/O blocks */
    static cilist io___170 = { 0, 9, 0, fmt_101, 0 };
    static cilist io___171 = { 0, 9, 0, fmt_102, 0 };



/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Computes the linearized Poincare map. This map is extracted from the */

/* decomposition of the Jacobian matrix as generated by BRBD. */


    /* Parameter adjustments */
    --ic;
    --ir;
    poin_dim1 = *ndim;
    poin_offset = poin_dim1 + 1;
    poin -= poin_offset;
    p1_dim1 = *ndim;
    p1_offset = p1_dim1 + 1;
    p1 -= p1_offset;
    p0_dim1 = *ndim;
    p0_offset = p0_dim1 + 1;
    p0 -= p0_offset;

    /* Function Body */
    i__1 = *ndim;
    for (i = 1; i <= i__1; ++i) {
	i__2 = *ndim;
	for (j = 1; j <= i__2; ++j) {
	    p0[i + j * p0_dim1] = -p0[i + j * p0_dim1];
/* L1: */
	}
/* L2: */
    }

    ge_(ndim, ndim, &p1[p1_offset], ndim, ndim, &poin[poin_offset], ndim, &p0[
	    p0_offset], &ir[1], &ic[1]);

    if (*iid > 2) {
	s_wsfe(&io___170);
	e_wsfe();
	i__1 = *ndim;
	for (i = 1; i <= i__1; ++i) {
	    s_wsfe(&io___171);
	    i__2 = *ndim;
	    for (j = 1; j <= i__2; ++j) {
		do_fio(&c__1, (char *)&poin[i + j * poin_dim1], (ftnlen)
			sizeof(doublereal));
	    }
	    e_wsfe();
/* L3: */
	}
    }


    return 0;
} /* poinc_ */


/*     ---------- ------ */
/* Subroutine */ int tpspbv_(ev, itp)
doublecomplex *ev;
integer *itp;
{
    /* System generated locals */
    integer i__1, i__2;
    doublereal d__1;
    doublecomplex z__1;

    /* Builtin functions */
    double z_abs(), d_imag(), sqrt(), asin();

    /* Local variables */
    static doublereal amin, d;
    static integer i;
    extern doublereal dreal_();
    static doublereal ad;
    static integer loc, loc1;
    static doublereal azm1;


/* Determines type of secondary periodic bifurcation. */

/* SGLE IMPLICIT REAL             (A-H,O-Z) */


/* SGLE COMPLEX  EV(NDIM) */

/* Find the eigenvalue closest to z=1. */

    /* Parameter adjustments */
    --ev;

    /* Function Body */
    loc = 1;
    amin = blrcn_1.rlarge;
    i__1 = blbcn_1.ndim;
    for (i = 1; i <= i__1; ++i) {
	i__2 = i;
	z__1.r = ev[i__2].r - blrcn_1.one, z__1.i = ev[i__2].i;
	azm1 = z_abs(&z__1);
/* SGLE    AZM1= CABS( EV(I) - ONE ) */
	if (azm1 > amin) {
	    goto L1;
	}
	amin = azm1;
	loc = i;
L1:
	;
    }

/* Find the eigenvalue closest to the unit circle */
/* (excluding the eigenvalue at z=1). */

    loc1 = 1;
    amin = blrcn_1.rlarge;
    i__1 = blbcn_1.ndim;
    for (i = 1; i <= i__1; ++i) {
	if (i == loc) {
	    goto L2;
	}
	d = z_abs(&ev[i]) - blrcn_1.one;
/* SGLE    D= CABS(EV(I)) - ONE */
	ad = abs(d);
/* SGLE    AD= ABS(D) */
	if (ad > amin) {
	    goto L2;
	}
	amin = ad;
	loc1 = i;
L2:
	;
    }

    if ((d__1 = d_imag(&ev[loc1]), abs(d__1)) > sqrt(bleps_1.epss)) {
/* SGLE IF( ABS(AIMAG(EV(LOC1))).GT. SQRT(EPSS))THEN */
/*       ** torus bifurcation */
	*itp = blitp_1.itpst * 10 + 8;
	blbcn_1.par[11] = asin((d_imag(&ev[loc1])));
/* SGLE   PAR(12)= ASIN(AIMAG(EV(LOC1))) */
    } else if (dreal_(&ev[loc1]) < -blrcn_1.half) {
/* SGLE ELSE IF( REAL(EV(LOC1)).LT.-HALF)THEN */
/*       ** period doubling */
	*itp = blitp_1.itpst * 10 + 7;
    } else {
/*       ** ordinary bifurcation or something else... */
	*itp = blitp_1.itpst * 10 + 6;
    }

    return 0;
} /* tpspbv_ */


/* ------------------------------------------------------------------- */
/* ------------------------------------------------------------------- */
/*                    Output (Boundary Value Problems) */
/* ------------------------------------------------------------------- */
/* ------------------------------------------------------------------- */

/*     ---------- ------ */
/* Subroutine */ int stplbv_(istop, itp, nitps, ntot, lab, ibr, m1u, ups, 
	udotps, tm, dtm, m1df)
integer *istop, *itp, *nitps, *ntot, *lab, *ibr, *m1u;
doublereal *ups, *udotps, *tm, *dtm;
integer *m1df;
{
    /* System generated locals */
    integer ups_dim1, ups_offset, udotps_dim1, udotps_offset, i__1;

    /* Builtin functions */
    double sqrt();

    /* Local variables */
    static integer itmp;
    static doublereal ulow[NAUTO];
    extern doublereal rnrm2_();
    static integer ntot1, i,iflag,i_1;
    static doublereal uhigh[NAUTO],u0[NAUTO],ubar[NAUTO];
    extern doublereal rintg_();
    static integer n2;
    extern /* Subroutine */ int wrtbv8_(), wrline_();
    extern doublereal rnrmsq_(), rmnups_(), rmxups_();
    static integer iab;
    static doublereal umx[7];
    static integer lab1, ibr1;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Writes the bifurcation diagram on unit 7 (Differential Equations) */
/* (Also controls the writing of complete solutions on unit 8). */
/* Every line written contains, in order, the following: */

/*  IBR    : The label of the branch. */
/*  NTOT   : The index of the point on the branch. */
/*           (Points are numbered consecutively along a branch). */
/*           If IPS=2 or 3, then the sign of NTOT indicates stability : */

/*            - = stable , + = unstable, or unknown. */
/*  ITP    : An integer indicating the type of point : */

/*             4  (  )  :   Output point (Every NPR steps along branch). 
*/
/*            -4  (UZ)  :   Output point (Zero of user function USZR). */
/*             5  (LP)  :   Limit point (fold). */
/*             6  (BP)  :   Bifurcation point. */
/*             7  (PD)  :   Period doubling bifurcation. */
/*             8  (TR)  :   Bifurcation to an invariant torus. */
/*             9  (EP)  :   End point of branch, normal termination. */
/*            -9  (MX)  :   End point of branch, abnormal termination. */

/*  LAB        : The label of a special point. */
/*  PAR(ICP(1)): The principal parameter. */
/*  A          : The L2-norm of the solution vector, or other measure of 
*/
/*               the solution (see the user-supplied parameter IPLT). */
/*  MAX U(*)   : The maxima of the first few solution components. */
/*  PAR(ICP(*)): Further free parameters (if any). */



    /* Parameter adjustments */
    --dtm;
    --tm;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;

    /* Function Body */
    ++(*ntot);

/* ITP is set to 4 every NPR steps along a branch of solns and the entire 
*/
/* solution is written on unit 8. */

    if (*ntot % blmax_1.npr == 0 && *itp % 10 == 0) {
	*itp = blitp_1.itpst * 10 + 4;
    }

/* Check whether limits of the bifurcation diagram have been reached : */

    iab = abs(blcde_1.iplt);
    if (iab == 0 || iab > blicn_1.ndm * 3) {
	blcrl_1.a = sqrt((rnrmsq_(&blicn_1.ndm, m1u, &ups[ups_offset], &dtm[1]
		)));
    }
/* SGLE  IF(IAB.EQ.0.OR.IAB.GT.3*NDM)A= SQRT(RNRMSQ(NDM,M1U,UPS,DTM)) */
    if (blcde_1.iplt > 0 && iab <= blicn_1.ndm) {
	blcrl_1.a = rmxups_(m1u, &iab, &ups[ups_offset]);
    }
    if (blcde_1.iplt > blicn_1.ndm && iab <= blicn_1.ndm << 1) {
	i__1 = iab - blicn_1.ndm;
	blcrl_1.a = rintg_(m1u, &i__1, &ups[ups_offset], &dtm[1]);
    }
    if (blcde_1.iplt > blicn_1.ndm << 1 && iab <= blicn_1.ndm * 3) {
	i__1 = iab - (blicn_1.ndm << 1);
	blcrl_1.a = rnrm2_(m1u, &i__1, &ups[ups_offset], &dtm[1]);
    }
    if (blcde_1.iplt < 0 && iab <= blicn_1.ndm) {
	blcrl_1.a = rmnups_(m1u, &iab, &ups[ups_offset]);
    }
/* Externel interupts */
    byeauto_(ntot, &iflag);
    if (*istop == 1) {
/*        ** Maximum number of iterations reached somewhere. */
	*itp = -9 - blitp_1.itpst * 10;
    } else {
	if (blbcn_1.par[blbcn_1.icp[0] - 1] < bllim_1.rl0 || blbcn_1.par[
		blbcn_1.icp[0] - 1] > bllim_1.rl1 || blcrl_1.a < bllim_1.a0 ||
		 blcrl_1.a > bllim_1.a1 || *ntot >= bllim_1.nmx || iflag==1) {
	    *istop = 1;
	    *itp = blitp_1.itpst * 10 + 9;
	}
    }

/* All special points receive label: */

    lab1 = 0;
    if (*itp % 10 != 0) {
	++(*lab);
	lab1 = *lab;
    }

/* Compute maxima of solution components. */

    n2 = blicn_1.ndm;
    if (n2 > 7) {
	n2 = 7;
    }
    i__1 = n2;
    for (i = 1; i <= i__1; ++i) {
	itmp = i;
	umx[i - 1] = rmxups_(m1u, &itmp, &ups[ups_offset]);
/* L1: */
    }

/* Determine stability, and write output on units 7 and 8. */

    ibr1 = *ibr;
    ntot1 = *ntot;
    if ((blbcn_1.ips == 2 || blbcn_1.ips == 3 || blbcn_1.ips == 6 || 
	    blbcn_1.ips == 12 || blbcn_1.ips == 13) && abs(blcde_1.isw) != 2) 
	    {
	ibr1 = -(*ibr);
	if (bldet_1.nins == blbcn_1.ndim) {
	    ntot1 = -(*ntot);
	}
    }
/*     WRLINE CALLED HERE !!!!!!! */

    i_1 = blicn_1.ndm;
    for (i = 1; i <= i_1; ++i) {
        u0[i-1] = ups[i * ups_dim1 + 1];
	uhigh[i - 1] = rmxups_(m1u, &i, &ups[ups_offset]);
	ulow[i - 1] = rmnups_(m1u, &i, &ups[ups_offset]);
        ubar[i-1] =rintg_(m1u,&i,&ups[ups_offset],&dtm[1]);  
/* L11: */
    }
/*     ADD BIFURCATION */

    addbif_(&ibr1, &ntot1, itp, &lab1, &blicn_1.nfpar,
	    &blcrl_1.a, uhigh, ulow, u0, ubar, &blicn_1.ndm);

    wrline_(&ibr1, &ntot1, itp, &lab1, &blcrl_1.a, umx);

/* Write plotting and restart data on unit 8. */

    if (blbcn_1.irs != 0 && *ntot == 1) {
	return 0;
    }
    if (*itp % 10 != 0) {
	wrtbv8_(itp, ntot, lab, ibr, m1u, &ups[ups_offset], &udotps[
		udotps_offset], &tm[1], &dtm[1]);
    }

    return 0;
} /* stplbv_ */


/*     ---------- ------ */
/* Subroutine */ int wrtbv8_(itp, ntot, lab, ibr, m1u, ups, udotps, tm, dtm)
integer *itp, *ntot, *lab, *ibr, *m1u;
doublereal *ups, *udotps, *tm, *dtm;
{
    /* Format strings */
    static char fmt_101[] = "(11i5,20i3)";
    static char fmt_102[] = "(4x,1p7e18.10)";

    /* System generated locals */
    integer ups_dim1, ups_offset, udotps_dim1, udotps_offset, i__1, i__2, 
	    i__3;

    /* Builtin functions */
    integer s_wsfe(), do_fio(), e_wsfe();

    /* Local variables */
    static integer ntpl;
    static doublereal time0, time1;
    static integer i, j, k;
    static doublereal t;
    static integer k1, k2;
    extern /* Subroutine */ int autim0_(), autim1_();
    static doublereal rn;
    static integer nrowpr, nar, nrd;

    /* Fortran I/O blocks */
    static cilist io___192 = { 0, 8, 0, fmt_101, 0 };
    static cilist io___199 = { 0, 8, 0, fmt_102, 0 };
    static cilist io___201 = { 0, 8, 0, fmt_102, 0 };
    static cilist io___202 = { 0, 8, 0, fmt_102, 0 };
    static cilist io___203 = { 0, 8, 0, fmt_102, 0 };
    static cilist io___204 = { 0, 8, 0, fmt_102, 0 };
    static cilist io___205 = { 0, 8, 0, fmt_102, 0 };



/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Writes plotting and restart data on unit 8, viz.: */
/* (1) data identifying the corresponding point on unit 7, */
/* (2) the complete solution, */
/* (3) the direction of the branch. */

/* Specifically the following is written: */

/*  IBR   : The index of the branch. */
/*  NTOT  : The index of the point. */
/*  ITP   : The type of point (see STPLBV above). */
/*  LAB   : The label of the point. */
/*  NFPAR : The number of free parameters used in the computation. */
/*  ISW   : The value of ISW used in the computation. */
/*  NTPL  : The number of points in the time interval [0,1] for which */
/*          solution values are wriiten. */
/*  NAR   : The number of values written per point. */
/*          (NAR=NDIM+1, since T and U(i), i=1,..,NDIM are written). */
/*  NROWPR: The number of lines printed following the identifying line */
/*          and before the next data set or the end of the file. */
/*          (Used for quickly skipping a data set when searching). */
/*  NTST  : The number of time intervals used in the discretization. */
/*  NCOL  : The number of collocation points used. */
/*   ICP  : The indices of the free parameters in PAR(.). */

/*  Following the above described identifying line there are NTPL lines */

/* containing : */
/*     T , U-1(T) , U-2(T) , ... , U-NDIM(T), */
/* where NDIM is the dimension of the system of differential equations. */


/* Following this is a line containing */
/*    RL-dot(i) , i=1,NFPAR, */

/* and following this are NTPL lines each containing */
/*    U-dot-1(T), U-dot-2(T), ... , U-dot-NDIM(T). */

/* Finally the parameter values PAR(i) , i=1,NPAR, are written. */

/*  Above, RL-dot(.) and U-dot(.) specify the direction of the branch. */



/* Set initial time (for timing of this subroutine). */

    /* Parameter adjustments */
    --dtm;
    --tm;
    udotps_dim1 = *m1u;
    udotps_offset = udotps_dim1 + 1;
    udotps -= udotps_offset;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;

    /* Function Body */
    autim0_(&time0);

/* Write information identifying the solution : */

    bldim_1.ntstp1 = blcde_1.ntst + 1;
    ntpl = blcde_1.ncol * blcde_1.ntst + 1;
    nar = blbcn_1.ndim + 1;
    nrd = blbcn_1.ndim / 7 + 2 + blbcn_1.ndim / 8;
    nrowpr = nrd * (blcde_1.ncol * blcde_1.ntst + 1) + blicn_1.nfpar / 8 + 1 
	    + blicn_1.npar / 8 + 1;
    s_wsfe(&io___192);
    do_fio(&c__1, (char *)&(*ibr), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*ntot), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*itp), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&(*lab), (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&blicn_1.nfpar, (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&blcde_1.isw, (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&ntpl, (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&nar, (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&nrowpr, (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&blcde_1.ntst, (ftnlen)sizeof(integer));
    do_fio(&c__1, (char *)&blcde_1.ncol, (ftnlen)sizeof(integer));
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	do_fio(&c__1, (char *)&blbcn_1.icp[i - 1], (ftnlen)sizeof(integer));
    }
    e_wsfe();

/* Write the entire solution on unit 8 : */

    i__1 = blcde_1.ntst;
    for (j = 1; j <= i__1; ++j) {
	rn = blrcn_1.one / blcde_1.ncol;
	i__2 = blcde_1.ncol;
	for (i = 1; i <= i__2; ++i) {
	    k1 = (i - 1) * blbcn_1.ndim + 1;
	    k2 = i * blbcn_1.ndim;
	    t = tm[j] + (i - 1) * rn * dtm[j];
	    s_wsfe(&io___199);
	    do_fio(&c__1, (char *)&t, (ftnlen)sizeof(doublereal));
	    i__3 = k2;
	    for (k = k1; k <= i__3; ++k) {
		do_fio(&c__1, (char *)&ups[j + k * ups_dim1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_wsfe();
/* L1: */
	}
/* L2: */
    }
    s_wsfe(&io___201);
    do_fio(&c__1, (char *)&tm[bldim_1.ntstp1], (ftnlen)sizeof(doublereal));
    i__1 = blbcn_1.ndim;
    for (i = 1; i <= i__1; ++i) {
	do_fio(&c__1, (char *)&ups[bldim_1.ntstp1 + i * ups_dim1], (ftnlen)
		sizeof(doublereal));
    }
    e_wsfe();

/* Store the direction of the branch: */

    s_wsfe(&io___202);
    i__1 = blicn_1.nfpar;
    for (i = 1; i <= i__1; ++i) {
	do_fio(&c__1, (char *)&blcrl_1.rldot[i - 1], (ftnlen)sizeof(
		doublereal));
    }
    e_wsfe();
    i__1 = blcde_1.ntst;
    for (j = 1; j <= i__1; ++j) {
	i__2 = blcde_1.ncol;
	for (i = 1; i <= i__2; ++i) {
	    k1 = (i - 1) * blbcn_1.ndim + 1;
	    k2 = i * blbcn_1.ndim;
	    s_wsfe(&io___203);
	    i__3 = k2;
	    for (k = k1; k <= i__3; ++k) {
		do_fio(&c__1, (char *)&udotps[j + k * udotps_dim1], (ftnlen)
			sizeof(doublereal));
	    }
	    e_wsfe();
/* L3: */
	}
/* L4: */
    }
    s_wsfe(&io___204);
    i__1 = blbcn_1.ndim;
    for (k = 1; k <= i__1; ++k) {
	do_fio(&c__1, (char *)&udotps[bldim_1.ntstp1 + k * udotps_dim1], (
		ftnlen)sizeof(doublereal));
    }
    e_wsfe();

/* Write the parameter values. */

    s_wsfe(&io___205);
    i__1 = blicn_1.npar;
    for (i = 1; i <= i__1; ++i) {
	do_fio(&c__1, (char *)&blbcn_1.par[i - 1], (ftnlen)sizeof(doublereal))
		;
    }
    e_wsfe();

/* Determine the time spent in this subroutine. */

    autim1_(&time1);
    bltim_1.twr8 = bltim_1.twr8 + time1 - time0;

    return 0;
} /* wrtbv8_ */


/*     ---------- ------ */
/* Subroutine */ int wrtbv9_(nitps, ibr, ntot, m1u, ups, tm, dtm)
integer *nitps, *ibr, *ntot, *m1u;
doublereal *ups, *tm, *dtm;
{
    /* Format strings */
    static char fmt_101[] = "(\002 BRANCH \002,i2,\002 N=\002,i4,1x,\002IT\
=\002,i2,1x,\002A=\002,1pe15.8,1x,\002 PAR= \002,1p5e15.8,3(/,47x,1p5e15.8))";

    static char fmt_103[] = "(\002 UPS :\002)";
    static char fmt_102[] = "(1x,1p7e18.10)";

    /* System generated locals */
    integer ups_dim1, ups_offset, i__1, i__2, i__3;

    /* Builtin functions */
    double sqrt();
    integer s_wsfe(), do_fio(), e_wsfe();

    /* Local variables */
    static integer i, j, k;
    static doublereal t;
    static integer k1, k2;
    static doublereal rn;
    extern doublereal rnrmsq_(), rmnups_(), rmxups_();
    static integer iab;

    /* Fortran I/O blocks */
    static cilist io___208 = { 0, 9, 0, fmt_101, 0 };
    static cilist io___210 = { 0, 9, 0, fmt_103, 0 };
    static cilist io___216 = { 0, 9, 0, fmt_102, 0 };
    static cilist io___218 = { 0, 9, 0, fmt_102, 0 };



/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Writes additional output on unit 9. */



    /* Parameter adjustments */
    --dtm;
    --tm;
    ups_dim1 = *m1u;
    ups_offset = ups_dim1 + 1;
    ups -= ups_offset;

    /* Function Body */
    iab = abs(blcde_1.iplt);
    if (iab == 0 || iab > blbcn_1.ndim) {
	blcrl_1.a = sqrt((rnrmsq_(&blicn_1.ndm, m1u, &ups[ups_offset], &dtm[1]
		)));
    }
/* SGLE  IF(IAB.EQ.0.OR.IAB.GT.NDIM)A= SQRT(RNRMSQ(NDM,M1U,UPS,DTM)) */
    if (blcde_1.iplt > 0 && iab <= blbcn_1.ndim) {
	blcrl_1.a = rmxups_(m1u, &iab, &ups[ups_offset]);
    }
    if (blcde_1.iplt < 0 && iab <= blbcn_1.ndim) {
	blcrl_1.a = rmnups_(m1u, &iab, &ups[ups_offset]);
    }
    if (blmax_1.iid >= 2) {
	s_wsfe(&io___208);
	do_fio(&c__1, (char *)&(*ibr), (ftnlen)sizeof(integer));
	i__1 = *ntot + 1;
	do_fio(&c__1, (char *)&i__1, (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&(*nitps), (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&blcrl_1.a, (ftnlen)sizeof(doublereal));
	i__2 = blicn_1.nfpar;
	for (i = 1; i <= i__2; ++i) {
	    do_fio(&c__1, (char *)&blcrl_1.rl[i - 1], (ftnlen)sizeof(
		    doublereal));
	}
	e_wsfe();
    }

    if (! (blmax_1.iid >= 5)) {
	return 0;
    }

    s_wsfe(&io___210);
    e_wsfe();
    i__1 = blcde_1.ntst;
    for (j = 1; j <= i__1; ++j) {
	rn = blrcn_1.one / blcde_1.ncol;
	i__2 = blcde_1.ncol;
	for (i = 1; i <= i__2; ++i) {
	    t = tm[j] + (i - 1) * rn * dtm[j];
	    k1 = (i - 1) * blbcn_1.ndim + 1;
	    k2 = i * blbcn_1.ndim;
	    s_wsfe(&io___216);
	    do_fio(&c__1, (char *)&t, (ftnlen)sizeof(doublereal));
	    i__3 = k2;
	    for (k = k1; k <= i__3; ++k) {
		do_fio(&c__1, (char *)&ups[j + k * ups_dim1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_wsfe();
/* L1: */
	}
/* L2: */
    }

    bldim_1.ntstp1 = blcde_1.ntst + 1;
    s_wsfe(&io___218);
    do_fio(&c__1, (char *)&tm[bldim_1.ntstp1], (ftnlen)sizeof(doublereal));
    i__1 = blbcn_1.ndim;
    for (i = 1; i <= i__1; ++i) {
	do_fio(&c__1, (char *)&ups[bldim_1.ntstp1 + i * ups_dim1], (ftnlen)
		sizeof(doublereal));
    }
    e_wsfe();


    return 0;
} /* wrtbv9_ */


/* ------------------------------------------------------------------- */
/* ------------------------------------------------------------------- */
/*                    Linear Equation Solver */
/* ------------------------------------------------------------------- */
/* ------------------------------------------------------------------- */

/*     ---------- ---- */
/* Subroutine */ int brbd_(na, nra, nca, ma1, ma2, a, ncb, mb1, mb2, b, nrc, 
	mc1, c, md1, d, mfa1, fa, fc, wkdr, idb, ifst, ir, ic, iwkdr, nllv)
integer *na, *nra, *nca, *ma1, *ma2;
doublereal *a;
integer *ncb, *mb1, *mb2;
doublereal *b;
integer *nrc, *mc1;
doublereal *c;
integer *md1;
doublereal *d;
integer *mfa1;
doublereal *fa, *fc, *wkdr;
integer *idb, *ifst, *ir, *ic, *iwkdr, *nllv;
{
    /* Format strings */
    static char fmt_101[] = "(\002 SUBROUTINE BRBD : IFST=\002,i1,1x,\002NLL\
V=\002,i2)";

    /* System generated locals */
    integer a_dim1, a_dim2, a_offset, b_dim1, b_dim2, b_offset, c_dim1, 
	    c_offset, d_dim1, d_offset, fa_dim1, fa_offset, i__1, i__2;

    /* Builtin functions */
    /* Subroutine */ int s_stop();
    integer s_wsfe(), do_fio(), e_wsfe();

    /* Local variables */
    static integer i, j, lrhse, lnext;
    extern /* Subroutine */ int print1_();
    static integer lb, lc, le;
    extern /* Subroutine */ int print3_(), print2_();
    static integer ne, ls;
    extern logical erbrbd_();
    static integer lt;
    extern /* Subroutine */ int reduce_(), dimrge_(), bcksub_(), infpar_(), 
	    conpar_(), redrhs_(), conrhs_(), copycp_();
    static integer la1, la2;
    extern /* Subroutine */ int cpyrhs_();
    static integer lfa, lic, llc, lir, lxe, nov, nap2, nov2;

    /* Fortran I/O blocks */
    static cilist io___221 = { 0, 9, 0, fmt_101, 0 };



/* SGLE IMPLICIT REAL             (A-H,O-Z) */

/* Solves linear systems with matrix profile: */

/*     ----------------------------------------------- */
/*     !XXXXXXXXXX                                !XX! */
/*     !XXXXXXXXXX                                !XX! */
/*     !XXXXXXXXXX                                !XX! */
/*     !XXXXXXXXXX                                !XX! */
/*     !XXXXXXXXXX                                !XX! */
/*     !XXXXXXXXXX                                !XX! */
/*     !XXXXXXXXXX                                !XX! */
/*     !XXXXXXXXXX                                !XX! */
/*     !        XXXXXXXXXX                        !XX! */
/*     !        XXXXXXXXXX                        !XX! */
/*     !        XXXXXXXXXX                        !XX! */
/*     !        XXXXXXXXXX                        !XX! */
/*     !        XXXXXXXXXX                        !XX! */
/*     !        XXXXXXXXXX                        !XX! */
/*     !        XXXXXXXXXX                        !XX! */
/*     !        XXXXXXXXXX                        !XX! */
/*     !                XXXXXXXXXX                !XX! */
/*     !                XXXXXXXXXX                !XX! */
/*     !                XXXXXXXXXX                !XX! */
/*     !                XXXXXXXXXX                !XX! */
/*     !                XXXXXXXXXX                !XX! */
/*     !                XXXXXXXXXX                !XX! */
/*     !                XXXXXXXXXX                !XX! */
/*     !                XXXXXXXXXX                !XX! */
/*     !                        XXXXXXXXXX        !XX! */
/*     !                        XXXXXXXXXX        !XX! */
/*     !                        XXXXXXXXXX        !XX! */
/*     !                        XXXXXXXXXX        !XX! */
/*     !                        XXXXXXXXXX        !XX! */
/*     !                        XXXXXXXXXX        !XX! */
/*     !                        XXXXXXXXXX        !XX! */
/*     !                        XXXXXXXXXX        !XX! */
/*     !                                XXXXXXXXXX!XX! */
/*     !                                XXXXXXXXXX!XX! */
/*     !                                XXXXXXXXXX!XX! */
/*     !                                XXXXXXXXXX!XX! */
/*     !                                XXXXXXXXXX!XX! */
/*     !                                XXXXXXXXXX!XX! */
/*     !                                XXXXXXXXXX!XX! */
/*     !                                XXXXXXXXXX!XX! */
/*     ----------------------------------------------- */
/*     !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX!XX! */
/*     !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX!XX! */
/*     !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX!XX! */
/*     !XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX!XX! */
/*     ----------------------------------------------- */


/* partioned as */


/*      --------- */
/*      !     ! !   !    !   !    ! */
/*      !  A  !B!   ! XA !   ! FA ! */
/*      !     ! ! . !    ! = !    !   . */
/*      !-----!-!   !----!   !----! */
/*      !  C  !D!   ! XC !   ! FC ! */
/*      !     ! !   ! XC !   ! FC ! */
/*      --------- */


/* Input parameters : */

/*   NA    number of blocks in A, */
/*   NRA   number of rows in each block of A, */
/*   NCA   number of columns in each block of A, */
/*   MA1   first dimension of A from DIMENSION statement, */
/*   MA2   second dimension of A from DIMENSION statement, */
/*   A     the matrix in the schematic representation above, */

/*   NCB   number of columns in each block of B, */
/*         (note that B is also three dimensional), */
/*   MB1   first dimension of B from DIMENSION statement, */
/*   MB2   second dimension of B from DIMENSION statement, */
/*   B     the matrix in the schema above, */

/*   NRC   the number of rows of the two dimensional matrix C, */
/*   MC1   the first dimension of C from DIMENSION statement, */
/*   C     the matrix C in the schema above, */

/*   MD1    the first dimension of D from DIMENSION statement, */
/*   D      the matrix D above, */

/*   MFA1  the first dimension of FA from DIMENSION statement, */
/*   FA     part of the right hand side vector, */
/*          (note that FA is also two dimensional), */
/*   FC     part of the right hand side vector. */

/*   WKDR: A one dimensional array used as workspace. */
/*          This array should be dimensioned at least */

/*   (4*NOV+NCB+NRC+1)*NOV*NA+(NRC+NOV)**2+2*(NRC+2*NOV)+NRC*NOV */

/*          with NOV defined by */

/*                  NA*(NCA-NRA)+NCB-NRC */
/*            NOV = -------------------- . */
/*                        NA-1 */

/*  IWKDR: Integer workspace array of dimension at least 3*NOV*(NA-1)+NA. 
*/

/*   IFST   = 1 on first call, */
/*          = 0  on subsequent calls with the same right hand side. */
/*              (WKDR,IWKDR should not be modified between such calls). */


/*   IDB   = 0    no debug output, */
/*         = 1,2  debug output on unit 9, */
/*         = 3    print residuals for test problem (see PRINT2), */

/*   IR, IC: Two integer arrays of dimension at least NRC+NOV. */

/*   NLLV : If NLLV>0 then the system is assumed to have a NLLV- */
/*          dimensional nullspace. */
/*          In this case a null vector will be returned. */
/*          If NLLVC = -1 then the system will be solved with zero right 
*/
/*          hand side, except for the last equation, for which the right 
*/
/*          hand side entry will be set to 1 (i.e., the last entry of FC 
*/
/*          will be set to 1, otherwise FA and FC are zero). */
/*          If the linear system is the same as in the preceding call */
/*          then IFST=0 may be used even if NLLV is nonzero. */

/* Returned values : */

/*   FA     Part of solution vector corresponding to XA in the diagram. */

/*   FC     Part of solution vector corresponding to XC in the diagram. */


/* Note : The number of columns of overlap for every two consecutive */
/*        blocks should be equal to the number NOV defined above. */




/* Check for consistency of data. */

    /* Parameter adjustments */
    --iwkdr;
    --ic;
    --ir;
    --wkdr;
    --fc;
    fa_dim1 = *mfa1;
    fa_offset = fa_dim1 + 1;
    fa -= fa_offset;
    d_dim1 = *md1;
    d_offset = d_dim1 + 1;
    d -= d_offset;
    c_dim1 = *mc1;
    c_offset = c_dim1 + 1;
    c -= c_offset;
    b_dim1 = *mb1;
    b_dim2 = *mb2;
    b_offset = b_dim1 * (b_dim2 + 1) + 1;
    b -= b_offset;
    a_dim1 = *ma1;
    a_dim2 = *ma2;
    a_offset = a_dim1 * (a_dim2 + 1) + 1;
    a -= a_offset;

    /* Function Body */
    if (erbrbd_(na, nra, nca, ncb, nrc)) {
	s_stop("", 0L);
    }
    nov = (*na * (*nca - *nra) + *ncb - *nrc) / (*na - 1);
    nov2 = nov << 1;

/* Print debug output. */

    if (*idb >= 2) {
	s_wsfe(&io___221);
	do_fio(&c__1, (char *)&(*ifst), (ftnlen)sizeof(integer));
	do_fio(&c__1, (char *)&(*nllv), (ftnlen)sizeof(integer));
	e_wsfe();
	print1_(&nov, na, nra, nca, ncb, nrc, ma1, ma2, &a[a_offset], mb1, 
		mb2, &b[b_offset], mc1, &c[c_offset], md1, &d[d_offset], mfa1,
		 &fa[fa_offset], &fc[1]);
    }

/* Eliminate local variables by "Condensation of Parameters". */

    if (*ifst != 0) {
	++blcnt_1.ndecom;
	conpar_(&nov, na, nra, nca, ma1, ma2, &a[a_offset], ncb, mb1, mb2, &b[
		b_offset], nrc, mc1, &c[c_offset], md1, &d[d_offset], &wkdr[1]
		, &iwkdr[1]);
/*       PRINT DEBUG OUTPUT */
	if (*idb >= 2) {
	    print1_(&nov, na, nra, nca, ncb, nrc, ma1, ma2, &a[a_offset], mb1,
		     mb2, &b[b_offset], mc1, &c[c_offset], md1, &d[d_offset], 
		    mfa1, &fa[fa_offset], &fc[1]);
	}
    }

/* Allocate workspace. */

    la1 = 1;
/* Computing 2nd power */
    i__1 = nov;
    la2 = la1 + i__1 * i__1 * *na;
/* Computing 2nd power */
    i__1 = nov;
    lb = la2 + i__1 * i__1 * *na;
    lc = lb + nov * *ncb * *na;
    lfa = lc + *nrc * nov * (*na + 1);
    ls = lfa + nov * (*na + 2);
/* Computing 2nd power */
    i__1 = nov;
    le = ls + i__1 * i__1 * *na;
/* Computing 2nd power */
    i__1 = *nrc + nov;
    lrhse = le + i__1 * i__1;
    lxe = lrhse + *nrc + nov;
    lt = lxe + nov + *nrc;
/* Computing 2nd power */
    i__1 = nov;
    lnext = lt + i__1 * i__1 * *na;

    lir = 1;
    lic = lir + (nov << 1) * (*na - 1);
    llc = lic + nov * (*na - 1);
    lnext = llc + *na;

    blbnd_1.nam1 = *na - 1;
    blbnd_1.nap1 = *na + 1;
    nap2 = *na + 2;
    blbnd_1.nxe = nov + *nrc;

/* Copy the reduced system generated by CONPAR into the workspace area */
/* for further processing by REDUCE. */

    if (*ifst != 0) {
	copycp_(na, &nov, nra, nca, ma1, ma2, &a[a_offset], ncb, mb1, mb2, &b[
		b_offset], nrc, mc1, &c[c_offset], &wkdr[la1], &wkdr[la2], &
		wkdr[lb], &wkdr[lc]);
/*       Reduction of the system. */
	reduce_(na, &nov, &nov2, ncb, nrc, md1, &d[d_offset], &wkdr[la1], &
		wkdr[la2], &wkdr[lb], &wkdr[lc], &wkdr[ls], &wkdr[lt], &iwkdr[
		lir], &iwkdr[lic]);
    }

/* Condensation of the right hand side following CONPAR. */

    if (*nllv == 0) {
	conrhs_(&nov, na, nra, nca, ma1, ma2, &a[a_offset], nrc, mc1, &c[
		c_offset], mfa1, &fa[fa_offset], &fc[1], &iwkdr[llc]);
/*       Copy the reduced right hand side into workspace. */
	cpyrhs_(na, nra, &nov, mfa1, &fa[fa_offset], &wkdr[lfa]);
/*       Print debug output. */
	if (*idb >= 2) {
	    print3_(na, &nov, &wkdr[lfa]);
	}
/*       Reduction of the right hand side following REDUCE. */
	redrhs_(na, &nov, &nov2, nrc, &wkdr[la1], &wkdr[la2], &wkdr[lc], &
		wkdr[lfa], &fc[1], &iwkdr[lir]);
    } else {
/*        Set right hand sides to zero. */
	i__1 = *na;
	for (i = 1; i <= i__1; ++i) {
	    i__2 = *nra;
	    for (j = 1; j <= i__2; ++j) {
		fa[i + j * fa_dim1] = blrcn_1.zero;
/* L1: */
	    }
/* L2: */
	}
	i__1 = *nrc;
	for (i = 1; i <= i__1; ++i) {
	    fc[i] = blrcn_1.zero;
/* L3: */
	}
	i__1 = ls - 1;
	for (i = lfa; i <= i__1; ++i) {
	    wkdr[i] = blrcn_1.zero;
/* L4: */
	}
    }

/* Print debug output */

    if (*idb >= 2) {
	print2_(idb, na, &nov, ncb, nrc, &wkdr[ls], &wkdr[la1], &wkdr[la2], &
		wkdr[lt], &wkdr[lb], &wkdr[lc], &wkdr[lfa]);
	print3_(na, &nov, &wkdr[lfa]);
    }

    ne = *nrc + nov;

/* Solve the system generated by REDUCE */
/* by Gauss elimination with complete pivoting. */

    dimrge_(na, &nov, nca, &wkdr[ls], &wkdr[la2], ncb, &wkdr[lb], &wkdr[lc], 
	    nrc, mc1, &c[c_offset], md1, &d[d_offset], &wkdr[lfa], &fc[1], &
	    ne, &wkdr[le], &wkdr[lrhse], &wkdr[lxe], idb, &ir[1], &ic[1], 
	    nllv);

/* Backsubstitution in the reduction process. */

    bcksub_(&nov, &nov2, na, ncb, &wkdr[ls], &wkdr[la1], &wkdr[la2], &wkdr[lt]
	    , &wkdr[lb], &wkdr[lfa], &wkdr[lxe], &fc[1], &iwkdr[lir], &iwkdr[
	    lic]);

/* Backsubstitution in the condensation of parameters process. */

    ++blcnt_1.nbcksb;
    infpar_(na, &nov, nra, nca, ma1, ma2, &a[a_offset], ncb, mb1, mb2, &b[
	    b_offset], mfa1, &fa[fa_offset], &bldim_1.ndrhs, &fc[1], &nap2, &
	    wkdr[lfa]);


    return 0;
} /* brbd_ */


/*     ------- -------- ------ */
logical erbrbd_(na, nra, nca, ncb, nrc)
integer *na, *nra, *nca, *ncb, *nrc;
{
    /* Format strings */
    static char fmt_101[] = "(\002 ERROR IN DATA IN SUBROUTINE -BRBD-\002,/\
,\002 (LINEAR EQUATION SOLVER)\002)";

    /* System generated locals */
    logical ret_val;

    /* Builtin functions */
    integer s_wsfe(), e_wsfe();

    /* Local variables */
    static integer id, in, nex, nov;

    /* Fortran I/O blocks */
    static cilist io___242 = { 0, 9, 0, fmt_101, 0 };
    static cilist io___245 = { 0, 9, 0, fmt_101, 0 };



/* Checks correctness of dimensions. */

    ret_val = FALSE_;

    in = *na * (*nca - *nra) + *ncb - *nrc;
    id = *na - 1;
    if (in % id != 0) {
	s_wsfe(&io___242);
	e_wsfe();
	ret_val = TRUE_;
	return ret_val;
    }

    nov = in / id;
    nex = *nca - (nov << 1);
    if (nex >= 0) {
	goto L1;
    }

    s_wsfe(&io___245);
    e_wsfe();
    ret_val = TRUE_;
    return ret_val;

L1:
    return ret_val;
} /* erbrbd_ */


/*     ---------- ------ */
/* Subroutine */ int copycp_(na, nov, nra, nca, ma1, ma2, a, ncb, mb1, mb2, b,
	 nrc, mc1, c, a1, a2, bc, cc)
integer *na, *nov, *nra, *nca, *ma1, *ma2;
doublereal *a;
integer *ncb, *mb1, *mb2;
doublereal *b;
integer *nrc, *mc1;
doublereal *c, *a1, *a2, *bc, *cc;
{
    /* System generated locals */
    integer a_dim1, a_dim2, a_offset, b_dim1, b_dim2, b_offset, c_dim1, 
	    c_offset, a1_dim1, a1_dim2, a1_offset, a2_dim1, a2_dim2, 
	    a2_offset, bc_dim1, bc_dim2, bc_offset, cc_dim1, cc_dim2, 
	    cc_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer i, ic, ir;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */



/* Copies the condensed system generated by CONPAR into the workspace. */

    /* Parameter adjustments */
    cc_dim1 = *nrc;
    cc_dim2 = *nov;
    cc_offset = cc_dim1 * (cc_dim2 + 1) + 1;
    cc -= cc_offset;
    bc_dim1 = *nov;
    bc_dim2 = *ncb;
    bc_offset = bc_dim1 * (bc_dim2 + 1) + 1;
    bc -= bc_offset;
    a2_dim1 = *nov;
    a2_dim2 = *nov;
    a2_offset = a2_dim1 * (a2_dim2 + 1) + 1;
    a2 -= a2_offset;
    a1_dim1 = *nov;
    a1_dim2 = *nov;
    a1_offset = a1_dim1 * (a1_dim2 + 1) + 1;
    a1 -= a1_offset;
    c_dim1 = *mc1;
    c_offset = c_dim1 + 1;
    c -= c_offset;
    b_dim1 = *mb1;
    b_dim2 = *mb2;
    b_offset = b_dim1 * (b_dim2 + 1) + 1;
    b -= b_offset;
    a_dim1 = *ma1;
    a_dim2 = *ma2;
    a_offset = a_dim1 * (a_dim2 + 1) + 1;
    a -= a_offset;

    /* Function Body */
    i__1 = *na;
    for (i = 1; i <= i__1; ++i) {

	i__2 = *nov;
	for (ir = 1; ir <= i__2; ++ir) {
	    i__3 = *nov;
	    for (ic = 1; ic <= i__3; ++ic) {
		a1[ir + (ic + i * a1_dim2) * a1_dim1] = a[i + (*nra - *nov + 
			ir + ic * a_dim2) * a_dim1];
		a2[ir + (ic + i * a2_dim2) * a2_dim1] = a[i + (*nra - *nov + 
			ir + (*nca - *nov + ic) * a_dim2) * a_dim1];
/* L1: */
	    }

	    i__3 = *ncb;
	    for (ic = 1; ic <= i__3; ++ic) {
		bc[ir + (ic + i * bc_dim2) * bc_dim1] = b[i + (*nra - *nov + 
			ir + ic * b_dim2) * b_dim1];
/* L2: */
	    }

/* L3: */
	}

/* L4: */
    }

    i__1 = blbnd_1.nap1;
    for (i = 1; i <= i__1; ++i) {
	i__2 = *nrc;
	for (ir = 1; ir <= i__2; ++ir) {
	    i__3 = *nov;
	    for (ic = 1; ic <= i__3; ++ic) {
		cc[ir + (ic + i * cc_dim2) * cc_dim1] = c[(i - 1) * (*nca - *
			nov) + ic + ir * c_dim1];
/* L5: */
	    }
/* L6: */
	}
/* L7: */
    }

    return 0;
} /* copycp_ */


/*     ---------- ------ */
/* Subroutine */ int reduce_(na, nov, nov2, ncb, nrc, md1, d, a1, a2, b, c, s,
	 t, ipr, ipc)
integer *na, *nov, *nov2, *ncb, *nrc, *md1;
doublereal *d, *a1, *a2, *b, *c, *s, *t;
integer *ipr, *ipc;
{
    /* System generated locals */
    integer a1_dim1, a1_dim2, a1_offset, a2_dim1, a2_dim2, a2_offset, s_dim1, 
	    s_dim2, s_offset, t_dim1, t_dim2, t_offset, ipr_dim1, ipr_offset, 
	    ipc_dim1, ipc_offset, b_dim1, b_dim2, b_offset, c_dim1, c_dim2, 
	    c_offset, d_dim1, d_offset, i__1, i__2, i__3, i__4;
    doublereal d__1;

    /* Local variables */
    static doublereal daba2, daba1;
    static integer itmp;
    static doublereal time0, time1, rmxa1, rmxa2;
    static integer novm1, i, k, i1, i2;
    extern /* Subroutine */ int autim0_(), autim1_();
    static integer ic, ir;
    static doublereal rm, tmp;
    static integer ica1, ica2, icp1, ira2, ira1;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */



/* Set initial time (For timing of this subroutine). */

    /* Parameter adjustments */
    ipc_dim1 = *nov;
    ipc_offset = ipc_dim1 + 1;
    ipc -= ipc_offset;
    ipr_dim1 = *nov2;
    ipr_offset = ipr_dim1 + 1;
    ipr -= ipr_offset;
    t_dim1 = *nov;
    t_dim2 = *nov;
    t_offset = t_dim1 * (t_dim2 + 1) + 1;
    t -= t_offset;
    s_dim1 = *nov;
    s_dim2 = *nov;
    s_offset = s_dim1 * (s_dim2 + 1) + 1;
    s -= s_offset;
    c_dim1 = *nrc;
    c_dim2 = *nov;
    c_offset = c_dim1 * (c_dim2 + 1) + 1;
    c -= c_offset;
    b_dim1 = *nov;
    b_dim2 = *ncb;
    b_offset = b_dim1 * (b_dim2 + 1) + 1;
    b -= b_offset;
    a2_dim1 = *nov;
    a2_dim2 = *nov;
    a2_offset = a2_dim1 * (a2_dim2 + 1) + 1;
    a2 -= a2_offset;
    a1_dim1 = *nov;
    a1_dim2 = *nov;
    a1_offset = a1_dim1 * (a1_dim2 + 1) + 1;
    a1 -= a1_offset;
    d_dim1 = *md1;
    d_offset = d_dim1 + 1;
    d -= d_offset;

    /* Function Body */
    autim0_(&time0);

    novm1 = *nov - 1;
    blbnd_1.nam1 = *na - 1;

/* CLEAR S AND T */

    i__1 = *na;
    for (i1 = 1; i1 <= i__1; ++i1) {
	i__2 = *nov;
	for (ir = 1; ir <= i__2; ++ir) {
	    i__3 = *nov;
	    for (ic = 1; ic <= i__3; ++ic) {
		t[ir + (ic + i1 * t_dim2) * t_dim1] = blrcn_1.zero;
		s[ir + (ic + i1 * s_dim2) * s_dim1] = blrcn_1.zero;
/* L1: */
	    }
/* L2: */
	}
/* L3: */
    }

/* Equate first block of S with first block of A1. */

    i__1 = *nov;
    for (ir = 1; ir <= i__1; ++ir) {
	i__2 = *nov;
	for (ic = 1; ic <= i__2; ++ic) {
	    s[ir + (ic + s_dim2) * s_dim1] = a1[ir + (ic + a1_dim2) * a1_dim1]
		    ;
/* L4: */
	}
/* L5: */
    }

/* Initialize pivot arrays. */

    i__1 = blbnd_1.nam1;
    for (i1 = 1; i1 <= i__1; ++i1) {
	i__2 = *nov;
	for (i = 1; i <= i__2; ++i) {
	    ipc[i + i1 * ipc_dim1] = i;
/* L6: */
	}
/* L7: */
    }

    i__1 = blbnd_1.nam1;
    for (i1 = 1; i1 <= i__1; ++i1) {

	i2 = i1 + 1;

	i__2 = *nov;
	for (ic = 1; ic <= i__2; ++ic) {
	    icp1 = ic + 1;

/*         Pivoting : */

	    rmxa2 = blrcn_1.zero;
	    i__3 = *nov;
	    for (ir = ic; ir <= i__3; ++ir) {
		i__4 = ic;
		for (k = ic; k <= i__4; ++k) {
		    daba2 = (d__1 = a2[ir + (k + i1 * a2_dim2) * a2_dim1], 
			    abs(d__1));
/* SGLE          DABA2= ABS(A2(IR,K,I1)) */
		    if (daba2 > rmxa2) {
			ira2 = ir;
			ica2 = k;
			rmxa2 = daba2;
		    }
/* L8: */
		}
/* L9: */
	    }

	    rmxa1 = blrcn_1.zero;
	    i__3 = *nov;
	    for (ir = 1; ir <= i__3; ++ir) {
		i__4 = ic;
		for (k = ic; k <= i__4; ++k) {
		    daba1 = (d__1 = a1[ir + (k + i2 * a1_dim2) * a1_dim1], 
			    abs(d__1));
/* SGLE          DABA1= ABS(A1(IR,K,I2)) */
		    if (daba1 > rmxa1) {
			ira1 = ir;
			ica1 = k;
			rmxa1 = daba1;
		    }
/* L10: */
		}
/* L11: */
	    }

	    if (rmxa2 > rmxa1) {
		ipr[ic + i1 * ipr_dim1] = ira2;
		itmp = ipc[ic + i1 * ipc_dim1];
		ipc[ic + i1 * ipc_dim1] = ipc[ica2 + i1 * ipc_dim1];
		ipc[ica2 + i1 * ipc_dim1] = itmp;
		i__3 = *nov;
		for (k = 1; k <= i__3; ++k) {
		    if (k >= ic) {
			tmp = a2[ic + (k + i1 * a2_dim2) * a2_dim1];
			a2[ic + (k + i1 * a2_dim2) * a2_dim1] = a2[ira2 + (k 
				+ i1 * a2_dim2) * a2_dim1];
			a2[ira2 + (k + i1 * a2_dim2) * a2_dim1] = tmp;
		    }
		    tmp = s[ic + (k + i1 * s_dim2) * s_dim1];
		    s[ic + (k + i1 * s_dim2) * s_dim1] = s[ira2 + (k + i1 * 
			    s_dim2) * s_dim1];
		    s[ira2 + (k + i1 * s_dim2) * s_dim1] = tmp;
		    tmp = t[ic + (k + i1 * t_dim2) * t_dim1];
		    t[ic + (k + i1 * t_dim2) * t_dim1] = t[ira2 + (k + i1 * 
			    t_dim2) * t_dim1];
		    t[ira2 + (k + i1 * t_dim2) * t_dim1] = tmp;
/* L12: */
		}
		i__3 = *ncb;
		for (k = 1; k <= i__3; ++k) {
		    tmp = b[ic + (k + i1 * b_dim2) * b_dim1];
		    b[ic + (k + i1 * b_dim2) * b_dim1] = b[ira2 + (k + i1 * 
			    b_dim2) * b_dim1];
		    b[ira2 + (k + i1 * b_dim2) * b_dim1] = tmp;
/* L13: */
		}

	    } else {

		ipr[ic + i1 * ipr_dim1] = ira1 + *nov;
		itmp = ipc[ic + i1 * ipc_dim1];
		ipc[ic + i1 * ipc_dim1] = ipc[ica1 + i1 * ipc_dim1];
		ipc[ica1 + i1 * ipc_dim1] = itmp;
		i__3 = *nov;
		for (k = 1; k <= i__3; ++k) {
		    if (k >= ic) {
			tmp = a2[ic + (k + i1 * a2_dim2) * a2_dim1];
			a2[ic + (k + i1 * a2_dim2) * a2_dim1] = a1[ira1 + (k 
				+ i2 * a1_dim2) * a1_dim1];
			a1[ira1 + (k + i2 * a1_dim2) * a1_dim1] = tmp;
		    }
		    tmp = s[ic + (k + i1 * s_dim2) * s_dim1];
		    s[ic + (k + i1 * s_dim2) * s_dim1] = s[ira1 + (k + i2 * 
			    s_dim2) * s_dim1];
		    s[ira1 + (k + i2 * s_dim2) * s_dim1] = tmp;
		    tmp = t[ic + (k + i1 * t_dim2) * t_dim1];
		    t[ic + (k + i1 * t_dim2) * t_dim1] = a2[ira1 + (k + i2 * 
			    a2_dim2) * a2_dim1];
		    a2[ira1 + (k + i2 * a2_dim2) * a2_dim1] = tmp;
/* L14: */
		}
		i__3 = *ncb;
		for (k = 1; k <= i__3; ++k) {
		    tmp = b[ic + (k + i1 * b_dim2) * b_dim1];
		    b[ic + (k + i1 * b_dim2) * b_dim1] = b[ira1 + (k + i2 * 
			    b_dim2) * b_dim1];
		    b[ira1 + (k + i2 * b_dim2) * b_dim1] = tmp;
/* L15: */
		}


	    }

/*         End of pivoting. */

	    if (ic != *nov) {
		i__3 = *nov;
		for (ir = icp1; ir <= i__3; ++ir) {
		    rm = a2[ir + (ic + i1 * a2_dim2) * a2_dim1] / a2[ic + (ic 
			    + i1 * a2_dim2) * a2_dim1];
		    a2[ir + (ic + i1 * a2_dim2) * a2_dim1] = rm;
		    i__4 = *nov;
		    for (k = icp1; k <= i__4; ++k) {
			a2[ir + (k + i1 * a2_dim2) * a2_dim1] -= rm * a2[ic + 
				(k + i1 * a2_dim2) * a2_dim1];
/* L16: */
		    }
		    i__4 = *nov;
		    for (k = 1; k <= i__4; ++k) {
			s[ir + (k + i1 * s_dim2) * s_dim1] -= rm * s[ic + (k 
				+ i1 * s_dim2) * s_dim1];
			t[ir + (k + i1 * t_dim2) * t_dim1] -= rm * t[ic + (k 
				+ i1 * t_dim2) * t_dim1];
/* L17: */
		    }
		    i__4 = *ncb;
		    for (k = 1; k <= i__4; ++k) {
			b[ir + (k + i1 * b_dim2) * b_dim1] -= rm * b[ic + (k 
				+ i1 * b_dim2) * b_dim1];
/* L18: */
		    }
/* L19: */
		}
	    }

	    i__3 = *nov;
	    for (ir = 1; ir <= i__3; ++ir) {
		rm = a1[ir + (ic + i2 * a1_dim2) * a1_dim1] / a2[ic + (ic + 
			i1 * a2_dim2) * a2_dim1];
		a1[ir + (ic + i2 * a1_dim2) * a1_dim1] = rm;
		if (icp1 > *nov) {
		    goto L21;
		}
		i__4 = *nov;
		for (k = icp1; k <= i__4; ++k) {
		    a1[ir + (k + i2 * a1_dim2) * a1_dim1] -= rm * a2[ic + (k 
			    + i1 * a2_dim2) * a2_dim1];
/* L20: */
		}
L21:
		i__4 = *nov;
		for (k = 1; k <= i__4; ++k) {
		    s[ir + (k + i2 * s_dim2) * s_dim1] -= rm * s[ic + (k + i1 
			    * s_dim2) * s_dim1];
		    a2[ir + (k + i2 * a2_dim2) * a2_dim1] -= rm * t[ic + (k + 
			    i1 * t_dim2) * t_dim1];
/* L22: */
		}
		i__4 = *ncb;
		for (k = 1; k <= i__4; ++k) {
		    b[ir + (k + i2 * b_dim2) * b_dim1] -= rm * b[ic + (k + i1 
			    * b_dim2) * b_dim1];
/* L23: */
		}
/* L24: */
	    }

/* L25: */
	}

	i__2 = *nov;
	for (ic = 1; ic <= i__2; ++ic) {
	    icp1 = ic + 1;
	    i__3 = *nrc;
	    for (ir = blcde_1.nbc + 1; ir <= i__3; ++ir) {
		rm = c[ir + (ic + i2 * c_dim2) * c_dim1] / a2[ic + (ic + i1 * 
			a2_dim2) * a2_dim1];
		c[ir + (ic + i2 * c_dim2) * c_dim1] = rm;
		if (icp1 > *nov) {
		    goto L27;
		}
		i__4 = *nov;
		for (k = icp1; k <= i__4; ++k) {
		    c[ir + (k + i2 * c_dim2) * c_dim1] -= rm * a2[ic + (k + 
			    i1 * a2_dim2) * a2_dim1];
/* L26: */
		}
L27:
		i__4 = *nov;
		for (k = 1; k <= i__4; ++k) {
		    c[ir + (k + c_dim2) * c_dim1] -= rm * s[ic + (k + i1 * 
			    s_dim2) * s_dim1];
		    c[ir + (k + (i2 + 1) * c_dim2) * c_dim1] -= rm * t[ic + (
			    k + i1 * t_dim2) * t_dim1];
/* L28: */
		}
		i__4 = *ncb;
		for (k = 1; k <= i__4; ++k) {
		    d[ir + k * d_dim1] -= rm * b[ic + (k + i1 * b_dim2) * 
			    b_dim1];
/* L29: */
		}
/* L30: */
	    }
/* L31: */
	}

/* L32: */
    }

/* Determine the time spent in this subroutine. */

    autim1_(&time1);
    bltim_1.treduc = bltim_1.treduc + time1 - time0;

    return 0;
} /* reduce_ */


/*     ---------- ------ */
/* Subroutine */ int cpyrhs_(na, nra, nov, mfa1, fa, fac)
integer *na, *nra, *nov, *mfa1;
doublereal *fa, *fac;
{
    /* System generated locals */
    integer fa_dim1, fa_offset, fac_dim1, fac_offset, i__1, i__2;

    /* Local variables */
    static integer i, ir;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */



    /* Parameter adjustments */
    fac_dim1 = *nov;
    fac_offset = fac_dim1 + 1;
    fac -= fac_offset;
    fa_dim1 = *mfa1;
    fa_offset = fa_dim1 + 1;
    fa -= fa_offset;

    /* Function Body */
    i__1 = *na;
    for (i = 1; i <= i__1; ++i) {
	i__2 = *nov;
	for (ir = 1; ir <= i__2; ++ir) {
	    fac[ir + i * fac_dim1] = fa[i + (*nra - *nov + ir) * fa_dim1];
/* L1: */
	}
/* L2: */
    }

    return 0;
} /* cpyrhs_ */


/*     ---------- ------ */
/* Subroutine */ int redrhs_(na, nov, nov2, nrc, a1, a2, c, fa, fc, ipr)
integer *na, *nov, *nov2, *nrc;
doublereal *a1, *a2, *c, *fa, *fc;
integer *ipr;
{
    /* System generated locals */
    integer a1_dim1, a1_dim2, a1_offset, a2_dim1, a2_dim2, a2_offset, c_dim1, 
	    c_dim2, c_offset, fa_dim1, fa_offset, ipr_dim1, ipr_offset, i__1, 
	    i__2, i__3;

    /* Local variables */
    static integer novm1, i1, i2, ic, ir;
    static doublereal rm, tmp;
    static integer icp1;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */



    /* Parameter adjustments */
    ipr_dim1 = *nov2;
    ipr_offset = ipr_dim1 + 1;
    ipr -= ipr_offset;
    --fc;
    fa_dim1 = *nov;
    fa_offset = fa_dim1 + 1;
    fa -= fa_offset;
    c_dim1 = *nrc;
    c_dim2 = *nov;
    c_offset = c_dim1 * (c_dim2 + 1) + 1;
    c -= c_offset;
    a2_dim1 = *nov;
    a2_dim2 = *nov;
    a2_offset = a2_dim1 * (a2_dim2 + 1) + 1;
    a2 -= a2_offset;
    a1_dim1 = *nov;
    a1_dim2 = *nov;
    a1_offset = a1_dim1 * (a1_dim2 + 1) + 1;
    a1 -= a1_offset;

    /* Function Body */
    novm1 = *nov - 1;

    i__1 = blbnd_1.nam1;
    for (i1 = 1; i1 <= i__1; ++i1) {
	i2 = i1 + 1;

	i__2 = *nov;
	for (ic = 1; ic <= i__2; ++ic) {
	    icp1 = ic + 1;

/* Interchanges. */

	    if (ipr[ic + i1 * ipr_dim1] <= *nov) {
		tmp = fa[ic + i1 * fa_dim1];
		fa[ic + i1 * fa_dim1] = fa[ipr[ic + i1 * ipr_dim1] + i1 * 
			fa_dim1];
		fa[ipr[ic + i1 * ipr_dim1] + i1 * fa_dim1] = tmp;
	    } else {
		tmp = fa[ic + i1 * fa_dim1];
		fa[ic + i1 * fa_dim1] = fa[ipr[ic + i1 * ipr_dim1] - *nov + 
			i2 * fa_dim1];
		fa[ipr[ic + i1 * ipr_dim1] - *nov + i2 * fa_dim1] = tmp;
	    }

/* End interchanges. */

	    if (ic != *nov) {
		i__3 = *nov;
		for (ir = icp1; ir <= i__3; ++ir) {
		    rm = a2[ir + (ic + i1 * a2_dim2) * a2_dim1];
		    fa[ir + i1 * fa_dim1] -= rm * fa[ic + i1 * fa_dim1];
/* L1: */
		}
	    }

	    i__3 = *nov;
	    for (ir = 1; ir <= i__3; ++ir) {
		rm = a1[ir + (ic + i2 * a1_dim2) * a1_dim1];
		fa[ir + i2 * fa_dim1] -= rm * fa[ic + i1 * fa_dim1];
/* L2: */
	    }
/* L3: */
	}

	i__2 = *nov;
	for (ic = 1; ic <= i__2; ++ic) {
	    i__3 = *nrc;
	    for (ir = 1; ir <= i__3; ++ir) {
		rm = c[ir + (ic + i2 * c_dim2) * c_dim1];
		fc[ir] -= rm * fa[ic + i1 * fa_dim1];
/* L4: */
	    }
/* L5: */
	}

/* L6: */
    }

    return 0;
} /* redrhs_ */


/*     ---------- ------ */
/* Subroutine */ int dimrge_(na, nov, nca, s, a2, ncb, b, cc, nrc, mc1, c, 
	md1, d, fa, fc, ne, e, rhse, xe, idb, ir, ic, nllv)
integer *na, *nov, *nca;
doublereal *s, *a2;
integer *ncb;
doublereal *b, *cc;
integer *nrc, *mc1;
doublereal *c;
integer *md1;
doublereal *d, *fa, *fc;
integer *ne;
doublereal *e, *rhse, *xe;
integer *idb, *ir, *ic, *nllv;
{
    /* Format strings */
    static char fmt_101[] = "(\002 REDUCED SYSTEM:\002)";
    static char fmt_102[] = "(1x,11e10.3)";
    static char fmt_103[] = "(\002 SOLUTION VECTOR :\002,/,10e10.3)";

    /* System generated locals */
    integer s_dim1, s_dim2, s_offset, a2_dim1, a2_dim2, a2_offset, b_dim1, 
	    b_dim2, b_offset, cc_dim1, cc_dim2, cc_offset, fa_dim1, fa_offset,
	     c_dim1, c_offset, d_dim1, d_offset, e_dim1, e_offset, i__1, i__2;


    /* Builtin functions */
    integer s_wsfe(), e_wsfe(), do_fio();

    /* Local variables */
    extern /* Subroutine */ int nlvc_();
    static integer i, j;
    extern /* Subroutine */ int ge_();

    /* Fortran I/O blocks */
    static cilist io___282 = { 0, 9, 0, fmt_101, 0 };
    static cilist io___283 = { 0, 9, 0, fmt_102, 0 };
    static cilist io___284 = { 0, 9, 0, fmt_103, 0 };



/* SGLE IMPLICIT REAL             (A-H,O-Z) */



    /* Parameter adjustments */
    --ic;
    --ir;
    --xe;
    --rhse;
    e_dim1 = *ne;
    e_offset = e_dim1 + 1;
    e -= e_offset;
    --fc;
    fa_dim1 = *nov;
    fa_offset = fa_dim1 + 1;
    fa -= fa_offset;
    d_dim1 = *md1;
    d_offset = d_dim1 + 1;
    d -= d_offset;
    c_dim1 = *mc1;
    c_offset = c_dim1 + 1;
    c -= c_offset;
    cc_dim1 = *nrc;
    cc_dim2 = *nov;
    cc_offset = cc_dim1 * (cc_dim2 + 1) + 1;
    cc -= cc_offset;
    b_dim1 = *nov;
    b_dim2 = *ncb;
    b_offset = b_dim1 * (b_dim2 + 1) + 1;
    b -= b_offset;
    a2_dim1 = *nov;
    a2_dim2 = *nov;
    a2_offset = a2_dim1 * (a2_dim2 + 1) + 1;
    a2 -= a2_offset;
    s_dim1 = *nov;
    s_dim2 = *nov;
    s_offset = s_dim1 * (s_dim2 + 1) + 1;
    s -= s_offset;

    /* Function Body */
    i__1 = *nov;
    for (i = 1; i <= i__1; ++i) {

	i__2 = *nov;
	for (j = 1; j <= i__2; ++j) {
	    e[i + j * e_dim1] = s[i + (j + *na * s_dim2) * s_dim1];
	    e[i + (*nov + j) * e_dim1] = a2[i + (j + *na * a2_dim2) * a2_dim1]
		    ;
/* L1: */
	}

	i__2 = *ncb;
	for (j = 1; j <= i__2; ++j) {
	    e[i + ((*nov << 1) + j) * e_dim1] = b[i + (j + *na * b_dim2) * 
		    b_dim1];
/* L2: */
	}

	rhse[i] = fa[i + *na * fa_dim1];

/* L3: */
    }

    i__1 = *nrc;
    for (i = 1; i <= i__1; ++i) {

	i__2 = *nov;
	for (j = 1; j <= i__2; ++j) {
	    e[*nov + i + j * e_dim1] = cc[i + (j + cc_dim2) * cc_dim1];
	    e[*nov + i + (*nov + j) * e_dim1] = cc[i + (j + (*na + 1) * 
		    cc_dim2) * cc_dim1];
/* L4: */
	}

	i__2 = *ncb;
	for (j = 1; j <= i__2; ++j) {
	    e[*nov + i + ((*nov << 1) + j) * e_dim1] = d[i + j * d_dim1];
/* L5: */
	}

	rhse[*nov + i] = fc[i];

/* L6: */
    }

    if (*idb != 0) {
	s_wsfe(&io___282);
	e_wsfe();
	i__1 = *ne;
	for (i = 1; i <= i__1; ++i) {
	    s_wsfe(&io___283);
	    i__2 = *ne;
	    for (j = 1; j <= i__2; ++j) {
		do_fio(&c__1, (char *)&e[i + j * e_dim1], (ftnlen)sizeof(
			doublereal));
	    }
	    do_fio(&c__1, (char *)&rhse[i], (ftnlen)sizeof(doublereal));
	    e_wsfe();
/* L7: */
	}
    }

    if (*nllv == 0) {
	ge_(ne, ne, &e[e_offset], &c__1, ne, &xe[1], ne, &rhse[1], &ir[1], &
		ic[1]);
    } else if (*nllv > 0) {
	i__1 = abs(*nllv);
	nlvc_(ne, ne, &i__1, &e[e_offset], &xe[1], &ir[1], &ic[1]);
    } else {
	i__1 = *ne - 1;
	for (i = 1; i <= i__1; ++i) {
	    rhse[i] = blrcn_1.zero;
/* L8: */
	}
	rhse[*ne] = blrcn_1.one;
	ge_(ne, ne, &e[e_offset], &c__1, ne, &xe[1], ne, &rhse[1], &ir[1], &
		ic[1]);
    }

    if (*idb != 0) {
	s_wsfe(&io___284);
	i__1 = *ne;
	for (i = 1; i <= i__1; ++i) {
	    do_fio(&c__1, (char *)&xe[i], (ftnlen)sizeof(doublereal));
	}
	e_wsfe();
    }

    i__1 = *nrc;
    for (i = 1; i <= i__1; ++i) {
	fc[i] = xe[*nov + i];
/* L9: */
    }


    return 0;
} /* dimrge_ */


/*     ---------- ---- */
/* Subroutine */ int bcksub_(nov, nov2, na, ncb, s, a1, a2, t, b, fa, xe, fc, 
	ipr, ipc)
integer *nov, *nov2, *na, *ncb;
doublereal *s, *a1, *a2, *t, *b, *fa, *xe, *fc;
integer *ipr, *ipc;
{
    /* System generated locals */
    integer s_dim1, s_dim2, s_offset, a1_dim1, a1_dim2, a1_offset, a2_dim1, 
	    a2_dim2, a2_offset, t_dim1, t_dim2, t_offset, b_dim1, b_dim2, 
	    b_offset, fa_dim1, fa_offset, ipr_dim1, ipr_offset, ipc_dim1, 
	    ipc_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer i, j, k, i1, i2, j1, k1;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */


/* Backsubstitution in the reduction process. */


    /* Parameter adjustments */
    ipc_dim1 = *nov;
    ipc_offset = ipc_dim1 + 1;
    ipc -= ipc_offset;
    ipr_dim1 = *nov2;
    ipr_offset = ipr_dim1 + 1;
    ipr -= ipr_offset;
    --fc;
    --xe;
    fa_dim1 = *nov;
    fa_offset = fa_dim1 + 1;
    fa -= fa_offset;
    b_dim1 = *nov;
    b_dim2 = *ncb;
    b_offset = b_dim1 * (b_dim2 + 1) + 1;
    b -= b_offset;
    t_dim1 = *nov;
    t_dim2 = *nov;
    t_offset = t_dim1 * (t_dim2 + 1) + 1;
    t -= t_offset;
    a2_dim1 = *nov;
    a2_dim2 = *nov;
    a2_offset = a2_dim1 * (a2_dim2 + 1) + 1;
    a2 -= a2_offset;
    a1_dim1 = *nov;
    a1_dim2 = *nov;
    a1_offset = a1_dim1 * (a1_dim2 + 1) + 1;
    a1 -= a1_offset;
    s_dim1 = *nov;
    s_dim2 = *nov;
    s_offset = s_dim1 * (s_dim2 + 1) + 1;
    s -= s_offset;

    /* Function Body */
    i__1 = *nov;
    for (i = 1; i <= i__1; ++i) {
	fa[i + (*na + 1) * fa_dim1] = xe[*nov + i];
/* L1: */
    }

    i__1 = blbnd_1.nam1;
    for (i = 1; i <= i__1; ++i) {
	i1 = *na - i;
	i2 = i1 + 1;
	i__2 = *nov;
	for (j1 = 1; j1 <= i__2; ++j1) {
	    j = *nov + 1 - j1;
	    fa[j + i2 * fa_dim1] = fa[j + i1 * fa_dim1];
	    i__3 = *nov;
	    for (k = 1; k <= i__3; ++k) {
		fa[j + i2 * fa_dim1] -= s[j + (k + i1 * s_dim2) * s_dim1] * 
			xe[k];
		fa[j + i2 * fa_dim1] -= t[j + (k + i1 * t_dim2) * t_dim1] * 
			fa[k + (i2 + 1) * fa_dim1];
/* L3: */
	    }
	    i__3 = *ncb;
	    for (k = 1; k <= i__3; ++k) {
		fa[j + i2 * fa_dim1] -= b[j + (k + i1 * b_dim2) * b_dim1] * 
			fc[*nov + k];
/* L4: */
	    }
	    if (j == *nov) {
		goto L6;
	    }
	    k1 = j + 1;
	    i__3 = *nov;
	    for (k = k1; k <= i__3; ++k) {
		fa[j + i2 * fa_dim1] -= a2[j + (k + i1 * a2_dim2) * a2_dim1] *
			 fa[k + i2 * fa_dim1];
/* L5: */
	    }
L6:
	    fa[j + i2 * fa_dim1] /= a2[j + (j + i1 * a2_dim2) * a2_dim1];
/* L7: */
	}
/* L8: */
    }

    i__1 = *nov;
    for (k = 1; k <= i__1; ++k) {
	fa[k + fa_dim1] = xe[k];
/* L9: */
    }

    return 0;
} /* bcksub_ */


/*     ---------- ------ */
/* Subroutine */ int print1_(nov, na, nra, nca, ncb, nrc, ma1, ma2, a, mb1, 
	mb2, b, mc1, c, md1, d, mfa1, fa, fc)
integer *nov, *na, *nra, *nca, *ncb, *nrc, *ma1, *ma2;
doublereal *a;
integer *mb1, *mb2;
doublereal *b;
integer *mc1;
doublereal *c;
integer *md1;
doublereal *d;
integer *mfa1;
doublereal *fa, *fc;
{
    /* Format strings */
    static char fmt_101[] = "(\002 A , B , FA (FULL DIMENSION) :\002)";
    static char fmt_102[] = "(\002 I=\002,i2)";
    static char fmt_103[] = "(1x,12e10.3)";
    static char fmt_104[] = "(\002 C (FULL DIMENSION) :\002)";
    static char fmt_106[] = "(\002 LAST NOV COLUMNS OF C :\002)";
    static char fmt_105[] = "(\002 D , FC\002)";

    /* System generated locals */
    integer a_dim1, a_dim2, a_offset, b_dim1, b_dim2, b_offset, c_dim1, 
	    c_offset, d_dim1, d_offset, fa_dim1, fa_offset, i__1, i__2, i__3, 
	    i__4;

    /* Builtin functions */
    integer s_wsfe(), e_wsfe(), do_fio();

    /* Local variables */
    static integer i, ic, ir, ic1, ic2;

    /* Fortran I/O blocks */
    static cilist io___292 = { 0, 9, 0, fmt_101, 0 };
    static cilist io___294 = { 0, 9, 0, fmt_102, 0 };
    static cilist io___296 = { 0, 9, 0, fmt_103, 0 };
    static cilist io___298 = { 0, 9, 0, fmt_104, 0 };
    static cilist io___299 = { 0, 9, 0, fmt_102, 0 };
    static cilist io___302 = { 0, 9, 0, fmt_103, 0 };
    static cilist io___303 = { 0, 9, 0, fmt_106, 0 };
    static cilist io___304 = { 0, 9, 0, fmt_103, 0 };
    static cilist io___305 = { 0, 9, 0, fmt_105, 0 };
    static cilist io___306 = { 0, 9, 0, fmt_103, 0 };



/* SGLE IMPLICIT REAL             (A-H,O-Z) */



    /* Parameter adjustments */
    --fc;
    fa_dim1 = *mfa1;
    fa_offset = fa_dim1 + 1;
    fa -= fa_offset;
    d_dim1 = *md1;
    d_offset = d_dim1 + 1;
    d -= d_offset;
    c_dim1 = *mc1;
    c_offset = c_dim1 + 1;
    c -= c_offset;
    b_dim1 = *mb1;
    b_dim2 = *mb2;
    b_offset = b_dim1 * (b_dim2 + 1) + 1;
    b -= b_offset;
    a_dim1 = *ma1;
    a_dim2 = *ma2;
    a_offset = a_dim1 * (a_dim2 + 1) + 1;
    a -= a_offset;

    /* Function Body */
    s_wsfe(&io___292);
    e_wsfe();
    i__1 = *na;
    for (i = 1; i <= i__1; ++i) {
	s_wsfe(&io___294);
	do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	e_wsfe();
	i__2 = *nra;
	for (ir = 1; ir <= i__2; ++ir) {
	    s_wsfe(&io___296);
	    i__3 = *nca;
	    for (ic = 1; ic <= i__3; ++ic) {
		do_fio(&c__1, (char *)&a[i + (ir + ic * a_dim2) * a_dim1], (
			ftnlen)sizeof(doublereal));
	    }
	    i__4 = *ncb;
	    for (ic = 1; ic <= i__4; ++ic) {
		do_fio(&c__1, (char *)&b[i + (ir + ic * b_dim2) * b_dim1], (
			ftnlen)sizeof(doublereal));
	    }
	    do_fio(&c__1, (char *)&fa[i + ir * fa_dim1], (ftnlen)sizeof(
		    doublereal));
	    e_wsfe();
/* L1: */
	}
/* L2: */
    }

    s_wsfe(&io___298);
    e_wsfe();

    i__1 = *na;
    for (i = 1; i <= i__1; ++i) {
	s_wsfe(&io___299);
	do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	e_wsfe();
	i__2 = *nrc;
	for (ir = 1; ir <= i__2; ++ir) {
	    ic1 = (i - 1) * (*nca - *nov) + 1;
	    ic2 = ic1 + *nca - *nov - 1;
	    s_wsfe(&io___302);
	    i__3 = ic2;
	    for (ic = ic1; ic <= i__3; ++ic) {
		do_fio(&c__1, (char *)&c[ic + ir * c_dim1], (ftnlen)sizeof(
			doublereal));
	    }
	    e_wsfe();
/* L3: */
	}
/* L4: */
    }

    s_wsfe(&io___303);
    e_wsfe();
    ic1 = *na * (*nca - *nov) + 1;
    ic2 = ic1 + *nov - 1;
    i__1 = *nrc;
    for (ir = 1; ir <= i__1; ++ir) {
	s_wsfe(&io___304);
	i__2 = ic2;
	for (ic = ic1; ic <= i__2; ++ic) {
	    do_fio(&c__1, (char *)&c[ic + ir * c_dim1], (ftnlen)sizeof(
		    doublereal));
	}
	e_wsfe();
/* L5: */
    }

    s_wsfe(&io___305);
    e_wsfe();

    i__1 = *nrc;
    for (ir = 1; ir <= i__1; ++ir) {
	s_wsfe(&io___306);
	i__2 = *ncb;
	for (ic = 1; ic <= i__2; ++ic) {
	    do_fio(&c__1, (char *)&d[ir + ic * d_dim1], (ftnlen)sizeof(
		    doublereal));
	}
	do_fio(&c__1, (char *)&fc[ir], (ftnlen)sizeof(doublereal));
	e_wsfe();
/* L6: */
    }


    return 0;
} /* print1_ */


/*     ---------- ------ */
/* Subroutine */ int print2_(idb, na, nov, ncb, nrc, s, a1, a2, t, b, c, f)
integer *idb, *na, *nov, *ncb, *nrc;
doublereal *s, *a1, *a2, *t, *b, *c, *f;
{
    /* Format strings */
    static char fmt_101[] = "(\002 A1 , A2 , B :\002)";
    static char fmt_104[] = "(\002 I=\002,i2)";
    static char fmt_102[] = "(1x,9e10.3)";
    static char fmt_105[] = "(\002 S AND T : \002)";
    static char fmt_103[] = "(\002 C :\002)";
    static char fmt_106[] = "(\002 RESIDUALS IN PRINT2\002)";
    static char fmt_107[] = "(1x,e10.3)";

    /* System generated locals */
    integer a1_dim1, a1_dim2, a1_offset, a2_dim1, a2_dim2, a2_offset, b_dim1, 
	    b_dim2, b_offset, s_dim1, s_dim2, s_offset, t_dim1, t_dim2, 
	    t_offset, f_dim1, f_offset, c_dim1, c_dim2, c_offset, i__1, i__2, 
	    i__3, i__4, i__5;

    /* Builtin functions */
    integer s_wsfe(), e_wsfe(), do_fio();

    /* Local variables */
    static integer i;
    static doublereal r;
    static integer i1, ic, ir;

    /* Fortran I/O blocks */
    static cilist io___307 = { 0, 9, 0, fmt_101, 0 };
    static cilist io___309 = { 0, 9, 0, fmt_104, 0 };
    static cilist io___311 = { 0, 9, 0, fmt_102, 0 };
    static cilist io___313 = { 0, 9, 0, fmt_105, 0 };
    static cilist io___314 = { 0, 9, 0, fmt_104, 0 };
    static cilist io___315 = { 0, 9, 0, fmt_102, 0 };
    static cilist io___316 = { 0, 9, 0, fmt_103, 0 };
    static cilist io___317 = { 0, 9, 0, fmt_104, 0 };
    static cilist io___318 = { 0, 9, 0, fmt_102, 0 };
    static cilist io___319 = { 0, 9, 0, fmt_106, 0 };
    static cilist io___322 = { 0, 9, 0, fmt_107, 0 };



/* SGLE IMPLICIT REAL             (A-H,O-Z) */



    /* Parameter adjustments */
    f_dim1 = *nov;
    f_offset = f_dim1 + 1;
    f -= f_offset;
    c_dim1 = *nrc;
    c_dim2 = *nov;
    c_offset = c_dim1 * (c_dim2 + 1) + 1;
    c -= c_offset;
    b_dim1 = *nov;
    b_dim2 = *ncb;
    b_offset = b_dim1 * (b_dim2 + 1) + 1;
    b -= b_offset;
    t_dim1 = *nov;
    t_dim2 = *nov;
    t_offset = t_dim1 * (t_dim2 + 1) + 1;
    t -= t_offset;
    a2_dim1 = *nov;
    a2_dim2 = *nov;
    a2_offset = a2_dim1 * (a2_dim2 + 1) + 1;
    a2 -= a2_offset;
    a1_dim1 = *nov;
    a1_dim2 = *nov;
    a1_offset = a1_dim1 * (a1_dim2 + 1) + 1;
    a1 -= a1_offset;
    s_dim1 = *nov;
    s_dim2 = *nov;
    s_offset = s_dim1 * (s_dim2 + 1) + 1;
    s -= s_offset;

    /* Function Body */
    s_wsfe(&io___307);
    e_wsfe();

    i__1 = *na;
    for (i = 1; i <= i__1; ++i) {
	s_wsfe(&io___309);
	do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	e_wsfe();
	i__2 = *nov;
	for (ir = 1; ir <= i__2; ++ir) {
	    s_wsfe(&io___311);
	    i__3 = *nov;
	    for (ic = 1; ic <= i__3; ++ic) {
		do_fio(&c__1, (char *)&a1[ir + (ic + i * a1_dim2) * a1_dim1], 
			(ftnlen)sizeof(doublereal));
	    }
	    i__4 = *nov;
	    for (ic = 1; ic <= i__4; ++ic) {
		do_fio(&c__1, (char *)&a2[ir + (ic + i * a2_dim2) * a2_dim1], 
			(ftnlen)sizeof(doublereal));
	    }
	    i__5 = *ncb;
	    for (ic = 1; ic <= i__5; ++ic) {
		do_fio(&c__1, (char *)&b[ir + (ic + i * b_dim2) * b_dim1], (
			ftnlen)sizeof(doublereal));
	    }
	    e_wsfe();
/* L1: */
	}
/* L2: */
    }

    s_wsfe(&io___313);
    e_wsfe();
    i__1 = *na;
    for (i = 1; i <= i__1; ++i) {
	s_wsfe(&io___314);
	do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	e_wsfe();
	i__2 = *nov;
	for (ir = 1; ir <= i__2; ++ir) {
	    s_wsfe(&io___315);
	    i__3 = *nov;
	    for (ic = 1; ic <= i__3; ++ic) {
		do_fio(&c__1, (char *)&s[ir + (ic + i * s_dim2) * s_dim1], (
			ftnlen)sizeof(doublereal));
	    }
	    i__4 = *nov;
	    for (ic = 1; ic <= i__4; ++ic) {
		do_fio(&c__1, (char *)&t[ir + (ic + i * t_dim2) * t_dim1], (
			ftnlen)sizeof(doublereal));
	    }
	    e_wsfe();
/* L3: */
	}
/* L4: */
    }

    s_wsfe(&io___316);
    e_wsfe();

    i__1 = *na;
    for (i = 1; i <= i__1; ++i) {
	s_wsfe(&io___317);
	do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	e_wsfe();
	i__2 = *nrc;
	for (ir = 1; ir <= i__2; ++ir) {
	    s_wsfe(&io___318);
	    i__3 = *nov;
	    for (ic = 1; ic <= i__3; ++ic) {
		do_fio(&c__1, (char *)&c[ir + (ic + i * c_dim2) * c_dim1], (
			ftnlen)sizeof(doublereal));
	    }
	    e_wsfe();
/* L5: */
	}
/* L6: */
    }

/* The following can be used if a linear system is given */
/* for which all all solution components are equal to 1. */
/* Compute residuals for test problem : */

    if (*idb == 3) {
	s_wsfe(&io___319);
	e_wsfe();
	i__1 = *na;
	for (i1 = 1; i1 <= i__1; ++i1) {
	    i__2 = *nov;
	    for (ir = 1; ir <= i__2; ++ir) {
		r = f[ir + i1 * f_dim1];
		i__3 = *nov;
		for (ic = 1; ic <= i__3; ++ic) {
		    r -= s[ir + (ic + i1 * s_dim2) * s_dim1];
		    r -= t[ir + (ic + i1 * t_dim2) * t_dim1];
		    if (ic >= ir) {
			r -= a2[ir + (ic + i1 * a2_dim2) * a2_dim1];
		    }
/* L7: */
		}
		i__3 = *ncb;
		for (ic = 1; ic <= i__3; ++ic) {
		    r -= b[ir + (ic + i1 * b_dim2) * b_dim1];
/* L8: */
		}
		s_wsfe(&io___322);
		do_fio(&c__1, (char *)&r, (ftnlen)sizeof(doublereal));
		e_wsfe();
/* L9: */
	    }
/* L10: */
	}
    }


    return 0;
} /* print2_ */


/*     ---------  ------ */
/* Subroutine */ int print3_(na, nov, fa)
integer *na, *nov;
doublereal *fa;
{
    /* Format strings */
    static char fmt_101[] = "(\002 FA : \002)";
    static char fmt_102[] = "(\002 I=\002,i3)";
    static char fmt_103[] = "(1x,9e10.3)";

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

    /* Builtin functions */
    integer s_wsfe(), e_wsfe(), do_fio();

    /* Local variables */
    static integer i, ir;

    /* Fortran I/O blocks */
    static cilist io___323 = { 0, 9, 0, fmt_101, 0 };
    static cilist io___325 = { 0, 9, 0, fmt_102, 0 };
    static cilist io___326 = { 0, 9, 0, fmt_103, 0 };



/* SGLE IMPLICIT REAL (A-H,O-Z) */


    /* Parameter adjustments */
    fa_dim1 = *nov;
    fa_offset = fa_dim1 + 1;
    fa -= fa_offset;

    /* Function Body */
    s_wsfe(&io___323);
    e_wsfe();
    i__1 = *na;
    for (i = 1; i <= i__1; ++i) {
	s_wsfe(&io___325);
	do_fio(&c__1, (char *)&i, (ftnlen)sizeof(integer));
	e_wsfe();
	s_wsfe(&io___326);
	i__2 = *nov;
	for (ir = 1; ir <= i__2; ++ir) {
	    do_fio(&c__1, (char *)&fa[ir + i * fa_dim1], (ftnlen)sizeof(
		    doublereal));
	}
	e_wsfe();
/* L1: */
    }


    return 0;
} /* print3_ */


/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */
/*    Vectorized Subroutines for the Linear Equation Solver BRBD */
/* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------- */

/*     ---------- ------ */
/* Subroutine */ int conpar_(nov, na, nra, nca, ma1, ma2, a, ncb, mb1, mb2, b,
	 nrc, mc1, c, md1, d, rm, lc)
integer *nov, *na, *nra, *nca, *ma1, *ma2;
doublereal *a;
integer *ncb, *mb1, *mb2;
doublereal *b;
integer *nrc, *mc1;
doublereal *c;
integer *md1;
doublereal *d, *rm;
integer *lc;
{
    /* System generated locals */
    integer a_dim1, a_dim2, a_offset, b_dim1, b_dim2, b_offset, c_dim1, 
	    c_offset, d_dim1, d_offset, i__1, i__2, i__3, i__4;

    /* Local variables */
    static doublereal time0, time1;
    static integer i, l, m1, m2;
    extern /* Subroutine */ int autim0_(), autim1_();
    static integer ic, ir, ir1, nex, irp, icp1;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */



/* Set initial time (for timing of this subroutine). */

    /* Parameter adjustments */
    --lc;
    --rm;
    d_dim1 = *md1;
    d_offset = d_dim1 + 1;
    d -= d_offset;
    c_dim1 = *mc1;
    c_offset = c_dim1 + 1;
    c -= c_offset;
    b_dim1 = *mb1;
    b_dim2 = *mb2;
    b_offset = b_dim1 * (b_dim2 + 1) + 1;
    b -= b_offset;
    a_dim1 = *ma1;
    a_dim2 = *ma2;
    a_offset = a_dim1 * (a_dim2 + 1) + 1;
    a -= a_offset;

    /* Function Body */
    autim0_(&time0);

    nex = *nca - (*nov << 1);
    if (nex == 0) {
	return 0;
    }

/*  Condensation of parameters  ( Elimination of "local" variables ). */

    m1 = *nov + 1;
    m2 = *nov + nex;

    i__1 = m2;
    for (ic = m1; ic <= i__1; ++ic) {
	ir1 = ic - *nov + 1;
	irp = ir1 - 1;
	icp1 = ic + 1;

	i__2 = *nra;
	for (ir = ir1; ir <= i__2; ++ir) {
	    i__3 = *na;
	    for (i = 1; i <= i__3; ++i) {
		rm[i] = a[i + (ir + ic * a_dim2) * a_dim1] / a[i + (irp + ic *
			 a_dim2) * a_dim1];
		a[i + (ir + ic * a_dim2) * a_dim1] = rm[i];
/* L10: */
	    }
	    i__3 = *nov;
	    for (l = 1; l <= i__3; ++l) {
		i__4 = *na;
		for (i = 1; i <= i__4; ++i) {
		    a[i + (ir + l * a_dim2) * a_dim1] -= rm[i] * a[i + (irp + 
			    l * a_dim2) * a_dim1];
/* L11: */
		}
/* L1: */
	    }
	    i__3 = *nca;
	    for (l = icp1; l <= i__3; ++l) {
		i__4 = *na;
		for (i = 1; i <= i__4; ++i) {
		    a[i + (ir + l * a_dim2) * a_dim1] -= rm[i] * a[i + (irp + 
			    l * a_dim2) * a_dim1];
/* L12: */
		}
/* L2: */
	    }
	    i__3 = *ncb;
	    for (l = 1; l <= i__3; ++l) {
		i__4 = *na;
		for (i = 1; i <= i__4; ++i) {
		    b[i + (ir + l * b_dim2) * b_dim1] -= rm[i] * b[i + (irp + 
			    l * b_dim2) * b_dim1];
/* L13: */
		}
/* L3: */
	    }
/* L4: */
	}

	i__2 = *nrc;
	for (ir = blcde_1.nbc + 1; ir <= i__2; ++ir) {
	    i__3 = *na;
	    for (i = 1; i <= i__3; ++i) {
		lc[i] = (i - 1) * (*nca - *nov) + ic;
/* L40: */
	    }
	    i__3 = *na;
	    for (i = 1; i <= i__3; ++i) {
		rm[i] = c[lc[i] + ir * c_dim1] / a[i + (irp + ic * a_dim2) * 
			a_dim1];
		c[lc[i] + ir * c_dim1] = rm[i];
/* L50: */
	    }
	    i__3 = *nov;
	    for (l = 1; l <= i__3; ++l) {
		i__4 = *na;
		for (i = 1; i <= i__4; ++i) {
		    lc[i] = (i - 1) * (*nca - *nov) + l;
/* L51: */
		}
		i__4 = *na;
		for (i = 1; i <= i__4; ++i) {
		    c[lc[i] + ir * c_dim1] -= rm[i] * a[i + (irp + l * a_dim2)
			     * a_dim1];
/* L52: */
		}
/* L5: */
	    }
	    i__3 = *nca;
	    for (l = icp1; l <= i__3; ++l) {
		i__4 = *na;
		for (i = 1; i <= i__4; ++i) {
		    lc[i] = (i - 1) * (*nca - *nov) + l;
/* L60: */
		}
		i__4 = *na;
		for (i = 1; i <= i__4; ++i) {
		    c[lc[i] + ir * c_dim1] -= rm[i] * a[i + (irp + l * a_dim2)
			     * a_dim1];
/* L61: */
		}
/* L6: */
	    }
	    i__3 = *ncb;
	    for (l = 1; l <= i__3; ++l) {
		i__4 = *na;
		for (i = 1; i <= i__4; ++i) {
		    d[ir + l * d_dim1] -= rm[i] * b[i + (irp + l * b_dim2) * 
			    b_dim1];
/* L71: */
		}
/* L7: */
	    }
/* L8: */
	}

/* L9: */
    }

/* Determine the time spent in this subroutine. */

    autim1_(&time1);
    bltim_1.tconpa = bltim_1.tconpa + time1 - time0;

    return 0;
} /* conpar_ */


/*     ---------- ------ */
/* Subroutine */ int conrhs_(nov, na, nra, nca, ma1, ma2, a, nrc, mc1, c, 
	mfa1, fa, fc, lc)
integer *nov, *na, *nra, *nca, *ma1, *ma2;
doublereal *a;
integer *nrc, *mc1;
doublereal *c;
integer *mfa1;
doublereal *fa, *fc;
integer *lc;
{
    /* System generated locals */
    integer a_dim1, a_dim2, a_offset, c_dim1, c_offset, fa_dim1, fa_offset, 
	    i__1, i__2, i__3;

    /* Local variables */
    static doublereal time0, time1;
    static integer i, m1, m2;
    extern /* Subroutine */ int autim0_(), autim1_();
    static integer ic, ir, ir1, nex, irp;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */



/* Set initial time (for timing of this subroutine). */

    /* Parameter adjustments */
    --lc;
    --fc;
    fa_dim1 = *mfa1;
    fa_offset = fa_dim1 + 1;
    fa -= fa_offset;
    c_dim1 = *mc1;
    c_offset = c_dim1 + 1;
    c -= c_offset;
    a_dim1 = *ma1;
    a_dim2 = *ma2;
    a_offset = a_dim1 * (a_dim2 + 1) + 1;
    a -= a_offset;

    /* Function Body */
    autim0_(&time0);

/* Condensation of the right hand side. */

    nex = *nca - (*nov << 1);
    if (nex == 0) {
	return 0;
    }

    m1 = *nov + 1;
    m2 = *nov + nex;

    i__1 = m2;
    for (ic = m1; ic <= i__1; ++ic) {
	ir1 = ic - *nov + 1;
	irp = ir1 - 1;

	i__2 = *nra;
	for (ir = ir1; ir <= i__2; ++ir) {
/* Note that  RM=A(I,IR,IC) is the multiplier. */
	    i__3 = *na;
	    for (i = 1; i <= i__3; ++i) {
		fa[i + ir * fa_dim1] -= a[i + (ir + ic * a_dim2) * a_dim1] * 
			fa[i + irp * fa_dim1];
/* L11: */
	    }
/* L1: */
	}

	i__2 = *nrc;
	for (ir = blcde_1.nbc + 1; ir <= i__2; ++ir) {
	    i__3 = *na;
	    for (i = 1; i <= i__3; ++i) {
		lc[i] = (i - 1) * (*nca - *nov) + ic;
/* L20: */
	    }
/* Note that  RM=C(LC(I),IR) is the multiplier. */
	    i__3 = *na;
	    for (i = 1; i <= i__3; ++i) {
		fc[ir] -= c[lc[i] + ir * c_dim1] * fa[i + irp * fa_dim1];
/* L21: */
	    }
/* L2: */
	}

/* L3: */
    }

/* Determine the time spent in this subroutine, */

    autim1_(&time1);
    bltim_1.tconrh = bltim_1.tconrh + time1 - time0;

    return 0;
} /* conrhs_ */


/*     ---------- ------ */
/* Subroutine */ int infpar_(na, nov, nra, nca, ma1, ma2, a, ncb, mb1, mb2, b,
	 mfa1, fa, nc, fc, nfadr, fadr)
integer *na, *nov, *nra, *nca, *ma1, *ma2;
doublereal *a;
integer *ncb, *mb1, *mb2;
doublereal *b;
integer *mfa1;
doublereal *fa;
integer *nc;
doublereal *fc;
integer *nfadr;
doublereal *fadr;
{
    /* System generated locals */
    integer a_dim1, a_dim2, a_offset, b_dim1, b_dim2, b_offset, fa_dim1, 
	    fa_offset, fadr_dim1, fadr_offset, i__1, i__2, i__3;

    /* Local variables */
    static integer nram;
    static doublereal time0, time1;
    static integer i, j, k, k1, k2;
    extern /* Subroutine */ int autim0_(), autim1_();
    static integer ir, ir1;


/* SGLE IMPLICIT REAL             (A-H,O-Z) */



/* Set initial time (for timing of this subroutine). */

    /* Parameter adjustments */
    fadr_dim1 = *nov;
    fadr_offset = fadr_dim1 + 1;
    fadr -= fadr_offset;
    --fc;
    fa_dim1 = *mfa1;
    fa_offset = fa_dim1 + 1;
    fa -= fa_offset;
    b_dim1 = *mb1;
    b_dim2 = *mb2;
    b_offset = b_dim1 * (b_dim2 + 1) + 1;
    b -= b_offset;
    a_dim1 = *ma1;
    a_dim2 = *ma2;
    a_offset = a_dim1 * (a_dim2 + 1) + 1;
    a -= a_offset;

    /* Function Body */
    autim0_(&time0);

/* Determine the "local" variables by backsubstitution. */

    nram = *nra - *nov;
    i__1 = nram;
    for (ir1 = 1; ir1 <= i__1; ++ir1) {
	ir = nram + 1 - ir1;
	i__2 = *nov;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *na;
	    for (i = 1; i <= i__3; ++i) {
		fa[i + ir * fa_dim1] -= a[i + (ir + k * a_dim2) * a_dim1] * 
			fadr[k + i * fadr_dim1];
/* L11: */
	    }
/* L1: */
	}
	i__2 = *nov;
	for (k = 1; k <= i__2; ++k) {
	    fadr[k + (*na + 2) * fadr_dim1] = fc[k];
/* L2: */
	}
	i__2 = *nov;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *na;
	    for (i = 1; i <= i__3; ++i) {
		fa[i + ir * fa_dim1] -= a[i + (ir + (*nca - *nov + k) * 
			a_dim2) * a_dim1] * fadr[k + (i + 1) * fadr_dim1];
/* L31: */
	    }
/* L3: */
	}
	i__2 = *ncb;
	for (k = 1; k <= i__2; ++k) {
	    i__3 = *na;
	    for (i = 1; i <= i__3; ++i) {
		fa[i + ir * fa_dim1] -= b[i + (ir + k * b_dim2) * b_dim1] * 
			fc[*nov + k];
/* L61: */
	    }
/* L6: */
	}
	if (ir1 == 1) {
	    goto L8;
	}
	k1 = *nca - *nov - ir1 + 2;
	k2 = *nca - *nov;
	i__2 = k2;
	for (k = k1; k <= i__2; ++k) {
	    i__3 = *na;
	    for (i = 1; i <= i__3; ++i) {
		fa[i + ir * fa_dim1] -= a[i + (ir + k * a_dim2) * a_dim1] * 
			fa[i + k * fa_dim1];
/* L71: */
	    }
/* L7: */
	}
L8:
	i__2 = *na;
	for (i = 1; i <= i__2; ++i) {
	    fa[i + (ir + *nov) * fa_dim1] = fa[i + ir * fa_dim1] / a[i + (ir 
		    + (*nca - *nov - ir1 + 1) * a_dim2) * a_dim1];
/* L91: */
	}
/* L9: */
    }

/* Copy the solution generated by DRBSUB (stored in FADR) */
/* into the final solution vector (stored in FA). */

    i__1 = *nov;
    for (j = 1; j <= i__1; ++j) {
	i__2 = *na;
	for (i = 1; i <= i__2; ++i) {
	    fa[i + j * fa_dim1] = fadr[j + i * fadr_dim1];
/* L21: */
	}
/* L22: */
    }

/* Determine the time spent in this subroutine. */

    autim1_(&time1);
    bltim_1.tinfpa = bltim_1.tinfpa + time1 - time0;

    return 0;
} /* infpar_ */