#include /* autlib3.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" /* Common Block Declarations */ typedef struct { int irot; int nrot[1000]; double torper; } ROTCHK; extern ROTCHK blrtn; struct { integer ndm, ndmp1, nrow, nclm, nrc, ncc, npar, nfpar, nbc0, nint0; } blicn_; #define blicn_1 blicn_ 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 { 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_; #define blwif_1 blwif_ 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 nmx, nuzr; doublereal rl0, rl1, a0, a1; } bllim_; #define bllim_1 bllim_ struct { doublereal ds, dsmin, dsmax; integer iads; } bldls_; #define bldls_1 bldls_ struct { doublereal u1zz[NAUTO], u2zz[NAUTO], f1zz[NAUTO], f2zz[NAUTO]; } bldif_; #define bldif_1 bldif_ /* Table of constant values */ static integer c__1 = 1; static integer c__0 = 0; static integer c__2 = 2; static integer c__3 = 3; /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Subroutines for the Continuation of Limit Points (Algebraic Problems) */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fnlp_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int fflp_(); static integer i, j; static doublereal ep, umx; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates the equations for the 2-par continuation of limit points. */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ fflp_(ndim, &u[1], uold, &icp[1], &par[1], &f[1], &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); if (*ijac == 0) { return 0; } /* Generate the Jacobian. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L1: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { blwif_1.u1xx[j - 1] = u[j]; blwif_1.u2xx[j - 1] = u[j]; /* L2: */ } blwif_1.u1xx[i - 1] -= ep; blwif_1.u2xx[i - 1] += ep; fflp_(ndim, blwif_1.u1xx, uold, &icp[1], &par[1], blwif_1.f1xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); fflp_(ndim, blwif_1.u2xx, uold, &icp[1], &par[1], blwif_1.f2xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[j + i * dfdu_dim1] = (blwif_1.f2xx[j - 1] - blwif_1.f1xx[j - 1]) / (ep * 2); /* L3: */ } /* L4: */ } par[icp[1]] += ep; fflp_(ndim, &u[1], uold, &icp[1], &par[1], blwif_1.f1xx, &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__1 = *ndim; for (j = 1; j <= i__1; ++j) { dfdp[j + icp[1] * dfdp_dim1] = (blwif_1.f1xx[j - 1] - f[j]) / ep; /* L5: */ } par[icp[1]] -= ep; return 0; } /* fnlp_ */ /* ---------- ---- */ /* Subroutine */ int fflp_(ndim1, u, uold, icp1, par1, f, ndm, dfdu, dfdp) integer *ndim1; doublereal *u, *uold; integer *icp1; doublereal *par1, *f; integer *ndm; doublereal *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ static integer ijac; extern /* Subroutine */ int fnds_(), funi_(); static integer i, j; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dfdp_dim1 = *ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --u; /* Function Body */ ijac = 1; blbcn_1.par[blbcn_1.icp[1] - 1] = u[blbcn_1.ndim]; if (blbcn_1.ips == -1) { fnds_(ndm, &u[1], uold, blbcn_1.icp, blbcn_1.par, &ijac, &f[1], &dfdu[ dfdu_offset], &dfdp[dfdp_offset]); } else { funi_(ndm, &u[1], uold, blbcn_1.icp, blbcn_1.par, &ijac, &f[1], &dfdu[ dfdu_offset], &dfdp[dfdp_offset]); } i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndm + i] = blrcn_1.zero; i__2 = *ndm; for (j = 1; j <= i__2; ++j) { f[*ndm + i] += dfdu[i + j * dfdu_dim1] * u[*ndm + j]; /* L1: */ } /* L2: */ } f[blbcn_1.ndim] = -blrcn_1.one; /* SGLE F(NDIM)=-ONE */ i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[blbcn_1.ndim] += u[*ndm + i] * u[*ndm + i]; /* L3: */ } return 0; } /* fflp_ */ /* ---------- ------ */ /* Subroutine */ int stpnlp_(ibr, u, ndm2, smat, dfdu, dfuxx, dfdp, v, f, ir, ic) integer *ibr; doublereal *u; integer *ndm2; doublereal *smat, *dfdu, *dfuxx, *dfdp, *v, *f; integer *ir, *ic; { /* System generated locals */ integer smat_dim1, smat_offset, dfdu_dim1, dfdu_offset, dfuxx_dim1, dfuxx_offset, dfdp_dim1, dfdp_offset, i__1; /* Local variables */ static integer ijac; extern /* Subroutine */ int fnds_(), funi_(), nlvc_(); static doublereal uold; static integer i; static logical found; extern /* Subroutine */ int nrmlz_(), readl3_(), findl3_(); static integer nfpar1, itp; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates starting data for the continuation of limit points. */ /* Parameter adjustments */ --ic; --ir; --f; --v; dfdp_dim1 = blicn_1.ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfuxx_dim1 = blicn_1.ndmp1; dfuxx_offset = dfuxx_dim1 + 1; dfuxx -= dfuxx_offset; dfdu_dim1 = blicn_1.ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; smat_dim1 = *ndm2; smat_offset = smat_dim1 + 1; smat -= smat_offset; --u; /* Function Body */ findl3_(&blbcn_1.irs, &itp, &nfpar1, &found); readl3_(&blbcn_1.ips, ibr, &u[1], blbcn_1.par); ijac = 1; if (blbcn_1.ips == -1) { fnds_(&blicn_1.ndm, &u[1], &uold, blbcn_1.icp, blbcn_1.par, &ijac, &f[ 1], &dfdu[dfdu_offset], &dfdp[dfdp_offset]); } else { funi_(&blicn_1.ndm, &u[1], &uold, blbcn_1.icp, blbcn_1.par, &ijac, &f[ 1], &dfdu[dfdu_offset], &dfdp[dfdp_offset]); } nlvc_(&blicn_1.ndm, &blicn_1.ndm, &c__1, &dfdu[dfdu_offset], &v[1], &ir[1] , &ic[1]); nrmlz_(&blicn_1.ndm, &v[1]); i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { u[blicn_1.ndm + i] = v[i]; /* L1: */ } u[blbcn_1.ndim] = blbcn_1.par[blbcn_1.icp[1] - 1]; return 0; } /* stpnlp_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Subroutines for the Optimization of Algebraic Systems */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fnc1_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1; /* Local variables */ extern /* Subroutine */ int fopi_(), funi_(); static integer i, j; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generate the equations for the continuation scheme used for */ /* the optimization of algebraic systems (one parameter). */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ par[icp[2]] = u[*ndim]; funi_(&blicn_1.ndm, &u[1], uold, &icp[1], &par[1], ijac, &f[1], &dfdu[ dfdu_offset], &dfdp[dfdp_offset]); /* Rearrange (Since dimensions in FNC1 and FUNI differ). */ if (*ijac != 0) { for (j = blicn_1.ndm; j >= 1; --j) { for (i = blicn_1.ndm; i >= 1; --i) { dfdu[i + j * dfdu_dim1] = dfdu[(j - 1) * blicn_1.ndm + i + dfdu_dim1]; /* L1: */ } /* L2: */ } for (j = blicn_1.npar; j >= 1; --j) { for (i = blicn_1.ndm; i >= 1; --i) { dfdp[i + j * dfdp_dim1] = dfdp[(j - 1) * blicn_1.ndm + i + dfdp_dim1]; /* L3: */ } /* L4: */ } } fopi_(&blicn_1.ndm, &u[1], &icp[1], &par[1], ijac, &f[*ndim], blwif_1.dduxx, blwif_1.ddpxx); f[*ndim] = par[icp[1]] - f[*ndim]; if (*ijac != 0) { i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { dfdu[*ndim + i * dfdu_dim1] = -blwif_1.dduxx[i - 1]; dfdu[i + *ndim * dfdu_dim1] = dfdp[i + icp[2] * dfdp_dim1]; dfdp[i + icp[1] * dfdp_dim1] = 0.; /* L5: */ } dfdu[*ndim + *ndim * dfdu_dim1] = -blwif_1.ddpxx[icp[2] - 1]; dfdp[*ndim + icp[1] * dfdp_dim1] = 1.; } return 0; } /* fnc1_ */ /* ---------- ------ */ /* Subroutine */ int stpnc1_(ibr, u, ndm2, smat, dfdu, dfuxx, dfdp, v, f, ir, ic) integer *ibr; doublereal *u; integer *ndm2; doublereal *smat, *dfdu, *dfuxx, *dfdp, *v, *f; integer *ir, *ic; { /* System generated locals */ integer smat_dim1, smat_offset, dfdu_dim1, dfdu_offset, dfuxx_dim1, dfuxx_offset, dfdp_dim1, dfdp_offset; /* Local variables */ extern /* Subroutine */ int fopi_(), stpnt_(); static doublereal fop; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generate starting data for optimization problems (one parameter). */ /* Parameter adjustments */ --ic; --ir; --f; --v; dfdp_dim1 = blicn_1.ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfuxx_dim1 = blicn_1.ndmp1; dfuxx_offset = dfuxx_dim1 + 1; dfuxx -= dfuxx_offset; dfdu_dim1 = blicn_1.ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; smat_dim1 = *ndm2; smat_offset = smat_dim1 + 1; smat -= smat_offset; --u; /* Function Body */ stpnt_(&blbcn_1.ndim, &u[1], blbcn_1.par); fopi_(&blicn_1.ndm, &u[1], blbcn_1.icp, blbcn_1.par, &c__0, &fop, &dfdu[ dfdu_offset], &dfdp[dfdp_offset]); blbcn_1.par[blbcn_1.icp[0] - 1] = fop; u[blbcn_1.ndim] = blbcn_1.par[blbcn_1.icp[1] - 1]; return 0; } /* stpnc1_ */ /* ---------- ---- */ /* Subroutine */ int fnc2_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ static integer i, j; static doublereal ep, umx; extern /* Subroutine */ int ffc2_(); /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generate the equations for the continuation scheme used for the */ /* optimization of algebraic systems (more than one parameter). */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ ffc2_(ndim, &u[1], uold, &icp[1], &par[1], &f[1], blwif_1.dfuxx, blwif_1.dfpxx); if (*ijac == 0) { return 0; } /* Generate the Jacobian. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L1: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { blwif_1.u1xx[j - 1] = u[j]; blwif_1.u2xx[j - 1] = u[j]; /* L2: */ } blwif_1.u1xx[i - 1] -= ep; blwif_1.u2xx[i - 1] += ep; ffc2_(ndim, blwif_1.u1xx, uold, &icp[1], &par[1], blwif_1.f1xx, blwif_1.dfuxx, blwif_1.dfpxx); ffc2_(ndim, blwif_1.u2xx, uold, &icp[1], &par[1], blwif_1.f2xx, blwif_1.dfuxx, blwif_1.dfpxx); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[j + i * dfdu_dim1] = (blwif_1.f2xx[j - 1] - blwif_1.f1xx[j - 1]) / (ep * 2); /* L3: */ } /* L4: */ } i__1 = *ndim; for (i = 1; i <= i__1; ++i) { dfdp[i + icp[1] * dfdp_dim1] = 0.; /* L5: */ } dfdp[*ndim + icp[1] * dfdp_dim1] = 1.; return 0; } /* fnc2_ */ /* ---------- ---- */ /* Subroutine */ int ffc2_(ndim, u, uold, icp, par, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par, *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ static integer ijac, icpm; extern /* Subroutine */ int fopi_(), funi_(); static integer i, j; static doublereal fop; static integer ndm2; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dfdp_dim1 = blicn_1.ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = blicn_1.ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ ijac = 1; i__1 = blicn_1.nfpar; for (i = 2; i <= i__1; ++i) { par[icp[i]] = u[(blicn_1.ndm << 1) + i]; /* L1: */ } funi_(&blicn_1.ndm, &u[1], uold, &icp[1], &par[1], &ijac, &f[1], &dfdu[ dfdu_offset], &dfdp[dfdp_offset]); fopi_(&blicn_1.ndm, &u[1], &icp[1], &par[1], &ijac, &fop, blwif_1.dduxx, blwif_1.ddpxx); i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { f[blicn_1.ndm + i] = blwif_1.dduxx[i - 1] * u[(blicn_1.ndm << 1) + 1]; i__2 = blicn_1.ndm; for (j = 1; j <= i__2; ++j) { f[blicn_1.ndm + i] += dfdu[j + i * dfdu_dim1] * u[blicn_1.ndm + j] ; /* L2: */ } /* L3: */ } ndm2 = blicn_1.ndm << 1; icpm = blicn_1.nfpar - 2; i__1 = icpm; for (i = 1; i <= i__1; ++i) { f[ndm2 + i] = blwif_1.ddpxx[icp[i + 1] - 1] * u[ndm2 + 1]; /* L4: */ } i__1 = icpm; for (i = 1; i <= i__1; ++i) { i__2 = blicn_1.ndm; for (j = 1; j <= i__2; ++j) { f[ndm2 + i] += u[blicn_1.ndm + j] * dfdp[j + icp[i + 1] * dfdp_dim1]; /* L5: */ } /* L6: */ } f[*ndim - 1] = u[ndm2 + 1] * u[ndm2 + 1] - 1; i__1 = blicn_1.ndm; for (j = 1; j <= i__1; ++j) { f[*ndim - 1] += u[blicn_1.ndm + j] * u[blicn_1.ndm + j]; /* L7: */ } f[*ndim] = par[icp[1]] - fop; return 0; } /* ffc2_ */ /* ---------- ------ */ /* Subroutine */ int stpnc2_(ibr, u, ndm2, smat, dfdu, dfu, dfdp, v, f, ir, ic) integer *ibr; doublereal *u; integer *ndm2; doublereal *smat, *dfdu, *dfu, *dfdp, *v, *f; integer *ir, *ic; { /* System generated locals */ integer smat_dim1, smat_offset, dfdu_dim1, dfdu_offset, dfu_dim1, dfu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ static integer ijac; extern /* Subroutine */ int fopi_(), funi_(), nlvc_(); static doublereal uold; static integer i, j; static logical found; extern /* Subroutine */ int nrmlz_(), readl3_(), findl3_(); static doublereal fop; static integer itp; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates starting data for the continuation equations for */ /* optimization of algebraic systems (More than one parameter). */ /* Parameter adjustments */ --ic; --ir; --f; --v; dfdp_dim1 = blicn_1.ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfu_dim1 = blicn_1.ndmp1; dfu_offset = dfu_dim1 + 1; dfu -= dfu_offset; dfdu_dim1 = blicn_1.ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; smat_dim1 = *ndm2; smat_offset = smat_dim1 + 1; smat -= smat_offset; --u; /* Function Body */ findl3_(&blbcn_1.irs, &itp, &blicn_1.nfpar, &found); ++blicn_1.nfpar; readl3_(&blbcn_1.ips, ibr, &u[1], blbcn_1.par); if (blicn_1.nfpar == 3) { ijac = 1; funi_(&blicn_1.ndm, &u[1], &uold, blbcn_1.icp, blbcn_1.par, &ijac, &f[ 1], &dfdu[dfdu_offset], &dfdp[dfdp_offset]); fopi_(&blicn_1.ndm, &u[1], blbcn_1.icp, blbcn_1.par, &ijac, &fop, blwif_1.dduxx, blwif_1.ddpxx); /* TRANSPOSE */ i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { i__2 = blicn_1.ndm; for (j = 1; j <= i__2; ++j) { dfu[i + j * dfu_dim1] = dfdu[j + i * dfdu_dim1]; /* L1: */ } /* L2: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { dfu[i + blicn_1.ndmp1 * dfu_dim1] = blwif_1.dduxx[i - 1]; dfu[blicn_1.ndmp1 + i * dfu_dim1] = dfdp[i + blbcn_1.icp[1] * dfdp_dim1]; /* L3: */ } dfu[blicn_1.ndmp1 + blicn_1.ndmp1 * dfu_dim1] = blwif_1.ddpxx[ blbcn_1.icp[1] - 1]; nlvc_(&blicn_1.ndmp1, &blicn_1.ndmp1, &c__1, &dfu[dfu_offset], &v[1], &ir[1], &ic[1]); nrmlz_(&blicn_1.ndmp1, &v[1]); i__1 = blicn_1.ndmp1; for (i = 1; i <= i__1; ++i) { u[blicn_1.ndm + i] = v[i]; /* L4: */ } blbcn_1.par[blbcn_1.icp[0] - 1] = fop; } i__1 = blicn_1.nfpar - 1; for (i = 1; i <= i__1; ++i) { u[blbcn_1.ndim - blicn_1.nfpar + 1 + i] = blbcn_1.par[blbcn_1.icp[i] - 1]; /* L5: */ } return 0; } /* stpnc2_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Subroutines for Discrete Dynamical Systems */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fnds_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1; /* Local variables */ extern /* Subroutine */ int funi_(); static integer i; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generate the equations for continuing fixed points. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ funi_(ndim, &u[1], uold, &icp[1], &par[1], ijac, &f[1], &dfdu[dfdu_offset] , &dfdp[dfdp_offset]); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { f[i] -= u[i]; /* L1: */ } if (*ijac == 0) { return 0; } i__1 = *ndim; for (i = 1; i <= i__1; ++i) { dfdu[i + i * dfdu_dim1] -= blrcn_1.one; /* L2: */ } return 0; } /* fnds_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Subroutines for the Continuation of Hopf Bifurcation Points (Maps) */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fnhd_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int ffhd_(); static integer i, j; static doublereal ep, umx; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates the equations for the 2-parameter continuation of Hopf */ /* bifurcation points for maps. */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --uold; --u; /* Function Body */ ffhd_(ndim, &u[1], &uold[1], &icp[1], &par[1], &f[1], &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); if (*ijac == 0) { return 0; } /* Generate the Jacobian. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L2: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { blwif_1.u1xx[j - 1] = u[j]; blwif_1.u2xx[j - 1] = u[j]; /* L3: */ } blwif_1.u1xx[i - 1] -= ep; blwif_1.u2xx[i - 1] += ep; ffhd_(ndim, blwif_1.u1xx, &uold[1], &icp[1], &par[1], blwif_1.f1xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); ffhd_(ndim, blwif_1.u2xx, &uold[1], &icp[1], &par[1], blwif_1.f2xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[j + i * dfdu_dim1] = (blwif_1.f2xx[j - 1] - blwif_1.f1xx[j - 1]) / (ep * 2); /* L4: */ } /* L5: */ } par[icp[1]] += ep; ffhd_(ndim, &u[1], &uold[1], &icp[1], &par[1], blwif_1.f1xx, &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__1 = *ndim; for (j = 1; j <= i__1; ++j) { dfdp[j + icp[1] * dfdp_dim1] = (blwif_1.f1xx[j - 1] - f[j]) / ep; /* L6: */ } par[icp[1]] -= ep; return 0; } /* fnhd_ */ /* ---------- ---- */ /* Subroutine */ int ffhd_(ndim, u, uold, icp, par, f, ndm, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par, *f; integer *ndm; doublereal *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Builtin functions */ double sin(), cos(); /* Local variables */ static integer ijac; extern /* Subroutine */ int funi_(); static integer i, j; static doublereal theta, c1, s1; static integer ndm2; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dfdp_dim1 = *ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --uold; --u; /* Function Body */ ndm2 = *ndm << 1; theta = u[*ndim - 1]; s1 = sin(theta); /* SGLE S1= SIN(THETA) */ c1 = cos(theta); /* SGLE C1= COS(THETA) */ par[icp[2]] = u[*ndim]; ijac = 1; funi_(ndm, &u[1], &uold[1], &icp[1], &par[1], &ijac, &f[1], &dfdu[ dfdu_offset], &dfdp[dfdp_offset]); i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[i] -= u[i]; dfdu[i + i * dfdu_dim1] -= c1; /* L1: */ } i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndm + i] = s1 * u[ndm2 + i]; f[ndm2 + i] = -s1 * u[*ndm + i]; i__2 = *ndm; for (j = 1; j <= i__2; ++j) { f[*ndm + i] += dfdu[i + j * dfdu_dim1] * u[*ndm + j]; f[ndm2 + i] += dfdu[i + j * dfdu_dim1] * u[ndm2 + j]; /* L2: */ } /* L3: */ } f[*ndim - 1] = -blrcn_1.one; i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndim - 1] = f[*ndim - 1] + u[*ndm + i] * u[*ndm + i] + u[ndm2 + i] * u[ndm2 + i]; /* L4: */ } f[*ndim] = blrcn_1.zero; i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndim] = f[*ndim] + uold[ndm2 + i] * u[*ndm + i] - uold[*ndm + i] * u[ndm2 + i]; /* L5: */ } return 0; } /* ffhd_ */ /* ---------- ------ */ /* Subroutine */ int stpnhd_(ibr, u, ndm2, smat, dfdu, dfuxx, dfdp, v, f, ir, ic) integer *ibr; doublereal *u; integer *ndm2; doublereal *smat, *dfdu, *dfuxx, *dfdp, *v, *f; integer *ir, *ic; { /* System generated locals */ integer smat_dim1, smat_offset, dfdu_dim1, dfdu_offset, dfuxx_dim1, dfuxx_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Builtin functions */ double sin(), cos(); /* Local variables */ static integer ijac; extern /* Subroutine */ int funi_(), nlvc_(); static doublereal uold; static integer i, j; static doublereal theta; static logical found; static doublereal c1; extern /* Subroutine */ int nrmlz_(); static doublereal s1; extern /* Subroutine */ int readl3_(), findl3_(); static integer nfpar1; extern doublereal pi_(); static integer itp; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates starting data for the continuation of Hopf bifurcation */ /* points for maps. */ /* Parameter adjustments */ --ic; --ir; --f; --v; dfdp_dim1 = blicn_1.ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfuxx_dim1 = blicn_1.ndmp1; dfuxx_offset = dfuxx_dim1 + 1; dfuxx -= dfuxx_offset; dfdu_dim1 = blicn_1.ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; smat_dim1 = *ndm2; smat_offset = smat_dim1 + 1; smat -= smat_offset; --u; /* Function Body */ findl3_(&blbcn_1.irs, &itp, &nfpar1, &found); readl3_(&blbcn_1.ips, ibr, &u[1], blbcn_1.par); theta = pi_(&blrcn_1.two) / blbcn_1.par[10]; s1 = sin(theta); /* SGLE S1= SIN(THETA) */ c1 = cos(theta); /* SGLE C1= COS(THETA) */ ijac = 1; funi_(&blicn_1.ndm, &u[1], &uold, blbcn_1.icp, blbcn_1.par, &ijac, &f[1], &dfdu[dfdu_offset], &dfdp[dfdp_offset]); *ndm2 = blicn_1.ndm << 1; i__1 = *ndm2; for (i = 1; i <= i__1; ++i) { i__2 = *ndm2; for (j = 1; j <= i__2; ++j) { smat[i + j * smat_dim1] = blrcn_1.zero; /* L1: */ } /* L2: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { smat[i + (blicn_1.ndm + i) * smat_dim1] = s1; /* L3: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { smat[blicn_1.ndm + i + i * smat_dim1] = -s1; /* L4: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { i__2 = blicn_1.ndm; for (j = 1; j <= i__2; ++j) { smat[i + j * smat_dim1] = dfdu[i + j * dfdu_dim1]; smat[blicn_1.ndm + i + (blicn_1.ndm + j) * smat_dim1] = dfdu[i + j * dfdu_dim1]; /* L5: */ } smat[i + i * smat_dim1] -= c1; smat[blicn_1.ndm + i + (blicn_1.ndm + i) * smat_dim1] -= c1; /* L6: */ } nlvc_(ndm2, ndm2, &c__2, &smat[smat_offset], &v[1], &ir[1], &ic[1]); nrmlz_(ndm2, &v[1]); i__1 = *ndm2; for (i = 1; i <= i__1; ++i) { u[blicn_1.ndm + i] = v[i]; /* L7: */ } u[blbcn_1.ndim - 1] = theta; u[blbcn_1.ndim] = blbcn_1.par[blbcn_1.icp[1] - 1]; return 0; } /* stpnhd_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Subroutines for the Continuation of Hopf Bifurcation Points (ODE) */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fnhb_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int ffhb_(); static integer i, j; static doublereal ep, umx; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates the equations for the 2-parameter continuation of Hopf */ /* bifurcation points in ODE. */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --uold; --u; /* Function Body */ ffhb_(ndim, &u[1], &uold[1], &icp[1], &par[1], &f[1], &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); if (*ijac == 0) { return 0; } /* Generate the Jacobian. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L2: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { blwif_1.u1xx[j - 1] = u[j]; blwif_1.u2xx[j - 1] = u[j]; /* L3: */ } blwif_1.u1xx[i - 1] -= ep; blwif_1.u2xx[i - 1] += ep; ffhb_(ndim, blwif_1.u1xx, &uold[1], &icp[1], &par[1], blwif_1.f1xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); ffhb_(ndim, blwif_1.u2xx, &uold[1], &icp[1], &par[1], blwif_1.f2xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[j + i * dfdu_dim1] = (blwif_1.f2xx[j - 1] - blwif_1.f1xx[j - 1]) / (ep * 2); /* L4: */ } /* L5: */ } par[icp[1]] += ep; ffhb_(ndim, &u[1], &uold[1], &icp[1], &par[1], blwif_1.f1xx, &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__1 = *ndim; for (j = 1; j <= i__1; ++j) { dfdp[j + icp[1] * dfdp_dim1] = (blwif_1.f1xx[j - 1] - f[j]) / ep; /* L6: */ } par[icp[1]] -= ep; return 0; } /* fnhb_ */ /* ---------- ---- */ /* Subroutine */ int ffhb_(ndim, u, uold, icp, par, f, ndm, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par, *f; integer *ndm; doublereal *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ static integer ijac; extern /* Subroutine */ int funi_(); static integer i, j; extern doublereal pi_(); static doublereal rom; static integer ndm2; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dfdp_dim1 = *ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --uold; --u; /* Function Body */ ndm2 = *ndm << 1; rom = u[*ndim - 1]; par[11] = rom * pi_(&blrcn_1.two); par[icp[2]] = u[*ndim]; ijac = 1; funi_(ndm, &u[1], &uold[1], &icp[1], &par[1], &ijac, &f[1], &dfdu[ dfdu_offset], &dfdp[dfdp_offset]); i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndm + i] = u[ndm2 + i]; f[ndm2 + i] = -u[*ndm + i]; i__2 = *ndm; for (j = 1; j <= i__2; ++j) { f[*ndm + i] += rom * dfdu[i + j * dfdu_dim1] * u[*ndm + j]; f[ndm2 + i] += rom * dfdu[i + j * dfdu_dim1] * u[ndm2 + j]; /* L1: */ } /* L2: */ } f[*ndim - 1] = -blrcn_1.one; i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndim - 1] = f[*ndim - 1] + u[*ndm + i] * u[*ndm + i] + u[ndm2 + i] * u[ndm2 + i]; /* L3: */ } f[*ndim] = blrcn_1.zero; i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndim] = f[*ndim] + uold[ndm2 + i] * (u[*ndm + i] - uold[*ndm + i]) - uold[*ndm + i] * (u[ndm2 + i] - uold[ndm2 + i]); /* L4: */ } return 0; } /* ffhb_ */ /* ---------- ------ */ /* Subroutine */ int stpnhb_(ibr, u, ndm2, smat, dfdu, dfuxx, dfdp, v, f, ir, ic) integer *ibr; doublereal *u; integer *ndm2; doublereal *smat, *dfdu, *dfuxx, *dfdp, *v, *f; integer *ir, *ic; { /* System generated locals */ integer smat_dim1, smat_offset, dfdu_dim1, dfdu_offset, dfuxx_dim1, dfuxx_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ static integer ijac; extern /* Subroutine */ int funi_(), nlvc_(); static doublereal uold; static integer i, j; static logical found; extern /* Subroutine */ int nrmlz_(), readl3_(), findl3_(); static integer nfpar1; extern doublereal pi_(); static doublereal rho; static integer itp; static doublereal rom; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates starting data for the 2-parameter continuation of */ /* Hopf bifurcation point (ODE). */ /* Parameter adjustments */ --ic; --ir; --f; --v; dfdp_dim1 = blicn_1.ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfuxx_dim1 = blicn_1.ndmp1; dfuxx_offset = dfuxx_dim1 + 1; dfuxx -= dfuxx_offset; dfdu_dim1 = blicn_1.ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; smat_dim1 = *ndm2; smat_offset = smat_dim1 + 1; smat -= smat_offset; --u; /* Function Body */ findl3_(&blbcn_1.irs, &itp, &nfpar1, &found); readl3_(&blbcn_1.ips, ibr, &u[1], blbcn_1.par); ijac = 1; rho = blbcn_1.par[10]; rom = rho / pi_(&blrcn_1.two); funi_(&blicn_1.ndm, &u[1], &uold, blbcn_1.icp, blbcn_1.par, &ijac, &f[1], &dfdu[dfdu_offset], &dfdp[dfdp_offset]); *ndm2 = blicn_1.ndm << 1; i__1 = *ndm2; for (i = 1; i <= i__1; ++i) { i__2 = *ndm2; for (j = 1; j <= i__2; ++j) { smat[i + j * smat_dim1] = blrcn_1.zero; /* L1: */ } /* L2: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { smat[i + (blicn_1.ndm + i) * smat_dim1] = blrcn_1.one; /* L3: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { smat[blicn_1.ndm + i + i * smat_dim1] = -blrcn_1.one; /* L4: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { i__2 = blicn_1.ndm; for (j = 1; j <= i__2; ++j) { smat[i + j * smat_dim1] = rom * dfdu[i + j * dfdu_dim1]; smat[blicn_1.ndm + i + (blicn_1.ndm + j) * smat_dim1] = rom * dfdu[i + j * dfdu_dim1]; /* L5: */ } /* L6: */ } nlvc_(ndm2, ndm2, &c__2, &smat[smat_offset], &v[1], &ir[1], &ic[1]); nrmlz_(ndm2, &v[1]); i__1 = *ndm2; for (i = 1; i <= i__1; ++i) { u[blicn_1.ndm + i] = v[i]; /* L7: */ } u[blbcn_1.ndim - 1] = rom; u[blbcn_1.ndim] = blbcn_1.par[blbcn_1.icp[1] - 1]; return 0; } /* stpnhb_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Subroutines for the Continuation of Hopf Bifurcation Points (Waves) */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fnhw_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int ffhw_(); static integer i, j; static doublereal ep, umx; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates the equations for the 2-parameter continuation of a */ /* bifurcation to a traveling wave. */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --uold; --u; /* Function Body */ ffhw_(ndim, &u[1], &uold[1], &icp[1], &par[1], &f[1], &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); if (*ijac == 0) { return 0; } /* Generate the Jacobian. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L2: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { blwif_1.u1xx[j - 1] = u[j]; blwif_1.u2xx[j - 1] = u[j]; /* L3: */ } blwif_1.u1xx[i - 1] -= ep; blwif_1.u2xx[i - 1] += ep; ffhw_(ndim, blwif_1.u1xx, &uold[1], &icp[1], &par[1], blwif_1.f1xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); ffhw_(ndim, blwif_1.u2xx, &uold[1], &icp[1], &par[1], blwif_1.f2xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[j + i * dfdu_dim1] = (blwif_1.f2xx[j - 1] - blwif_1.f1xx[j - 1]) / (ep * 2); /* L4: */ } /* L5: */ } par[icp[1]] += ep; ffhw_(ndim, &u[1], &uold[1], &icp[1], &par[1], blwif_1.f1xx, &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__1 = *ndim; for (j = 1; j <= i__1; ++j) { dfdp[j + icp[1] * dfdp_dim1] = (blwif_1.f1xx[j - 1] - f[j]) / ep; /* L6: */ } par[icp[1]] -= ep; return 0; } /* fnhw_ */ /* ---------- ---- */ /* Subroutine */ int ffhw_(ndim, u, uold, icp, par, f, ndm, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par, *f; integer *ndm; doublereal *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ static integer ijac; extern /* Subroutine */ int fnws_(); static integer i, j; static doublereal rom; static integer ndm2; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dfdp_dim1 = *ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --uold; --u; /* Function Body */ ndm2 = *ndm << 1; rom = u[*ndim - 1]; par[icp[2]] = u[*ndim]; ijac = 1; fnws_(ndm, &u[1], &uold[1], &icp[1], &par[1], &ijac, &f[1], &dfdu[ dfdu_offset], &dfdp[dfdp_offset]); i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndm + i] = u[ndm2 + i]; f[ndm2 + i] = -u[*ndm + i]; i__2 = *ndm; for (j = 1; j <= i__2; ++j) { f[*ndm + i] += rom * dfdu[i + j * dfdu_dim1] * u[*ndm + j]; f[ndm2 + i] += rom * dfdu[i + j * dfdu_dim1] * u[ndm2 + j]; /* L1: */ } /* L2: */ } f[*ndim - 1] = -blrcn_1.one; i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndim - 1] = f[*ndim - 1] + u[*ndm + i] * u[*ndm + i] + u[ndm2 + i] * u[ndm2 + i]; /* L3: */ } f[*ndim] = blrcn_1.zero; i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndim] = f[*ndim] + uold[ndm2 + i] * (u[*ndm + i] - uold[*ndm + i]) - uold[*ndm + i] * (u[ndm2 + i] - uold[ndm2 + i]); /* L4: */ } return 0; } /* ffhw_ */ /* ---------- ------ */ /* Subroutine */ int stpnhw_(ibr, u, ndm2, smat, dfdu, dfuxx, dfdp, v, f, ir, ic) integer *ibr; doublereal *u; integer *ndm2; doublereal *smat, *dfdu, *dfuxx, *dfdp, *v, *f; integer *ir, *ic; { /* System generated locals */ integer smat_dim1, smat_offset, dfdu_dim1, dfdu_offset, dfuxx_dim1, dfuxx_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ static integer ijac; extern /* Subroutine */ int nlvc_(); static doublereal uold; extern /* Subroutine */ int fnws_(); static integer i, j; static logical found; extern /* Subroutine */ int nrmlz_(), readl3_(), findl3_(); static integer nfpar1; extern doublereal pi_(); static doublereal rho; static integer itp; static doublereal rom; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates starting data for the continuation of a bifurcation to a */ /* traveling wave. */ /* Parameter adjustments */ --ic; --ir; --f; --v; dfdp_dim1 = blicn_1.ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfuxx_dim1 = blicn_1.ndmp1; dfuxx_offset = dfuxx_dim1 + 1; dfuxx -= dfuxx_offset; dfdu_dim1 = blicn_1.ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; smat_dim1 = *ndm2; smat_offset = smat_dim1 + 1; smat -= smat_offset; --u; /* Function Body */ findl3_(&blbcn_1.irs, &itp, &nfpar1, &found); readl3_(&blbcn_1.ips, ibr, &u[1], blbcn_1.par); ijac = 1; rho = blbcn_1.par[10]; rom = rho / pi_(&blrcn_1.two); fnws_(&blicn_1.ndm, &u[1], &uold, blbcn_1.icp, blbcn_1.par, &ijac, &f[1], &dfdu[dfdu_offset], &dfdp[dfdp_offset]); *ndm2 = blicn_1.ndm << 1; i__1 = *ndm2; for (i = 1; i <= i__1; ++i) { i__2 = *ndm2; for (j = 1; j <= i__2; ++j) { smat[i + j * smat_dim1] = blrcn_1.zero; /* L1: */ } /* L2: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { smat[i + (blicn_1.ndm + i) * smat_dim1] = blrcn_1.one; /* L3: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { smat[blicn_1.ndm + i + i * smat_dim1] = -blrcn_1.one; /* L4: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { i__2 = blicn_1.ndm; for (j = 1; j <= i__2; ++j) { smat[i + j * smat_dim1] = rom * dfdu[i + j * dfdu_dim1]; smat[blicn_1.ndm + i + (blicn_1.ndm + j) * smat_dim1] = rom * dfdu[i + j * dfdu_dim1]; /* L5: */ } /* L6: */ } nlvc_(ndm2, ndm2, &c__2, &smat[smat_offset], &v[1], &ir[1], &ic[1]); nrmlz_(ndm2, &v[1]); i__1 = *ndm2; for (i = 1; i <= i__1; ++i) { u[blicn_1.ndm + i] = v[i]; /* L7: */ } u[blbcn_1.ndim - 1] = rom; u[blbcn_1.ndim] = blbcn_1.par[blbcn_1.icp[1] - 1]; return 0; } /* stpnhw_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Subroutines for the Continuation of Periodic Limit Points */ /* and Period Doubling Bifurcations */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fnpl_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int ffpl_(); static integer i, j; static doublereal ep, umx; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Subroutines for the 2-parameter continuation of of limit points */ /* on branches of periodic solutions. */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ ffpl_(ndim, &u[1], uold, &icp[1], &par[1], &f[1], &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); if (*ijac == 0) { return 0; } /* Generate the Jacobian. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L1: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { blwif_1.u1xx[j - 1] = u[j]; blwif_1.u2xx[j - 1] = u[j]; /* L2: */ } blwif_1.u1xx[i - 1] -= ep; blwif_1.u2xx[i - 1] += ep; ffpl_(ndim, blwif_1.u1xx, uold, &icp[1], &par[1], blwif_1.f1xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); ffpl_(ndim, blwif_1.u2xx, uold, &icp[1], &par[1], blwif_1.f2xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[j + i * dfdu_dim1] = (blwif_1.f2xx[j - 1] - blwif_1.f1xx[j - 1]) / (ep * 2); /* L3: */ } /* L4: */ } i__1 = blicn_1.nfpar; for (i = 1; i <= i__1; ++i) { par[icp[i]] += ep; ffpl_(ndim, &u[1], uold, &icp[1], &par[1], blwif_1.f1xx, &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdp[j + icp[i] * dfdp_dim1] = (blwif_1.f1xx[j - 1] - f[j]) / ep; /* L5: */ } par[icp[i]] -= ep; /* L6: */ } return 0; } /* fnpl_ */ /* ---------- ---- */ /* Subroutine */ int ffpl_(ndim, u, uold, icp, par, f, ndm, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par, *f; integer *ndm; doublereal *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ static integer ijac; extern /* Subroutine */ int funi_(); static integer i, j; static doublereal c1, c2; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dfdp_dim1 = *ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ ijac = 1; c1 = par[icp[3]]; /* SGLE C1=PAR(ICP(3)) */ c2 = par[icp[4]]; /* SGLE C2=PAR(ICP(4)) */ funi_(ndm, &u[1], uold, &icp[1], &par[1], &ijac, &f[1], &dfdu[dfdu_offset] , &dfdp[dfdp_offset]); i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndm + i] = blrcn_1.zero; i__2 = *ndm; for (j = 1; j <= i__2; ++j) { f[*ndm + i] += dfdu[i + j * dfdu_dim1] * u[*ndm + j]; /* L1: */ } f[*ndm + i] = c1 * f[*ndm + i] + c2 * f[i]; f[i] = c1 * f[i]; /* L2: */ } return 0; } /* ffpl_ */ /* ---------- ---- */ /* Subroutine */ int bcpl_(ndim, par, icp, nbc, u0, u1, f, ijac, dbc) integer *ndim; doublereal *par; integer *icp, *nbc; doublereal *u0, *u1, *f; integer *ijac; doublereal *dbc; { /* System generated locals */ integer dbc_dim1, dbc_offset, i__1, i__2; /* Local variables */ static integer i, j, nn; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dbc_dim1 = *nbc; dbc_offset = dbc_dim1 + 1; dbc -= dbc_offset; --f; --u1; --u0; --icp; --par; /* Function Body */ i__1 = *ndim; for (i = 1; i <= i__1; ++i) { f[i] = u0[i] - u1[i]; /* L1: */ } if (blrtn.irot != 0) { for (i= 1; i <= i__1; ++i) { if (blrtn.nrot[i- 1] != 0) { f[i] += blrtn.torper* blrtn.nrot[i- 1]; } } } if (*ijac == 0) { return 0; } nn = (*ndim << 1) + blicn_1.npar; i__1 = *nbc; for (i = 1; i <= i__1; ++i) { i__2 = nn; for (j = 1; j <= i__2; ++j) { dbc[i + j * dbc_dim1] = blrcn_1.zero; /* L2: */ } /* L3: */ } i__1 = *ndim; for (i = 1; i <= i__1; ++i) { dbc[i + i * dbc_dim1] = blrcn_1.one; dbc[i + (*ndim + i) * dbc_dim1] = -blrcn_1.one; /* L4: */ } return 0; } /* bcpl_ */ /* ---------- ---- */ /* Subroutine */ int bcpd_(ndim, par, icp, nbc, u0, u1, f, ijac, dbc) integer *ndim; doublereal *par; integer *icp, *nbc; doublereal *u0, *u1, *f; integer *ijac; doublereal *dbc; { /* System generated locals */ integer dbc_dim1, dbc_offset, i__1, i__2; /* Local variables */ static integer i, j, nn; /* Generate boundary conditions for the 2-parameter continuation */ /* of period doubling bifurcations. */ /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dbc_dim1 = *nbc; dbc_offset = dbc_dim1 + 1; dbc -= dbc_offset; --f; --u1; --u0; --icp; --par; /* Function Body */ i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { f[i] = u0[i] - u1[i]; f[blicn_1.ndm + i] = u0[blicn_1.ndm + i] + u1[blicn_1.ndm + i]; /* L1: */ } if (blrtn.irot != 0) { for (i= 1; i <= i__1; ++i) { if (blrtn.nrot[i- 1] != 0) { f[i] += blrtn.torper* blrtn.nrot[i- 1]; } } } if (*ijac == 0) { return 0; } nn = (*ndim << 1) + blicn_1.npar; i__1 = *nbc; for (i = 1; i <= i__1; ++i) { i__2 = nn; for (j = 1; j <= i__2; ++j) { dbc[i + j * dbc_dim1] = blrcn_1.zero; /* L2: */ } /* L3: */ } i__1 = *ndim; for (i = 1; i <= i__1; ++i) { dbc[i + i * dbc_dim1] = blrcn_1.one; if (i <= blicn_1.ndm) { dbc[i + (*ndim + i) * dbc_dim1] = -blrcn_1.one; } else { dbc[i + (*ndim + i) * dbc_dim1] = blrcn_1.one; } /* L4: */ } return 0; } /* bcpd_ */ /* ---------- ---- */ /* Subroutine */ int icpl_(ndim, par, icp, nint, u, uold, udot, upold, f, ijac, dint) integer *ndim; doublereal *par; integer *icp, *nint; doublereal *u, *uold, *udot, *upold, *f; integer *ijac; doublereal *dint; { /* System generated locals */ integer dint_dim1, dint_offset, i__1, i__2; doublereal d__1; /* Local variables */ static integer i, j, nn; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dint_dim1 = *nint; dint_offset = dint_dim1 + 1; dint -= dint_offset; --f; --upold; --udot; --uold; --u; --icp; --par; /* Function Body */ f[1] = blrcn_1.zero; f[2] = blrcn_1.zero; /* Computing 2nd power */ d__1 = par[icp[4]]; f[3] = d__1 * d__1 - blrcn_1.one; i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { if (blrtn.nrot[i - 1] == 0) { f[1] += u[i] * upold[i]; f[2] += u[blicn_1.ndm + i] * upold[i];} f[3] += u[blicn_1.ndm + i] * u[blicn_1.ndm + i]; /* L1: */ } if (*ijac == 0) { return 0; } nn = *ndim + blicn_1.npar; i__1 = *nint; for (i = 1; i <= i__1; ++i) { i__2 = nn; for (j = 1; j <= i__2; ++j) { dint[i + j * dint_dim1] = blrcn_1.zero; /* L2: */ } /* L3: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { if (blrtn.nrot[i- 1] == 0) { dint[i * dint_dim1 + 1] = upold[i]; dint[(blicn_1.ndm + i) * dint_dim1 + 2] = upold[i]; } else { dint[i * dint_dim1 + 1] = 0; dint[(blicn_1.ndm + i) * dint_dim1 + 2] = 0; } dint[(blicn_1.ndm + i) * dint_dim1 + 3] = blrcn_1.two * u[blicn_1.ndm + i]; /* L4: */ } dint[(*ndim + icp[4]) * dint_dim1 + 3] = blrcn_1.two * par[icp[4]]; return 0; } /* icpl_ */ /* ---------- ------ */ /* Subroutine */ int stpnpl_(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]; static integer ntpl1, nrsp1, ntot1, i, j, k; static logical found; static integer icprs[20], k1, k2; extern /* Subroutine */ int findl3_(); static integer nfpar1, nskip1, lab1, nar1, itp1, isw1; /* Fortran I/O blocks */ static cilist io___113 = { 0, 3, 0, 0, 0 }; static cilist io___126 = { 0, 3, 0, fmt_101, 0 }; static cilist io___129 = { 0, 3, 0, fmt_101, 0 }; static cilist io___130 = { 0, 3, 0, fmt_101, 0 }; static cilist io___131 = { 0, 3, 0, fmt_101, 0 }; static cilist io___132 = { 0, 3, 0, fmt_101, 0 }; static cilist io___133 = { 0, 3, 0, fmt_101, 0 }; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates starting data for the 2-parameter continuation of limit */ /* points on a branch of periodic solutions. */ /* 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___113); 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 *)&lab1, (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)); do_lio(&c__3, &c__1, (char *)&icprs[0], (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&icprs[1], (ftnlen)sizeof(integer)); e_rsle(); nrsp1 = *ntstrs + 1; 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 = k1 + blicn_1.ndm - 1; s_rsfe(&io___126); 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(); /* L1: */ } tm[j] = temp[0]; /* L2: */ } s_rsfe(&io___129); do_fio(&c__1, (char *)&tm[nrsp1], (ftnlen)sizeof(doublereal)); i__1 = blicn_1.ndm; for (k = 1; k <= i__1; ++k) { do_fio(&c__1, (char *)&ups[nrsp1 + k * ups_dim1], (ftnlen)sizeof( doublereal)); } e_rsfe(); s_rsfe(&io___130); do_fio(&c__1, (char *)&blcrl_1.rldot[0], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&blcrl_1.rldot[1], (ftnlen)sizeof(doublereal)); e_rsfe(); /* Read U-dot (derivative with respect to arclength). */ 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 + blicn_1.ndm + 1; k2 = i * blbcn_1.ndim; s_rsfe(&io___131); 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(); /* L3: */ } /* L4: */ } k1 = blicn_1.ndm + 1; s_rsfe(&io___132); i__1 = blbcn_1.ndim; for (k = k1; k <= i__1; ++k) { do_fio(&c__1, (char *)&ups[nrsp1 + k * ups_dim1], (ftnlen)sizeof( doublereal)); } e_rsfe(); /* Read the parameter values. */ s_rsfe(&io___133); 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(); blbcn_1.par[blbcn_1.icp[3] - 1] = blcrl_1.rldot[1]; 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: */ } *nodir = 1; return 0; } /* stpnpl_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Subroutines for the Continuation of Limit Points for BVP. */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fnbl_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int ffbl_(); static integer i, j; static doublereal ep, umx; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates the equations for the 2-parameter continuation */ /* of limit points (BVP). */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ ffbl_(ndim, &u[1], uold, &icp[1], &par[1], &f[1], blwif_1.dfuxx, blwif_1.dfpxx); if (*ijac == 0) { return 0; } /* Generate the Jacobian. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L1: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { blwif_1.u1xx[j - 1] = u[j]; blwif_1.u2xx[j - 1] = u[j]; /* L2: */ } blwif_1.u1xx[i - 1] -= ep; blwif_1.u2xx[i - 1] += ep; ffbl_(ndim, blwif_1.u1xx, uold, &icp[1], &par[1], blwif_1.f1xx, blwif_1.dfuxx, blwif_1.dfpxx); ffbl_(ndim, blwif_1.u2xx, uold, &icp[1], &par[1], blwif_1.f2xx, blwif_1.dfuxx, blwif_1.dfpxx); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[j + i * dfdu_dim1] = (blwif_1.f2xx[j - 1] - blwif_1.f1xx[j - 1]) / (ep * 2); /* L3: */ } /* L4: */ } i__1 = blicn_1.nfpar; for (i = 1; i <= i__1; ++i) { par[icp[i]] += ep; ffbl_(ndim, &u[1], uold, &icp[1], &par[1], blwif_1.f1xx, blwif_1.dfuxx, blwif_1.dfpxx); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdp[j + icp[i] * dfdp_dim1] = (blwif_1.f1xx[j - 1] - f[j]) / ep; /* L5: */ } par[icp[i]] -= ep; /* L6: */ } return 0; } /* fnbl_ */ /* ---------- ---- */ /* Subroutine */ int ffbl_(ndim, u, uold, icp, par, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par, *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ static integer ijac; extern /* Subroutine */ int funi_(); static integer nfpx, i, j; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dfdp_dim1 = blicn_1.ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = blicn_1.ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ ijac = 1; funi_(&blicn_1.ndm, &u[1], uold, &icp[1], &par[1], &ijac, &f[1], &dfdu[ dfdu_offset], &dfdp[dfdp_offset]); nfpx = blicn_1.nfpar / 2 - 1; i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { f[blicn_1.ndm + i] = blrcn_1.zero; i__2 = blicn_1.ndm; for (j = 1; j <= i__2; ++j) { f[blicn_1.ndm + i] += dfdu[i + j * dfdu_dim1] * u[blicn_1.ndm + j] ; /* L1: */ } if (nfpx > 0) { i__2 = nfpx; for (j = 1; j <= i__2; ++j) { f[blicn_1.ndm + i] += dfdp[i + icp[j + 1] * dfdp_dim1] * par[ icp[blicn_1.nfpar - nfpx + j]]; /* L2: */ } } /* L3: */ } return 0; } /* ffbl_ */ /* ---------- ---- */ /* Subroutine */ int bcbl_(ndim, par, icp, nbc, u0, u1, f, ijac, dbc) integer *ndim; doublereal *par; integer *icp, *nbc; doublereal *u0, *u1, *f; integer *ijac; doublereal *dbc; { /* System generated locals */ integer dbc_dim1, dbc_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int fbbl_(); static integer i, j; static doublereal ep, umx; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates the boundary conditions for the 2-parameter continuation */ /* of limit points (BVP). */ /* Generate the function. */ /* Parameter adjustments */ dbc_dim1 = *nbc; dbc_offset = dbc_dim1 + 1; dbc -= dbc_offset; --f; --u1; --u0; --icp; --par; /* Function Body */ fbbl_(ndim, &par[1], &icp[1], nbc, &u0[1], &u1[1], &f[1], blwif_1.dfuxx); if (*ijac == 0) { return 0; } /* Derivatives with respect to U0. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u0[i], abs(d__1)) > umx) { umx = (d__2 = u0[i], abs(d__2)); } /* SGLE IF( ABS(U0(I)).GT.UMX)UMX= ABS(U0(I)) */ /* L1: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { blwif_1.u1xx[j - 1] = u0[j]; blwif_1.u2xx[j - 1] = u0[j]; /* L2: */ } blwif_1.u1xx[i - 1] -= ep; blwif_1.u2xx[i - 1] += ep; fbbl_(ndim, &par[1], &icp[1], nbc, blwif_1.u1xx, &u1[1], blwif_1.f1xx, blwif_1.dfuxx); fbbl_(ndim, &par[1], &icp[1], nbc, blwif_1.u2xx, &u1[1], blwif_1.f2xx, blwif_1.dfuxx); i__2 = *nbc; for (j = 1; j <= i__2; ++j) { dbc[j + i * dbc_dim1] = (blwif_1.f2xx[j - 1] - blwif_1.f1xx[j - 1] ) / (ep * 2); /* L3: */ } /* L4: */ } /* Derivatives with respect to U1. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u1[i], abs(d__1)) > umx) { umx = (d__2 = u1[i], abs(d__2)); } /* SGLE IF( ABS(U1(I)).GT.UMX)UMX= ABS(U1(I)) */ /* L5: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { blwif_1.u1xx[j - 1] = u1[j]; blwif_1.u2xx[j - 1] = u1[j]; /* L6: */ } blwif_1.u1xx[i - 1] -= ep; blwif_1.u2xx[i - 1] += ep; fbbl_(ndim, &par[1], &icp[1], nbc, &u0[1], blwif_1.u1xx, blwif_1.f1xx, blwif_1.dfuxx); fbbl_(ndim, &par[1], &icp[1], nbc, &u0[1], blwif_1.u2xx, blwif_1.f2xx, blwif_1.dfuxx); i__2 = *nbc; for (j = 1; j <= i__2; ++j) { dbc[j + (*ndim + i) * dbc_dim1] = (blwif_1.f2xx[j - 1] - blwif_1.f1xx[j - 1]) / (ep * 2); /* L7: */ } /* L8: */ } i__1 = blicn_1.nfpar; for (i = 1; i <= i__1; ++i) { par[icp[i]] += ep; fbbl_(ndim, &par[1], &icp[1], nbc, &u0[1], &u1[1], blwif_1.f2xx, blwif_1.dfuxx); i__2 = *nbc; for (j = 1; j <= i__2; ++j) { dbc[j + ((*ndim << 1) + icp[i]) * dbc_dim1] = (blwif_1.f2xx[j - 1] - f[j]) / ep; /* L9: */ } par[icp[i]] -= ep; /* L10: */ } return 0; } /* bcbl_ */ /* ---------- ---- */ /* Subroutine */ int fbbl_(ndim, par, icp, nbc, u0, u1, f, dbc) integer *ndim; doublereal *par; integer *icp, *nbc; doublereal *u0, *u1, *f, *dbc; { /* System generated locals */ integer dbc_dim1, dbc_offset, i__1, i__2; /* Local variables */ static integer ijac; extern /* Subroutine */ int bcni_(); static integer nfpx, i, j; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dbc_dim1 = blicn_1.nbc0; dbc_offset = dbc_dim1 + 1; dbc -= dbc_offset; --f; --u1; --u0; --icp; --par; /* Function Body */ ijac = 1; nfpx = blicn_1.nfpar / 2 - 1; bcni_(&blicn_1.ndm, &par[1], &icp[1], &blicn_1.nbc0, &u0[1], &u1[1], &f[1] , &ijac, &dbc[dbc_offset]); i__1 = blicn_1.nbc0; for (i = 1; i <= i__1; ++i) { f[blicn_1.nbc0 + i] = blrcn_1.zero; i__2 = blicn_1.ndm; for (j = 1; j <= i__2; ++j) { f[blicn_1.nbc0 + i] += dbc[i + j * dbc_dim1] * u0[blicn_1.ndm + j] ; f[blicn_1.nbc0 + i] += dbc[i + (blicn_1.ndm + j) * dbc_dim1] * u1[ blicn_1.ndm + j]; /* L1: */ } if (nfpx != 0) { i__2 = nfpx; for (j = 1; j <= i__2; ++j) { f[blicn_1.nbc0 + i] += dbc[i + (*ndim + icp[j + 1]) * dbc_dim1] * par[icp[blicn_1.nfpar - nfpx + j]]; /* L2: */ } } /* L3: */ } return 0; } /* fbbl_ */ /* ---------- ---- */ /* Subroutine */ int icbl_(ndim, par, icp, nint, u, uold, udot, upold, f, ijac, dint) integer *ndim; doublereal *par; integer *icp, *nint; doublereal *u, *uold, *udot, *upold, *f; integer *ijac; doublereal *dint; { /* System generated locals */ integer dint_dim1, dint_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int fibl_(); static integer i, j; static doublereal ep, umx; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates integral conditions for the 2-parameter continuation of */ /* limit points (BVP). */ /* Generate the function. */ /* Parameter adjustments */ dint_dim1 = *nint; dint_offset = dint_dim1 + 1; dint -= dint_offset; --f; --upold; --udot; --uold; --u; --icp; --par; /* Function Body */ fibl_(ndim, &par[1], &icp[1], nint, &u[1], &uold[1], &udot[1], &upold[1], &f[1], blwif_1.dfuxx); if (*ijac == 0) { return 0; } /* Generate the Jacobian. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L1: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { blwif_1.u1xx[j - 1] = u[j]; blwif_1.u2xx[j - 1] = u[j]; /* L2: */ } blwif_1.u1xx[i - 1] -= ep; blwif_1.u2xx[i - 1] += ep; fibl_(ndim, &par[1], &icp[1], nint, blwif_1.u1xx, &uold[1], &udot[1], &upold[1], blwif_1.f1xx, blwif_1.dfuxx); fibl_(ndim, &par[1], &icp[1], nint, blwif_1.u2xx, &uold[1], &udot[1], &upold[1], blwif_1.f2xx, blwif_1.dfuxx); i__2 = *nint; for (j = 1; j <= i__2; ++j) { dint[j + i * dint_dim1] = (blwif_1.f2xx[j - 1] - blwif_1.f1xx[j - 1]) / (ep * 2); /* L3: */ } /* L4: */ } i__1 = blicn_1.nfpar; for (i = 1; i <= i__1; ++i) { par[icp[i]] += ep; fibl_(ndim, &par[1], &icp[1], nint, &u[1], &uold[1], &udot[1], &upold[ 1], blwif_1.f1xx, blwif_1.dfuxx); i__2 = *nint; for (j = 1; j <= i__2; ++j) { dint[j + (*ndim + icp[i]) * dint_dim1] = (blwif_1.f1xx[j - 1] - f[ j]) / ep; /* L5: */ } par[icp[i]] -= ep; /* L6: */ } return 0; } /* icbl_ */ /* ---------- ---- */ /* Subroutine */ int fibl_(ndim, par, icp, nint, u, uold, udot, upold, f, dint) integer *ndim; doublereal *par; integer *icp, *nint; doublereal *u, *uold, *udot, *upold, *f, *dint; { /* System generated locals */ integer dint_dim1, dint_offset, i__1, i__2; doublereal d__1; /* Local variables */ static integer ijac; extern /* Subroutine */ int icni_(); static integer nfpx, i, j; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dint_dim1 = blicn_1.nint0; dint_offset = dint_dim1 + 1; dint -= dint_offset; --f; --upold; --udot; --uold; --u; --icp; --par; /* Function Body */ if (blicn_1.nint0 > 0) { nfpx = blicn_1.nfpar / 2 - 1; ijac = 1; icni_(&blicn_1.ndm, &par[1], &icp[1], &blicn_1.nint0, &u[1], &uold[1], &udot[1], &upold[1], &f[1], &ijac, &dint[dint_offset]); i__1 = blicn_1.nint0; for (i = 1; i <= i__1; ++i) { f[blicn_1.nint0 + i] = blrcn_1.zero; i__2 = blicn_1.ndm; for (j = 1; j <= i__2; ++j) { f[blicn_1.nint0 + i] += dint[i + j * dint_dim1] * u[ blicn_1.ndm + j]; /* L1: */ } if (nfpx != 0) { i__2 = nfpx; for (j = 1; j <= i__2; ++j) { f[blicn_1.nint0 + i] += dint[i + (blicn_1.ndm + icp[j + 1] ) * dint_dim1] * par[icp[blicn_1.nfpar - nfpx + j] ]; /* L2: */ } } /* L3: */ } } f[*nint] = -blrcn_1.one; i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { f[*nint] += u[blicn_1.ndm + i] * u[blicn_1.ndm + i]; /* L4: */ } if (nfpx != 0) { i__1 = nfpx; for (i = 1; i <= i__1; ++i) { /* Computing 2nd power */ d__1 = par[icp[blicn_1.nfpar - nfpx + i]]; f[*nint] += d__1 * d__1; /* L5: */ } } return 0; } /* fibl_ */ /* ---------- ------ */ /* Subroutine */ int stpnbl_(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]; static integer nfpx, ntpl1, nrsp1, ntot1, i, j, k; static logical found; static integer icprs[20], k1, k2; extern /* Subroutine */ int findl3_(); static integer nfpar0, nfpar1, nskip1, lab1, nar1, itp1, isw1; /* Fortran I/O blocks */ static cilist io___161 = { 0, 3, 0, 0, 0 }; static cilist io___174 = { 0, 3, 0, fmt_101, 0 }; static cilist io___177 = { 0, 3, 0, fmt_101, 0 }; static cilist io___179 = { 0, 3, 0, fmt_101, 0 }; static cilist io___180 = { 0, 3, 0, fmt_101, 0 }; static cilist io___181 = { 0, 3, 0, fmt_101, 0 }; static cilist io___182 = { 0, 3, 0, fmt_101, 0 }; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates starting data for the 2-parameter continuation of limit */ /* points (BVP). */ /* 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___161); 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 *)&lab1, (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)); do_lio(&c__3, &c__1, (char *)&icprs[0], (ftnlen)sizeof(integer)); e_rsle(); nrsp1 = *ntstrs + 1; 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 = k1 + blicn_1.ndm - 1; s_rsfe(&io___174); 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(); /* L1: */ } tm[j] = temp[0]; /* L2: */ } s_rsfe(&io___177); do_fio(&c__1, (char *)&tm[nrsp1], (ftnlen)sizeof(doublereal)); i__1 = blicn_1.ndm; for (k = 1; k <= i__1; ++k) { do_fio(&c__1, (char *)&ups[nrsp1 + k * ups_dim1], (ftnlen)sizeof( doublereal)); } e_rsfe(); nfpar0 = blicn_1.nfpar / 2; s_rsfe(&io___179); i__1 = nfpar0; 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 (Derivative with respect to arclength). */ 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 + blicn_1.ndm + 1; k2 = i * blbcn_1.ndim; s_rsfe(&io___180); 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(); /* L3: */ } /* L4: */ } k1 = blicn_1.ndm + 1; s_rsfe(&io___181); i__1 = blbcn_1.ndim; for (k = k1; k <= i__1; ++k) { do_fio(&c__1, (char *)&ups[nrsp1 + k * ups_dim1], (ftnlen)sizeof( doublereal)); } e_rsfe(); /* Read the parameter values. */ s_rsfe(&io___182); 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(); nfpx = blicn_1.nfpar / 2 - 1; if (nfpx > 0) { i__1 = nfpx; for (i = 1; i <= i__1; ++i) { blbcn_1.par[blbcn_1.icp[nfpar0 + 1 + i - 1] - 1] = blcrl_1.rldot[ i]; /* L5: */ } } 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]; /* L6: */ } *nodir = 1; return 0; } /* stpnbl_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Subroutines for the Continuation of Bifurcations to Tori */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fntr_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int fftr_(); static integer i, j; static doublereal ep, umx; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates the equations for the 2-parameter continuation of */ /* bifurcations to tori. */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ fftr_(ndim, &u[1], uold, &icp[1], &par[1], &f[1], &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); if (*ijac == 0) { return 0; } /* Generate the Jacobian. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L1: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { blwif_1.u1xx[j - 1] = u[j]; blwif_1.u2xx[j - 1] = u[j]; /* L2: */ } blwif_1.u1xx[i - 1] -= ep; blwif_1.u2xx[i - 1] += ep; fftr_(ndim, blwif_1.u1xx, uold, &icp[1], &par[1], blwif_1.f1xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); fftr_(ndim, blwif_1.u2xx, uold, &icp[1], &par[1], blwif_1.f2xx, & blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[j + i * dfdu_dim1] = (blwif_1.f2xx[j - 1] - blwif_1.f1xx[j - 1]) / (ep * 2); /* L3: */ } /* L4: */ } i__1 = blicn_1.nfpar; for (i = 1; i <= i__1; ++i) { par[icp[i]] += ep; fftr_(ndim, &u[1], uold, &icp[1], &par[1], blwif_1.f1xx, &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdp[j + icp[i] * dfdp_dim1] = (blwif_1.f1xx[j - 1] - f[j]) / ep; /* L5: */ } par[icp[i]] -= ep; /* L6: */ } return 0; } /* fntr_ */ /* ---------- ---- */ /* Subroutine */ int fftr_(ndim, u, uold, icp, par, f, ndm, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par, *f; integer *ndm; doublereal *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ static integer ijac; extern /* Subroutine */ int funi_(); static integer i, j; static doublereal c1; static integer ndm2; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dfdp_dim1 = *ndm; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndm; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ ijac = 1; c1 = par[11]; funi_(ndm, &u[1], uold, &icp[1], &par[1], &ijac, &f[1], &dfdu[dfdu_offset] , &dfdp[dfdp_offset]); ndm2 = *ndm << 1; i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndm + i] = blrcn_1.zero; f[ndm2 + i] = blrcn_1.zero; i__2 = *ndm; for (j = 1; j <= i__2; ++j) { f[*ndm + i] += dfdu[i + j * dfdu_dim1] * u[*ndm + j]; f[ndm2 + i] += dfdu[i + j * dfdu_dim1] * u[ndm2 + j]; /* L1: */ } f[*ndm + i] = c1 * f[*ndm + i]; f[ndm2 + i] = c1 * f[ndm2 + i]; f[i] = c1 * f[i]; /* L2: */ } return 0; } /* fftr_ */ /* ---------- ---- */ /* Subroutine */ int bctr_(ndim, par, icp, nbc, u0, u1, f, ijac, dbc) integer *ndim; doublereal *par; integer *icp, *nbc; doublereal *u0, *u1, *f; integer *ijac; doublereal *dbc; { /* System generated locals */ integer dbc_dim1, dbc_offset, i__1, i__2; /* Builtin functions */ double sin(), cos(); /* Local variables */ static integer i, j; static doublereal cs; static integer nn; static doublereal ss, sig; static integer ndm2; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dbc_dim1 = *nbc; dbc_offset = dbc_dim1 + 1; dbc -= dbc_offset; --f; --u1; --u0; --icp; --par; /* Function Body */ ndm2 = blicn_1.ndm << 1; sig = par[12]; ss = sin(sig); /* SGLE SS= SIN(SIG) */ cs = cos(sig); /* SGLE CS= COS(SIG) */ i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { f[i] = u0[i] - u1[i]; f[blicn_1.ndm + i] = u1[blicn_1.ndm + i] - cs * u0[blicn_1.ndm + i] + ss * u0[ndm2 + i]; f[ndm2 + i] = u1[ndm2 + i] - cs * u0[ndm2 + i] - ss * u0[blicn_1.ndm + i]; /* L1: */ } if (blrtn.irot != 0) { for (i= 1; i <= i__1; ++i) { if (blrtn.nrot[i- 1] != 0) { f[i] += blrtn.torper* blrtn.nrot[i- 1]; } } } if (*ijac == 0) { return 0; } nn = (*ndim << 1) + blicn_1.npar; i__1 = *nbc; for (i = 1; i <= i__1; ++i) { i__2 = nn; for (j = 1; j <= i__2; ++j) { dbc[i + j * dbc_dim1] = blrcn_1.zero; /* L2: */ } /* L3: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { dbc[i + i * dbc_dim1] = 1.; dbc[i + (*ndim + i) * dbc_dim1] = -1.; dbc[blicn_1.ndm + i + (blicn_1.ndm + i) * dbc_dim1] = -cs; dbc[blicn_1.ndm + i + (ndm2 + i) * dbc_dim1] = ss; dbc[blicn_1.ndm + i + (*ndim + blicn_1.ndm + i) * dbc_dim1] = 1.; dbc[blicn_1.ndm + i + ((*ndim << 1) + icp[4]) * dbc_dim1] = cs * u0[ ndm2 + i] + ss * u0[blicn_1.ndm + i]; dbc[ndm2 + i + (blicn_1.ndm + i) * dbc_dim1] = -ss; dbc[ndm2 + i + (ndm2 + i) * dbc_dim1] = -cs; dbc[ndm2 + i + (*ndim + ndm2 + i) * dbc_dim1] = 1.; dbc[ndm2 + i + ((*ndim << 1) + icp[4]) * dbc_dim1] = ss * u0[ndm2 + i] - cs * u0[blicn_1.ndm + i]; /* L4: */ } return 0; } /* bctr_ */ /* ---------- ---- */ /* Subroutine */ int ictr_(ndim, par, icp, nint, u, uold, udot, upold, f, ijac, dint) integer *ndim; doublereal *par; integer *icp, *nint; doublereal *u, *uold, *udot, *upold, *f; integer *ijac; doublereal *dint; { /* System generated locals */ integer dint_dim1, dint_offset, i__1, i__2; /* Local variables */ static integer i, j, nn, ndm2; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dint_dim1 = *nint; dint_offset = dint_dim1 + 1; dint -= dint_offset; --f; --upold; --udot; --uold; --u; --icp; --par; /* Function Body */ f[1] = blrcn_1.zero; f[2] = blrcn_1.zero; f[3] = -par[13]; ndm2 = blicn_1.ndm << 1; i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { if (blrtn.nrot[i- 1] == 0) { f[1] += u[i] * upold[i]; } f[2] = f[2] + u[blicn_1.ndm + i] * uold[ndm2 + i] - u[ndm2 + i] * uold[blicn_1.ndm + i]; f[3] = f[3] + u[blicn_1.ndm + i] * u[blicn_1.ndm + i] + u[ndm2 + i] * u[ndm2 + i]; /* L1: */ } if (*ijac == 0) { return 0; } nn = *ndim + blicn_1.npar; i__1 = *nint; for (i = 1; i <= i__1; ++i) { i__2 = nn; for (j = 1; j <= i__2; ++j) { dint[i + j * dint_dim1] = blrcn_1.zero; /* L2: */ } /* L3: */ } i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { if (blrtn.nrot[i- 1] == 0) { dint[i * dint_dim1 + 1] = upold[i]; }else {dint[i * dint_dim1 + 1] = 0.;} dint[(blicn_1.ndm + i) * dint_dim1 + 2] = uold[ndm2 + i]; dint[(ndm2 + i) * dint_dim1 + 2] = -uold[blicn_1.ndm + i]; dint[(blicn_1.ndm + i) * dint_dim1 + 3] = blrcn_1.two * u[blicn_1.ndm + i]; dint[(ndm2 + i) * dint_dim1 + 3] = blrcn_1.two * u[ndm2 + i]; /* L4: */ } dint[(*ndim + 13) * dint_dim1 + 3] = -blrcn_1.one; return 0; } /* ictr_ */ /* ---------- ------ */ /* Subroutine */ int stpntr_(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(); double sin(); /* Local variables */ static doublereal temp[7]; static integer ntpl1, nrsp1, ntot1, i, j, k; static logical found; static integer icprs[20], k1, k2, k3; extern /* Subroutine */ int findl3_(); static integer nfpar1, nskip1, k2p1, lab1, nar1, itp1, isw1; /* Fortran I/O blocks */ static cilist io___207 = { 0, 3, 0, 0, 0 }; static cilist io___220 = { 0, 3, 0, fmt_101, 0 }; static cilist io___225 = { 0, 3, 0, fmt_101, 0 }; static cilist io___226 = { 0, 3, 0, fmt_101, 0 }; static cilist io___227 = { 0, 3, 0, fmt_101, 0 }; static cilist io___228 = { 0, 3, 0, fmt_101, 0 }; static cilist io___229 = { 0, 3, 0, fmt_101, 0 }; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates starting data for the 2-parameter continuation of torus */ /* bifurcations. */ /* 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___207); 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 *)&lab1, (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)); do_lio(&c__3, &c__1, (char *)&icprs[0], (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&icprs[1], (ftnlen)sizeof(integer)); e_rsle(); nrsp1 = *ntstrs + 1; 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 = k1 + blicn_1.ndm - 1; s_rsfe(&io___220); 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(); k2p1 = k2 + 1; k3 = k2 + blicn_1.ndm; i__3 = k3; for (k = k2p1; k <= i__3; ++k) { ups[j + k * ups_dim1] = blrcn_1.rsmall * sin(temp[i - 1]); /* SGLE UPS(J,K)=RSMALL* SIN(TEMP(I)) */ ups[j + (k + blicn_1.ndm) * ups_dim1] = blrcn_1.zero; /* L1: */ } /* L2: */ } tm[j] = temp[0]; /* L3: */ } s_rsfe(&io___225); do_fio(&c__1, (char *)&tm[nrsp1], (ftnlen)sizeof(doublereal)); i__1 = blicn_1.ndm; for (k = 1; k <= i__1; ++k) { do_fio(&c__1, (char *)&ups[nrsp1 + k * ups_dim1], (ftnlen)sizeof( doublereal)); } e_rsfe(); i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { ups[nrsp1 + (blicn_1.ndm + i) * ups_dim1] = blrcn_1.zero; ups[nrsp1 + ((blicn_1.ndm << 1) + i) * ups_dim1] = blrcn_1.zero; /* L4: */ } s_rsfe(&io___226); do_fio(&c__1, (char *)&blcrl_1.rldot[0], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&blcrl_1.rldot[2], (ftnlen)sizeof(doublereal)); e_rsfe(); blcrl_1.rldot[1] = blrcn_1.zero; blcrl_1.rldot[3] = blrcn_1.zero; /* Read the direction vector and initialize the starting direction. */ 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 = k1 + blicn_1.ndm - 1; s_rsfe(&io___227); 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(); k2p1 = k2 + 1; k3 = k2 + blicn_1.ndm; i__3 = k3; for (k = k2p1; k <= i__3; ++k) { udotps[j + k * udotps_dim1] = blrcn_1.zero; udotps[j + (k + blicn_1.ndm) * udotps_dim1] = blrcn_1.zero; /* L5: */ } /* L6: */ } /* L7: */ } s_rsfe(&io___228); i__1 = blicn_1.ndm; for (k = 1; k <= i__1; ++k) { do_fio(&c__1, (char *)&udotps[nrsp1 + k * udotps_dim1], (ftnlen) sizeof(doublereal)); } e_rsfe(); i__1 = blicn_1.ndm; for (i = 1; i <= i__1; ++i) { udotps[nrsp1 + (blicn_1.ndm + i) * udotps_dim1] = blrcn_1.zero; udotps[nrsp1 + ((blicn_1.ndm << 1) + i) * udotps_dim1] = blrcn_1.zero; /* L8: */ } /* Read the parameter values. */ s_rsfe(&io___229); 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(); blbcn_1.par[12] = blrcn_1.zero; 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]; /* L9: */ } *nodir = 0; return 0; } /* stpntr_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Periodic Solutions and Fixed Period Orbits */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fnps_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ extern /* Subroutine */ int funi_(); static integer i, j; static doublereal rho; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates the equations for the continuation of periodic orbits. */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ funi_(ndim, &u[1], uold, &icp[1], &par[1], ijac, &f[1], &dfdu[dfdu_offset] , &dfdp[dfdp_offset]); rho = par[icp[2]]; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { dfdp[i + icp[2] * dfdp_dim1] = f[i]; f[i] = rho * dfdp[i + icp[2] * dfdp_dim1]; /* L1: */ } if (*ijac == 0) { return 0; } /* Generate the Jacobian. */ i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[i + j * dfdu_dim1] = rho * dfdu[i + j * dfdu_dim1]; /* L2: */ } dfdp[i + icp[1] * dfdp_dim1] = rho * dfdp[i + icp[1] * dfdp_dim1]; /* L3: */ } return 0; } /* fnps_ */ /* ---------- ---- */ /* Subroutine */ int fnfp_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ extern /* Subroutine */ int funi_(); static integer i, j; static doublereal rho; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates the equations for the 2-parameter continuation of */ /* periodic orbits of fixed period. */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ rho = par[icp[3]]; funi_(ndim, &u[1], uold, &icp[1], &par[1], ijac, &f[1], &dfdu[dfdu_offset] , &dfdp[dfdp_offset]); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { f[i] = rho * f[i]; /* L1: */ } if (*ijac == 0) { return 0; } /* Generate the Jacobian. */ i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[i + j * dfdu_dim1] = rho * dfdu[i + j * dfdu_dim1]; /* L2: */ } for (j = 1; j <= 2; ++j) { dfdp[i + icp[j] * dfdp_dim1] = rho * dfdp[i + icp[j] * dfdp_dim1]; /* L3: */ } /* L4: */ } return 0; } /* fnfp_ */ /* ---------- ---- */ /* Subroutine */ int bcps_(ndim, par, icp, nbc, u0, u1, f, ijac, dbc) integer *ndim; doublereal *par; integer *icp, *nbc; doublereal *u0, *u1, *f; integer *ijac; doublereal *dbc; { /* System generated locals */ integer dbc_dim1, dbc_offset, i__1, i__2; /* Local variables */ static integer i, j, nn; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dbc_dim1 = *nbc; dbc_offset = dbc_dim1 + 1; dbc -= dbc_offset; --f; --u1; --u0; --icp; --par; /* Function Body */ i__1 = *ndim; for (i = 1; i <= i__1; ++i) { f[i] = u0[i] - u1[i]; /* L1: */ } if (blrtn.irot != 0) { i__1 = *ndim; for (i= 1; i <= i__1; ++i) { if (blrtn.nrot[i- 1] != 0) { f[i] += blrtn.torper* blrtn.nrot[i- 1]; } } } if (*ijac == 0) { return 0; } nn = (*ndim << 1) + blicn_1.npar; i__1 = *nbc; for (i = 1; i <= i__1; ++i) { i__2 = nn; for (j = 1; j <= i__2; ++j) { dbc[i + j * dbc_dim1] = blrcn_1.zero; /* L2: */ } /* L3: */ } i__1 = *ndim; for (i = 1; i <= i__1; ++i) { dbc[i + i * dbc_dim1] = blrcn_1.one; dbc[i + (*ndim + i) * dbc_dim1] = -blrcn_1.one; /* L4: */ } return 0; } /* bcps_ */ /* ---------- ---- */ /* Subroutine */ int icps_(ndim, par, icp, nint, u, uold, udot, upold, f, ijac, dint) integer *ndim; doublereal *par; integer *icp, *nint; doublereal *u, *uold, *udot, *upold, *f; integer *ijac; doublereal *dint; { /* System generated locals */ integer dint_dim1, dint_offset, i__1, i__2; /* Local variables */ static integer i, j, nn; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dint_dim1 = *nint; dint_offset = dint_dim1 + 1; dint -= dint_offset; --f; --upold; --udot; --uold; --u; --icp; --par; /* Function Body */ f[1] = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if(blrtn.nrot[i-1]==0){ f[1] += u[i] * upold[i]; } /* L1: */ } if (*ijac == 0) { return 0; } nn = *ndim + blicn_1.npar; i__1 = *nint; for (i = 1; i <= i__1; ++i) { i__2 = nn; for (j = 1; j <= i__2; ++j) { dint[i + j * dint_dim1] = blrcn_1.zero; /* L2: */ } /* L3: */ } i__1 = *ndim; for (i = 1; i <= i__1; ++i) { dint[i * dint_dim1 + 1] = upold[i]; /* L4: */ } return 0; } /* icps_ */ /* ---------- ----- */ /* Subroutine */ int pdble_(ndim, ntst, ncol, m1u, ups, udotps, tm, rho) integer *ndim, *ntst, *ncol, *m1u; doublereal *ups, *udotps, *tm, *rho; { /* System generated locals */ integer ups_dim1, ups_offset, udotps_dim1, udotps_offset, i__1, i__2, i__3; /* Local variables */ static integer i, j, i1, i2; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Preprocesses restart data for switching branches at a period doubling */ /* bifurcation. */ /* Parameter adjustments */ --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 */ *rho = blrcn_1.two * *rho; if (blrtn.irot != 0) { blrtn.torper *= 2.; } i__1 = *ntst; for (i = 1; i <= i__1; ++i) { tm[i] = blrcn_1.half * tm[i]; tm[*ntst + i] = blrcn_1.half + tm[i]; /* L1: */ } tm[(*ntst << 1) + 1] = blrcn_1.one; i__1 = *ntst + 1; for (j = 1; j <= i__1; ++j) { i__2 = *ndim; for (i1 = 1; i1 <= i__2; ++i1) { i__3 = *ncol; for (i2 = 1; i2 <= i__3; ++i2) { i = (i2 - 1) * *ndim + i1; ups[*ntst + j + i * ups_dim1] = ups[*ntst + 1 + i1 * ups_dim1] + ups[j + i * ups_dim1] - ups[i1 * ups_dim1 + 1]; udotps[*ntst + j + i * udotps_dim1] = udotps[*ntst + 1 + i1 * udotps_dim1] + udotps[j + i * udotps_dim1] - udotps[ i1 * udotps_dim1 + 1]; /* L2: */ } } /* L3: */ } *ntst <<= 1; return 0; } /* pdble_ */ /* ---------- ------ */ /* Subroutine */ int stpnps_(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; /* Builtin functions */ double sin(), cos(); /* Local variables */ static integer ijac; extern /* Subroutine */ int funi_(), nlvc_(); static doublereal uold, c; static integer i, j, k; static doublereal s, t, rimhb; static logical found; static integer k1; extern /* Subroutine */ int nrmlz_(), readl3_(), findl3_(); static integer nfpar1; static doublereal dt; extern doublereal pi_(); extern /* Subroutine */ int scalebb_(), msh_(); static doublereal rho; static integer itp; static doublereal tpi; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates starting data for the continuation of a branch of periodic */ /* solutions from a Hopf bifurcation point. */ /* 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, &itp, &nfpar1, &found); readl3_(&blbcn_1.ips, ibr, &u[1], blbcn_1.par); 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]; /* L1: */ } rho = blbcn_1.par[10]; tpi = pi_(&blrcn_1.two); rimhb = tpi / rho; *ntstrs = blcde_1.ntst; *ncolrs = blcde_1.ncol; i__1 = *ndim2; for (i = 1; i <= i__1; ++i) { i__2 = *ndim2; for (j = 1; j <= i__2; ++j) { smat[i + j * smat_dim1] = blrcn_1.zero; /* L2: */ } /* L3: */ } i__1 = blbcn_1.ndim; for (i = 1; i <= i__1; ++i) { smat[i + i * smat_dim1] = -rimhb; smat[blbcn_1.ndim + i + (blbcn_1.ndim + i) * smat_dim1] = rimhb; /* L4: */ } ijac = 1; funi_(&blbcn_1.ndim, &u[1], &uold, blbcn_1.icp, blbcn_1.par, &ijac, &f[1], &dfdu[dfdu_offset], &dfdp[dfdp_offset]); i__1 = blbcn_1.ndim; for (i = 1; i <= i__1; ++i) { i__2 = blbcn_1.ndim; for (j = 1; j <= i__2; ++j) { smat[i + (blbcn_1.ndim + j) * smat_dim1] = dfdu[i + j * dfdu_dim1] ; smat[blbcn_1.ndim + i + j * smat_dim1] = dfdu[i + j * dfdu_dim1]; /* L5: */ } /* L6: */ } nlvc_(ndim2, ndim2, &c__2, &smat[smat_offset], &rnllv[1], &ir[1], &ic[1]); nrmlz_(ndim2, &rnllv[1]); /* Generate the (initially uniform) mesh. */ msh_(&tm[1]); dt = blrcn_1.one / blcde_1.ntst; i__1 = blcde_1.ntst + 1; for (j = 1; j <= i__1; ++j) { t = tm[j]; s = sin(tpi * t); /* SGLE S=SIN(TPI*T) */ c = cos(tpi * t); /* SGLE C=COS(TPI*T) */ i__2 = blbcn_1.ndim; for (k = 1; k <= i__2; ++k) { udotps[j + k * udotps_dim1] = s * rnllv[k] + c * rnllv[ blbcn_1.ndim + k]; upoldp[j + k * upoldp_dim1] = c * rnllv[k] - s * rnllv[ blbcn_1.ndim + k]; ups[j + k * ups_dim1] = u[k]; /* L7: */ } /* L8: */ } i__1 = blcde_1.ncol - 1; for (i = 1; i <= i__1; ++i) { i__2 = blcde_1.ntst; for (j = 1; j <= i__2; ++j) { t = tm[j] + i * (tm[j + 1] - tm[j]) / blcde_1.ncol; s = sin(tpi * t); /* SGLE S=SIN(TPI*T) */ c = cos(tpi * t); /* SGLE C=COS(TPI*T) */ i__3 = blbcn_1.ndim; for (k = 1; k <= i__3; ++k) { k1 = i * blbcn_1.ndim + k; udotps[j + k1 * udotps_dim1] = s * rnllv[k] + c * rnllv[ blbcn_1.ndim + k]; upoldp[j + k1 * upoldp_dim1] = c * rnllv[k] - s * rnllv[ blbcn_1.ndim + k]; ups[j + k1 * ups_dim1] = u[k]; /* L9: */ } /* L10: */ } /* L11: */ } blcrl_1.rldot[0] = blrcn_1.zero; blcrl_1.rldot[1] = blrcn_1.zero; i__1 = blcde_1.ntst; for (i = 1; i <= i__1; ++i) { dtm[i] = dt; /* L12: */ } scalebb_(m1u, &udotps[udotps_offset], blcrl_1.rldot, &dtm[1]); *nodir = -1; return 0; } /* stpnps_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Travelling Wave Solutions to Diffusive Systems */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fnws_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset; /* Local variables */ extern /* Subroutine */ int ffws_(); static integer ndm2; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Sets up equations for the continuation of spatially homogeneous */ /* solutions to parabolic systems, for the purpose of finding */ /* bifurcations to travelling wave solutions. */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ ndm2 = blicn_1.ndm / 2; ffws_(ndim, &u[1], uold, &icp[1], &blicn_1.nfpar, &par[1], ijac, &f[1], & dfdu[dfdu_offset], &dfdp[dfdp_offset], &ndm2, blwif_1.dfuxx, blwif_1.dfpxx); return 0; } /* fnws_ */ /* ---------- ---- */ /* Subroutine */ int ffws_(ndim, u, uold, icp, nfpar, par, ijac, f, dfdu, dfdp, ndm, dfuxx, dfpxx) integer *ndim; doublereal *u, *uold; integer *icp, *nfpar; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; integer *ndm; doublereal *dfuxx, *dfpxx; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, dfuxx_dim1, dfuxx_offset, dfpxx_dim1, dfpxx_offset, i__1, i__2; /* Local variables */ extern /* Subroutine */ int funi_(); static doublereal c; static integer i, j; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dfpxx_dim1 = *ndm; dfpxx_offset = dfpxx_dim1 + 1; dfpxx -= dfpxx_offset; dfuxx_dim1 = *ndm; dfuxx_offset = dfuxx_dim1 + 1; dfuxx -= dfuxx_offset; dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ c = par[10]; funi_(ndm, &u[1], uold, &icp[1], &par[1], ijac, &f[1], &dfuxx[ dfuxx_offset], &dfpxx[dfpxx_offset]); i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[*ndm + i] = -(c * u[*ndm + i] + f[i]) / par[i + 14]; f[i] = u[*ndm + i]; /* L1: */ } if (*ijac == 0) { return 0; } i__1 = *ndm; for (i = 1; i <= i__1; ++i) { i__2 = *ndm; for (j = 1; j <= i__2; ++j) { dfdu[i + j * dfdu_dim1] = blrcn_1.zero; dfdu[i + (j + *ndm) * dfdu_dim1] = blrcn_1.zero; dfdu[i + *ndm + j * dfdu_dim1] = -dfuxx[i + j * dfuxx_dim1] / par[ i + 14]; dfdu[i + *ndm + (j + *ndm) * dfdu_dim1] = blrcn_1.zero; /* L2: */ } dfdu[i + (i + *ndm) * dfdu_dim1] = blrcn_1.one; dfdu[i + *ndm + (i + *ndm) * dfdu_dim1] = -c / par[i + 14]; if (icp[1] < 10) { dfdp[i + icp[1] * dfdp_dim1] = blrcn_1.zero; dfdp[i + *ndm + icp[1] * dfdp_dim1] = -dfpxx[i + icp[1] * dfpxx_dim1] / par[i + 14]; } if (*nfpar > 1 && icp[2] < 10) { dfdp[i + icp[2] * dfdp_dim1] = blrcn_1.zero; dfdp[i + *ndm + icp[2] * dfdp_dim1] = -dfpxx[i + icp[2] * dfpxx_dim1] / par[i + 14]; } /* L3: */ } /* Derivative with respect to the wave speed. */ i__1 = *ndm; for (i = 1; i <= i__1; ++i) { dfdp[i + dfdp_dim1 * 10] = blrcn_1.zero; dfdp[i + *ndm + dfdp_dim1 * 10] = -u[*ndm + i] / par[i + 14]; /* L4: */ } /* Derivatives with respect to the diffusion coefficients. */ i__1 = *ndm; for (j = 1; j <= i__1; ++j) { i__2 = *ndm; for (i = 1; i <= i__2; ++i) { dfdp[i + (j + 14) * dfdp_dim1] = blrcn_1.zero; dfdp[i + *ndm + (j + 14) * dfdp_dim1] = blrcn_1.zero; /* L5: */ } dfdp[j + *ndm + (j + 14) * dfdp_dim1] = -f[j + *ndm] / par[j + 14]; /* L6: */ } return 0; } /* ffws_ */ /* ---------- ---- */ /* Subroutine */ int fnwp_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ extern /* Subroutine */ int fnws_(); static integer i, j; static doublereal rho; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Equations for the continuation of traveling waves. */ /* Generate the function and Jacobian. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ fnws_(ndim, &u[1], uold, &icp[1], &par[1], ijac, &f[1], &dfdu[dfdu_offset] , &dfdp[dfdp_offset]); rho = par[icp[2]]; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { dfdp[i + icp[2] * dfdp_dim1] = f[i]; f[i] = rho * f[i]; /* L1: */ } if (*ijac == 0) { return 0; } i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[i + j * dfdu_dim1] = rho * dfdu[i + j * dfdu_dim1]; /* L2: */ } } i__2 = *ndim; for (i = 1; i <= i__2; ++i) { dfdp[i + icp[1] * dfdp_dim1] = rho * dfdp[i + icp[1] * dfdp_dim1]; /* L3: */ } return 0; } /* fnwp_ */ /* ---------- ---- */ /* Subroutine */ int fnwf_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; /* Local variables */ extern /* Subroutine */ int fnws_(); static integer i, j; static doublereal rho; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates the equations for the two parameter continuation */ /* of wavetrains of fixed wave length. */ /* Generate the function and Jacobian. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --u; /* Function Body */ fnws_(ndim, &u[1], uold, &icp[1], &par[1], ijac, &f[1], &dfdu[dfdu_offset] , &dfdp[dfdp_offset]); rho = par[icp[3]]; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { f[i] = rho * f[i]; /* L1: */ } if (*ijac == 0) { return 0; } i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[i + j * dfdu_dim1] = rho * dfdu[i + j * dfdu_dim1]; /* L2: */ } } i__2 = *ndim; for (i = 1; i <= i__2; ++i) { for (j = 1; j <= 2; ++j) { dfdp[i + icp[j] * dfdp_dim1] = rho * dfdp[i + icp[j] * dfdp_dim1]; /* L3: */ } } return 0; } /* fnwf_ */ /* ---------- ------ */ /* Subroutine */ int stpnwp_(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; /* Builtin functions */ double sin(), cos(); /* Local variables */ static integer ijac; extern /* Subroutine */ int nlvc_(); static doublereal uold; extern /* Subroutine */ int fnws_(); static doublereal c; static integer i, j, k; static doublereal s, t, rimhb; static logical found; static integer k1; extern /* Subroutine */ int nrmlz_(), readl3_(), findl3_(); static integer nfpar1; static doublereal dt; extern doublereal pi_(); extern /* Subroutine */ int scalebb_(), msh_(); static doublereal rho; static integer itp; static doublereal tpi; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates starting data for the continuation of a branch of periodic */ /* solutions starting from a Hopf bifurcation point (Waves). */ /* 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, &itp, &nfpar1, &found); readl3_(&blbcn_1.ips, ibr, &u[1], blbcn_1.par); 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]; /* L1: */ } rho = blbcn_1.par[10]; tpi = pi_(&blrcn_1.two); rimhb = tpi / rho; *ntstrs = blcde_1.ntst; *ncolrs = blcde_1.ncol; i__1 = *ndim2; for (i = 1; i <= i__1; ++i) { i__2 = *ndim2; for (j = 1; j <= i__2; ++j) { smat[i + j * smat_dim1] = blrcn_1.zero; /* L2: */ } /* L3: */ } i__1 = blbcn_1.ndim; for (i = 1; i <= i__1; ++i) { smat[i + i * smat_dim1] = -rimhb; smat[blbcn_1.ndim + i + (blbcn_1.ndim + i) * smat_dim1] = rimhb; /* L4: */ } ijac = 1; fnws_(&blbcn_1.ndim, &u[1], &uold, blbcn_1.icp, blbcn_1.par, &ijac, &f[1], &dfdu[dfdu_offset], &dfdp[dfdp_offset]); i__1 = blbcn_1.ndim; for (i = 1; i <= i__1; ++i) { i__2 = blbcn_1.ndim; for (j = 1; j <= i__2; ++j) { smat[i + (blbcn_1.ndim + j) * smat_dim1] = dfdu[i + j * dfdu_dim1] ; smat[blbcn_1.ndim + i + j * smat_dim1] = dfdu[i + j * dfdu_dim1]; /* L5: */ } /* L6: */ } nlvc_(ndim2, ndim2, &c__2, &smat[smat_offset], &rnllv[1], &ir[1], &ic[1]); nrmlz_(ndim2, &rnllv[1]); /* Generate the (initially uniform) mesh. */ msh_(&tm[1]); dt = blrcn_1.one / blcde_1.ntst; i__1 = blcde_1.ntst + 1; for (j = 1; j <= i__1; ++j) { t = tm[j]; s = sin(tpi * t); /* SGLE S=SIN(TPI*T) */ c = cos(tpi * t); /* SGLE C=COS(TPI*T) */ i__2 = blbcn_1.ndim; for (k = 1; k <= i__2; ++k) { udotps[j + k * udotps_dim1] = s * rnllv[k] + c * rnllv[ blbcn_1.ndim + k]; upoldp[j + k * upoldp_dim1] = c * rnllv[k] - s * rnllv[ blbcn_1.ndim + k]; ups[j + k * ups_dim1] = u[k]; /* L7: */ } /* L8: */ } i__1 = blcde_1.ncol - 1; for (i = 1; i <= i__1; ++i) { i__2 = blcde_1.ntst; for (j = 1; j <= i__2; ++j) { t = tm[j] + i * (tm[j + 1] - tm[j]) / blcde_1.ncol; s = sin(tpi * t); /* SGLE S=SIN(TPI*T) */ c = cos(tpi * t); /* SGLE C=COS(TPI*T) */ i__3 = blbcn_1.ndim; for (k = 1; k <= i__3; ++k) { k1 = i * blbcn_1.ndim + k; udotps[j + k1 * udotps_dim1] = s * rnllv[k] + c * rnllv[ blbcn_1.ndim + k]; upoldp[j + k1 * upoldp_dim1] = c * rnllv[k] - s * rnllv[ blbcn_1.ndim + k]; ups[j + k1 * ups_dim1] = u[k]; /* L9: */ } /* L10: */ } /* L11: */ } blcrl_1.rldot[0] = blrcn_1.zero; blcrl_1.rldot[1] = blrcn_1.zero; i__1 = blcde_1.ntst; for (i = 1; i <= i__1; ++i) { dtm[i] = dt; /* L12: */ } scalebb_(m1u, &udotps[udotps_offset], blcrl_1.rldot, &dtm[1]); *nodir = -1; return 0; } /* stpnwp_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Diffusive Systems : Time Evolution */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int fnpe_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset; /* Local variables */ extern /* Subroutine */ int ffpe_(); /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Generates the equations for taking one time step (Implicit Euler). */ /* Generate the function and Jacobian. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --uold; --u; /* Function Body */ ffpe_(ndim, &u[1], &uold[1], &icp[1], &par[1], ijac, &f[1], &dfdu[ dfdu_offset], &dfdp[dfdp_offset], &blicn_1.ndm, blwif_1.dfuxx, blwif_1.dfpxx); return 0; } /* fnpe_ */ /* ---------- ---- */ /* Subroutine */ int ffpe_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp, ndm, dfuxx, dfpxx) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; integer *ndm; doublereal *dfuxx, *dfpxx; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, dfuxx_dim1, dfuxx_offset, dfpxx_dim1, dfpxx_offset, i__1, i__2; doublereal d__1; /* Local variables */ extern /* Subroutine */ int funi_(); static integer i, j; static doublereal t, dt, rho; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Parameter adjustments */ dfpxx_dim1 = *ndm; dfpxx_offset = dfpxx_dim1 + 1; dfpxx -= dfpxx_offset; dfuxx_dim1 = *ndm; dfuxx_offset = dfuxx_dim1 + 1; dfuxx -= dfuxx_offset; dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --uold; --u; /* Function Body */ rho = par[11]; t = par[icp[1]]; dt = t - blcrl_1.rlold[0]; if (abs(dt) < bldls_1.dsmin) { dt = bldls_1.ds; } /* SGLE IF( ABS(DT).LT.DSMIN)DT=DS */ funi_(ndm, &u[1], &uold[1], &icp[1], &par[1], ijac, &f[*ndm + 1], &dfuxx[ dfuxx_offset], &dfpxx[dfpxx_offset]); i__1 = *ndm; for (i = 1; i <= i__1; ++i) { f[i] = rho * u[*ndm + i]; f[*ndm + i] = rho * ((u[i] - uold[i]) / dt - f[*ndm + i]) / par[i + 14]; /* L1: */ } if (*ijac == 0) { return 0; } i__1 = *ndm; for (i = 1; i <= i__1; ++i) { i__2 = *ndm; for (j = 1; j <= i__2; ++j) { dfdu[i + j * dfdu_dim1] = blrcn_1.zero; dfdu[i + (j + *ndm) * dfdu_dim1] = blrcn_1.zero; dfdu[i + *ndm + j * dfdu_dim1] = -rho * dfuxx[i + j * dfuxx_dim1] / par[i + 14]; dfdu[i + *ndm + (j + *ndm) * dfdu_dim1] = blrcn_1.zero; /* L2: */ } dfdu[i + (i + *ndm) * dfdu_dim1] = rho; dfdu[i + *ndm + i * dfdu_dim1] += rho / (dt * par[i + 14]); dfdp[i + icp[1] * dfdp_dim1] = blrcn_1.zero; /* Computing 2nd power */ d__1 = dt; dfdp[i + *ndm + icp[1] * dfdp_dim1] = -rho * (u[i] - uold[i]) / (d__1 * d__1 * par[i + 14]); /* L3: */ } return 0; } /* ffpe_ */ /* Subroutine */ int icpe_(ndim, par, icp, nint, u, uold, udot, upold, f, ijac, dint) integer *ndim; real *par; integer *icp, *nint; real *u, *uold, *udot, *upold, *f; integer *ijac; real *dint; { /* ---------- ---- */ /* Dummy integral condition subroutine for parabolic systems. */ return 0; } /* icpe_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Routines for Interface with User Supplied Routines */ /* (To generate Jacobian by differencing, if not supplied analytically) */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* ---------- ---- */ /* Subroutine */ int funi_(ndim, u, uold, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u, *uold; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer dfdu_dim1, dfdu_offset, dfdp_dim1, dfdp_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int func_(); static integer i, j; static doublereal ep, umx; /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Interface subroutine to user supplied FUNC. */ /* Generate the function. */ /* Parameter adjustments */ dfdp_dim1 = *ndim; dfdp_offset = dfdp_dim1 + 1; dfdp -= dfdp_offset; dfdu_dim1 = *ndim; dfdu_offset = dfdu_dim1 + 1; dfdu -= dfdu_offset; --f; --par; --icp; --uold; --u; /* Function Body */ func_(ndim, &u[1], &icp[1], &par[1], ijac, &f[1], &dfdu[dfdu_offset], & dfdp[dfdp_offset]); if (blmax_1.jac == 1 || *ijac == 0) { return 0; } /* Generate the Jacobian by differencing. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L1: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { bldif_1.u1zz[j - 1] = u[j]; bldif_1.u2zz[j - 1] = u[j]; /* L2: */ } bldif_1.u1zz[i - 1] -= ep; bldif_1.u2zz[i - 1] += ep; func_(ndim, bldif_1.u1zz, &icp[1], &par[1], &c__0, bldif_1.f1zz, & dfdu[dfdu_offset], &dfdp[dfdp_offset]); func_(ndim, bldif_1.u2zz, &icp[1], &par[1], &c__0, bldif_1.f2zz, & dfdu[dfdu_offset], &dfdp[dfdp_offset]); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdu[j + i * dfdu_dim1] = (bldif_1.f2zz[j - 1] - bldif_1.f1zz[j - 1]) / (ep * 2); /* L3: */ } /* L4: */ } i__1 = blicn_1.nfpar; for (i = 1; i <= i__1; ++i) { ep = blrcn_1.hmach * (blrcn_1.one + (d__1 = par[icp[i]], abs(d__1))); /* SGLE EP=HMACH*( ONE + ABS(PAR(ICP(I))) ) */ par[icp[i]] += ep; func_(ndim, &u[1], &icp[1], &par[1], &c__0, bldif_1.f1zz, &dfdu[ dfdu_offset], &dfdp[dfdp_offset]); i__2 = *ndim; for (j = 1; j <= i__2; ++j) { dfdp[j + icp[i] * dfdp_dim1] = (bldif_1.f1zz[j - 1] - f[j]) / ep; /* L5: */ } par[icp[i]] -= ep; /* L6: */ } return 0; } /* funi_ */ /* Subroutine */ int bcni_(ndim, par, icp, nbc, u0, u1, f, ijac, dbc) integer *ndim; doublereal *par; integer *icp, *nbc; doublereal *u0, *u1, *f; integer *ijac; doublereal *dbc; { /* System generated locals */ integer dbc_dim1, dbc_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int bcnd_(); static integer i, j; static doublereal ep, umx; /* ---------- ---- */ /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Interface subroutine to the user supplied BCND. */ /* Generate the function. */ /* Parameter adjustments */ dbc_dim1 = *nbc; dbc_offset = dbc_dim1 + 1; dbc -= dbc_offset; --f; --u1; --u0; --icp; --par; /* Function Body */ bcnd_(ndim, &par[1], &icp[1], nbc, &u0[1], &u1[1], &f[1], ijac, &dbc[ dbc_offset]); if (blmax_1.jac == 1 || *ijac == 0) { return 0; } /* Generate the Jacobian by differencing. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u0[i], abs(d__1)) > umx) { umx = (d__2 = u0[i], abs(d__2)); } /* SGLE IF( ABS(U0(I)).GT.UMX)UMX= ABS(U0(I)) */ /* L1: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { bldif_1.u1zz[j - 1] = u0[j]; bldif_1.u2zz[j - 1] = u0[j]; /* L2: */ } bldif_1.u1zz[i - 1] -= ep; bldif_1.u2zz[i - 1] += ep; bcnd_(ndim, &par[1], &icp[1], nbc, bldif_1.u1zz, &u1[1], bldif_1.f1zz, &c__0, &dbc[dbc_offset]); bcnd_(ndim, &par[1], &icp[1], nbc, bldif_1.u2zz, &u1[1], bldif_1.f2zz, &c__0, &dbc[dbc_offset]); i__2 = *nbc; for (j = 1; j <= i__2; ++j) { dbc[j + i * dbc_dim1] = (bldif_1.f2zz[j - 1] - bldif_1.f1zz[j - 1] ) / (ep * 2); /* L3: */ } /* L4: */ } umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u1[i], abs(d__1)) > umx) { umx = (d__2 = u1[i], abs(d__2)); } /* SGLE IF( ABS(U1(I)).GT.UMX)UMX= ABS(U1(I)) */ /* L5: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { bldif_1.u1zz[j - 1] = u1[j]; bldif_1.u2zz[j - 1] = u1[j]; /* L6: */ } bldif_1.u1zz[i - 1] -= ep; bldif_1.u2zz[i - 1] += ep; bcnd_(ndim, &par[1], &icp[1], nbc, &u0[1], bldif_1.u1zz, bldif_1.f1zz, &c__0, &dbc[dbc_offset]); bcnd_(ndim, &par[1], &icp[1], nbc, &u0[1], bldif_1.u2zz, bldif_1.f2zz, &c__0, &dbc[dbc_offset]); i__2 = *nbc; for (j = 1; j <= i__2; ++j) { dbc[j + (*ndim + i) * dbc_dim1] = (bldif_1.f2zz[j - 1] - bldif_1.f1zz[j - 1]) / (ep * 2); /* L7: */ } /* L8: */ } i__1 = blicn_1.nfpar; for (i = 1; i <= i__1; ++i) { ep = blrcn_1.hmach * (blrcn_1.one + (d__1 = par[icp[i]], abs(d__1))); /* SGLE EP=HMACH*( ONE + ABS(PAR(ICP(I))) ) */ par[icp[i]] += ep; bcnd_(ndim, &par[1], &icp[1], nbc, &u0[1], &u1[1], bldif_1.f1zz, & c__0, &dbc[dbc_offset]); i__2 = *nbc; for (j = 1; j <= i__2; ++j) { dbc[j + ((*ndim << 1) + icp[i]) * dbc_dim1] = (bldif_1.f1zz[j - 1] - f[j]) / ep; /* L9: */ } par[icp[i]] -= ep; /* L10: */ } return 0; } /* bcni_ */ /* Subroutine */ int icni_(ndim, par, icp, nint, u, uold, udot, upold, f, ijac, dint) integer *ndim; doublereal *par; integer *icp, *nint; doublereal *u, *uold, *udot, *upold, *f; integer *ijac; doublereal *dint; { /* System generated locals */ integer dint_dim1, dint_offset, i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int icnd_(); static integer i, j; static doublereal ep, umx; /* ---------- ---- */ /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Interface subroutine to user supplied ICND. */ /* Generate the integrand. */ /* Parameter adjustments */ dint_dim1 = *nint; dint_offset = dint_dim1 + 1; dint -= dint_offset; --f; --upold; --udot; --uold; --u; --icp; --par; /* Function Body */ icnd_(ndim, &par[1], &icp[1], nint, &u[1], &uold[1], &udot[1], &upold[1], &f[1], ijac, &dint[dint_offset]); if (blmax_1.jac == 1 || *ijac == 0) { return 0; } /* Generate the Jacobian by differencing. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L1: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { bldif_1.u1zz[j - 1] = u[j]; bldif_1.u2zz[j - 1] = u[j]; /* L2: */ } bldif_1.u1zz[i - 1] -= ep; bldif_1.u2zz[i - 1] += ep; icnd_(ndim, &par[1], &icp[1], nint, bldif_1.u1zz, &uold[1], &udot[1], &upold[1], bldif_1.f1zz, &c__0, &dint[dint_offset]); icnd_(ndim, &par[1], &icp[1], nint, bldif_1.u2zz, &uold[1], &udot[1], &upold[1], bldif_1.f2zz, &c__0, &dint[dint_offset]); i__2 = *nint; for (j = 1; j <= i__2; ++j) { dint[j + i * dint_dim1] = (bldif_1.f2zz[j - 1] - bldif_1.f1zz[j - 1]) / (ep * 2); /* L3: */ } /* L4: */ } i__1 = blicn_1.nfpar; for (i = 1; i <= i__1; ++i) { ep = blrcn_1.hmach * (blrcn_1.one + (d__1 = par[icp[i]], abs(d__1))); /* SGLE EP=HMACH*( ONE + ABS(PAR(ICP(I))) ) */ par[icp[i]] += ep; icnd_(ndim, &par[1], &icp[1], nint, &u[1], &uold[1], &udot[1], &upold[ 1], bldif_1.f1zz, &c__0, &dint[dint_offset]); i__2 = *nint; for (j = 1; j <= i__2; ++j) { dint[j + (*ndim + icp[i]) * dint_dim1] = (bldif_1.f1zz[j - 1] - f[ j]) / ep; /* L5: */ } par[icp[i]] -= ep; /* L6: */ } return 0; } /* icni_ */ /* Subroutine */ int fopi_(ndim, u, icp, par, ijac, f, dfdu, dfdp) integer *ndim; doublereal *u; integer *icp; doublereal *par; integer *ijac; doublereal *f, *dfdu, *dfdp; { /* System generated locals */ integer i__1, i__2; doublereal d__1, d__2; /* Local variables */ extern /* Subroutine */ int fopt_(); static integer i, j; static doublereal f1, f2, ep, umx; /* ---------- ---- */ /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Interface subroutine to user supplied FOPT. */ /* Generate the objective function. */ /* Parameter adjustments */ --dfdp; --dfdu; --par; --icp; --u; /* Function Body */ fopt_(ndim, &u[1], &icp[1], &par[1], ijac, f, &dfdu[1], &dfdp[1]); if (blmax_1.jac == 1 || *ijac == 0) { return 0; } /* Generate the Jacobian by differencing. */ umx = blrcn_1.zero; i__1 = *ndim; for (i = 1; i <= i__1; ++i) { if ((d__1 = u[i], abs(d__1)) > umx) { umx = (d__2 = u[i], abs(d__2)); } /* SGLE IF( ABS(U(I)).GT.UMX)UMX= ABS(U(I)) */ /* L1: */ } ep = blrcn_1.hmach * (blrcn_1.one + umx); i__1 = *ndim; for (i = 1; i <= i__1; ++i) { i__2 = *ndim; for (j = 1; j <= i__2; ++j) { bldif_1.u1zz[j - 1] = u[j]; bldif_1.u2zz[j - 1] = u[j]; /* L2: */ } bldif_1.u1zz[i - 1] -= ep; bldif_1.u2zz[i - 1] += ep; fopt_(ndim, bldif_1.u1zz, &icp[1], &par[1], &c__0, &f1, &dfdu[1], & dfdp[1]); fopt_(ndim, bldif_1.u2zz, &icp[1], &par[1], &c__0, &f2, &dfdu[1], & dfdp[1]); dfdu[i] = (f2 - f1) / (ep * 2); /* L4: */ } i__1 = blicn_1.nfpar; for (i = 1; i <= i__1; ++i) { ep = blrcn_1.hmach * (blrcn_1.one + (d__1 = par[icp[i]], abs(d__1))); /* SGLE EP=HMACH*( ONE + ABS(PAR(ICP(I))) ) */ par[icp[i]] += ep; fopt_(ndim, &u[1], &icp[1], &par[1], &c__0, &f1, &dfdu[1], &dfdp[1]); dfdp[icp[i]] = (f1 - *f) / ep; par[icp[i]] -= ep; /* L6: */ } return 0; } /* fopi_ */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Installation Dependent Subroutines for Timing AUTO */ /* ----------------------------------------------------------------------- */ /* ----------------------------------------------------------------------- */ /* Subroutine */ int autim0_(t) doublereal *t; { extern doublereal etime_(); static real timaray[2]; /* ---------- ------ */ /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Set initial time for measuring CPU time used. */ *t = etime_(timaray); return 0; } /* autim0_ */ /* Subroutine */ int autim1_(t) doublereal *t; { extern doublereal etime_(); static real timaray[2]; /* ---------- ------ */ /* SGLE IMPLICIT REAL (A-H,O-Z) */ /* Set final time for measuring CPU time used. */ *t = etime_(timaray); return 0; } /* autim1_ */