/****************************************************************************** This collection o c-function for the stiff integration has been converted from F77 by mean of f2c and then made indipendent from the libf2c. Further it has been included a C-function jacket that overrides the F77 calling convention (by address) by a more suitable c-calling convention. Functions in this file: Radau5Integrator = the jacket to function radau5 radau5_ = this check parameter and allocate the needed memory and call the core integrator radcor_ = the core integrator, is this one that make count contr5rad = provide the collocation solution between two steps of integration Those functions compose the linear algebra nedded by the previous functions dec_, dech_, decc_, dechc_, decb_, decbc_ sol_, solh_, dolc_, solhc_, solb_, solbc_ elmhes_, decomr_, decomc_ slvrad_ estrad_ Here it is how must be defined the external functions: Rigth hand side of the system void RHS(long int n,double t,double *X,double *DXDT); n = order of the system t = time where to evaluate the RHS X[0 .. n-1] = vector of state variables where to evaluate the RHS DXDT[0 .. n-1] = vector of the values of the RHS Jacobian of the RHS void JAC(long int n, double t, double *X, double *DFDX, long int rowDFDX); n = order of the system t = time where to evaluate the Jacobian X[0 .. n-1] = vector of the state variable where to evaluate the RHS DFDX[0 .. rowDFDX,0 .. n-1] = matrix of the values of the Jacobian of the RHS rowDFDX = the colum-lenght of the matrix DFDX, it is furnished by the calling program. If the jacobian is supposed to be full colDFDX is equal to n and the partial derivatives are stored in DFDX as DFDX[i,j] = PARTIAL F(i+1) / PARTIAL X(j+1) else if the jacobian is taken as banded the partial derivatives are stored diagonal-wise as DFDY[i-j+mujac+1,j] = PARTIAL F(i+1) / PARTIAL Y(j+1) where mujac is the upper diagonal that is non zero and rowDFDX=mljac+mujac+1 where mljac is lowest diagonal that is non zero Mass maxtrix of the system void MAS(long int n, double t, double *M, long int rowM); n = order of the system M[0 .. rowM,0 .. n-1] = mass matrix of the system rowM = the colum-lenght of the matrix M, it is furnished by the calling program. If the matrix M is supposed to be full colM is equal to n and it is stored as is M[i,j] = element(i+1,j+1,M) else if the matrix M is taken as banded its elements are stored diagonal-wise as M[i-j+mumas+1,j] = element(i+1,j+1,M) where mumas is the upper diagonal that is non zero and rowM=mlmas+mumas+1 where mlmas is lowest diagonal that is non zero Every step function int EVRSTP(long int n,double told,double tnew,double *X,double *CV); n = order of the system told = time of the old step tnew = time of the new step X[0 .. n-1] = vector of the state variable at tnew CV = vector of the coefficent needed by the function contr5rad, see below returned values 0 everithing ok 1 stopped by user 2 something wrong in some evaluation This function is called at every end of step. By mean of this function it is possible to check zeros of user function along the trajectories, printout the solution, check if the user want to stop the integration ... During calls to EVRSTP, a continuous solution for the interval [told,tnew] is available through the function Xi(t) = (double) contr5rad(long int index,double t,double *CV) which provides an approximation to the i-th component of the solution at the time t. The value of t should lie in the interval [told,tnew]. Note that the index is the real index of the i-th vatiable (1,2,3 ...) !!!!!!!!!!!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!!!!!!!!!!!!!! DO NOT CHANGES THE ENTRIES OF THE VECTOR CV ******************************************************************************/ /********************* BEGIN of Radau5Integrator ****************************/ /****************************************************************************** int Radau5Integrator( long int n, int (*RHS)(), double tin, double tend,double *tcur, double *X0, double h, double rtol, double atol, int (*JAC)(), long int ijac, long int mljac, long int mujac, int (*MAS)(), long int imas, long int mlmas, long int mumas, double *rpar, long int *ipar, int (*EVRSTP)() ) Numerical solution of a stiff (or differential algebraic) system of first order ordinary differential equations of the kind M*X'=F(t,X) The system can be (linearly) implicit (the mass-matrix M is not the identity matrix) or explicit (M=I). The method used is an implicit runge-kutta method (radau IIa) of order 5 with step size control and continuous output. c.f. section iv.8 authors: e. hairer and g. wanner universite de geneve, dept. de mathematiques ch-1211 geneve 24, switzerland e-mail: hairer@divsun.unige.ch, wanner@divsun.unige.ch adapted in C by Oscar De Feo with the help of f2c translator odefeo@circhp.epfl.ch this code is part of the book: e. hairer and g. wanner, solving ordinary differential equations ii. stiff and differential-algebraic problems. springer series in computational mathematics 14, springer-verlag 1991, second edition 1996. version of September 10, 1997 input parameters ---------------- - n Dimension of the system - RHS Name of the subroutine computing the value of F(t,Y), see above for its definition - tin Initial t value - tend Final t value - X0[0 .. n-1] Vector of the state variable at tin - h Initial step size guess; for stiff equations with initial transient, h=1.d0/(norm of f'), usually 1.d-3 or 1.d-5, is good. This choice is not very important, the step size is quickly adapted. if h=0.d0, the code puts h=1.d-6. - rtol,atol Relative and absolute error tolerances. The code keep the local error of every state variable Xi below the value rtol*|Xi|+atol. - JAC Name of the subroutine which computes the partial derivatives of F(t,X) with respect to X. This routine is only called if ijac=1; supply a dummy subroutine for ijac=0. See above for its definition. - ijac Switch for the computation of the jacobian: ijac=0: jacobian is computed internally by finite differences, subroutine "jac" is never called. ijac=1: jacobian is supplied by subroutine JAC. - mljac Switch for the banded structure of the jacobian: mljac=n: jacobian is a full matrix. the linear algebra is done by full-matrix gauss-elimination. 0<=mljac= number of non-zero diagonals below the main diagonal). - mujac Upper bandwith of jacobian matrix (>= number of non-zero diagonals above the main diagonal). Need not be defined if mljac=n. - MAS Name of the subroutine computing the mass-matrix M. If imas=0, this matrix is assumed to be the identity matrix and needs not to be defined;supply a dummy subroutine in this case. if imas=1, matrix M is supplied by subroutine MAS. See above for its definition. - imas Gives information on the mass-matrix: imas=0: m is supposed to be the identity matrix, MAS is never called. imas=1: mass-matrix is supplied. - mlmas Switch for the banded structure of the mass-matrix: mlmas=n: the full matrix case. the linear algebra is done by full-matrix gauss-elimination. 0<=mlmas= number of non-zero diagonals below the main diagonal). mlmas is supposed to be lower ot equal to mljac. - mumas Upper bandwith of mass-matrix (>= number of non-zero diagonals above the main diagonal). Need not be defined if mlmas=n. mumas is supposed to be lower or equal mujac. - rpar (must be an array of 20 elements) Array of real advanced parameter. rpar(1), rpar(2),.., rpar(20) serve as parameters for the code. For standard use of the code rpar(1),..,rpar(20) must be set to zero before calling. see below for a more sophisticated use. - ipar (must be an array of 20 elements) Array of integer advanced parameter. ipar(1),ipar(2),...,ipar(20) serve as parameters for the code. For standard use, set ipar(1),.., ipar(20) to zero before calling. - evrstp Name of the subroutine that is called at every step during the integration. See above for its definition. returned value -------------- The returned value idid has the following meanings idid= 1 computation succesful but interrupted by user idid= 0 computation successful idid=-1x input is not consistent idid=-11 round factor incosistent with inner constants idid=-12 tolerances are too small idid=-13 number of newton iteration cannot be negative idid=-14 the number of indexed variables is inconsistent with the dimension n idid=-15 incosistent input for parameter m1 and/or m2 idid=-16 unacceptable value for step size control safety factor idid=-17 unacceptable value for the jacobian computation cost idid=-18 inconsistent newton stopping criterion and relative tolerance idid=-19 inconsistent changing step size criterion parameter idid=-110 inconsistent max and min step size ratio idid=-111 the bandwidth of mas cannot be smaller than thos of jac idid=-112 the hessenberg option can be applied only for ode with full jacobian idid=-2 step size becomes too small idid=-3 matrix is repeatedly singular idid=-4 computation of an utf failed sophisticated setting of parameters ----------------------------------- Several parameters of the code are tuned at priori to make it work, they may be changed by setting rpar(1),... as well as ipar(1),... different from zero. For zero input, the code chooses default values: - ipar(1) If ipar(1) is not 0, the code transforms the jacobian matrix to hessenberg form. This is particularly advantageous for large systems with full jacobian. It does not work for banded jacobian (mljac 1. The function-subroutine should be written such that the index 1,2,3 variables appear in this order. In estimating the error the index 2 variables are multiplied by h, the index 3 variables by h**2. - ipar(5) dimension of the index 1 variables (must be > 0). For ode's this equals the dimension of the system. Default ipar(5)=n. - ipar(6) Dimension of the index 2 variables. default ipar(6)=0. - ipar(7) Dimension of the index 3 variables. default ipar(7)=0. - ipar(8) Switch for step size strategy. If ipar(8) is 1 mod. predictive controller (gustafsson) If ipar(8) is 2 classical step size control The default value (for ipar(8)=0) is ipar(8)=1. The choice ipar(8)=1 seems to produce safer results; for simple problems, the choice ipar(8)=2 produces often slightly faster runs If the differential system has the special structure that y(i)' = y(i+m2) for i=1,...,m1, with m1 a multiple of m2, a substantial gain in computertime can be achieved by setting the parameters ipar(9) and ipar(10). e.g., for second order systems p'=v, v'=g(p,v), where p and v are vectors of dimension n/2, one has to put m1=m2=n/2. - ipar(9) The value of m1. Default m1=0. - ipar(10) The value of m2. Default m2=m1. For m1>0 some of the input parameters have different meanings: JAC: only the elements of the non-trivial part of the jacobian have to be stored if (mljac=(n-m1)) the jacobian is supposed to be full DFDX[i,j] = partial f(i+m1+1) / partial x(j+1) for i=0,n-m1-1 and j=0,n-1. else, the jacobian is banded ( m1 = m2 * mm ) DFDX(i-j+mujac+1,j+k*m2) = partial f(i+m1+1)/partial x(j+k*m2+1) for i=0,mljac+mujac and j=0,m2-1 and k=0,mm. mljac: mljac=n-m1: if the non-trivial part of the jacobian is full 0<=mljac /*F77 Equivalent constants*/ #define TRUE_ (1) #define FALSE_ (0) /*F77 Equivalent functions*/ #undef min #undef max #undef copysign #define min(a,b) ((a) <= (b) ? (a) : (b)) #define max(a,b) ((a) >= (b) ? (a) : (b)) #define copysign(a,b) ((b)>=0?fabs((a)):-fabs((a))) /*Definition of the F77 equivalent type*/ typedef long int integer; typedef double doublereal; typedef long int logical; typedef /* Unknown type subroutine*/ int (*U_fp)(); typedef /* Subroutine */ int (*S_fp)(); /* Common Block of F77 Declarations */ struct { integer nn, nn2, nn3, nn4; doublereal xsol, hsol, c2m1, c1m1; } conra5_; #define conra5_1 conra5_ struct { integer mle, mue, mbjac, mbb, mdiag, mdiff, mbdiag; } linal_; #define linal_1 linal_ /* Table of constant values */ static doublereal c_b2 = .5; static doublereal c_b4 = 81.; static doublereal c_b5 = .33333333333333331; static doublereal c_b6 = 9.; static doublereal c_b16 = 1.; static doublereal c_b26 = .8; static doublereal c_b28 = .25; static integer c__1 = 1; /************************** BEGIN of radau5_ ********************************/ /* Subroutine */ int radau5_(integer *n, voidsub fcn, doublereal *x, doublereal * xend, doublereal *y, doublereal *h__, doublereal *rtol, doublereal * atol, voidsub jac, integer *ijac, integer *mljac, integer *mujac, voidsub mas, integer *imas, integer *mlmas, integer *mumas, doublereal *rpar, integer *ipar, intfunc evrstp, integer *idid) { /* System generated locals */ doublereal d__1, d__2, d__3, d__4; /* Local variables */ doublereal facl; integer ndec, njac; doublereal facr, safe; integer ijob, nfcn; logical pred; doublereal hmax; integer nmax; doublereal thet, expm; integer nsol; doublereal quot; integer iee2i, iee2r, ieip1, ieip2, nind1, nind2, nind3; doublereal quot1, quot2; integer iejac, ldjac; logical jband; integer iecon, iemas, ldmas, ieiph; doublereal fnewt; integer nstep, m1, m2; doublereal tolst; integer ldmas2, iescal, naccpt; extern /* Subroutine */ int radcor_(integer *, voidsub, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, voidsub, integer *, integer *, integer *, voidsub, integer *, integer *, integer *, doublereal *, doublereal * , doublereal *, doublereal *, doublereal *, doublereal *, integer *, integer *, logical *, integer *, integer *, integer *, logical *, doublereal *, doublereal *, integer *, integer *, integer *, logical *, logical *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, integer *, integer *, doublereal *, integer *, integer *, integer *, integer *, integer *, integer *, integer *, U_fp, integer *); integer nrejct; logical implct; integer istore; logical startn; doublereal uround; integer nm1, nit, iee1, ief1, lde1, ief2, ief3, iey0, iez1, iez2, iez3; /* *** *** *** *** *** *** *** */ /* SETTING THE PARAMETERS */ /* *** *** *** *** *** *** *** */ /* Parameter adjustments */ --y; --rpar; --ipar; /* Function Body */ nfcn = 0; njac = 0; nstep = 0; naccpt = 0; nrejct = 0; ndec = 0; nsol = 0; /* -------- UROUND SMALLEST NUMBER SATISFYING 1.0D0+UROUND>1.0D0 */ if (rpar[1] == 0.) { uround = 1e-16; } else { uround = rpar[1]; if (uround <= 1e-19 || uround >= 1.) { /* return(-11); */ *idid = -11; return 0; } } /* -------- CHECK AND CHANGE THE TOLERANCES */ expm = .66666666666666663; if (*atol <= 0. || *rtol <= uround * 10.) { /* return(-12); */ *idid = -12; return 0; } else { quot = *atol / *rtol; *rtol = pow(*rtol, expm) * .1; *atol = *rtol * quot; } /* -------- NIT MAXIMAL NUMBER OF NEWTON ITERATIONS */ if (ipar[3] == 0) { nit = 7; } else { nit = ipar[3]; if (nit <= 0) { /* return(-13); */ *idid = -13; return 0; } } /* -------- STARTN SWITCH FOR STARTING VALUES OF NEWTON ITERATIONS */ if (ipar[4] == 0) { startn = FALSE_; } else { startn = TRUE_; } /* -------- PARAMETER FOR DIFFERENTIAL-ALGEBRAIC COMPONENTS */ nind1 = ipar[5]; nind2 = ipar[6]; nind3 = ipar[7]; if (nind1 == 0) { nind1 = *n; } if (nind1 + nind2 + nind3 != *n) { /* return(-14); */ *idid = -14; return 0; } /* -------- PRED STEP SIZE CONTROL */ if (ipar[8] <= 1) { pred = TRUE_; } else { pred = FALSE_; } /* -------- PARAMETER FOR SECOND ORDER EQUATIONS */ m1 = ipar[9]; m2 = ipar[10]; nm1 = *n - m1; if (m1 == 0) { m2 = *n; } if (m2 == 0) { m2 = m1; } if (m1 < 0 || m2 < 0 || m1 + m2 > *n) { /* return(-15); */ *idid = -15; return 0; } /* --------- SAFE SAFETY FACTOR IN STEP SIZE PREDICTION */ if (rpar[2] == 0.) { safe = .9; } else { safe = rpar[2]; if (safe <= .001 || safe >= 1.) { /* return(-16); */ *idid = -16; return 0; } } /* ------ THET DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED; */ if (rpar[3] == 0.) { thet = .001; } else { thet = rpar[3]; if (thet >= 1.) { /* return(-17); */ *idid = -17; return 0; } } /* --- FNEWT STOPPING CRIERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1. */ tolst = *rtol; if (rpar[4] == 0.) { /* Computing MAX */ /* Computing MIN */ d__3 = .03, d__4 = pow(tolst,c_b2); d__1 = uround * 10 / tolst, d__2 = min(d__3,d__4); fnewt = max(d__1,d__2); } else { fnewt = rpar[4]; if (fnewt <= uround / tolst) { /* return(-18); */ *idid = -18; return 0; } } /* --- QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST. */ if (rpar[5] == 0.) { quot1 = 1.; } else { quot1 = rpar[5]; } if (rpar[6] == 0.) { quot2 = 1.2; } else { quot2 = rpar[6]; } if (quot1 > 1. || quot2 < 1.) { /* return(-19); */ *idid = -19; return 0; } /* -------- MAXIMAL STEP SIZE */ if (rpar[7] == 0.) { hmax = *xend - *x; } else { hmax = rpar[7]; } /* ------- FACL,FACR PARAMETERS FOR STEP SIZE SELECTION */ if (rpar[8] == 0.) { facl = 5.; } else { facl = 1. / rpar[8]; } if (rpar[9] == 0.) { facr = .125; } else { facr = 1. / rpar[9]; } if (facl < 1. || facr > 1.) { /* return(-110); */ *idid = -110; return 0; } /* *** *** *** *** *** *** *** *** *** *** *** *** *** */ /* COMPUTATION OF ARRAY ENTRIES */ /* *** *** *** *** *** *** *** *** *** *** *** *** *** */ /* ---- IMPLICIT, BANDED OR NOT ? */ implct = *imas != 0; jband = *mljac < nm1; /* -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS --- */ /* -- JACOBIAN AND MATRICES E1, E2 */ if (jband) { ldjac = *mljac + *mujac + 1; lde1 = *mljac + ldjac; } else { *mljac = nm1; *mujac = nm1; ldjac = nm1; lde1 = nm1; } /* -- MASS MATRIX */ if (implct) { if (*mlmas != nm1) { ldmas = *mlmas + *mumas + 1; if (jband) { ijob = 4; } else { ijob = 3; } } else { ldmas = nm1; ijob = 5; } /* ------ BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF "JAC" */ if (*mlmas > *mljac || *mumas > *mujac) { /* return(-111); */ *idid = -111; return 0; } } else { ldmas = 0; if (jband) { ijob = 2; } else { ijob = 1; if (*n > 2 && ipar[1] != 0) { ijob = 7; } } } ldmas2 = max(1,ldmas); /* ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN */ if ((implct || jband) && ijob == 7) { /* return(-112); */ *idid = -112; return 0; } /* ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- */ iez1 = 1; iez2 = iez1 + *n; iez3 = iez2 + *n; iey0 = iez3 + *n; iescal = iey0 + *n; ief1 = iescal + *n; ief2 = ief1 + *n; ief3 = ief2 + *n; iecon = ief3 + *n; iejac = iecon + (*n << 2); iemas = iejac + *n * ldjac; iee1 = iemas + nm1 * ldmas; iee2r = iee1 + nm1 * lde1; iee2i = iee2r + nm1 * lde1; /* ------ TOTAL STORAGE REQUIREMENT ----------- */ istore = iee2i + nm1 * lde1; /*Now allocates the memory in dynamic way*/ work=(doublereal *)MemNew(istore*sizeof(doublereal)); /* ------- ENTRY POINTS FOR INTEGER WORKSPACE ----- */ ieip1 = 1; ieip2 = ieip1 + nm1; ieiph = ieip2 + nm1; /* --------- TOTAL REQUIREMENT --------------- */ istore = ieiph + nm1; /*Now allocates the memory in dynamic way*/ iwork=(integer *)MemNew(istore*sizeof(integer)); /* -------- CALL TO CORE INTEGRATOR ------------ */ radcor_(n, (voidsub)fcn, x, xend, &y[1], &hmax, h__, rtol, atol, (voidsub)jac, ijac, mljac, mujac, (voidsub)mas, mlmas, mumas, &nmax, &uround, & safe, &thet, &fnewt, "1, "2, &nit, &ijob, &startn, &nind1, &nind2, &nind3, &pred, &facl, &facr, &m1, &m2, &nm1, &implct, & jband, &ldjac, &lde1, &ldmas2, &work[iez1 - 1], &work[iez2 - 1], & work[iez3 - 1], &work[iey0 - 1], &work[iescal - 1], &work[ief1 - 1], &work[ief2 - 1], &work[ief3 - 1], &work[iejac - 1], &work[ iee1 - 1], &work[iee2r - 1], &work[iee2i - 1], &work[iemas - 1], & iwork[ieip1 - 1], &iwork[ieip2 - 1], &iwork[ieiph - 1], &work[ iecon - 1], &nfcn, &njac, &nstep, &naccpt, &nrejct, &ndec, &nsol, (intfunc)evrstp, idid); ipar[14] = nfcn; ipar[15] = njac; ipar[16] = nstep; ipar[17] = naccpt; ipar[18] = nrejct; ipar[19] = ndec; ipar[20] = nsol; /* -------- RESTORE TOLERANCES */ expm = 1. / expm; quot = *atol / *rtol; d__1 = *rtol * 10.; *rtol = pow(d__1,expm); *atol = *rtol * quot; /* ----------- RETURN ----------- */ iwork=MemFree(iwork); work=MemFree(work); return 0; } /* radau5_ */ /*************************** END of radau5_ *********************************/ /************************** BEGIN of radcor_ ********************************/ /* Subroutine */ int radcor_(integer *n, voidsub fcn, doublereal *x, doublereal * xend, doublereal *y, doublereal *hmax, doublereal *h__, doublereal * rtol, doublereal *atol, voidsub jac, integer *ijac, integer *mljac, integer *mujac, voidsub mas, integer *mlmas, integer *mumas, integer * nmax, doublereal *uround, doublereal *safe, doublereal *thet, doublereal *fnewt, doublereal *quot1, doublereal *quot2, integer *nit, integer *ijob, logical *startn, integer *nind1, integer *nind2, integer *nind3, logical *pred, doublereal *facl, doublereal *facr, integer *m1, integer *m2, integer *nm1, logical *implct, logical * banded, integer *ldjac, integer *lde1, integer *ldmas, doublereal *z1, doublereal *z2, doublereal *z3, doublereal *y0, doublereal *scal, doublereal *f1, doublereal *f2, doublereal *f3, doublereal *fjac, doublereal *e1, doublereal *e2r, doublereal *e2i, doublereal *fmas, integer *ip1, integer *ip2, integer *iphes, doublereal *cont, integer *nfcn, integer *njac, integer *nstep, integer *naccpt, integer * nrejct, integer *ndec, integer *nsol, intfunc evrstp, integer *idid) { /* System generated locals */ integer fjac_dim1, fjac_offset, fmas_dim1, fmas_offset, e1_dim1, e1_offset, e2r_dim1, e2r_offset, e2i_dim1, e2i_offset, i__1, i__2, i__3, i__4; doublereal d__1, d__2, d__3, d__4; /* Local variables */ doublereal c1mc2, beta; integer lbeg; doublereal alph, hold; integer lend; doublereal delt, hnew; logical last; doublereal xold, dyno, dyth, hopt; integer newt; doublereal quot; integer i__; doublereal hhfac; integer k, j, l; doublereal betan, alphn, denom, theta, ysafe, hmaxn; integer nsing; doublereal c1, c2; logical first; integer j1, irtrn, nrsol, nsolu, n2; doublereal xosol; integer n3; doublereal u1, a1, a2, a3, qnewt, acont3; logical index1, index2, index3; doublereal ak; logical caljac; integer md; doublereal t11, t12, t13, t21, t22, t23, t31, faccon; logical calhes; integer mm; extern /* Subroutine */ int decomc_(integer *, doublereal *, integer *, doublereal *, integer *, integer *, integer *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, integer *, integer *, integer *); doublereal erracc, qt; integer mujacj; extern /* Subroutine */ int decomr_(integer *, doublereal *, integer *, doublereal *, integer *, integer *, integer *, integer *, integer *, integer *, doublereal *, doublereal *, integer *, integer *, integer *, integer *, logical *, integer *); logical reject; doublereal facgus; integer mujacp; extern /* Subroutine */ int estrad_(integer *, doublereal *, integer *, integer *, integer *, doublereal *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, voidsub, integer *, doublereal *, doublereal *, integer *, doublereal *, integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, integer *, doublereal *, doublereal *, logical *, logical *, doublereal *); doublereal dd1, dd2, dd3, posneg, ak1, ak2, ak3; extern /* Subroutine */ int slvrad_(integer *, doublereal *, integer *, integer *, integer *, doublereal *, integer *, integer *, integer *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, integer *, integer *, integer *, integer *); doublereal dynold, thqold, f1i, c1q, c2q, c3q, f2i, f3i, z1i, z2i, z3i, sq6, fac, ti11, cno; integer lrc; doublereal ti12, ti13, ti21, ti22, ti23, ti31, ti32, ti33; integer ier; doublereal xph, thq, err, fac1, cfac, hacc; /* ---------------------------------------------------------- */ /* CORE INTEGRATOR FOR RADAU5 */ /* PARAMETERS SAME AS IN RADAU5 WITH WORKSPACE ADDED */ /* ---------------------------------------------------------- */ /* DECLARATIONS */ /* ---------------------------------------------------------- */ /* *** *** *** *** *** *** *** */ /* INITIALISATIONS */ /* *** *** *** *** *** *** *** */ /* --------- DUPLIFY N FOR COMMON BLOCK CONT ----- */ /* Parameter adjustments */ --cont; --f3; --f2; --f1; --scal; --y0; --z3; --z2; --z1; --y; --iphes; --ip2; --ip1; fjac_dim1 = *ldjac; fjac_offset = fjac_dim1 + 1; fjac -= fjac_offset; e2i_dim1 = *lde1; e2i_offset = e2i_dim1 + 1; e2i -= e2i_offset; e2r_dim1 = *lde1; e2r_offset = e2r_dim1 + 1; e2r -= e2r_offset; e1_dim1 = *lde1; e1_offset = e1_dim1 + 1; e1 -= e1_offset; fmas_dim1 = *ldmas; fmas_offset = fmas_dim1 + 1; fmas -= fmas_offset; /* Function Body */ conra5_1.nn = *n; conra5_1.nn2 = *n << 1; conra5_1.nn3 = *n * 3; lrc = *n << 2; /* -------- CHECK THE INDEX OF THE PROBLEM ----- */ index1 = *nind1 != 0; index2 = *nind2 != 0; index3 = *nind3 != 0; /* ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ---------- */ if (*implct) { (*mas)(*nm1, &fmas[fmas_offset], *ldmas); } /* ---------- CONSTANTS --------- */ sq6 = sqrt(6.); c1 = (4. - sq6) / 10.; c2 = (sq6 + 4.) / 10.; conra5_1.c1m1 = c1 - 1.; conra5_1.c2m1 = c2 - 1.; c1mc2 = c1 - c2; dd1 = -(sq6 * 7. + 13.) / 3.; dd2 = (sq6 * 7. - 13.) / 3.; dd3 = -.33333333333333331; u1 = (pow(c_b4,c_b5) + 6. - pow(c_b6,c_b5)) / 30.; alph = (12. - pow(c_b4,c_b5) + pow(c_b6,c_b5)) / 60.; beta = (pow(c_b4,c_b5) + pow(c_b6,c_b5)) * sqrt(3.) / 60.; /* Computing 2nd power */ d__1 = alph; /* Computing 2nd power */ d__2 = beta; cno = d__1 * d__1 + d__2 * d__2; u1 = 1. / u1; alph /= cno; beta /= cno; t11 = .091232394870892942792; t12 = -.14125529502095420843; t13 = -.030029194105147424492; t21 = .24171793270710701896; t22 = .20412935229379993199; t23 = .38294211275726193779; t31 = .96604818261509293619; ti11 = 4.325579890063155351; ti12 = .33919925181580986954; ti13 = .54177053993587487119; ti21 = -4.1787185915519047273; ti22 = -.32768282076106238708; ti23 = .47662355450055045196; ti31 = -.50287263494578687595; ti32 = 2.5719269498556054292; ti33 = -.59603920482822492497; if (*m1 > 0) { *ijob += 10; } d__1 = *xend - *x; posneg = copysign(c_b16,d__1); /* Computing MIN */ d__2 = fabs(*hmax), d__3 = (d__1 = *xend - *x, fabs(d__1)); hmaxn = min(d__2,d__3); if (fabs(*h__) <= *uround * 10.) { *h__ = 1e-6; } /* Computing MIN */ d__1 = fabs(*h__); *h__ = min(d__1,hmaxn); *h__ = copysign(*h__,posneg); hold = *h__; reject = FALSE_; first = TRUE_; last = FALSE_; if ((*x + *h__ * 1.0001 - *xend) * posneg >= 0.) { *h__ = *xend - *x; last = TRUE_; } faccon = 1.; cfac = *safe * ((*nit << 1) + 1); nsing = 0; xold = *x; irtrn = 0; nrsol = 1; xosol = xold; conra5_1.xsol = *x; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { cont[i__] = y[i__]; } nsolu = *n; conra5_1.hsol = hold; /* The call has been modified to be C-like */ irtrn=(*evrstp)(nsolu, xosol, conra5_1.xsol, &y[1], &cont[1]); /* End of the modfication */ if (irtrn == 2) { goto L178; } if (irtrn == 1) { goto L179; } linal_1.mle = *mljac; linal_1.mue = *mujac; linal_1.mbjac = *mljac + *mujac + 1; linal_1.mbb = *mlmas + *mumas + 1; linal_1.mdiag = linal_1.mle + linal_1.mue + 1; linal_1.mdiff = linal_1.mle + linal_1.mue - *mumas; linal_1.mbdiag = *mumas + 1; n2 = *n << 1; n3 = *n * 3; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { scal[i__] = *atol + *rtol * (d__1 = y[i__], fabs(d__1)); } hhfac = *h__; (*fcn)(*n, *x, &y[1], &y0[1]); ++(*nfcn); /* --- BASIC INTEGRATION STEP */ L10: /* *** *** *** *** *** *** *** */ /* COMPUTATION OF THE JACOBIAN */ /* *** *** *** *** *** *** *** */ ++(*njac); if (*ijac == 0) { /* --- COMPUTE JACOBIAN MATRIX NUMERICALLY */ if (*banded) { /* --- JACOBIAN IS BANDED */ mujacp = *mujac + 1; md = min(linal_1.mbjac,*m2); i__1 = *m1 / *m2 + 1; for (mm = 1; mm <= i__1; ++mm) { i__2 = md; for (k = 1; k <= i__2; ++k) { j = k + (mm - 1) * *m2; L12: f1[j] = y[j]; /* Computing MAX */ d__2 = 1e-5, d__3 = (d__1 = y[j], fabs(d__1)); f2[j] = sqrt(*uround * max(d__2,d__3)); y[j] += f2[j]; j += md; if (j <= mm * *m2) { goto L12; } (*fcn)(*n, *x, &y[1], &cont[1]); j = k + (mm - 1) * *m2; j1 = k; /* Computing MAX */ i__3 = 1, i__4 = j1 - *mujac; lbeg = max(i__3,i__4) + *m1; L14: /* Computing MIN */ i__3 = *m2, i__4 = j1 + *mljac; lend = min(i__3,i__4) + *m1; y[j] = f1[j]; mujacj = mujacp - j1 - *m1; i__3 = lend; for (l = lbeg; l <= i__3; ++l) { fjac[l + mujacj + j * fjac_dim1] = (cont[l] - y0[l]) / f2[j]; } j += md; j1 += md; lbeg = lend + 1; if (j <= mm * *m2) { goto L14; } } } } else { /* --- JACOBIAN IS FULL */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { ysafe = y[i__]; /* Computing MAX */ d__1 = 1e-5, d__2 = fabs(ysafe); delt = sqrt(*uround * max(d__1,d__2)); y[i__] = ysafe + delt; (*fcn)(*n, *x, &y[1], &cont[1]); i__2 = *n; for (j = *m1 + 1; j <= i__2; ++j) { fjac[j - *m1 + i__ * fjac_dim1] = (cont[j] - y0[j]) / delt; } y[i__] = ysafe; } } } else { /* --- COMPUTE JACOBIAN MATRIX ANALYTICALLY */ (*jac)(*n, *x, &y[1], &fjac[fjac_offset], *ldjac); } caljac = TRUE_; calhes = TRUE_; L20: /* --- COMPUTE THE MATRICES E1 AND E2 AND THEIR DECOMPOSITIONS */ fac1 = u1 / *h__; alphn = alph / *h__; betan = beta / *h__; decomr_(n, &fjac[fjac_offset], ldjac, &fmas[fmas_offset], ldmas, mlmas, mumas, m1, m2, nm1, &fac1, &e1[e1_offset], lde1, &ip1[1], &ier, ijob, &calhes, &iphes[1]); if (ier != 0) { goto L78; } decomc_(n, &fjac[fjac_offset], ldjac, &fmas[fmas_offset], ldmas, mlmas, mumas, m1, m2, nm1, &alphn, &betan, &e2r[e2r_offset], &e2i[ e2i_offset], lde1, &ip2[1], &ier, ijob); if (ier != 0) { goto L78; } ++(*ndec); L30: ++(*nstep); if (fabs(*h__) * .1 <= fabs(*x) * *uround) { goto L177; } if (index2) { i__1 = *nind1 + *nind2; for (i__ = *nind1 + 1; i__ <= i__1; ++i__) { scal[i__] /= hhfac; } } if (index3) { i__1 = *nind1 + *nind2 + *nind3; for (i__ = *nind1 + *nind2 + 1; i__ <= i__1; ++i__) { scal[i__] /= hhfac * hhfac; } } xph = *x + *h__; /* *** *** *** *** *** *** *** */ /* STARTING VALUES FOR NEWTON ITERATION */ /* *** *** *** *** *** *** *** */ if (first || *startn) { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { z1[i__] = 0.; z2[i__] = 0.; z3[i__] = 0.; f1[i__] = 0.; f2[i__] = 0.; f3[i__] = 0.; } } else { c3q = *h__ / hold; c1q = c1 * c3q; c2q = c2 * c3q; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { ak1 = cont[i__ + *n]; ak2 = cont[i__ + n2]; ak3 = cont[i__ + n3]; z1i = c1q * (ak1 + (c1q - conra5_1.c2m1) * (ak2 + (c1q - conra5_1.c1m1) * ak3)); z2i = c2q * (ak1 + (c2q - conra5_1.c2m1) * (ak2 + (c2q - conra5_1.c1m1) * ak3)); z3i = c3q * (ak1 + (c3q - conra5_1.c2m1) * (ak2 + (c3q - conra5_1.c1m1) * ak3)); z1[i__] = z1i; z2[i__] = z2i; z3[i__] = z3i; f1[i__] = ti11 * z1i + ti12 * z2i + ti13 * z3i; f2[i__] = ti21 * z1i + ti22 * z2i + ti23 * z3i; f3[i__] = ti31 * z1i + ti32 * z2i + ti33 * z3i; } } /* *** *** *** *** *** *** *** */ /* LOOP FOR THE SIMPLIFIED NEWTON ITERATION */ /* *** *** *** *** *** *** *** */ newt = 0; d__1 = max(faccon,*uround); faccon = pow(d__1,c_b26); theta = fabs(*thet); L40: if (newt >= *nit) { goto L78; } /* --- COMPUTE THE RIGHT-HAND SIDE */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { cont[i__] = y[i__] + z1[i__]; } d__1 = *x + c1 * *h__; (*fcn)(*n, d__1, &cont[1], &z1[1]); i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { cont[i__] = y[i__] + z2[i__]; } d__1 = *x + c2 * *h__; (*fcn)(*n, d__1, &cont[1], &z2[1]); i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { cont[i__] = y[i__] + z3[i__]; } (*fcn)(*n, xph, &cont[1], &z3[1]); *nfcn += 3; /* --- SOLVE THE LINEAR SYSTEMS */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { a1 = z1[i__]; a2 = z2[i__]; a3 = z3[i__]; z1[i__] = ti11 * a1 + ti12 * a2 + ti13 * a3; z2[i__] = ti21 * a1 + ti22 * a2 + ti23 * a3; z3[i__] = ti31 * a1 + ti32 * a2 + ti33 * a3; } slvrad_(n, &fjac[fjac_offset], ldjac, mljac, mujac, &fmas[fmas_offset], ldmas, mlmas, mumas, m1, m2, nm1, &fac1, &alphn, &betan, &e1[ e1_offset], &e2r[e2r_offset], &e2i[e2i_offset], lde1, &z1[1], &z2[ 1], &z3[1], &f1[1], &f2[1], &f3[1], &cont[1], &ip1[1], &ip2[1], & iphes[1], &ier, ijob); ++(*nsol); ++newt; dyno = 0.; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { denom = scal[i__]; /* Computing 2nd power */ d__1 = z1[i__] / denom; /* Computing 2nd power */ d__2 = z2[i__] / denom; /* Computing 2nd power */ d__3 = z3[i__] / denom; dyno = dyno + d__1 * d__1 + d__2 * d__2 + d__3 * d__3; } dyno = sqrt(dyno / n3); /* --- BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE */ if (newt > 1 && newt < *nit) { thq = dyno / dynold; if (newt == 2) { theta = thq; } else { theta = sqrt(thq * thqold); } thqold = thq; if (theta < .99) { faccon = theta / (1. - theta); i__1 = *nit - 1 - newt; dyth = faccon * dyno * pow(theta,(double)(i__1)) / *fnewt; if (dyth >= 1.) { /* Computing MAX */ d__1 = 1e-4, d__2 = min(20.,dyth); qnewt = max(d__1,d__2); d__1 = -1. / (*nit + 4. - 1 - newt); hhfac = pow(qnewt,d__1) * .8; *h__ = hhfac * *h__; reject = TRUE_; last = FALSE_; if (caljac) { goto L20; } goto L10; } } else { goto L78; } } dynold = max(dyno,*uround); i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { f1i = f1[i__] + z1[i__]; f2i = f2[i__] + z2[i__]; f3i = f3[i__] + z3[i__]; f1[i__] = f1i; f2[i__] = f2i; f3[i__] = f3i; z1[i__] = t11 * f1i + t12 * f2i + t13 * f3i; z2[i__] = t21 * f1i + t22 * f2i + t23 * f3i; z3[i__] = t31 * f1i + f2i; } if (faccon * dyno > *fnewt) { goto L40; } /* --- ERROR ESTIMATION */ estrad_(n, &fjac[fjac_offset], ldjac, mljac, mujac, &fmas[fmas_offset], ldmas, mlmas, mumas, h__, &dd1, &dd2, &dd3, (voidsub)fcn, nfcn, &y0[ 1], &y[1], ijob, x, m1, m2, nm1, &e1[e1_offset], lde1, &z1[1], & z2[1], &z3[1], &cont[1], &f1[1], &f2[1], &ip1[1], &iphes[1], & scal[1], &err, &first, &reject, &fac1); /* --- COMPUTATION OF HNEW */ /* --- WE REQUIRE .2<=HNEW/H<=8. */ /* Computing MIN */ d__1 = *safe, d__2 = cfac / (newt + (*nit << 1)); fac = min(d__1,d__2); /* Computing MAX */ /* Computing MIN */ d__3 = *facl, d__4 = pow(err,c_b28) / fac; d__1 = *facr, d__2 = min(d__3,d__4); quot = max(d__1,d__2); hnew = *h__ / quot; /* *** *** *** *** *** *** *** */ /* IS THE ERROR SMALL ENOUGH ? */ /* *** *** *** *** *** *** *** */ if (err < 1.) { /* --- STEP IS ACCEPTED */ first = FALSE_; ++(*naccpt); if (*pred) { /* --- PREDICTIVE CONTROLLER OF GUSTAFSSON */ if (*naccpt > 1) { /* Computing 2nd power */ d__2 = err; d__1 = d__2 * d__2 / erracc; facgus = hacc / *h__ * pow(d__1,c_b28) / *safe; /* Computing MAX */ d__1 = *facr, d__2 = min(*facl,facgus); facgus = max(d__1,d__2); quot = max(quot,facgus); hnew = *h__ / quot; } hacc = *h__; erracc = max(.01,err); } xold = *x; hold = *h__; *x = xph; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { y[i__] += z3[i__]; z2i = z2[i__]; z1i = z1[i__]; cont[i__ + *n] = (z2i - z3[i__]) / conra5_1.c2m1; ak = (z1i - z2i) / c1mc2; acont3 = z1i / c1; acont3 = (ak - acont3) / c2; cont[i__ + n2] = (ak - cont[i__ + *n]) / conra5_1.c1m1; cont[i__ + n3] = cont[i__ + n2] - acont3; } i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { scal[i__] = *atol + *rtol * (d__1 = y[i__], fabs(d__1)); } nrsol = *naccpt + 1; conra5_1.xsol = *x; xosol = xold; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { cont[i__] = y[i__]; } nsolu = *n; conra5_1.hsol = hold; /* The call has been modified to be C-like */ irtrn=(*evrstp)(nsolu, xosol, conra5_1.xsol, &y[1], &cont[1]); /* End of the modfication */ if (irtrn == 2) { goto L178; } if (irtrn == 1) { goto L179; } caljac = FALSE_; if (last) { *h__ = hopt; *idid = 0; return 0; } (*fcn)(*n, *x, &y[1], &y0[1]); ++(*nfcn); /* Computing MIN */ d__1 = fabs(hnew); hnew = posneg * min(d__1,hmaxn); hopt = hnew; hopt = min(*h__,hnew); if (reject) { /* Computing MIN */ d__1 = fabs(hnew), d__2 = fabs(*h__); hnew = posneg * min(d__1,d__2); } reject = FALSE_; if ((*x + hnew / *quot1 - *xend) * posneg >= 0.) { *h__ = *xend - *x; last = TRUE_; } else { qt = hnew / *h__; hhfac = *h__; if (theta <= *thet && qt >= *quot1 && qt <= *quot2) { goto L30; } *h__ = hnew; } hhfac = *h__; if (theta <= *thet) { goto L20; } goto L10; } else { /* --- STEP IS REJECTED */ reject = TRUE_; last = FALSE_; if (first) { *h__ *= .1; hhfac = .1; } else { hhfac = hnew / *h__; *h__ = hnew; } if (*naccpt >= 1) { ++(*nrejct); } if (caljac) { goto L20; } goto L10; } /* --- UNEXPECTED STEP-REJECTION */ L78: if (ier != 0) { ++nsing; if (nsing >= 5) { goto L176; } } *h__ *= .5; hhfac = .5; reject = TRUE_; last = FALSE_; if (caljac) { goto L20; } goto L10; /* --- FAIL EXIT */ L176: /* MATRIX IS REPEATEDLY SINGULAR */ /* return(-3); */ *idid = -3; return 0; L177: /* STEP SIZE BECOMES TOO SMALL */ /* return(-2) */ *idid = -2; return 0; L178: /* UNDEFINED USER TEST FUNCTION */ /* return(-4) */ *idid = -4; return 0; /* --- EXIT CAUSED BY EVRSTP */ L179: /* USER STOPPED */ /* return(1); */ *idid = 1; return 0; } /* radcor_ */ /*************************** END of radcor_ *********************************/ /************************** BEGIN of contr5rad ********************************/ /****************************************************************************** WARNING It has been modified such that the index of the solution (i), and the current time x are passed by value and not by address, becouse in such a way is easyest to use this function. ******************************************************************************/ doublereal contr5rad(integer isol, doublereal x, doublereal *cont) { /* System generated locals */ doublereal ret_val; /* Local variables */ doublereal s; /* ---------------------------------------------------------- */ /* THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT. IT PROVIDES AN */ /* APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT X. */ /* IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR */ /* THE LAST SUCCESSFULLY COMPUTED STEP (BY RADAU5). */ /* ---------------------------------------------------------- */ /* Parameter adjustments */ --cont; /* Function Body */ s = (x - conra5_1.xsol) / conra5_1.hsol; ret_val = cont[isol] + s * (cont[isol + conra5_1.nn] + (s - conra5_1.c2m1) * (cont[isol + conra5_1.nn2] + (s - conra5_1.c1m1) * cont[isol + conra5_1.nn3])); return ret_val; } /* contr5rad */ /*************************** END of contr5rad *********************************/ /****************************************************************************** !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! WARNING BEGIN OF LINEAR ALGEBRA FUNCTIONS NOTHING UNDER THIS LINE HAS BEEN TOUCHED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ******************************************************************************/ /* Subroutine */ int dec_(integer *n, integer *ndim, doublereal *a, integer * ip, integer *ier) { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3; doublereal d__1, d__2; /* Local variables */ static integer i__, j, k, m; static doublereal t; static integer nm1, kp1; /* VERSION REAL DOUBLE PRECISION */ /* ----------------------------------------------------------------------- */ /* MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION. */ /* INPUT.. */ /* N = ORDER OF MATRIX. */ /* NDIM = DECLARED DIMENSION OF ARRAY A . */ /* A = MATRIX TO BE TRIANGULARIZED. */ /* OUTPUT.. */ /* A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U . */ /* A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. */ /* IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW. */ /* IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O . */ /* IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE */ /* SINGULAR AT STAGE K. */ /* USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. */ /* DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N). */ /* IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. */ /* REFERENCE.. */ /* C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, */ /* C.A.C.M. 15 (1972), P. 274. */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; a_dim1 = *ndim; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ *ier = 0; ip[*n] = 1; if (*n == 1) { goto L70; } nm1 = *n - 1; i__1 = nm1; for (k = 1; k <= i__1; ++k) { kp1 = k + 1; m = k; i__2 = *n; for (i__ = kp1; i__ <= i__2; ++i__) { if ((d__1 = a[i__ + k * a_dim1], fabs(d__1)) > (d__2 = a[m + k * a_dim1], fabs(d__2))) { m = i__; } /* L10: */ } ip[k] = m; t = a[m + k * a_dim1]; if (m == k) { goto L20; } ip[*n] = -ip[*n]; a[m + k * a_dim1] = a[k + k * a_dim1]; a[k + k * a_dim1] = t; L20: if (t == 0.) { goto L80; } t = 1. / t; i__2 = *n; for (i__ = kp1; i__ <= i__2; ++i__) { /* L30: */ a[i__ + k * a_dim1] = -a[i__ + k * a_dim1] * t; } i__2 = *n; for (j = kp1; j <= i__2; ++j) { t = a[m + j * a_dim1]; a[m + j * a_dim1] = a[k + j * a_dim1]; a[k + j * a_dim1] = t; if (t == 0.) { goto L45; } i__3 = *n; for (i__ = kp1; i__ <= i__3; ++i__) { /* L40: */ a[i__ + j * a_dim1] += a[i__ + k * a_dim1] * t; } L45: /* L50: */ ; } /* L60: */ } L70: k = *n; if (a[*n + *n * a_dim1] == 0.) { goto L80; } return 0; L80: *ier = k; ip[*n] = 0; return 0; /* ----------------------- END OF SUBROUTINE DEC ------------------------- */ } /* dec_ */ /* Subroutine */ int sol_(integer *n, integer *ndim, doublereal *a, doublereal *b, integer *ip) { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2; /* Local variables */ static integer i__, k, m; static doublereal t; static integer kb, km1, nm1, kp1; /* VERSION REAL DOUBLE PRECISION */ /* ----------------------------------------------------------------------- */ /* SOLUTION OF LINEAR SYSTEM, A*X = B . */ /* INPUT.. */ /* N = ORDER OF MATRIX. */ /* NDIM = DECLARED DIMENSION OF ARRAY A . */ /* A = TRIANGULARIZED MATRIX OBTAINED FROM DEC. */ /* B = RIGHT HAND SIDE VECTOR. */ /* IP = PIVOT VECTOR OBTAINED FROM DEC. */ /* DO NOT USE IF DEC HAS SET IER .NE. 0. */ /* OUTPUT.. */ /* B = SOLUTION VECTOR, X . */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; --b; a_dim1 = *ndim; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ if (*n == 1) { goto L50; } nm1 = *n - 1; i__1 = nm1; for (k = 1; k <= i__1; ++k) { kp1 = k + 1; m = ip[k]; t = b[m]; b[m] = b[k]; b[k] = t; i__2 = *n; for (i__ = kp1; i__ <= i__2; ++i__) { /* L10: */ b[i__] += a[i__ + k * a_dim1] * t; } /* L20: */ } i__1 = nm1; for (kb = 1; kb <= i__1; ++kb) { km1 = *n - kb; k = km1 + 1; b[k] /= a[k + k * a_dim1]; t = -b[k]; i__2 = km1; for (i__ = 1; i__ <= i__2; ++i__) { /* L30: */ b[i__] += a[i__ + k * a_dim1] * t; } /* L40: */ } L50: b[1] /= a[a_dim1 + 1]; return 0; /* ----------------------- END OF SUBROUTINE SOL ------------------------- */ } /* sol_ */ /* Subroutine */ int dech_(integer *n, integer *ndim, doublereal *a, integer * lb, integer *ip, integer *ier) { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3; doublereal d__1, d__2; /* Local variables */ static integer i__, j, k, m; static doublereal t; static integer na, nm1, kp1; /* VERSION REAL DOUBLE PRECISION */ /* ----------------------------------------------------------------------- */ /* MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A HESSENBERG */ /* MATRIX WITH LOWER BANDWIDTH LB */ /* INPUT.. */ /* N = ORDER OF MATRIX A. */ /* NDIM = DECLARED DIMENSION OF ARRAY A . */ /* A = MATRIX TO BE TRIANGULARIZED. */ /* LB = LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED, LB.GE.1). */ /* OUTPUT.. */ /* A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U . */ /* A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. */ /* IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW. */ /* IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O . */ /* IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE */ /* SINGULAR AT STAGE K. */ /* USE SOLH TO OBTAIN SOLUTION OF LINEAR SYSTEM. */ /* DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N). */ /* IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. */ /* REFERENCE.. */ /* THIS IS A SLIGHT MODIFICATION OF */ /* C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, */ /* C.A.C.M. 15 (1972), P. 274. */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; a_dim1 = *ndim; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ *ier = 0; ip[*n] = 1; if (*n == 1) { goto L70; } nm1 = *n - 1; i__1 = nm1; for (k = 1; k <= i__1; ++k) { kp1 = k + 1; m = k; /* Computing MIN */ i__2 = *n, i__3 = *lb + k; na = min(i__2,i__3); i__2 = na; for (i__ = kp1; i__ <= i__2; ++i__) { if ((d__1 = a[i__ + k * a_dim1], fabs(d__1)) > (d__2 = a[m + k * a_dim1], fabs(d__2))) { m = i__; } /* L10: */ } ip[k] = m; t = a[m + k * a_dim1]; if (m == k) { goto L20; } ip[*n] = -ip[*n]; a[m + k * a_dim1] = a[k + k * a_dim1]; a[k + k * a_dim1] = t; L20: if (t == 0.) { goto L80; } t = 1. / t; i__2 = na; for (i__ = kp1; i__ <= i__2; ++i__) { /* L30: */ a[i__ + k * a_dim1] = -a[i__ + k * a_dim1] * t; } i__2 = *n; for (j = kp1; j <= i__2; ++j) { t = a[m + j * a_dim1]; a[m + j * a_dim1] = a[k + j * a_dim1]; a[k + j * a_dim1] = t; if (t == 0.) { goto L45; } i__3 = na; for (i__ = kp1; i__ <= i__3; ++i__) { /* L40: */ a[i__ + j * a_dim1] += a[i__ + k * a_dim1] * t; } L45: /* L50: */ ; } /* L60: */ } L70: k = *n; if (a[*n + *n * a_dim1] == 0.) { goto L80; } return 0; L80: *ier = k; ip[*n] = 0; return 0; /* ----------------------- END OF SUBROUTINE DECH ------------------------ */ } /* dech_ */ /* Subroutine */ int solh_(integer *n, integer *ndim, doublereal *a, integer * lb, doublereal *b, integer *ip) { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3; /* Local variables */ static integer i__, k, m; static doublereal t; static integer kb, na, km1, nm1, kp1; /* VERSION REAL DOUBLE PRECISION */ /* ----------------------------------------------------------------------- */ /* SOLUTION OF LINEAR SYSTEM, A*X = B . */ /* INPUT.. */ /* N = ORDER OF MATRIX A. */ /* NDIM = DECLARED DIMENSION OF ARRAY A . */ /* A = TRIANGULARIZED MATRIX OBTAINED FROM DECH. */ /* LB = LOWER BANDWIDTH OF A. */ /* B = RIGHT HAND SIDE VECTOR. */ /* IP = PIVOT VECTOR OBTAINED FROM DEC. */ /* DO NOT USE IF DECH HAS SET IER .NE. 0. */ /* OUTPUT.. */ /* B = SOLUTION VECTOR, X . */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; --b; a_dim1 = *ndim; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ if (*n == 1) { goto L50; } nm1 = *n - 1; i__1 = nm1; for (k = 1; k <= i__1; ++k) { kp1 = k + 1; m = ip[k]; t = b[m]; b[m] = b[k]; b[k] = t; /* Computing MIN */ i__2 = *n, i__3 = *lb + k; na = min(i__2,i__3); i__2 = na; for (i__ = kp1; i__ <= i__2; ++i__) { /* L10: */ b[i__] += a[i__ + k * a_dim1] * t; } /* L20: */ } i__1 = nm1; for (kb = 1; kb <= i__1; ++kb) { km1 = *n - kb; k = km1 + 1; b[k] /= a[k + k * a_dim1]; t = -b[k]; i__2 = km1; for (i__ = 1; i__ <= i__2; ++i__) { /* L30: */ b[i__] += a[i__ + k * a_dim1] * t; } /* L40: */ } L50: b[1] /= a[a_dim1 + 1]; return 0; /* ----------------------- END OF SUBROUTINE SOLH ------------------------ */ } /* solh_ */ /* Subroutine */ int decc_(integer *n, integer *ndim, doublereal *ar, doublereal *ai, integer *ip, integer *ier) { /* System generated locals */ integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3; doublereal d__1, d__2, d__3, d__4; /* Local variables */ static integer i__, j, k, m; static doublereal prodi, prodr, ti, tr; static integer nm1, kp1; static doublereal den; /* VERSION COMPLEX DOUBLE PRECISION */ /* ----------------------------------------------------------------------- */ /* MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION */ /* ------ MODIFICATION FOR COMPLEX MATRICES -------- */ /* INPUT.. */ /* N = ORDER OF MATRIX. */ /* NDIM = DECLARED DIMENSION OF ARRAYS AR AND AI . */ /* (AR, AI) = MATRIX TO BE TRIANGULARIZED. */ /* OUTPUT.. */ /* AR(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; REAL PART. */ /* AI(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; IMAGINARY PART. */ /* AR(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. */ /* REAL PART. */ /* AI(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. */ /* IMAGINARY PART. */ /* IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW. */ /* IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O . */ /* IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE */ /* SINGULAR AT STAGE K. */ /* USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. */ /* IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. */ /* REFERENCE.. */ /* C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, */ /* C.A.C.M. 15 (1972), P. 274. */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; ai_dim1 = *ndim; ai_offset = ai_dim1 + 1; ai -= ai_offset; ar_dim1 = *ndim; ar_offset = ar_dim1 + 1; ar -= ar_offset; /* Function Body */ *ier = 0; ip[*n] = 1; if (*n == 1) { goto L70; } nm1 = *n - 1; i__1 = nm1; for (k = 1; k <= i__1; ++k) { kp1 = k + 1; m = k; i__2 = *n; for (i__ = kp1; i__ <= i__2; ++i__) { if ((d__1 = ar[i__ + k * ar_dim1], fabs(d__1)) + (d__2 = ai[i__ + k * ai_dim1], fabs(d__2)) > (d__3 = ar[m + k * ar_dim1], fabs(d__3)) + (d__4 = ai[m + k * ai_dim1], fabs(d__4))) { m = i__; } /* L10: */ } ip[k] = m; tr = ar[m + k * ar_dim1]; ti = ai[m + k * ai_dim1]; if (m == k) { goto L20; } ip[*n] = -ip[*n]; ar[m + k * ar_dim1] = ar[k + k * ar_dim1]; ai[m + k * ai_dim1] = ai[k + k * ai_dim1]; ar[k + k * ar_dim1] = tr; ai[k + k * ai_dim1] = ti; L20: if (fabs(tr) + fabs(ti) == 0.) { goto L80; } den = tr * tr + ti * ti; tr /= den; ti = -ti / den; i__2 = *n; for (i__ = kp1; i__ <= i__2; ++i__) { prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; ar[i__ + k * ar_dim1] = -prodr; ai[i__ + k * ai_dim1] = -prodi; /* L30: */ } i__2 = *n; for (j = kp1; j <= i__2; ++j) { tr = ar[m + j * ar_dim1]; ti = ai[m + j * ai_dim1]; ar[m + j * ar_dim1] = ar[k + j * ar_dim1]; ai[m + j * ai_dim1] = ai[k + j * ai_dim1]; ar[k + j * ar_dim1] = tr; ai[k + j * ai_dim1] = ti; if (fabs(tr) + fabs(ti) == 0.) { goto L48; } if (ti == 0.) { i__3 = *n; for (i__ = kp1; i__ <= i__3; ++i__) { prodr = ar[i__ + k * ar_dim1] * tr; prodi = ai[i__ + k * ai_dim1] * tr; ar[i__ + j * ar_dim1] += prodr; ai[i__ + j * ai_dim1] += prodi; /* L40: */ } goto L48; } if (tr == 0.) { i__3 = *n; for (i__ = kp1; i__ <= i__3; ++i__) { prodr = -ai[i__ + k * ai_dim1] * ti; prodi = ar[i__ + k * ar_dim1] * ti; ar[i__ + j * ar_dim1] += prodr; ai[i__ + j * ai_dim1] += prodi; /* L45: */ } goto L48; } i__3 = *n; for (i__ = kp1; i__ <= i__3; ++i__) { prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; ar[i__ + j * ar_dim1] += prodr; ai[i__ + j * ai_dim1] += prodi; /* L47: */ } L48: /* L50: */ ; } /* L60: */ } L70: k = *n; if ((d__1 = ar[*n + *n * ar_dim1], fabs(d__1)) + (d__2 = ai[*n + *n * ai_dim1], fabs(d__2)) == 0.) { goto L80; } return 0; L80: *ier = k; ip[*n] = 0; return 0; /* ----------------------- END OF SUBROUTINE DECC ------------------------ */ } /* decc_ */ /* Subroutine */ int solc_(integer *n, integer *ndim, doublereal *ar, doublereal *ai, doublereal *br, doublereal *bi, integer *ip) { /* System generated locals */ integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2; /* Local variables */ static integer i__, k, m; static doublereal prodi, prodr; static integer kb; static doublereal ti, tr; static integer km1, nm1, kp1; static doublereal den; /* VERSION COMPLEX DOUBLE PRECISION */ /* ----------------------------------------------------------------------- */ /* SOLUTION OF LINEAR SYSTEM, A*X = B . */ /* INPUT.. */ /* N = ORDER OF MATRIX. */ /* NDIM = DECLARED DIMENSION OF ARRAYS AR AND AI. */ /* (AR,AI) = TRIANGULARIZED MATRIX OBTAINED FROM DEC. */ /* (BR,BI) = RIGHT HAND SIDE VECTOR. */ /* IP = PIVOT VECTOR OBTAINED FROM DEC. */ /* DO NOT USE IF DEC HAS SET IER .NE. 0. */ /* OUTPUT.. */ /* (BR,BI) = SOLUTION VECTOR, X . */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; --bi; --br; ai_dim1 = *ndim; ai_offset = ai_dim1 + 1; ai -= ai_offset; ar_dim1 = *ndim; ar_offset = ar_dim1 + 1; ar -= ar_offset; /* Function Body */ if (*n == 1) { goto L50; } nm1 = *n - 1; i__1 = nm1; for (k = 1; k <= i__1; ++k) { kp1 = k + 1; m = ip[k]; tr = br[m]; ti = bi[m]; br[m] = br[k]; bi[m] = bi[k]; br[k] = tr; bi[k] = ti; i__2 = *n; for (i__ = kp1; i__ <= i__2; ++i__) { prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; br[i__] += prodr; bi[i__] += prodi; /* L10: */ } /* L20: */ } i__1 = nm1; for (kb = 1; kb <= i__1; ++kb) { km1 = *n - kb; k = km1 + 1; den = ar[k + k * ar_dim1] * ar[k + k * ar_dim1] + ai[k + k * ai_dim1] * ai[k + k * ai_dim1]; prodr = br[k] * ar[k + k * ar_dim1] + bi[k] * ai[k + k * ai_dim1]; prodi = bi[k] * ar[k + k * ar_dim1] - br[k] * ai[k + k * ai_dim1]; br[k] = prodr / den; bi[k] = prodi / den; tr = -br[k]; ti = -bi[k]; i__2 = km1; for (i__ = 1; i__ <= i__2; ++i__) { prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; br[i__] += prodr; bi[i__] += prodi; /* L30: */ } /* L40: */ } L50: den = ar[ar_dim1 + 1] * ar[ar_dim1 + 1] + ai[ai_dim1 + 1] * ai[ai_dim1 + 1]; prodr = br[1] * ar[ar_dim1 + 1] + bi[1] * ai[ai_dim1 + 1]; prodi = bi[1] * ar[ar_dim1 + 1] - br[1] * ai[ai_dim1 + 1]; br[1] = prodr / den; bi[1] = prodi / den; return 0; /* ----------------------- END OF SUBROUTINE SOLC ------------------------ */ } /* solc_ */ /* Subroutine */ int dechc_(integer *n, integer *ndim, doublereal *ar, doublereal *ai, integer *lb, integer *ip, integer *ier) { /* System generated locals */ integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3; doublereal d__1, d__2, d__3, d__4; /* Local variables */ static integer i__, j, k, m; static doublereal prodi, prodr; static integer na; static doublereal ti, tr; static integer nm1, kp1; static doublereal den; /* VERSION COMPLEX DOUBLE PRECISION */ /* ----------------------------------------------------------------------- */ /* MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION */ /* ------ MODIFICATION FOR COMPLEX MATRICES -------- */ /* INPUT.. */ /* N = ORDER OF MATRIX. */ /* NDIM = DECLARED DIMENSION OF ARRAYS AR AND AI . */ /* (AR, AI) = MATRIX TO BE TRIANGULARIZED. */ /* OUTPUT.. */ /* AR(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; REAL PART. */ /* AI(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; IMAGINARY PART. */ /* AR(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. */ /* REAL PART. */ /* AI(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L. */ /* IMAGINARY PART. */ /* LB = LOWER BANDWIDTH OF A (DIAGONAL NOT COUNTED), LB.GE.1. */ /* IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW. */ /* IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O . */ /* IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE */ /* SINGULAR AT STAGE K. */ /* USE SOL TO OBTAIN SOLUTION OF LINEAR SYSTEM. */ /* IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO. */ /* REFERENCE.. */ /* C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, */ /* C.A.C.M. 15 (1972), P. 274. */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; ai_dim1 = *ndim; ai_offset = ai_dim1 + 1; ai -= ai_offset; ar_dim1 = *ndim; ar_offset = ar_dim1 + 1; ar -= ar_offset; /* Function Body */ *ier = 0; ip[*n] = 1; if (*lb == 0) { goto L70; } if (*n == 1) { goto L70; } nm1 = *n - 1; i__1 = nm1; for (k = 1; k <= i__1; ++k) { kp1 = k + 1; m = k; /* Computing MIN */ i__2 = *n, i__3 = *lb + k; na = min(i__2,i__3); i__2 = na; for (i__ = kp1; i__ <= i__2; ++i__) { if ((d__1 = ar[i__ + k * ar_dim1], fabs(d__1)) + (d__2 = ai[i__ + k * ai_dim1], fabs(d__2)) > (d__3 = ar[m + k * ar_dim1], fabs(d__3)) + (d__4 = ai[m + k * ai_dim1], fabs(d__4))) { m = i__; } /* L10: */ } ip[k] = m; tr = ar[m + k * ar_dim1]; ti = ai[m + k * ai_dim1]; if (m == k) { goto L20; } ip[*n] = -ip[*n]; ar[m + k * ar_dim1] = ar[k + k * ar_dim1]; ai[m + k * ai_dim1] = ai[k + k * ai_dim1]; ar[k + k * ar_dim1] = tr; ai[k + k * ai_dim1] = ti; L20: if (fabs(tr) + fabs(ti) == 0.) { goto L80; } den = tr * tr + ti * ti; tr /= den; ti = -ti / den; i__2 = na; for (i__ = kp1; i__ <= i__2; ++i__) { prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; ar[i__ + k * ar_dim1] = -prodr; ai[i__ + k * ai_dim1] = -prodi; /* L30: */ } i__2 = *n; for (j = kp1; j <= i__2; ++j) { tr = ar[m + j * ar_dim1]; ti = ai[m + j * ai_dim1]; ar[m + j * ar_dim1] = ar[k + j * ar_dim1]; ai[m + j * ai_dim1] = ai[k + j * ai_dim1]; ar[k + j * ar_dim1] = tr; ai[k + j * ai_dim1] = ti; if (fabs(tr) + fabs(ti) == 0.) { goto L48; } if (ti == 0.) { i__3 = na; for (i__ = kp1; i__ <= i__3; ++i__) { prodr = ar[i__ + k * ar_dim1] * tr; prodi = ai[i__ + k * ai_dim1] * tr; ar[i__ + j * ar_dim1] += prodr; ai[i__ + j * ai_dim1] += prodi; /* L40: */ } goto L48; } if (tr == 0.) { i__3 = na; for (i__ = kp1; i__ <= i__3; ++i__) { prodr = -ai[i__ + k * ai_dim1] * ti; prodi = ar[i__ + k * ar_dim1] * ti; ar[i__ + j * ar_dim1] += prodr; ai[i__ + j * ai_dim1] += prodi; /* L45: */ } goto L48; } i__3 = na; for (i__ = kp1; i__ <= i__3; ++i__) { prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; ar[i__ + j * ar_dim1] += prodr; ai[i__ + j * ai_dim1] += prodi; /* L47: */ } L48: /* L50: */ ; } /* L60: */ } L70: k = *n; if ((d__1 = ar[*n + *n * ar_dim1], fabs(d__1)) + (d__2 = ai[*n + *n * ai_dim1], fabs(d__2)) == 0.) { goto L80; } return 0; L80: *ier = k; ip[*n] = 0; return 0; /* ----------------------- END OF SUBROUTINE DECHC ----------------------- */ } /* dechc_ */ /* Subroutine */ int solhc_(integer *n, integer *ndim, doublereal *ar, doublereal *ai, integer *lb, doublereal *br, doublereal *bi, integer * ip) { /* System generated locals */ integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3, i__4; /* Local variables */ static integer i__, k, m; static doublereal prodi, prodr; static integer kb; static doublereal ti, tr; static integer km1, nm1, kp1; static doublereal den; /* VERSION COMPLEX DOUBLE PRECISION */ /* ----------------------------------------------------------------------- */ /* SOLUTION OF LINEAR SYSTEM, A*X = B . */ /* INPUT.. */ /* N = ORDER OF MATRIX. */ /* NDIM = DECLARED DIMENSION OF ARRAYS AR AND AI. */ /* (AR,AI) = TRIANGULARIZED MATRIX OBTAINED FROM DEC. */ /* (BR,BI) = RIGHT HAND SIDE VECTOR. */ /* LB = LOWER BANDWIDTH OF A. */ /* IP = PIVOT VECTOR OBTAINED FROM DEC. */ /* DO NOT USE IF DEC HAS SET IER .NE. 0. */ /* OUTPUT.. */ /* (BR,BI) = SOLUTION VECTOR, X . */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; --bi; --br; ai_dim1 = *ndim; ai_offset = ai_dim1 + 1; ai -= ai_offset; ar_dim1 = *ndim; ar_offset = ar_dim1 + 1; ar -= ar_offset; /* Function Body */ if (*n == 1) { goto L50; } nm1 = *n - 1; if (*lb == 0) { goto L25; } i__1 = nm1; for (k = 1; k <= i__1; ++k) { kp1 = k + 1; m = ip[k]; tr = br[m]; ti = bi[m]; br[m] = br[k]; bi[m] = bi[k]; br[k] = tr; bi[k] = ti; /* Computing MIN */ i__3 = *n, i__4 = *lb + k; i__2 = min(i__3,i__4); for (i__ = kp1; i__ <= i__2; ++i__) { prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; br[i__] += prodr; bi[i__] += prodi; /* L10: */ } /* L20: */ } L25: i__1 = nm1; for (kb = 1; kb <= i__1; ++kb) { km1 = *n - kb; k = km1 + 1; den = ar[k + k * ar_dim1] * ar[k + k * ar_dim1] + ai[k + k * ai_dim1] * ai[k + k * ai_dim1]; prodr = br[k] * ar[k + k * ar_dim1] + bi[k] * ai[k + k * ai_dim1]; prodi = bi[k] * ar[k + k * ar_dim1] - br[k] * ai[k + k * ai_dim1]; br[k] = prodr / den; bi[k] = prodi / den; tr = -br[k]; ti = -bi[k]; i__2 = km1; for (i__ = 1; i__ <= i__2; ++i__) { prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; br[i__] += prodr; bi[i__] += prodi; /* L30: */ } /* L40: */ } L50: den = ar[ar_dim1 + 1] * ar[ar_dim1 + 1] + ai[ai_dim1 + 1] * ai[ai_dim1 + 1]; prodr = br[1] * ar[ar_dim1 + 1] + bi[1] * ai[ai_dim1 + 1]; prodi = bi[1] * ar[ar_dim1 + 1] - br[1] * ai[ai_dim1 + 1]; br[1] = prodr / den; bi[1] = prodi / den; return 0; /* ----------------------- END OF SUBROUTINE SOLHC ----------------------- */ } /* solhc_ */ /* Subroutine */ int decb_(integer *n, integer *ndim, doublereal *a, integer * ml, integer *mu, integer *ip, integer *ier) { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3, i__4; doublereal d__1, d__2; /* Local variables */ static integer i__, j, k, m; static doublereal t; static integer md, jk, mm, ju, md1, nm1, kp1, mdl, ijk; /* ----------------------------------------------------------------------- */ /* MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED */ /* MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU */ /* INPUT.. */ /* N ORDER OF THE ORIGINAL MATRIX A. */ /* NDIM DECLARED DIMENSION OF ARRAY A. */ /* A CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS */ /* OF THE MATRIX ARE STORED IN THE COLUMNS OF A AND */ /* THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS */ /* ML+1 THROUGH 2*ML+MU+1 OF A. */ /* ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). */ /* MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). */ /* OUTPUT.. */ /* A AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND */ /* THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. */ /* IP INDEX VECTOR OF PIVOT INDICES. */ /* IP(N) (-1)**(NUMBER OF INTERCHANGES) OR O . */ /* IER = 0 IF MATRIX A IS NONSINGULAR, OR = K IF FOUND TO BE */ /* SINGULAR AT STAGE K. */ /* USE SOLB TO OBTAIN SOLUTION OF LINEAR SYSTEM. */ /* DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N) WITH MD=ML+MU+1. */ /* IF IP(N)=O, A IS SINGULAR, SOLB WILL DIVIDE BY ZERO. */ /* REFERENCE.. */ /* THIS IS A MODIFICATION OF */ /* C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, */ /* C.A.C.M. 15 (1972), P. 274. */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; a_dim1 = *ndim; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ *ier = 0; ip[*n] = 1; md = *ml + *mu + 1; md1 = md + 1; ju = 0; if (*ml == 0) { goto L70; } if (*n == 1) { goto L70; } if (*n < *mu + 2) { goto L7; } i__1 = *n; for (j = *mu + 2; j <= i__1; ++j) { i__2 = *ml; for (i__ = 1; i__ <= i__2; ++i__) { /* L5: */ a[i__ + j * a_dim1] = 0.; } } L7: nm1 = *n - 1; i__2 = nm1; for (k = 1; k <= i__2; ++k) { kp1 = k + 1; m = md; /* Computing MIN */ i__1 = *ml, i__3 = *n - k; mdl = min(i__1,i__3) + md; i__1 = mdl; for (i__ = md1; i__ <= i__1; ++i__) { if ((d__1 = a[i__ + k * a_dim1], fabs(d__1)) > (d__2 = a[m + k * a_dim1], fabs(d__2))) { m = i__; } /* L10: */ } ip[k] = m + k - md; t = a[m + k * a_dim1]; if (m == md) { goto L20; } ip[*n] = -ip[*n]; a[m + k * a_dim1] = a[md + k * a_dim1]; a[md + k * a_dim1] = t; L20: if (t == 0.) { goto L80; } t = 1. / t; i__1 = mdl; for (i__ = md1; i__ <= i__1; ++i__) { /* L30: */ a[i__ + k * a_dim1] = -a[i__ + k * a_dim1] * t; } /* Computing MIN */ /* Computing MAX */ i__3 = ju, i__4 = *mu + ip[k]; i__1 = max(i__3,i__4); ju = min(i__1,*n); mm = md; if (ju < kp1) { goto L55; } i__1 = ju; for (j = kp1; j <= i__1; ++j) { --m; --mm; t = a[m + j * a_dim1]; if (m == mm) { goto L35; } a[m + j * a_dim1] = a[mm + j * a_dim1]; a[mm + j * a_dim1] = t; L35: if (t == 0.) { goto L45; } jk = j - k; i__3 = mdl; for (i__ = md1; i__ <= i__3; ++i__) { ijk = i__ - jk; /* L40: */ a[ijk + j * a_dim1] += a[i__ + k * a_dim1] * t; } L45: /* L50: */ ; } L55: /* L60: */ ; } L70: k = *n; if (a[md + *n * a_dim1] == 0.) { goto L80; } return 0; L80: *ier = k; ip[*n] = 0; return 0; /* ----------------------- END OF SUBROUTINE DECB ------------------------ */ } /* decb_ */ /* Subroutine */ int solb_(integer *n, integer *ndim, doublereal *a, integer * ml, integer *mu, doublereal *b, integer *ip) { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3; /* Local variables */ static integer i__, k, m; static doublereal t; static integer kb, md, lm, md1, nm1, imd, kmd, mdl, mdm; /* ----------------------------------------------------------------------- */ /* SOLUTION OF LINEAR SYSTEM, A*X = B . */ /* INPUT.. */ /* N ORDER OF MATRIX A. */ /* NDIM DECLARED DIMENSION OF ARRAY A . */ /* A TRIANGULARIZED MATRIX OBTAINED FROM DECB. */ /* ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). */ /* MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). */ /* B RIGHT HAND SIDE VECTOR. */ /* IP PIVOT VECTOR OBTAINED FROM DECB. */ /* DO NOT USE IF DECB HAS SET IER .NE. 0. */ /* OUTPUT.. */ /* B SOLUTION VECTOR, X . */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; --b; a_dim1 = *ndim; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ md = *ml + *mu + 1; md1 = md + 1; mdm = md - 1; nm1 = *n - 1; if (*ml == 0) { goto L25; } if (*n == 1) { goto L50; } i__1 = nm1; for (k = 1; k <= i__1; ++k) { m = ip[k]; t = b[m]; b[m] = b[k]; b[k] = t; /* Computing MIN */ i__2 = *ml, i__3 = *n - k; mdl = min(i__2,i__3) + md; i__2 = mdl; for (i__ = md1; i__ <= i__2; ++i__) { imd = i__ + k - md; /* L10: */ b[imd] += a[i__ + k * a_dim1] * t; } /* L20: */ } L25: i__1 = nm1; for (kb = 1; kb <= i__1; ++kb) { k = *n + 1 - kb; b[k] /= a[md + k * a_dim1]; t = -b[k]; kmd = md - k; /* Computing MAX */ i__2 = 1, i__3 = kmd + 1; lm = max(i__2,i__3); i__2 = mdm; for (i__ = lm; i__ <= i__2; ++i__) { imd = i__ - kmd; /* L30: */ b[imd] += a[i__ + k * a_dim1] * t; } /* L40: */ } L50: b[1] /= a[md + a_dim1]; return 0; /* ----------------------- END OF SUBROUTINE SOLB ------------------------ */ } /* solb_ */ /* Subroutine */ int decbc_(integer *n, integer *ndim, doublereal *ar, doublereal *ai, integer *ml, integer *mu, integer *ip, integer *ier) { /* System generated locals */ integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3, i__4; doublereal d__1, d__2, d__3, d__4; /* Local variables */ static integer i__, j, k, m; static doublereal prodi, prodr; static integer md, jk, mm; static doublereal ti; static integer ju; static doublereal tr; static integer md1, nm1, kp1; static doublereal den; static integer mdl, ijk; /* ----------------------------------------------------------------------- */ /* MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED COMPLEX */ /* MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU */ /* INPUT.. */ /* N ORDER OF THE ORIGINAL MATRIX A. */ /* NDIM DECLARED DIMENSION OF ARRAY A. */ /* AR, AI CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS */ /* OF THE MATRIX ARE STORED IN THE COLUMNS OF AR (REAL */ /* PART) AND AI (IMAGINARY PART) AND */ /* THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS */ /* ML+1 THROUGH 2*ML+MU+1 OF AR AND AI. */ /* ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). */ /* MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). */ /* OUTPUT.. */ /* AR, AI AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND */ /* THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. */ /* IP INDEX VECTOR OF PIVOT INDICES. */ /* IP(N) (-1)**(NUMBER OF INTERCHANGES) OR O . */ /* IER = 0 IF MATRIX A IS NONSINGULAR, OR = K IF FOUND TO BE */ /* SINGULAR AT STAGE K. */ /* USE SOLBC TO OBTAIN SOLUTION OF LINEAR SYSTEM. */ /* DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N) WITH MD=ML+MU+1. */ /* IF IP(N)=O, A IS SINGULAR, SOLBC WILL DIVIDE BY ZERO. */ /* REFERENCE.. */ /* THIS IS A MODIFICATION OF */ /* C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER, */ /* C.A.C.M. 15 (1972), P. 274. */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; ai_dim1 = *ndim; ai_offset = ai_dim1 + 1; ai -= ai_offset; ar_dim1 = *ndim; ar_offset = ar_dim1 + 1; ar -= ar_offset; /* Function Body */ *ier = 0; ip[*n] = 1; md = *ml + *mu + 1; md1 = md + 1; ju = 0; if (*ml == 0) { goto L70; } if (*n == 1) { goto L70; } if (*n < *mu + 2) { goto L7; } i__1 = *n; for (j = *mu + 2; j <= i__1; ++j) { i__2 = *ml; for (i__ = 1; i__ <= i__2; ++i__) { ar[i__ + j * ar_dim1] = 0.; ai[i__ + j * ai_dim1] = 0.; /* L5: */ } } L7: nm1 = *n - 1; i__2 = nm1; for (k = 1; k <= i__2; ++k) { kp1 = k + 1; m = md; /* Computing MIN */ i__1 = *ml, i__3 = *n - k; mdl = min(i__1,i__3) + md; i__1 = mdl; for (i__ = md1; i__ <= i__1; ++i__) { if ((d__1 = ar[i__ + k * ar_dim1], fabs(d__1)) + (d__2 = ai[i__ + k * ai_dim1], fabs(d__2)) > (d__3 = ar[m + k * ar_dim1], fabs(d__3)) + (d__4 = ai[m + k * ai_dim1], fabs(d__4))) { m = i__; } /* L10: */ } ip[k] = m + k - md; tr = ar[m + k * ar_dim1]; ti = ai[m + k * ai_dim1]; if (m == md) { goto L20; } ip[*n] = -ip[*n]; ar[m + k * ar_dim1] = ar[md + k * ar_dim1]; ai[m + k * ai_dim1] = ai[md + k * ai_dim1]; ar[md + k * ar_dim1] = tr; ai[md + k * ai_dim1] = ti; L20: if (fabs(tr) + fabs(ti) == 0.) { goto L80; } den = tr * tr + ti * ti; tr /= den; ti = -ti / den; i__1 = mdl; for (i__ = md1; i__ <= i__1; ++i__) { prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; ar[i__ + k * ar_dim1] = -prodr; ai[i__ + k * ai_dim1] = -prodi; /* L30: */ } /* Computing MIN */ /* Computing MAX */ i__3 = ju, i__4 = *mu + ip[k]; i__1 = max(i__3,i__4); ju = min(i__1,*n); mm = md; if (ju < kp1) { goto L55; } i__1 = ju; for (j = kp1; j <= i__1; ++j) { --m; --mm; tr = ar[m + j * ar_dim1]; ti = ai[m + j * ai_dim1]; if (m == mm) { goto L35; } ar[m + j * ar_dim1] = ar[mm + j * ar_dim1]; ai[m + j * ai_dim1] = ai[mm + j * ai_dim1]; ar[mm + j * ar_dim1] = tr; ai[mm + j * ai_dim1] = ti; L35: if (fabs(tr) + fabs(ti) == 0.) { goto L48; } jk = j - k; if (ti == 0.) { i__3 = mdl; for (i__ = md1; i__ <= i__3; ++i__) { ijk = i__ - jk; prodr = ar[i__ + k * ar_dim1] * tr; prodi = ai[i__ + k * ai_dim1] * tr; ar[ijk + j * ar_dim1] += prodr; ai[ijk + j * ai_dim1] += prodi; /* L40: */ } goto L48; } if (tr == 0.) { i__3 = mdl; for (i__ = md1; i__ <= i__3; ++i__) { ijk = i__ - jk; prodr = -ai[i__ + k * ai_dim1] * ti; prodi = ar[i__ + k * ar_dim1] * ti; ar[ijk + j * ar_dim1] += prodr; ai[ijk + j * ai_dim1] += prodi; /* L45: */ } goto L48; } i__3 = mdl; for (i__ = md1; i__ <= i__3; ++i__) { ijk = i__ - jk; prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; ar[ijk + j * ar_dim1] += prodr; ai[ijk + j * ai_dim1] += prodi; /* L47: */ } L48: /* L50: */ ; } L55: /* L60: */ ; } L70: k = *n; if ((d__1 = ar[md + *n * ar_dim1], fabs(d__1)) + (d__2 = ai[md + *n * ai_dim1], fabs(d__2)) == 0.) { goto L80; } return 0; L80: *ier = k; ip[*n] = 0; return 0; /*----------------------- END OF SUBROUTINE DECBC ----------------------- -*/ } /* decbc_ */ /* Subroutine */ int solbc_(integer *n, integer *ndim, doublereal *ar, doublereal *ai, integer *ml, integer *mu, doublereal *br, doublereal * bi, integer *ip) { /* System generated locals */ integer ar_dim1, ar_offset, ai_dim1, ai_offset, i__1, i__2, i__3; /* Local variables */ static integer i__, k, m; static doublereal prodi, prodr; static integer kb, md, lm; static doublereal ti, tr; static integer md1, nm1; static doublereal den; static integer imd, kmd, mdl, mdm; /* ----------------------------------------------------------------------- */ /* SOLUTION OF LINEAR SYSTEM, A*X = B , */ /* VERSION BANDED AND COMPLEX-DOUBLE PRECISION. */ /* INPUT.. */ /* N ORDER OF MATRIX A. */ /* NDIM DECLARED DIMENSION OF ARRAY A . */ /* AR, AI TRIANGULARIZED MATRIX OBTAINED FROM DECB (REAL AND IMAG. PART) .*/ /* ML LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). */ /* MU UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED). */ /* BR, BI RIGHT HAND SIDE VECTOR (REAL AND IMAG. PART). */ /* IP PIVOT VECTOR OBTAINED FROM DECBC. */ /* DO NOT USE IF DECB HAS SET IER .NE. 0. */ /* OUTPUT.. */ /* BR, BI SOLUTION VECTOR, X (REAL AND IMAG. PART). */ /* ----------------------------------------------------------------------- */ /* Parameter adjustments */ --ip; --bi; --br; ai_dim1 = *ndim; ai_offset = ai_dim1 + 1; ai -= ai_offset; ar_dim1 = *ndim; ar_offset = ar_dim1 + 1; ar -= ar_offset; /* Function Body */ md = *ml + *mu + 1; md1 = md + 1; mdm = md - 1; nm1 = *n - 1; if (*ml == 0) { goto L25; } if (*n == 1) { goto L50; } i__1 = nm1; for (k = 1; k <= i__1; ++k) { m = ip[k]; tr = br[m]; ti = bi[m]; br[m] = br[k]; bi[m] = bi[k]; br[k] = tr; bi[k] = ti; /* Computing MIN */ i__2 = *ml, i__3 = *n - k; mdl = min(i__2,i__3) + md; i__2 = mdl; for (i__ = md1; i__ <= i__2; ++i__) { imd = i__ + k - md; prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; br[imd] += prodr; bi[imd] += prodi; /* L10: */ } /* L20: */ } L25: i__1 = nm1; for (kb = 1; kb <= i__1; ++kb) { k = *n + 1 - kb; den = ar[md + k * ar_dim1] * ar[md + k * ar_dim1] + ai[md + k * ai_dim1] * ai[md + k * ai_dim1]; prodr = br[k] * ar[md + k * ar_dim1] + bi[k] * ai[md + k * ai_dim1]; prodi = bi[k] * ar[md + k * ar_dim1] - br[k] * ai[md + k * ai_dim1]; br[k] = prodr / den; bi[k] = prodi / den; tr = -br[k]; ti = -bi[k]; kmd = md - k; /* Computing MAX */ i__2 = 1, i__3 = kmd + 1; lm = max(i__2,i__3); i__2 = mdm; for (i__ = lm; i__ <= i__2; ++i__) { imd = i__ - kmd; prodr = ar[i__ + k * ar_dim1] * tr - ai[i__ + k * ai_dim1] * ti; prodi = ai[i__ + k * ai_dim1] * tr + ar[i__ + k * ar_dim1] * ti; br[imd] += prodr; bi[imd] += prodi; /* L30: */ } /* L40: */ } den = ar[md + ar_dim1] * ar[md + ar_dim1] + ai[md + ai_dim1] * ai[md + ai_dim1]; prodr = br[1] * ar[md + ar_dim1] + bi[1] * ai[md + ai_dim1]; prodi = bi[1] * ar[md + ar_dim1] - br[1] * ai[md + ai_dim1]; br[1] = prodr / den; bi[1] = prodi / den; L50: return 0; /*----------------------- END OF SUBROUTINE SOLBC ----------------------- -*/ } /* solbc_ */ /* Subroutine */ int elmhes_(integer *nm, integer *n, integer *low, integer * igh, doublereal *a, integer *int__) { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3; doublereal d__1; /* Local variables */ static integer i__, j, m; static doublereal x, y; static integer la, mm1, kp1, mp1; /* this subroutine is a translation of the algol procedure elmhes, */ /* num. math. 12, 349-368(1968) by martin and wilkinson. */ /* handbook for auto. comp., vol.ii-linear algebra, 339-358(1971). */ /* given a real general matrix, this subroutine */ /* reduces a submatrix situated in rows and columns */ /* low through igh to upper hessenberg form by */ /* stabilized elementary similarity transformations. */ /* on input: */ /* nm must be set to the row dimension of two-dimensional */ /* array parameters as declared in the calling program */ /* dimension statement; */ /* n is the order of the matrix; */ /* low and igh are integers determined by the balancing */ /* subroutine balanc. if balanc has not been used, */ /* set low=1, igh=n; */ /* a contains the input matrix. */ /* on output: */ /* a contains the hessenberg matrix. the multipliers */ /* which were used in the reduction are stored in the */ /* remaining triangle under the hessenberg matrix; */ /* int contains information on the rows and columns */ /* interchanged in the reduction. */ /* only elements low through igh are used. */ /* questions and comments should be directed to b. s. garbow, */ /* applied mathematics division, argonne national laboratory */ /* ------------------------------------------------------------------ */ /* Parameter adjustments */ a_dim1 = *nm; a_offset = a_dim1 + 1; a -= a_offset; --int__; /* Function Body */ la = *igh - 1; kp1 = *low + 1; if (la < kp1) { goto L200; } i__1 = la; for (m = kp1; m <= i__1; ++m) { mm1 = m - 1; x = 0.; i__ = m; i__2 = *igh; for (j = m; j <= i__2; ++j) { if ((d__1 = a[j + mm1 * a_dim1], fabs(d__1)) <= fabs(x)) { goto L100; } x = a[j + mm1 * a_dim1]; i__ = j; L100: ; } int__[m] = i__; if (i__ == m) { goto L130; } /* :::::::::: interchange rows and columns of a :::::::::: */ i__2 = *n; for (j = mm1; j <= i__2; ++j) { y = a[i__ + j * a_dim1]; a[i__ + j * a_dim1] = a[m + j * a_dim1]; a[m + j * a_dim1] = y; /* L110: */ } i__2 = *igh; for (j = 1; j <= i__2; ++j) { y = a[j + i__ * a_dim1]; a[j + i__ * a_dim1] = a[j + m * a_dim1]; a[j + m * a_dim1] = y; /* L120: */ } /* :::::::::: end interchange :::::::::: */ L130: if (x == 0.) { goto L180; } mp1 = m + 1; i__2 = *igh; for (i__ = mp1; i__ <= i__2; ++i__) { y = a[i__ + mm1 * a_dim1]; if (y == 0.) { goto L160; } y /= x; a[i__ + mm1 * a_dim1] = y; i__3 = *n; for (j = m; j <= i__3; ++j) { /* L140: */ a[i__ + j * a_dim1] -= y * a[m + j * a_dim1]; } i__3 = *igh; for (j = 1; j <= i__3; ++j) { /* L150: */ a[j + m * a_dim1] += y * a[j + i__ * a_dim1]; } L160: ; } L180: ; } L200: return 0; /* :::::::::: last card of elmhes :::::::::: */ } /* elmhes_ */ /* ****************************************** */ /* VERSION OF SEPTEMBER 18, 1995 */ /* ****************************************** */ /* Subroutine */ int decomr_(integer *n, doublereal *fjac, integer *ldjac, doublereal *fmas, integer *ldmas, integer *mlmas, integer *mumas, integer *m1, integer *m2, integer *nm1, doublereal *fac1, doublereal * e1, integer *lde1, integer *ip1, integer *ier, integer *ijob, logical *calhes, integer *iphes) { /* System generated locals */ integer fjac_dim1, fjac_offset, fmas_dim1, fmas_offset, e1_dim1, e1_offset, i__1, i__2, i__3, i__4, i__5, i__6; /* Local variables */ extern /* Subroutine */ int dech_(integer *, integer *, doublereal *, integer *, integer *, integer *); static integer i__, j, k, j1, ib, mm; extern /* Subroutine */ int elmhes_(integer *, integer *, integer *, integer *, doublereal *, integer *); static integer jm1; extern /* Subroutine */ int dec_(integer *, integer *, doublereal *, integer *, integer *); static doublereal sum; extern /* Subroutine */ int decb_(integer *, integer *, doublereal *, integer *, integer *, integer *, integer *); /* Parameter adjustments */ --iphes; fjac_dim1 = *ldjac; fjac_offset = fjac_dim1 + 1; fjac -= fjac_offset; --ip1; fmas_dim1 = *ldmas; fmas_offset = fmas_dim1 + 1; fmas -= fmas_offset; e1_dim1 = *lde1; e1_offset = e1_dim1 + 1; e1 -= e1_offset; /* Function Body */ switch (*ijob) { case 1: goto L1; case 2: goto L2; case 3: goto L3; case 4: goto L4; case 5: goto L5; case 6: goto L6; case 7: goto L7; case 8: goto L55; case 9: goto L55; case 10: goto L55; case 11: goto L11; case 12: goto L12; case 13: goto L13; case 14: goto L14; case 15: goto L15; } /* ----------------------------------------------------------- */ L1: /* --- B=IDENTITY, JACOBIAN A FULL MATRIX */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { e1[i__ + j * e1_dim1] = -fjac[i__ + j * fjac_dim1]; } e1[j + j * e1_dim1] += *fac1; } dec_(n, lde1, &e1[e1_offset], &ip1[1], ier); return 0; /* ----------------------------------------------------------- */ L11: /* --- B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *nm1; for (j = 1; j <= i__1; ++j) { jm1 = j + *m1; i__2 = *nm1; for (i__ = 1; i__ <= i__2; ++i__) { e1[i__ + j * e1_dim1] = -fjac[i__ + jm1 * fjac_dim1]; } e1[j + j * e1_dim1] += *fac1; } L45: mm = *m1 / *m2; i__1 = *m2; for (j = 1; j <= i__1; ++j) { i__2 = *nm1; for (i__ = 1; i__ <= i__2; ++i__) { sum = 0.; i__3 = mm - 1; for (k = 0; k <= i__3; ++k) { sum = (sum + fjac[i__ + (j + k * *m2) * fjac_dim1]) / *fac1; } e1[i__ + j * e1_dim1] -= sum; } } dec_(nm1, lde1, &e1[e1_offset], &ip1[1], ier); return 0; /* ----------------------------------------------------------- */ L2: /* --- B=IDENTITY, JACOBIAN A BANDED MATRIX */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = linal_1.mbjac; for (i__ = 1; i__ <= i__2; ++i__) { e1[i__ + linal_1.mle + j * e1_dim1] = -fjac[i__ + j * fjac_dim1]; } e1[linal_1.mdiag + j * e1_dim1] += *fac1; } decb_(n, lde1, &e1[e1_offset], &linal_1.mle, &linal_1.mue, &ip1[1], ier); return 0; /* ----------------------------------------------------------- */ L12: /* --- B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER */ i__1 = *nm1; for (j = 1; j <= i__1; ++j) { jm1 = j + *m1; i__2 = linal_1.mbjac; for (i__ = 1; i__ <= i__2; ++i__) { e1[i__ + linal_1.mle + j * e1_dim1] = -fjac[i__ + jm1 * fjac_dim1] ; } e1[linal_1.mdiag + j * e1_dim1] += *fac1; } L46: mm = *m1 / *m2; i__1 = *m2; for (j = 1; j <= i__1; ++j) { i__2 = linal_1.mbjac; for (i__ = 1; i__ <= i__2; ++i__) { sum = 0.; i__3 = mm - 1; for (k = 0; k <= i__3; ++k) { sum = (sum + fjac[i__ + (j + k * *m2) * fjac_dim1]) / *fac1; } e1[i__ + linal_1.mle + j * e1_dim1] -= sum; } } decb_(nm1, lde1, &e1[e1_offset], &linal_1.mle, &linal_1.mue, &ip1[1], ier) ; return 0; /* ----------------------------------------------------------- */ L3: /* --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { e1[i__ + j * e1_dim1] = -fjac[i__ + j * fjac_dim1]; } /* Computing MAX */ i__2 = 1, i__3 = j - *mumas; /* Computing MIN */ i__5 = *n, i__6 = j + *mlmas; i__4 = min(i__5,i__6); for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) { e1[i__ + j * e1_dim1] += *fac1 * fmas[i__ - j + linal_1.mbdiag + j * fmas_dim1]; } } dec_(n, lde1, &e1[e1_offset], &ip1[1], ier); return 0; /* ----------------------------------------------------------- */ L13: /* --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *nm1; for (j = 1; j <= i__1; ++j) { jm1 = j + *m1; i__4 = *nm1; for (i__ = 1; i__ <= i__4; ++i__) { e1[i__ + j * e1_dim1] = -fjac[i__ + jm1 * fjac_dim1]; } /* Computing MAX */ i__4 = 1, i__2 = j - *mumas; /* Computing MIN */ i__5 = *nm1, i__6 = j + *mlmas; i__3 = min(i__5,i__6); for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) { e1[i__ + j * e1_dim1] += *fac1 * fmas[i__ - j + linal_1.mbdiag + j * fmas_dim1]; } } goto L45; /* ----------------------------------------------------------- */ L4: /* --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__3 = linal_1.mbjac; for (i__ = 1; i__ <= i__3; ++i__) { e1[i__ + linal_1.mle + j * e1_dim1] = -fjac[i__ + j * fjac_dim1]; } i__3 = linal_1.mbb; for (i__ = 1; i__ <= i__3; ++i__) { ib = i__ + linal_1.mdiff; e1[ib + j * e1_dim1] += *fac1 * fmas[i__ + j * fmas_dim1]; } } decb_(n, lde1, &e1[e1_offset], &linal_1.mle, &linal_1.mue, &ip1[1], ier); return 0; /* ----------------------------------------------------------- */ L14: /* --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER */ i__1 = *nm1; for (j = 1; j <= i__1; ++j) { jm1 = j + *m1; i__3 = linal_1.mbjac; for (i__ = 1; i__ <= i__3; ++i__) { e1[i__ + linal_1.mle + j * e1_dim1] = -fjac[i__ + jm1 * fjac_dim1] ; } i__3 = linal_1.mbb; for (i__ = 1; i__ <= i__3; ++i__) { ib = i__ + linal_1.mdiff; e1[ib + j * e1_dim1] += *fac1 * fmas[i__ + j * fmas_dim1]; } } goto L46; /* ----------------------------------------------------------- */ L5: /* --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__3 = *n; for (i__ = 1; i__ <= i__3; ++i__) { e1[i__ + j * e1_dim1] = fmas[i__ + j * fmas_dim1] * *fac1 - fjac[ i__ + j * fjac_dim1]; } } dec_(n, lde1, &e1[e1_offset], &ip1[1], ier); return 0; /* ----------------------------------------------------------- */ L15: /* --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *nm1; for (j = 1; j <= i__1; ++j) { jm1 = j + *m1; i__3 = *nm1; for (i__ = 1; i__ <= i__3; ++i__) { e1[i__ + j * e1_dim1] = fmas[i__ + j * fmas_dim1] * *fac1 - fjac[ i__ + jm1 * fjac_dim1]; } } goto L45; /* ----------------------------------------------------------- */ L6: /* --- B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX */ /* --- THIS OPTION IS NOT PROVIDED */ return 0; /* ----------------------------------------------------------- */ L7: /* --- B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION */ if (*calhes) { elmhes_(ldjac, n, &c__1, n, &fjac[fjac_offset], &iphes[1]); } *calhes = FALSE_; i__1 = *n - 1; for (j = 1; j <= i__1; ++j) { j1 = j + 1; e1[j1 + j * e1_dim1] = -fjac[j1 + j * fjac_dim1]; } i__1 = *n; for (j = 1; j <= i__1; ++j) { i__3 = j; for (i__ = 1; i__ <= i__3; ++i__) { e1[i__ + j * e1_dim1] = -fjac[i__ + j * fjac_dim1]; } e1[j + j * e1_dim1] += *fac1; } dech_(n, lde1, &e1[e1_offset], &c__1, &ip1[1], ier); return 0; /* ----------------------------------------------------------- */ L55: return 0; } /* decomr_ */ /* END OF SUBROUTINE DECOMR */ /* *********************************************************** */ /* Subroutine */ int decomc_(integer *n, doublereal *fjac, integer *ldjac, doublereal *fmas, integer *ldmas, integer *mlmas, integer *mumas, integer *m1, integer *m2, integer *nm1, doublereal *alphn, doublereal *betan, doublereal *e2r, doublereal *e2i, integer *lde1, integer *ip2, integer *ier, integer *ijob) { /* System generated locals */ integer fjac_dim1, fjac_offset, fmas_dim1, fmas_offset, e2r_dim1, e2r_offset, e2i_dim1, e2i_offset, i__1, i__2, i__3, i__4, i__5, i__6; doublereal d__1, d__2; /* Local variables */ static doublereal ffma, abno; static integer imle; static doublereal sumi, sumr, sums; extern /* Subroutine */ int decbc_(integer *, integer *, doublereal *, doublereal *, integer *, integer *, integer *, integer *), dechc_( integer *, integer *, doublereal *, doublereal *, integer *, integer *, integer *); static integer i__, j, k, j1; static doublereal bb; static integer ib, mm, jm1; static doublereal bet, alp; extern /* Subroutine */ int decc_(integer *, integer *, doublereal *, doublereal *, integer *, integer *); /* Parameter adjustments */ fjac_dim1 = *ldjac; fjac_offset = fjac_dim1 + 1; fjac -= fjac_offset; --ip2; fmas_dim1 = *ldmas; fmas_offset = fmas_dim1 + 1; fmas -= fmas_offset; e2i_dim1 = *lde1; e2i_offset = e2i_dim1 + 1; e2i -= e2i_offset; e2r_dim1 = *lde1; e2r_offset = e2r_dim1 + 1; e2r -= e2r_offset; /* Function Body */ switch (*ijob) { case 1: goto L1; case 2: goto L2; case 3: goto L3; case 4: goto L4; case 5: goto L5; case 6: goto L6; case 7: goto L7; case 8: goto L55; case 9: goto L55; case 10: goto L55; case 11: goto L11; case 12: goto L12; case 13: goto L13; case 14: goto L14; case 15: goto L15; } /* ----------------------------------------------------------- */ L1: /* --- B=IDENTITY, JACOBIAN A FULL MATRIX */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { e2r[i__ + j * e2r_dim1] = -fjac[i__ + j * fjac_dim1]; e2i[i__ + j * e2i_dim1] = 0.; } e2r[j + j * e2r_dim1] += *alphn; e2i[j + j * e2i_dim1] = *betan; } decc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &ip2[1], ier); return 0; /* ----------------------------------------------------------- */ L11: /* --- B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *nm1; for (j = 1; j <= i__1; ++j) { jm1 = j + *m1; i__2 = *nm1; for (i__ = 1; i__ <= i__2; ++i__) { e2r[i__ + j * e2r_dim1] = -fjac[i__ + jm1 * fjac_dim1]; e2i[i__ + j * e2i_dim1] = 0.; } e2r[j + j * e2r_dim1] += *alphn; e2i[j + j * e2i_dim1] = *betan; } L45: mm = *m1 / *m2; /* Computing 2nd power */ d__1 = *alphn; /* Computing 2nd power */ d__2 = *betan; abno = d__1 * d__1 + d__2 * d__2; alp = *alphn / abno; bet = *betan / abno; i__1 = *m2; for (j = 1; j <= i__1; ++j) { i__2 = *nm1; for (i__ = 1; i__ <= i__2; ++i__) { sumr = 0.; sumi = 0.; i__3 = mm - 1; for (k = 0; k <= i__3; ++k) { sums = sumr + fjac[i__ + (j + k * *m2) * fjac_dim1]; sumr = sums * alp + sumi * bet; sumi = sumi * alp - sums * bet; } e2r[i__ + j * e2r_dim1] -= sumr; e2i[i__ + j * e2i_dim1] -= sumi; } } decc_(nm1, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &ip2[1], ier); return 0; /* ----------------------------------------------------------- */ L2: /* --- B=IDENTITY, JACOBIAN A BANDED MATRIX */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = linal_1.mbjac; for (i__ = 1; i__ <= i__2; ++i__) { imle = i__ + linal_1.mle; e2r[imle + j * e2r_dim1] = -fjac[i__ + j * fjac_dim1]; e2i[imle + j * e2i_dim1] = 0.; } e2r[linal_1.mdiag + j * e2r_dim1] += *alphn; e2i[linal_1.mdiag + j * e2i_dim1] = *betan; } decbc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &linal_1.mle, & linal_1.mue, &ip2[1], ier); return 0; /* ----------------------------------------------------------- */ L12: /* --- B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER */ i__1 = *nm1; for (j = 1; j <= i__1; ++j) { jm1 = j + *m1; i__2 = linal_1.mbjac; for (i__ = 1; i__ <= i__2; ++i__) { e2r[i__ + linal_1.mle + j * e2r_dim1] = -fjac[i__ + jm1 * fjac_dim1]; e2i[i__ + linal_1.mle + j * e2i_dim1] = 0.; } e2r[linal_1.mdiag + j * e2r_dim1] += *alphn; e2i[linal_1.mdiag + j * e2i_dim1] += *betan; } L46: mm = *m1 / *m2; /* Computing 2nd power */ d__1 = *alphn; /* Computing 2nd power */ d__2 = *betan; abno = d__1 * d__1 + d__2 * d__2; alp = *alphn / abno; bet = *betan / abno; i__1 = *m2; for (j = 1; j <= i__1; ++j) { i__2 = linal_1.mbjac; for (i__ = 1; i__ <= i__2; ++i__) { sumr = 0.; sumi = 0.; i__3 = mm - 1; for (k = 0; k <= i__3; ++k) { sums = sumr + fjac[i__ + (j + k * *m2) * fjac_dim1]; sumr = sums * alp + sumi * bet; sumi = sumi * alp - sums * bet; } imle = i__ + linal_1.mle; e2r[imle + j * e2r_dim1] -= sumr; e2i[imle + j * e2i_dim1] -= sumi; } } decbc_(nm1, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &linal_1.mle, & linal_1.mue, &ip2[1], ier); return 0; /* ----------------------------------------------------------- */ L3: /* --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { e2r[i__ + j * e2r_dim1] = -fjac[i__ + j * fjac_dim1]; e2i[i__ + j * e2i_dim1] = 0.; } } i__1 = *n; for (j = 1; j <= i__1; ++j) { /* Computing MAX */ i__2 = 1, i__3 = j - *mumas; /* Computing MIN */ i__5 = *n, i__6 = j + *mlmas; i__4 = min(i__5,i__6); for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) { bb = fmas[i__ - j + linal_1.mbdiag + j * fmas_dim1]; e2r[i__ + j * e2r_dim1] += *alphn * bb; e2i[i__ + j * e2i_dim1] = *betan * bb; } } decc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &ip2[1], ier); return 0; /* ----------------------------------------------------------- */ L13: /* --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *nm1; for (j = 1; j <= i__1; ++j) { jm1 = j + *m1; i__4 = *nm1; for (i__ = 1; i__ <= i__4; ++i__) { e2r[i__ + j * e2r_dim1] = -fjac[i__ + jm1 * fjac_dim1]; e2i[i__ + j * e2i_dim1] = 0.; } /* Computing MAX */ i__4 = 1, i__2 = j - *mumas; /* Computing MIN */ i__5 = *nm1, i__6 = j + *mlmas; i__3 = min(i__5,i__6); for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) { ffma = fmas[i__ - j + linal_1.mbdiag + j * fmas_dim1]; e2r[i__ + j * e2r_dim1] += *alphn * ffma; e2i[i__ + j * e2i_dim1] += *betan * ffma; } } goto L45; /* ----------------------------------------------------------- */ L4: /* --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__3 = linal_1.mbjac; for (i__ = 1; i__ <= i__3; ++i__) { imle = i__ + linal_1.mle; e2r[imle + j * e2r_dim1] = -fjac[i__ + j * fjac_dim1]; e2i[imle + j * e2i_dim1] = 0.; } /* Computing MAX */ i__3 = 1, i__4 = *mumas + 2 - j; /* Computing MIN */ i__5 = linal_1.mbb, i__6 = *mumas + 1 - j + *n; i__2 = min(i__5,i__6); for (i__ = max(i__3,i__4); i__ <= i__2; ++i__) { ib = i__ + linal_1.mdiff; bb = fmas[i__ + j * fmas_dim1]; e2r[ib + j * e2r_dim1] += *alphn * bb; e2i[ib + j * e2i_dim1] = *betan * bb; } } decbc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &linal_1.mle, & linal_1.mue, &ip2[1], ier); return 0; /* ----------------------------------------------------------- */ L14: /* --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER */ i__1 = *nm1; for (j = 1; j <= i__1; ++j) { jm1 = j + *m1; i__2 = linal_1.mbjac; for (i__ = 1; i__ <= i__2; ++i__) { e2r[i__ + linal_1.mle + j * e2r_dim1] = -fjac[i__ + jm1 * fjac_dim1]; e2i[i__ + linal_1.mle + j * e2i_dim1] = 0.; } i__2 = linal_1.mbb; for (i__ = 1; i__ <= i__2; ++i__) { ib = i__ + linal_1.mdiff; ffma = fmas[i__ + j * fmas_dim1]; e2r[ib + j * e2r_dim1] += *alphn * ffma; e2i[ib + j * e2i_dim1] += *betan * ffma; } } goto L46; /* ----------------------------------------------------------- */ L5: /* --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX */ i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { bb = fmas[i__ + j * fmas_dim1]; e2r[i__ + j * e2r_dim1] = bb * *alphn - fjac[i__ + j * fjac_dim1]; e2i[i__ + j * e2i_dim1] = bb * *betan; } } decc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &ip2[1], ier); return 0; /* ----------------------------------------------------------- */ L15: /* --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *nm1; for (j = 1; j <= i__1; ++j) { jm1 = j + *m1; i__2 = *nm1; for (i__ = 1; i__ <= i__2; ++i__) { e2r[i__ + j * e2r_dim1] = *alphn * fmas[i__ + j * fmas_dim1] - fjac[i__ + jm1 * fjac_dim1]; e2i[i__ + j * e2i_dim1] = *betan * fmas[i__ + j * fmas_dim1]; } } goto L45; /* ----------------------------------------------------------- */ L6: /* --- B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX */ /* --- THIS OPTION IS NOT PROVIDED */ return 0; /* ----------------------------------------------------------- */ L7: /* --- B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION */ i__1 = *n - 1; for (j = 1; j <= i__1; ++j) { j1 = j + 1; e2r[j1 + j * e2r_dim1] = -fjac[j1 + j * fjac_dim1]; e2i[j1 + j * e2i_dim1] = 0.; } i__1 = *n; for (j = 1; j <= i__1; ++j) { i__2 = j; for (i__ = 1; i__ <= i__2; ++i__) { e2i[i__ + j * e2i_dim1] = 0.; e2r[i__ + j * e2r_dim1] = -fjac[i__ + j * fjac_dim1]; } e2r[j + j * e2r_dim1] += *alphn; e2i[j + j * e2i_dim1] = *betan; } dechc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &c__1, &ip2[1], ier); return 0; /* ----------------------------------------------------------- */ L55: return 0; } /* decomc_ */ /* END OF SUBROUTINE DECOMC */ /* *********************************************************** */ /* Subroutine */ int slvrad_(integer *n, doublereal *fjac, integer *ldjac, integer *mljac, integer *mujac, doublereal *fmas, integer *ldmas, integer *mlmas, integer *mumas, integer *m1, integer *m2, integer * nm1, doublereal *fac1, doublereal *alphn, doublereal *betan, doublereal *e1, doublereal *e2r, doublereal *e2i, integer *lde1, doublereal *z1, doublereal *z2, doublereal *z3, doublereal *f1, doublereal *f2, doublereal *f3, doublereal *cont, integer *ip1, integer *ip2, integer *iphes, integer *ier, integer *ijob) { /* System generated locals */ integer fjac_dim1, fjac_offset, fmas_dim1, fmas_offset, e1_dim1, e1_offset, e2r_dim1, e2r_offset, e2i_dim1, e2i_offset, i__1, i__2, i__3, i__4, i__5, i__6; doublereal d__1, d__2; /* Local variables */ static doublereal ffja, abno; extern /* Subroutine */ int solb_(integer *, integer *, doublereal *, integer *, integer *, doublereal *, integer *), solc_(integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, integer *), solh_(integer *, integer *, doublereal *, integer *, doublereal *, integer *); static doublereal sumh, e1imp; static integer i__, j, k; extern /* Subroutine */ int solbc_(integer *, integer *, doublereal *, doublereal *, integer *, integer *, doublereal *, doublereal *, integer *); static doublereal zsafe; extern /* Subroutine */ int solhc_(integer *, integer *, doublereal *, doublereal *, integer *, doublereal *, doublereal *, integer *); static doublereal s1, s2, s3, bb; static integer mm, mp, j1b, j2b, im1, jm1, mp1; static doublereal z2i, z3i; static integer jkm, mpi; extern /* Subroutine */ int sol_(integer *, integer *, doublereal *, doublereal *, integer *); static doublereal sum1, sum2, sum3; /* Parameter adjustments */ --iphes; --f3; --f2; --f1; --z3; --z2; --z1; fjac_dim1 = *ldjac; fjac_offset = fjac_dim1 + 1; fjac -= fjac_offset; --ip2; --ip1; fmas_dim1 = *ldmas; fmas_offset = fmas_dim1 + 1; fmas -= fmas_offset; e2i_dim1 = *lde1; e2i_offset = e2i_dim1 + 1; e2i -= e2i_offset; e2r_dim1 = *lde1; e2r_offset = e2r_dim1 + 1; e2r -= e2r_offset; e1_dim1 = *lde1; e1_offset = e1_dim1 + 1; e1 -= e1_offset; /* Function Body */ switch (*ijob) { case 1: goto L1; case 2: goto L2; case 3: goto L3; case 4: goto L4; case 5: goto L5; case 6: goto L6; case 7: goto L7; case 8: goto L55; case 9: goto L55; case 10: goto L55; case 11: goto L11; case 12: goto L12; case 13: goto L13; case 14: goto L13; case 15: goto L15; } /* ----------------------------------------------------------- */ L1: /* --- B=IDENTITY, JACOBIAN A FULL MATRIX */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { s2 = -f2[i__]; s3 = -f3[i__]; z1[i__] -= f1[i__] * *fac1; z2[i__] = z2[i__] + s2 * *alphn - s3 * *betan; z3[i__] = z3[i__] + s3 * *alphn + s2 * *betan; } sol_(n, lde1, &e1[e1_offset], &z1[1], &ip1[1]); solc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &z2[1], &z3[1], &ip2[1] ); return 0; /* ----------------------------------------------------------- */ L11: /* --- B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { s2 = -f2[i__]; s3 = -f3[i__]; z1[i__] -= f1[i__] * *fac1; z2[i__] = z2[i__] + s2 * *alphn - s3 * *betan; z3[i__] = z3[i__] + s3 * *alphn + s2 * *betan; } L48: /* Computing 2nd power */ d__1 = *alphn; /* Computing 2nd power */ d__2 = *betan; abno = d__1 * d__1 + d__2 * d__2; mm = *m1 / *m2; i__1 = *m2; for (j = 1; j <= i__1; ++j) { sum1 = 0.; sum2 = 0.; sum3 = 0.; for (k = mm - 1; k >= 0; --k) { jkm = j + k * *m2; sum1 = (z1[jkm] + sum1) / *fac1; sumh = (z2[jkm] + sum2) / abno; sum3 = (z3[jkm] + sum3) / abno; sum2 = sumh * *alphn + sum3 * *betan; sum3 = sum3 * *alphn - sumh * *betan; i__2 = *nm1; for (i__ = 1; i__ <= i__2; ++i__) { im1 = i__ + *m1; z1[im1] += fjac[i__ + jkm * fjac_dim1] * sum1; z2[im1] += fjac[i__ + jkm * fjac_dim1] * sum2; z3[im1] += fjac[i__ + jkm * fjac_dim1] * sum3; } } } sol_(nm1, lde1, &e1[e1_offset], &z1[*m1 + 1], &ip1[1]); solc_(nm1, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &z2[*m1 + 1], &z3[* m1 + 1], &ip2[1]); L49: for (i__ = *m1; i__ >= 1; --i__) { mpi = *m2 + i__; z1[i__] = (z1[i__] + z1[mpi]) / *fac1; z2i = z2[i__] + z2[mpi]; z3i = z3[i__] + z3[mpi]; z3[i__] = (z3i * *alphn - z2i * *betan) / abno; z2[i__] = (z2i * *alphn + z3i * *betan) / abno; } return 0; /* ----------------------------------------------------------- */ L2: /* --- B=IDENTITY, JACOBIAN A BANDED MATRIX */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { s2 = -f2[i__]; s3 = -f3[i__]; z1[i__] -= f1[i__] * *fac1; z2[i__] = z2[i__] + s2 * *alphn - s3 * *betan; z3[i__] = z3[i__] + s3 * *alphn + s2 * *betan; } solb_(n, lde1, &e1[e1_offset], &linal_1.mle, &linal_1.mue, &z1[1], &ip1[1] ); solbc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &linal_1.mle, & linal_1.mue, &z2[1], &z3[1], &ip2[1]); return 0; /* ----------------------------------------------------------- */ L12: /* --- B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { s2 = -f2[i__]; s3 = -f3[i__]; z1[i__] -= f1[i__] * *fac1; z2[i__] = z2[i__] + s2 * *alphn - s3 * *betan; z3[i__] = z3[i__] + s3 * *alphn + s2 * *betan; } L45: /* Computing 2nd power */ d__1 = *alphn; /* Computing 2nd power */ d__2 = *betan; abno = d__1 * d__1 + d__2 * d__2; mm = *m1 / *m2; i__1 = *m2; for (j = 1; j <= i__1; ++j) { sum1 = 0.; sum2 = 0.; sum3 = 0.; for (k = mm - 1; k >= 0; --k) { jkm = j + k * *m2; sum1 = (z1[jkm] + sum1) / *fac1; sumh = (z2[jkm] + sum2) / abno; sum3 = (z3[jkm] + sum3) / abno; sum2 = sumh * *alphn + sum3 * *betan; sum3 = sum3 * *alphn - sumh * *betan; /* Computing MAX */ i__2 = 1, i__3 = j - *mujac; /* Computing MIN */ i__5 = *nm1, i__6 = j + *mljac; i__4 = min(i__5,i__6); for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) { im1 = i__ + *m1; ffja = fjac[i__ + *mujac + 1 - j + jkm * fjac_dim1]; z1[im1] += ffja * sum1; z2[im1] += ffja * sum2; z3[im1] += ffja * sum3; } } } solb_(nm1, lde1, &e1[e1_offset], &linal_1.mle, &linal_1.mue, &z1[*m1 + 1], &ip1[1]); solbc_(nm1, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &linal_1.mle, & linal_1.mue, &z2[*m1 + 1], &z3[*m1 + 1], &ip2[1]); goto L49; /* ----------------------------------------------------------- */ L3: /* --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { s1 = 0.; s2 = 0.; s3 = 0.; /* Computing MAX */ i__4 = 1, i__2 = i__ - *mlmas; /* Computing MIN */ i__5 = *n, i__6 = i__ + *mumas; i__3 = min(i__5,i__6); for (j = max(i__4,i__2); j <= i__3; ++j) { bb = fmas[i__ - j + linal_1.mbdiag + j * fmas_dim1]; s1 -= bb * f1[j]; s2 -= bb * f2[j]; s3 -= bb * f3[j]; } z1[i__] += s1 * *fac1; z2[i__] = z2[i__] + s2 * *alphn - s3 * *betan; z3[i__] = z3[i__] + s3 * *alphn + s2 * *betan; } sol_(n, lde1, &e1[e1_offset], &z1[1], &ip1[1]); solc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &z2[1], &z3[1], &ip2[1] ); return 0; /* ----------------------------------------------------------- */ L13: /* --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *m1; for (i__ = 1; i__ <= i__1; ++i__) { s2 = -f2[i__]; s3 = -f3[i__]; z1[i__] -= f1[i__] * *fac1; z2[i__] = z2[i__] + s2 * *alphn - s3 * *betan; z3[i__] = z3[i__] + s3 * *alphn + s2 * *betan; } i__1 = *nm1; for (i__ = 1; i__ <= i__1; ++i__) { im1 = i__ + *m1; s1 = 0.; s2 = 0.; s3 = 0.; /* Computing MAX */ i__3 = 1, i__4 = i__ - *mlmas; j1b = max(i__3,i__4); /* Computing MIN */ i__3 = *nm1, i__4 = i__ + *mumas; j2b = min(i__3,i__4); i__3 = j2b; for (j = j1b; j <= i__3; ++j) { jm1 = j + *m1; bb = fmas[i__ - j + linal_1.mbdiag + j * fmas_dim1]; s1 -= bb * f1[jm1]; s2 -= bb * f2[jm1]; s3 -= bb * f3[jm1]; } z1[im1] += s1 * *fac1; z2[im1] = z2[im1] + s2 * *alphn - s3 * *betan; z3[im1] = z3[im1] + s3 * *alphn + s2 * *betan; } if (*ijob == 14) { goto L45; } goto L48; /* ----------------------------------------------------------- */ L4: /* --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { s1 = 0.; s2 = 0.; s3 = 0.; /* Computing MAX */ i__3 = 1, i__4 = i__ - *mlmas; /* Computing MIN */ i__5 = *n, i__6 = i__ + *mumas; i__2 = min(i__5,i__6); for (j = max(i__3,i__4); j <= i__2; ++j) { bb = fmas[i__ - j + linal_1.mbdiag + j * fmas_dim1]; s1 -= bb * f1[j]; s2 -= bb * f2[j]; s3 -= bb * f3[j]; } z1[i__] += s1 * *fac1; z2[i__] = z2[i__] + s2 * *alphn - s3 * *betan; z3[i__] = z3[i__] + s3 * *alphn + s2 * *betan; } solb_(n, lde1, &e1[e1_offset], &linal_1.mle, &linal_1.mue, &z1[1], &ip1[1] ); solbc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &linal_1.mle, & linal_1.mue, &z2[1], &z3[1], &ip2[1]); return 0; /* ----------------------------------------------------------- */ L5: /* --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { s1 = 0.; s2 = 0.; s3 = 0.; i__2 = *n; for (j = 1; j <= i__2; ++j) { bb = fmas[i__ + j * fmas_dim1]; s1 -= bb * f1[j]; s2 -= bb * f2[j]; s3 -= bb * f3[j]; } z1[i__] += s1 * *fac1; z2[i__] = z2[i__] + s2 * *alphn - s3 * *betan; z3[i__] = z3[i__] + s3 * *alphn + s2 * *betan; } sol_(n, lde1, &e1[e1_offset], &z1[1], &ip1[1]); solc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &z2[1], &z3[1], &ip2[1] ); return 0; /* ----------------------------------------------------------- */ L15: /* --- B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *m1; for (i__ = 1; i__ <= i__1; ++i__) { s2 = -f2[i__]; s3 = -f3[i__]; z1[i__] -= f1[i__] * *fac1; z2[i__] = z2[i__] + s2 * *alphn - s3 * *betan; z3[i__] = z3[i__] + s3 * *alphn + s2 * *betan; } i__1 = *nm1; for (i__ = 1; i__ <= i__1; ++i__) { im1 = i__ + *m1; s1 = 0.; s2 = 0.; s3 = 0.; i__2 = *nm1; for (j = 1; j <= i__2; ++j) { jm1 = j + *m1; bb = fmas[i__ + j * fmas_dim1]; s1 -= bb * f1[jm1]; s2 -= bb * f2[jm1]; s3 -= bb * f3[jm1]; } z1[im1] += s1 * *fac1; z2[im1] = z2[im1] + s2 * *alphn - s3 * *betan; z3[im1] = z3[im1] + s3 * *alphn + s2 * *betan; } goto L48; /* ----------------------------------------------------------- */ L6: /* --- B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX */ /* --- THIS OPTION IS NOT PROVIDED */ return 0; /* ----------------------------------------------------------- */ L7: /* --- B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { s2 = -f2[i__]; s3 = -f3[i__]; z1[i__] -= f1[i__] * *fac1; z2[i__] = z2[i__] + s2 * *alphn - s3 * *betan; z3[i__] = z3[i__] + s3 * *alphn + s2 * *betan; } for (mm = *n - 2; mm >= 1; --mm) { mp = *n - mm; mp1 = mp - 1; i__ = iphes[mp]; if (i__ == mp) { goto L746; } zsafe = z1[mp]; z1[mp] = z1[i__]; z1[i__] = zsafe; zsafe = z2[mp]; z2[mp] = z2[i__]; z2[i__] = zsafe; zsafe = z3[mp]; z3[mp] = z3[i__]; z3[i__] = zsafe; L746: i__1 = *n; for (i__ = mp + 1; i__ <= i__1; ++i__) { e1imp = fjac[i__ + mp1 * fjac_dim1]; z1[i__] -= e1imp * z1[mp]; z2[i__] -= e1imp * z2[mp]; z3[i__] -= e1imp * z3[mp]; } } solh_(n, lde1, &e1[e1_offset], &c__1, &z1[1], &ip1[1]); solhc_(n, lde1, &e2r[e2r_offset], &e2i[e2i_offset], &c__1, &z2[1], &z3[1], &ip2[1]); i__1 = *n - 2; for (mm = 1; mm <= i__1; ++mm) { mp = *n - mm; mp1 = mp - 1; i__2 = *n; for (i__ = mp + 1; i__ <= i__2; ++i__) { e1imp = fjac[i__ + mp1 * fjac_dim1]; z1[i__] += e1imp * z1[mp]; z2[i__] += e1imp * z2[mp]; z3[i__] += e1imp * z3[mp]; } i__ = iphes[mp]; if (i__ == mp) { goto L750; } zsafe = z1[mp]; z1[mp] = z1[i__]; z1[i__] = zsafe; zsafe = z2[mp]; z2[mp] = z2[i__]; z2[i__] = zsafe; zsafe = z3[mp]; z3[mp] = z3[i__]; z3[i__] = zsafe; L750: ; } return 0; /* ----------------------------------------------------------- */ L55: return 0; } /* slvrad_ */ /* END OF SUBROUTINE SLVRAD */ /* *********************************************************** */ /* Subroutine */ int estrad_(integer *n, doublereal *fjac, integer *ldjac, integer *mljac, integer *mujac, doublereal *fmas, integer *ldmas, integer *mlmas, integer *mumas, doublereal *h__, doublereal *dd1, doublereal *dd2, doublereal *dd3, voidsub fcn, integer *nfcn, doublereal *y0, doublereal *y, integer *ijob, doublereal *x, integer *m1, integer *m2, integer *nm1, doublereal *e1, integer *lde1, doublereal * z1, doublereal *z2, doublereal *z3, doublereal *cont, doublereal *f1, doublereal *f2, integer *ip1, integer *iphes, doublereal *scal, doublereal *err, logical *first, logical *reject, doublereal *fac1) { /* System generated locals */ integer fjac_dim1, fjac_offset, fmas_dim1, fmas_offset, e1_dim1, e1_offset, i__1, i__2, i__3, i__4, i__5, i__6; doublereal d__1; /* Local variables */ extern /* Subroutine */ int solb_(integer *, integer *, doublereal *, integer *, integer *, doublereal *, integer *), solh_(integer *, integer *, doublereal *, integer *, doublereal *, integer *); static integer i__, j, k; static doublereal zsafe; static integer mm, mp, im1; extern /* Subroutine */ int sol_(integer *, integer *, doublereal *, doublereal *, integer *); static doublereal sum, hee1, hee2, hee3, sum1; /* Parameter adjustments */ --scal; --iphes; --f2; --f1; --cont; --z3; --z2; --z1; --y; --y0; fjac_dim1 = *ldjac; fjac_offset = fjac_dim1 + 1; fjac -= fjac_offset; --ip1; fmas_dim1 = *ldmas; fmas_offset = fmas_dim1 + 1; fmas -= fmas_offset; e1_dim1 = *lde1; e1_offset = e1_dim1 + 1; e1 -= e1_offset; /* Function Body */ hee1 = *dd1 / *h__; hee2 = *dd2 / *h__; hee3 = *dd3 / *h__; switch (*ijob) { case 1: goto L1; case 2: goto L2; case 3: goto L3; case 4: goto L4; case 5: goto L5; case 6: goto L6; case 7: goto L7; case 8: goto L55; case 9: goto L55; case 10: goto L55; case 11: goto L11; case 12: goto L12; case 13: goto L13; case 14: goto L14; case 15: goto L15; } L1: /* ------ B=IDENTITY, JACOBIAN A FULL MATRIX */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { f2[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; cont[i__] = f2[i__] + y0[i__]; } sol_(n, lde1, &e1[e1_offset], &cont[1], &ip1[1]); goto L77; L11: /* ------ B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { f2[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; cont[i__] = f2[i__] + y0[i__]; } L48: mm = *m1 / *m2; i__1 = *m2; for (j = 1; j <= i__1; ++j) { sum1 = 0.; for (k = mm - 1; k >= 0; --k) { sum1 = (cont[j + k * *m2] + sum1) / *fac1; i__2 = *nm1; for (i__ = 1; i__ <= i__2; ++i__) { im1 = i__ + *m1; cont[im1] += fjac[i__ + (j + k * *m2) * fjac_dim1] * sum1; } } } sol_(nm1, lde1, &e1[e1_offset], &cont[*m1 + 1], &ip1[1]); for (i__ = *m1; i__ >= 1; --i__) { cont[i__] = (cont[i__] + cont[*m2 + i__]) / *fac1; } goto L77; L2: /* ------ B=IDENTITY, JACOBIAN A BANDED MATRIX */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { f2[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; cont[i__] = f2[i__] + y0[i__]; } solb_(n, lde1, &e1[e1_offset], &linal_1.mle, &linal_1.mue, &cont[1], &ip1[ 1]); goto L77; L12: /* ------ B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { f2[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; cont[i__] = f2[i__] + y0[i__]; } L45: mm = *m1 / *m2; i__1 = *m2; for (j = 1; j <= i__1; ++j) { sum1 = 0.; for (k = mm - 1; k >= 0; --k) { sum1 = (cont[j + k * *m2] + sum1) / *fac1; /* Computing MAX */ i__2 = 1, i__3 = j - *mujac; /* Computing MIN */ i__5 = *nm1, i__6 = j + *mljac; i__4 = min(i__5,i__6); for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) { im1 = i__ + *m1; cont[im1] += fjac[i__ + *mujac + 1 - j + (j + k * *m2) * fjac_dim1] * sum1; } } } solb_(nm1, lde1, &e1[e1_offset], &linal_1.mle, &linal_1.mue, &cont[*m1 + 1], &ip1[1]); for (i__ = *m1; i__ >= 1; --i__) { cont[i__] = (cont[i__] + cont[*m2 + i__]) / *fac1; } goto L77; L3: /* ------ B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { f1[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; } i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { sum = 0.; /* Computing MAX */ i__4 = 1, i__2 = i__ - *mlmas; /* Computing MIN */ i__5 = *n, i__6 = i__ + *mumas; i__3 = min(i__5,i__6); for (j = max(i__4,i__2); j <= i__3; ++j) { sum += fmas[i__ - j + linal_1.mbdiag + j * fmas_dim1] * f1[j]; } f2[i__] = sum; cont[i__] = sum + y0[i__]; } sol_(n, lde1, &e1[e1_offset], &cont[1], &ip1[1]); goto L77; L13: /* ------ B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *m1; for (i__ = 1; i__ <= i__1; ++i__) { f2[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; cont[i__] = f2[i__] + y0[i__]; } i__1 = *n; for (i__ = *m1 + 1; i__ <= i__1; ++i__) { f1[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; } i__1 = *nm1; for (i__ = 1; i__ <= i__1; ++i__) { sum = 0.; /* Computing MAX */ i__3 = 1, i__4 = i__ - *mlmas; /* Computing MIN */ i__5 = *nm1, i__6 = i__ + *mumas; i__2 = min(i__5,i__6); for (j = max(i__3,i__4); j <= i__2; ++j) { sum += fmas[i__ - j + linal_1.mbdiag + j * fmas_dim1] * f1[j + * m1]; } im1 = i__ + *m1; f2[im1] = sum; cont[im1] = sum + y0[im1]; } goto L48; L4: /* ------ B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { f1[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; } i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { sum = 0.; /* Computing MAX */ i__2 = 1, i__3 = i__ - *mlmas; /* Computing MIN */ i__5 = *n, i__6 = i__ + *mumas; i__4 = min(i__5,i__6); for (j = max(i__2,i__3); j <= i__4; ++j) { sum += fmas[i__ - j + linal_1.mbdiag + j * fmas_dim1] * f1[j]; } f2[i__] = sum; cont[i__] = sum + y0[i__]; } solb_(n, lde1, &e1[e1_offset], &linal_1.mle, &linal_1.mue, &cont[1], &ip1[ 1]); goto L77; L14: /* ------ B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER */ i__1 = *m1; for (i__ = 1; i__ <= i__1; ++i__) { f2[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; cont[i__] = f2[i__] + y0[i__]; } i__1 = *n; for (i__ = *m1 + 1; i__ <= i__1; ++i__) { f1[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; } i__1 = *nm1; for (i__ = 1; i__ <= i__1; ++i__) { sum = 0.; /* Computing MAX */ i__4 = 1, i__2 = i__ - *mlmas; /* Computing MIN */ i__5 = *nm1, i__6 = i__ + *mumas; i__3 = min(i__5,i__6); for (j = max(i__4,i__2); j <= i__3; ++j) { sum += fmas[i__ - j + linal_1.mbdiag + j * fmas_dim1] * f1[j + * m1]; } im1 = i__ + *m1; f2[im1] = sum; cont[im1] = sum + y0[im1]; } goto L45; L5: /* ------ B IS A FULL MATRIX, JACOBIAN A FULL MATRIX */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { f1[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; } i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { sum = 0.; i__3 = *n; for (j = 1; j <= i__3; ++j) { sum += fmas[i__ + j * fmas_dim1] * f1[j]; } f2[i__] = sum; cont[i__] = sum + y0[i__]; } sol_(n, lde1, &e1[e1_offset], &cont[1], &ip1[1]); goto L77; L15: /* ------ B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER */ i__1 = *m1; for (i__ = 1; i__ <= i__1; ++i__) { f2[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; cont[i__] = f2[i__] + y0[i__]; } i__1 = *n; for (i__ = *m1 + 1; i__ <= i__1; ++i__) { f1[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; } i__1 = *nm1; for (i__ = 1; i__ <= i__1; ++i__) { sum = 0.; i__3 = *nm1; for (j = 1; j <= i__3; ++j) { sum += fmas[i__ + j * fmas_dim1] * f1[j + *m1]; } im1 = i__ + *m1; f2[im1] = sum; cont[im1] = sum + y0[im1]; } goto L48; L6: /* ------ B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX */ /* ------ THIS OPTION IS NOT PROVIDED */ return 0; L7: /* ------ B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION */ i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { f2[i__] = hee1 * z1[i__] + hee2 * z2[i__] + hee3 * z3[i__]; cont[i__] = f2[i__] + y0[i__]; } for (mm = *n - 2; mm >= 1; --mm) { mp = *n - mm; i__ = iphes[mp]; if (i__ == mp) { goto L310; } zsafe = cont[mp]; cont[mp] = cont[i__]; cont[i__] = zsafe; L310: i__1 = *n; for (i__ = mp + 1; i__ <= i__1; ++i__) { cont[i__] -= fjac[i__ + (mp - 1) * fjac_dim1] * cont[mp]; } } solh_(n, lde1, &e1[e1_offset], &c__1, &cont[1], &ip1[1]); i__1 = *n - 2; for (mm = 1; mm <= i__1; ++mm) { mp = *n - mm; i__3 = *n; for (i__ = mp + 1; i__ <= i__3; ++i__) { cont[i__] += fjac[i__ + (mp - 1) * fjac_dim1] * cont[mp]; } i__ = iphes[mp]; if (i__ == mp) { goto L440; } zsafe = cont[mp]; cont[mp] = cont[i__]; cont[i__] = zsafe; L440: ; } /* -------------------------------------- */ L77: *err = 0.; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { /* Computing 2nd power */ d__1 = cont[i__] / scal[i__]; *err += d__1 * d__1; } /* Computing MAX */ d__1 = sqrt(*err / *n); *err = max(d__1,1e-10); if (*err < 1.) { return 0; } if (*first || *reject) { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { cont[i__] = y[i__] + cont[i__]; } (*fcn)(*n, *x, &cont[1], &f1[1]); ++(*nfcn); i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { cont[i__] = f1[i__] + f2[i__]; } switch (*ijob) { case 1: goto L31; case 2: goto L32; case 3: goto L31; case 4: goto L32; case 5: goto L31; case 6: goto L32; case 7: goto L33; case 8: goto L55; case 9: goto L55; case 10: goto L55; case 11: goto L41; case 12: goto L42; case 13: goto L41; case 14: goto L42; case 15: goto L41; } /* ------ FULL MATRIX OPTION */ L31: sol_(n, lde1, &e1[e1_offset], &cont[1], &ip1[1]); goto L88; /* ------ FULL MATRIX OPTION, SECOND ORDER */ L41: i__1 = *m2; for (j = 1; j <= i__1; ++j) { sum1 = 0.; for (k = mm - 1; k >= 0; --k) { sum1 = (cont[j + k * *m2] + sum1) / *fac1; i__3 = *nm1; for (i__ = 1; i__ <= i__3; ++i__) { im1 = i__ + *m1; cont[im1] += fjac[i__ + (j + k * *m2) * fjac_dim1] * sum1; } } } sol_(nm1, lde1, &e1[e1_offset], &cont[*m1 + 1], &ip1[1]); for (i__ = *m1; i__ >= 1; --i__) { cont[i__] = (cont[i__] + cont[*m2 + i__]) / *fac1; } goto L88; /* ------ BANDED MATRIX OPTION */ L32: solb_(n, lde1, &e1[e1_offset], &linal_1.mle, &linal_1.mue, &cont[1], & ip1[1]); goto L88; /* ------ BANDED MATRIX OPTION, SECOND ORDER */ L42: i__1 = *m2; for (j = 1; j <= i__1; ++j) { sum1 = 0.; for (k = mm - 1; k >= 0; --k) { sum1 = (cont[j + k * *m2] + sum1) / *fac1; /* Computing MAX */ i__3 = 1, i__4 = j - *mujac; /* Computing MIN */ i__5 = *nm1, i__6 = j + *mljac; i__2 = min(i__5,i__6); for (i__ = max(i__3,i__4); i__ <= i__2; ++i__) { im1 = i__ + *m1; cont[im1] += fjac[i__ + *mujac + 1 - j + (j + k * *m2) * fjac_dim1] * sum1; } } } solb_(nm1, lde1, &e1[e1_offset], &linal_1.mle, &linal_1.mue, &cont[* m1 + 1], &ip1[1]); for (i__ = *m1; i__ >= 1; --i__) { cont[i__] = (cont[i__] + cont[*m2 + i__]) / *fac1; } goto L88; /* ------ HESSENBERG MATRIX OPTION */ L33: for (mm = *n - 2; mm >= 1; --mm) { mp = *n - mm; i__ = iphes[mp]; if (i__ == mp) { goto L510; } zsafe = cont[mp]; cont[mp] = cont[i__]; cont[i__] = zsafe; L510: i__1 = *n; for (i__ = mp + 1; i__ <= i__1; ++i__) { cont[i__] -= fjac[i__ + (mp - 1) * fjac_dim1] * cont[mp]; } } solh_(n, lde1, &e1[e1_offset], &c__1, &cont[1], &ip1[1]); i__1 = *n - 2; for (mm = 1; mm <= i__1; ++mm) { mp = *n - mm; i__2 = *n; for (i__ = mp + 1; i__ <= i__2; ++i__) { cont[i__] += fjac[i__ + (mp - 1) * fjac_dim1] * cont[mp]; } i__ = iphes[mp]; if (i__ == mp) { goto L640; } zsafe = cont[mp]; cont[mp] = cont[i__]; cont[i__] = zsafe; L640: ; } /* ----------------------------------- */ L88: *err = 0.; i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { /* Computing 2nd power */ d__1 = cont[i__] / scal[i__]; *err += d__1 * d__1; } /* Computing MAX */ d__1 = sqrt(*err / *n); *err = max(d__1,1e-10); } return 0; /* ----------------------------------------------------------- */ L55: return 0; } /* estrad_ */