SUBROUTINE BVPSOL(FCN,BC,IVPSOL,N,M,T,X,EPS,IOPT,INFO,IRW,RW, & IIW,IW) c*********************************************************************72 c cc bvpsol() solves highly nonlinear two point boundary value problems. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) EXTERNAL FCN,BC,IVPSOL INTEGER N,M DOUBLE PRECISION T(M) DOUBLE PRECISION X(N,M) DOUBLE PRECISION EPS INTEGER INFO,IOPT(5) INTEGER IRW DOUBLE PRECISION RW(IRW) INTEGER IIW INTEGER IW(IIW) C C ------------------------------------------------------------ C C* Title C C (B)oundary (V)alue (P)roblem (Sol)lver for highly nonlinear C two point boundary value problems using a local linear C solver (condensing algorithm) or a global sparse linear solver C for the solution of the arising C linear subproblems. C C* Written by P. Deuflhard, G.Bader, L. Weimann C* Purpose Solution of nonlinear two-point boundary value C problems. C* Method Local and Global Nonlinear two-point Boundary Value C Problems solver (Multiple shooting approach) C* Category I1b2a - Differential and integral equations C Two point boundary value problems C* Keywords Nonlinear boundary value problems, Multiple C shooting, Newton methods C* Version 1.2 C* Revision February 2002 C* Latest Change January 2004 C* Library CodeLib C* Code Fortran 77, Double Precision C* Environment Standard Fortran 77 environment on PC's, C workstations and hosts. C* Copyright (c) Konrad-Zuse-Zentrum fuer C Informationstechnik Berlin (ZIB) C Takustrasse 7, D-14195 Berlin-Dahlem C phone : + 49/30/84185-0 C fax : + 49/30/84185-125 C* Contact Lutz Weimann C ZIB, Division Scientific Computing, C Department Numerical Analysis and Modelling C phone : + 49/30/84185-185 C fax : + 49/30/84185-107 C e-mail: weimann@zib.de C C* References: C C /1/ R.Bulirsch: C Die Mehrzielmethode zur numerischen Loesung von C nichtlinearen Randwertproblemen und Aufgaben der C optimalen Steuerung. C Carl-Cranz-Gesellschaft: Tech.Rep. (Oct.1971) C C /2/ J.Stoer, R.Bulirsch: C Einfuehrung in die Numerische Mathematik II. C Berlin, Heidelberg, New York: Springer (1st Ed. 1973) C C /3/ P.Deuflhard: C A Modified Newton Method for the Solution of C Ill-Conditioned Systems of Nonlinear Equations with C Application to Multiple Shooting. C Numer. Math. 22, 289-315 (1974) C C /4/ P.Deuflhard: C Recent Advances in Multiple Shooting Techniques. C (Survey Article including further References) C In: I.Gladwell, D.K.Sayers (Ed.): Computational C Techniques for Ordinary Differential Equations. C Section 10, P.217-272. C London, New York: Academic Press (1980) C C /5/ P.Deuflhard, G.Bader: C Multiple Shooting Techniques Revisited. C Univ. Heidelberg, SFB 123, Tech. Rep. 163 (1982) C C /6/ P. Deuflhard: C Newton Methods for Nonlinear Problems. - C Affine Invariance and Adaptive Algorithms. C Series Computational Mathematics 35, Springer (2004) C C --------------------------------------------------------------- C C* Licence C You may use or modify this code for your own non commercial C purposes for an unlimited time. C In any case you should not deliver this code without a special C permission of ZIB. C In case you intend to use the code commercially, we oblige you C to sign an according licence agreement with ZIB. C C* Warranty C This code has been tested up to a certain level. Defects and C weaknesses, which may be included in the code, do not establish C any warranties by ZIB. ZIB does not take over any liabilities C which may follow from aquisition or application of this code. C C* Software status C This code is under partial care of ZIB and belongs to ZIB software C class 2. C C ------------------------------------------------------------ C C External subroutines (to be supplied by the user) C ================================================= C C FCN(N,T,Y,DY) Right-hand side of system of C first-order differential equations C N Input: Number of first order ODE's C T Input: Actual position in the C interval ³ A , B ü C Y(N) Input: Values at T C DY(N) Output: Derivatives at T C C BC(YA,YB,R) Two-point boundary conditions at ( A C = T(1), B = T(M)) C YA(N) Input: Values at A = T(1) C YB(N) Input: Values at B = T(M) C R(N) Output: Values of C boundary conditions function C C IVPSOL(N,FCN,T,Y,TEND,TOL,HMAX,H,KFLAG) C Initial value problem (IVP) C integrator C N Number of first-order ODE's C FCN Right-hand side of the ODE's system C ( see above ) C T Input: Starting point of integration C T.LT.TEND C Output: Achieved final point of C integration C Y(N) Input and Output: Values at T C TEND Input: Prescribed final point of C integration C TOL Input: Prescribed relative precision C (>0) C HMAX Input: Maximum permitted stepsize C H Input: Initial stepsize guess C Output: Stepsize proposal for next C integration step ( H.EQ.0 , if C IVPSOL fails to proceed ) C KFLAG Input: Print parameter C Output: Error flag ( KFLAG.LT.0 C indicates an error ) . C For further details, see IVPSOL . C C Input parameters (* marks inout parameters) C =========================================== C C N Number of first-order ordinary C differential equations. C M Number of Shooting nodes. C =2 Single Shooting C >2 Multiple Shooting C T(M) Single/Multiple Shooting Nodes C ( T(1)= A , T(M)= B ) C * X(N,M) Start data for Newton iteration. C EPS Required relative precision of C solution. C IOPT(5) Options array. The options have the C meanings as described below: C IOPT(1) - MAXIT Maximum permitted number of C iteration steps. C IOPT(2) - NONLIN Boundary value problem C classification by user: C 0 Linear boundary value problem. C 1 Nonlinear boundary value C problem. Good initial data C available. C 2 Highly nonlinear boundary C value problem. Only bad C initial data available. Small C initial damping factor in C Gauss Newton method. C 3 Highly nonlinear boundary C value problem. Only bad C initial data available. Small C initial damping factor in C Gauss Newton method. C Additionally initial rank C reduction to separable linear C boundary conditions. C IOPT(3) - IBVPSL Solution method switch: C 0 use local linear solver with C condensing algorithm C 1 use global sparse linear solver C (formerly realized as "BVPSOG") C IOPT(4) - IPRINT Print parameter: C -1 No print C 0 Print initial data, iterative C values of level functions, C solution data (or final data, C respectively) C +1 Additionally print iterates C T(J),X(I,J), I = 1,...,N , J C = 1,...,M C IOPT(5) - LUPRI Print output unit: C The unit where the output, controlled C by IOPT(4), is written to. C If set to 0, output unit will be 6. C NRW Dimension of real workspace RW C If IOPT(3) = 0: C NRW.GE.N*N*(M+5)+10*M*N+10*N+M C If IOPT(3) = 1: C NRW.GE.N*N*(M+1)+12*N*M+4*N+M-1+LICN C with: C NZ = N*N*(M+1)+N*(M-1) C LICN = 2*NZ C LIRN = DMIN1(DMAX1(1.5*NZ,NZ+4*M*N), C LICN) C RW(NRW) Real workspace C C NIW Dimension of integer workspace IW C If IOPT(3) = 0: C NIW.GE.2*N*N+4*N C If IOPT(3) = 1: C NIW.GE.2*N*N+3*N+8*M*N+8*M*N C +2*NZ+LICN+LIRN C IW(NIW) Integer workspace C C C C Output parameters: C ================== C C X(N,M) Solution data ( or final data, C respectively ) C INFO Information output parameter C >0 Number of iterations performed C to obtain the solution C <0 BVPSOL termination C -1 If IOPT(3) = 0: C Iteration stops at stationary C point C If IOPT(3) = 1: C Gaussian elimination failed due C to singular Jacobian C -2 Iteration stops after ITMAX C iteration steps ( as indicated C by option ITMAX=IOPT(1) ) C -3 Integrator failed to complete C the trajectory computation C -4 Gauss Newton method failed to C converge C -5 Given initial values C inconsistent with separable C linear boundary conditions C -6 If IOPT(3) = 0: C Iterative refinement failed to C converge C If IOPT(3) = 1: C Termination since Multiple C Shooting condition or C condition of Jacobian is too bad C -7 Reliable relative accuracy C greater than 1.0D-2 C -8 Condensing algorithm for C linear block system fails, use C global linear solver in C boundary value problem routine, C i.e. set IOPT(3)=1 C -9 If IOPT(3) = 1: C Sparse linear solver failed, C possibly workspace is C too small C -10 Real or integer work-space C exhausted C -11 If IOPT(3) = 0: C Rank reduction failed - C resulting rank is zero C C ------------------------------------------------------------ C C* End Prologue C: SMALL = squareroot of "smallest positive machine number C divided by relative machine precision" DOUBLE PRECISION SMALL PARAMETER (SMALL=4.94D-32) INTEGER M1,NM,NM1,NN,NRW,NIW,ITMAX,NONLIN,IBVPSL DOUBLE PRECISION RELDIF,TOL,XTHR,V C: Begin C ------------------------------------------------------------ C 1 Internal parameters C Standard values fixed below C Monitor output unit LUMON = IOPT(5) IF (LUMON.LE.0.OR.LUMON.GE.100) LUMON=6 C Scaling threshold XTHR = SMALL C Prescribed relative precision for numerical integration TOL = EPS*1.0D-2 C Prescribed relative deviation for numerical differentiation RELDIF = DSQRT(TOL) ITMAX = IOPT(1) NONLIN = IOPT(2) IBVPSL = IOPT(3) INFO = IOPT(4) IF(INFO.GE.0)THEN C Print BVPSOL heading lines 1 FORMAT('B V P S O L',2X,5('*'),2X,'V e r s i o n',2 * X,'1 . 2',1X,3('*'),//,1X,'Newton',1('-'),'Method ','for ', * 'the ','solution ','of ','boundary ','value ','problems',/ * /) write ( lumon, '(a)' ) ' ' WRITE(LUMON,1) ENDIF IF (IBVPSL.EQ.0) THEN C Starting value for pseudo - rank of sensitivity matrix E IRANK = N C Initial preparations M1 = M-1 NN = N*N NM = N*M NM1 = N*M1 C: WorkSpace: IW L4=1 L5=L4+N L6=L5+N L7=L6+N L8=L7+N L9=L8+N*N L10=L9+N*N NIW=L10-1 C. End WorkSpace at NIW C: WorkSpace: RW L11=1 L12=L11+N*N*M1 L13=L12+N*N L14=L13+N*N L15=L14+N*N L16=L15+N*N L17=L16+N*N L18=L17+NM L19=L18+NM L20=L19+NM L21=L20+NM L22=L21+NM L23=L22+NM L24=L23+NM1 L25=L24+NM1 L26=L25+NM1 L27=L26+NM1 L28=L27+N L29=L28+N L30=L29+N L31=L30+N L32=L31+N L33=L32+N L34=L33+N L35=L34+N L36=L35+N L37=L36+N L38=L37+N L39=L38+N L40=L39+N L41=L40+M L42=L41+N L43=L42+N*N NRW=L43-1 C. End WorkSpace at NRW C ------------------------------------------------------------ C 2 Check for sufficient real/integer workspace IF (INFO.GE.0) THEN 2 FORMAT('Minimal ','required ','work-space ',':',/, * 'Real ','array ','RW(',I4,')',/,'Integer ','array ', * 'IW(',I4,')') write ( lumon, '(a)' ) ' ' WRITE(LUMON,2)NRW,NIW ENDIF IF(NRW.LE.IRW.AND.NIW.LE.IIW)THEN CALL BVPL(FCN,BC,IVPSOL,N,M,M1,NM,NM1,T,X,EPS,TOL,RELDIF, * NONLIN,IRANK,ITMAX,INFO,XTHR,IW(L4),IW(L5),IW(L6),IW(L7), * IW(L8),IW(L9),RW(L11),RW(L12),RW(L13),RW(L14),RW(L15),RW( * L16),RW(L17),RW(L18),RW(L19),RW(L20),RW(L21),RW(L22),RW( * L23),RW(L24),RW(L25),RW(L26),RW(L27),RW(L28),RW(L29),RW( * L30),RW(L31),RW(L32),RW(L33),RW(L34),RW(L35),RW(L36),RW( * L37),RW(L38),RW(L39),RW(L40),RW(L41),RW(L42),LUMON) ELSE C Fail exit work-space exhausted IF(INFO.GE.0.AND.NRW.GT.IRW)THEN 3 FORMAT('Error: ','real ','work ','- ','space ', * 'exhausted',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,3) ENDIF IF(INFO.GE.0.AND.NIW.GT.IIW)THEN 4 FORMAT('Error: ','integer ','work ','- ','space ', * 'exhausted',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,4) ENDIF INFO = -10 ENDIF ELSE IF (IBVPSL.EQ.1) THEN C Initial preparations M1 = M-1 NN = N*N NM = N*M NM1 = N*M1 NMX8 = 8*NM NZ = N*N*(M+1)+N*(M-1) LICNQ = 2*NZ LIRNQ = MAX0(3*NZ/2,NZ+4*M*N) LIRNQ = MIN0(LIRNQ,LICNQ) V = DBLE(LIRNQ)/DBLE(LICNQ) NRW = N*N*(M+1)+12*M*N+4*N+M-1 NI2W = 2*NZ+8*M*N NIW = NMX8+2*N*N+3*N L1 = IRW-NRW L2 = IDINT(DBLE(IIW-NIW-NI2W)/(V+1.0D0)) LICN = MIN0(L1,L2) LIRN = IDINT(V*DBLE(LICN)) LISNQ = LICNQ+LIRNQ NKEEP = NMX8 MINI2W = NI2W + LISNQ C: WorkSpace: I2W L4=1 L5=L4+LIRN L6=L5+LICN L7=L6+NZ L8=L7+NZ L9=L8+NKEEP NI2W=L9-1 C. End WorkSpace at NI2W C: WorkSpace: IW L10=L9 L11=L10+N L12=L11+N L13=L12+N L14=L13+N*N L15=L14+N*N L16=L15+NMX8 NIW=L16-1 C. End WorkSpace at NIW C: WorkSpace: RW L17=1 L18=L17+N*N*M1 L19=L18+N*N L20=L19+N*N L21=L20+LICN L22=L21+NM L23=L22+NM L24=L23+NM L25=L24+NM L26=L25+NM L27=L26+NM L28=L27+NM1 L29=L28+NM1 L30=L29+NM1 L31=L30+NM1 L32=L31+N L33=L32+N L34=L33+N L35=L34+N L36=L35+NM L37=L36+NM L38=L37+N L39=L38+N L40=L39+N L41=L40+N L42=L41+M1 NRW=L42-1 C. End WorkSpace at NRW C ------------------------------------------------------------ C 2 Check for sufficient real/integer workspace IF (INFO.GE.0) THEN 5 FORMAT('Minimal ','required ','work-space ',':',/, * 'Real ','array ','RW( ',I5,')',/, * 'Integer ','array ','IW( ',I5,')') write ( lumon, '(a)' ) ' ' WRITE(LUMON,5)NRW-LICN+LICNQ,NIW-LIRN-LICN+LISNQ ENDIF IF(NRW.LE.IRW.AND.NIW.LE.IIW * .AND.LICNQ.LE.LICN+1.AND.LIRNQ.LE.LIRN+1)THEN CALL BVPG(FCN,BC,IVPSOL,N,M,M1,NM,NM1,NMX8,NZ,LICN,LIRN, * LISNQ,NKEEP,T,X,EPS,TOL,RELDIF,NONLIN,ITMAX,INFO,XTHR,IW( * L10),IW(L11),IW(L12),IW(L13),IW(L14),IW(L15),RW(L17),RW( * L18),RW(L19),RW(L20),RW(L21),RW(L22),RW(L23),RW(L24),RW( * L25),RW(L26),RW(L27),RW(L28),RW(L29),RW(L30),RW(L31),RW( * L32),RW(L33),RW(L34),RW(L35),RW(L36),RW(L37),RW(L38),RW( * L39),RW(L40),RW(L41),IW(L6),IW(L7),IW(L5),IW(L4),IW( * L8),LUMON) ELSE C Fail exit work-space exhausted IF(INFO.GE.0.AND.NRW.GT.IRW)THEN 6 FORMAT('Error: ',A,'work ','- ', * 'space ','exhausted',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,6) 'real ' ENDIF IF(INFO.GE.0.AND.NIW.GT.IIW)THEN WRITE(LUMON,6) 'integer ' ENDIF INFO = -10 ENDIF ELSE WRITE(LUMON,*) ' Invalid option IOPT(3) set - must be 0 or 1' ENDIF RETURN END SUBROUTINE BVPL(FCN,BC,IVPSOL,N,M,M1,NM,NM1,T,X,EPS,TOL, *RELDIF,NONLIN,IRANK,ITMAX,INFO,XTHR,IROW,ICOL,ICOLB,PIVOT,IA, *IB,G,A,B,BG,E,QE,DX,DDX,DXQ,DXQA,XA,XW,XU,HH,DHH,HHA,D,DE,R, *DR,RA,U,DU,QU,X1,XM,T1,T2,DX1,RF,US,EH,LUMON) c*********************************************************************72 c cc BVPL() carries out the local linear solver approach. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) EXTERNAL FCN,BC,IVPSOL INTEGER N,M,M1,NM,NM1,LUMON DOUBLE PRECISION T(M),X(NM) DOUBLE PRECISION EPS DOUBLE PRECISION TOL,RELDIF INTEGER NONLIN,ITMAX,IRANK INTEGER INFO DOUBLE PRECISION XTHR INTEGER IROW(N),ICOL(N),ICOLB(N),PIVOT(N) INTEGER IA(N,N),IB(N,N) DOUBLE PRECISION G(N,N,M1) DOUBLE PRECISION A(N,N),B(N,N),BG(N,N),E(N,N),QE(N,N) DOUBLE PRECISION DX(NM),DDX(NM),DXQ(NM),DXQA(NM),XA(NM),XW( *NM),XU(NM1),HH(NM1),DHH(NM1),HHA(NM1),D(N),DE(N),R(N),DR(N), *RA(N),U(N),DU(N),QU(N),X1(N),XM(N),T1(N),T2(N),DX1(N),RF(M), *US(N) DOUBLE PRECISION EH(N,N) C C Additional dimensional integer variables C ======================================== C C M1 M-1 C NM N*M C NM1 N*(M-1) C C Internal real arrays (workspace) : C ================================== C C G(N,N,M1) (N,N) -Wronskian Matrices G(1),...,G(M-1) C . C A(N,N) Wronskian Matrix on left boundary C dBC/dX(X(1,...,N),T(1)). C B(N,N) Wronskian Matrix on right boundary C dBC/dX(X((N-1)*M+1,...,N*M),T(M)). C BG(N,N) Workspace for subroutine RHS1. Holds C subsequently the matrices B,B*G(M-1),..., C B*G(M-1)*...*G(2)during computation of C the right hand side of the condensed C linear system. C E(N,N) Sensitivity Matrix of the boundary value C problem: E = A+B*G(M-1)*...*G(1). C EH(N,N) Holds a copy of the row- and column C scaled, but not decomposed sensitivity C matrix E needed for iterative refinement C computations. C QE(N,N) Workspace for DECCON and SOLCON to hold C the updating part of the pseudoinverse C of E in case of it's rank defeciency. C D(N) Diagonal elements of upper triangular C matrix of QR-decomposition computed by C DECCON . C DDX(NM) Workspace for subroutine BLSOLI - gets C the iterative refinement vector of DXQ . C DE(N) Holds row scaling factors for the C sensitivity matrix. C DHH(NM1) Holds the recursive refinement of HH C computed in BLSOLI . C DR(N) Workspace for subroutine BLSOLI to hold C the boundary residual C BC(DXQ(1,...,N),DXQ((M-1)*N+1,...,M*N))+ C (A*DXQ(1,...,N))+B*DXQ((M-1)*N+1,...,M*N) C . C DU(N) Workspace for subroutine BLSOLI to hold C the right hand side of the linear system C E*T2 = U-E*DX1 solved to compute the C iterative refinement for DX1 . C DX(NM) Actual newton correction. C DXQ(NM) Simplified Newton correction J(k-1)*X(k) C with the Jacobian J(k) and the iterate C vector X(k) at the k-th iterate. C DXQA(NM) Previous simplified Newton correction C J(k-2)*X(k-1). C DX1(N) Workspace to receive the solution output C of SOLCON within computation of DXQ and C it's iterative refinement. C HH(NM1) Elements (J-1)*N+1 to J*N are holding C the values C Y(T(J+1),X((J-1)*N+1,...,J*N))-X(J*N+1, C ...,(J+1)*N) C ( with the trajectory Y in C ³ T(J),T(J+1)ü , J = 1,...,M-1 ). C HHA(NM1) Holds the previous values of HH . C QU(N) Savespace to hold U(N)for restoring it C in the case of rank reduction. C R(N) Value of the boundary condition function C BC for the current iterate. C RA(N) Previous values of R . C RF(M) Workspace for subroutine BLSOLI - R(J) C gets the maximum norm of the C components (J-1)*N,...,J*N of the C iterative refinement vector to DXQ for C determination of the sweep index for C subsequent iterations. C T1(N) Workspace used for miscellaneous C purposes temporarely. C T2(N) Workspace used for miscellaneous C purposes temporarely. C U(N) Holds the right hand side of the C condensed linear system computed by C subroutine BLRHS1 . C US(N) Workspace for subroutine BLSOLI to save C U(N). C XA(NM) Previous Newton iterate. C XU(NM1) Elements (J-1)*N+1 to J*N are holding C the values Y(T(J+1),X((J-1)*N+1,...,J*N)) C of the trajectory in the interval ³ T(J), C T(J+1)ü , (for J = 1,...,M-1 ). C XW(NM) Scaling factors for iteration vector. C X1(N) Components of the iteration vector C corresponding to the left boundary C A = T(1). C XM(N) Components of the iteration vector C corresponding to the right boundary C B = T(M). C C Internal integer arrays (workspace) C =================================== C C IROW(N) Row permutations of sensitivity matrix. C ICOL(N) Column permutations of matrix A C (left boundary). C ICOLB(N) Column permutations of matrix B C (right boundary). C PIVOT(N) Workspace for DECCON and SOLCON to hold C permutations performed during C QR-decomposition. C IA(N,N) Reflects the sparse structure of matrix C A by values 0, 1. C IB(N,N) Reflects the sparse structure of matrix C B by values 0, 1. C C Internal real variables C ======================= C C COND Gets the condition of the sensitivity C matrix determined by DECCON . Not C further used. C CONDE Maximum permitted subcondition of C sensitivity matrix E . C COND1 Condition of boundary conditions part of C the decomposed sensitivity matrix E . C COND2 Condition of the remaining part of the C sensitivity matrix E . C CONV Scaled maximum norm of DXQ computed by C subroutine BLLVLS . Used for convergence C test. C CONVA Holds the previous value of CONV . C DEL Becomes the Euklid's normsquare of the C defect belonging to the condensed linear C system in case of rank defeciency. Used C for computation of damping factor C predictor. C EPSMIN Smallest reasonable permitted accuracy C EPS that can be prescribed by the user. C FC Actual Gauss Newton iteration damping C factor. C FCA Previous Gauss Newton iteration damping C factor. C FCDNM Used to compute the denominator of the C damping factor FC during computation of C it's predictor, corrector and C aposteriori estimate (in the case of C performing a Rank1 update) . C FCH Temporarely used for storing the new FC C when computing aposteriori estimate. C FCMIN Minimum permitted relaxation factor. If C FC becomes smaller than this value, one C of the following may occur: C a. Recomputation of the sensitivity C matrix by means of difference C approximation (instead of Rank1 C update), if Rank1 - update C previously was used C b. Rank reduction of sensitivity C matrix E , if difference C approximation was used previously C and Rank(E).NE.0 C c. Fail exit otherwise C FCMINH DSQRT(FCMIN). Used for rank C determination computations. C FCMIN2 FCMIN**2 . Used for FC-predictor C computation. C FCNUM Gets the numerator of the aposteriori C estimate of FC . C FCNUMK Gets the numerator of the corrector C computation of FC . C FCNUMP Gets the numerator of the predictor C computation of FC . C H Actual integrator stepsize. C HMAX Maximum permitted integrator stepsize. C Set to the length of the integration C interval, e.g. the distance of the C effected Shooting points. C HSAVE Stepsize saved across the call of the C integrator. C HSTART Start stepsize for integration used by C subroutines BLFCNI and BLDERG . C MUE Temporary value used during computation C of damping factors predictor. C REDH Multi purpose reduction factor. (???) C RELDIF Relative deviation for numerical C differentation. C SENS1 Sensitivity of boundary conditions part C of the decomposed sensitivity matrix E . C SENS2 Sensitivity of the remaining part of the C matrix E . C SIGMA Decision parameter for Jacobian Rank1 C updates (SIGMA.GT.1) . Rank1 updates are C inhibited, if SIGMA.GT.1/FCMIN is set. C SKAP Used to compute and print out the C incompatibility factor of the nonlinear C boundary value (e.g. least squares) C problem. C SUMF Standard level of the current iterate, C e.g. Norm2(F(X))**2 C with the nonlinear model function F on C which Newton iteration is performed, C arising from the Multiple Shooting C approach. C SUMX Natural level of the current iterate, C e.g. Norm2(DX) C with the Newton correcture DX C (see above). C SUMXA Natural level of the previous iterate. C TFAIL Used to get and print out in case of an C integrator failure the last reached T C value as a proposal for insertion of a C new Shooting point. C TOL Prescribed relative precision for C numerical integration. C TOLH Temporary used for computation of TOL C (may be obmitted|). C TOLMIN Lower bound value for TOL . C XTHR Threshold for scaling. C TJ Used by BLFCNI to hold T(J). C TJ1 Used by BLFCNI to hold T(J+1). C CORR Used by BLSOLI to compute RF(J), J = 1, C ...,M . C EPX1 Used by BLSOLI to get the maximum norm C of DX1 corrected by iterative C refinement. C EPDX1 Used by BLSOLI to get the maximum norm C of the iterative refinement correcture C for DX1 . C EPX1H Relative accuracy of DX1 after iterative C refinement (= EPDX1/EPX1 ). If EPX1H.GT. C 1/2 , refinement is considered to be C failed. C EPH Used by BLSOLI as tolerance for the C iterative refinement vectors maximum C norm. If it exceeds at index J the C tolerance, refinement will be performed C next time for the components associated C to Shooting points T(J),...,T(M). C SIGDEL Used by BLSOLI to compute the required C integrator accuracy from the multiple C shooting condition. C SIGDLH Used by BLSOLI to compute the multiple C shooting condition C Max( RF(J+1)/RF(J) ) , ( J = 1,...,M-1 C ). See above for RF . C C Internal integer variables C ========================== C C IC Permutated index. Used by BLSOLI . C ICA Temporarely used during search for C separable linear boundary conditions at C T(1). C ICB Temporarely used during search for C separable linear boundary conditions at C T(M). C IH Temporarely used during search for C separable linear boundary conditions at C T(1) and T(M). C IR Temporary storage for a row permutation C index. C IRANK Rank of sensitivity matrix E of current C iteration step. C IRANKA Rank of sensitivity matrix E of previous C iteration step. C IRANKB Rank of the decomposed sensitivity C matrix part belonging to the boundary C conditions. C IREPET Parameter of subroutines DECCON and C SOLCON indicating the mode of operation C of these routines: C >=0 Full QR-decomposition required C <0 Rank reduction required C IRKMAX Holds the maximum applied rank of all C previous iterates. Must be necessary C .EQ. NE for convergence. C IS Additional DO loop index. C ISUM Used for determination of sparse C structure of matrices A and B as C nonzeros location counter. C ITER Iteration count. C JA Previous sweep index. Used in subroutine C BLSOLI . C JIN Current sweep index. Used in subroutine C BLSOLI . C JJ Used as "reverse DO loop" index: C JJ = IUPB-J in a loop like DO J = 1,IUPB C ... C JN New sweep index determined in subroutine C BLSOLI . C JRED Damping factor reduction count during an C iterate. C KC Temporary storage for a column C permutation index. C KFLAG Gets the subintervall number of the C failure from subroutine BLDERG C if the integrator failed. C KOUNT Trajectory evaluations count. C KPRINT Print parameter - copy of input C parameter INFO . C LEVEL Flow control parameter needed by C subroutine BLSOLI : C 0 indicates computation of Gauss C Newton correcture, C 1 indicates computation of C simplified Gauss Newton correcture C (after computation of the C preliminary new iterate) C NB Number of separable boundary conditions C at T(M) C NE Number of not separable boundary C conditions at T(1) (and number of rows C and columns of sensitivity matrix) C NEW Count of subsequent performed Rank1 C (Broyden) updates. C NY Iterative refinement sweep count for an C iterate ( used in subroutine BLSOLI ) . C NYMAX Highest allowed iterative refinement C sweep index. C: End Parameter C: EPMACH = relative machine precision DOUBLE PRECISION EPMACH,SMALL DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) DOUBLE PRECISION HALF PARAMETER (HALF=0.5D0) DOUBLE PRECISION REDH PARAMETER (REDH=1.0D-2) DOUBLE PRECISION ONE PARAMETER (ONE=1.0D0) DOUBLE PRECISION TWO PARAMETER (TWO=2.0D0) DOUBLE PRECISION EIGHT PARAMETER (EIGHT=8.0D0) DOUBLE PRECISION TEN PARAMETER (TEN=1.0D1) DOUBLE PRECISION FCMIN PARAMETER (FCMIN=1.0D-2) DOUBLE PRECISION FCMIN2 PARAMETER (FCMIN2=1.0D-4) DOUBLE PRECISION FCNLIN PARAMETER (FCNLIN=1.0D-2) DOUBLE PRECISION SIGMA PARAMETER (SIGMA=2.0D0) INTEGER I,ICA,ICB,IH,IR,IRANKA,IRANKB,IREPET,IRKMAX,IS, *ISUM,ITER,I0,J,JJ,JN,JRED,J1,K,KC,KFLAG, *KOUNT,KPRINT,L,LEVEL,NB,NE,NEW,NY,NYMAX DOUBLE PRECISION COND,COND1,COND2,CONV,CONVA,DEL,EPH, *EPSMIN,FC,FCA,FCDNM,FCH,FCMINH,FCNMP2,FCNUM,FCNUMK,FCNUMP, *HSTART,MUE,S,SENS1,SENS2,SIGDEL,SIGDLH,SKAP, *SMALIN,ST,SUMF,SUMX,SUMXA,TFAIL,TH,TOLH,TOLMIN, *EPX1H,CONDE LOGICAL DIFAPP,FCOMPT,JACRFR,NEXT,REDUCT INTEGER L1,L2 DOUBLE PRECISION S1 C: Begin C: Begin of Segment Bvpsol.Body C ---------------------------------------------------------- C 1 Initialization C ---------------------------------------------------------- C 1.1 Internal parameters C Standard values fixed below C Minimum relative precision of integrator ( to be adapted ) CALL ZIBCONST(EPMACH,SMALL) TOLMIN = EPMACH*TEN*TEN C Maximum permitted number of iterative refinements sweeps C Maximum permitted sub - condition number of senitivity C matrix E CONDE = ONE/EPMACH NYMAX = M-1 C FCMINH = DSQRT(FCMIN) C ---------------------------------------------------------- C 1.1.1 Common parameters C Starting value of relaxation factor (FCMIN.LE.FC.LE.1.0) IF(NONLIN.LE.1)THEN C for linear or mildly nonlinear problems FC = ONE ELSE C for highly nonlinear problems FC = FCNLIN ENDIF C Starting value for pseudo - rank of matrix A IRANK = N C Minimum reasonable value for EPS EPSMIN = DSQRT(TEN*EPMACH) IF(EPS.LT.EPSMIN) EPS = EPSMIN C ---------------------------------------------------------- C 1.2 Initial preparations IF(FC.LT.FCMIN) FC = FCMIN IF(FC.GT.ONE) FC = ONE KPRINT = INFO ITER = 0 KOUNT = 0 IREPET = 0 INFO = -1000 FCA = FC CONV = ONE JACRFR = .FALSE. C: Begin SetVec.Vec DO 5 L1=1,NM XA(L1)=X(L1) 5 CONTINUE C. End SetVec.Vec IF(TOL.LE.ZERO) TOL = EPS/TEN IF(TOL.LT.TOLMIN) TOL = TOLMIN DIFAPP = .TRUE. HSTART =(T(2)-T(1))*REDH SENS1 = ZERO COND1 = ONE SENS2 = ZERO COND2 = ONE IRKMAX = 0 IRANKB = 0 SIGDLH = ZERO SUMF = ONE C: Mat IA = Scalar (Rows 1,N ; Cols 1,N) L1 = 0 DO 6 L2=1,N DO 6 L44=1,N IA(L2,L44)=L1 6 CONTINUE C. End SetIntMat.S C: Mat IB = Scalar (Rows 1,N ; Cols 1,N) L1 = 0 DO 7 L2=1,N DO 7 L44=1,N IB(L2,L44)=L1 7 CONTINUE C. End SetIntMat.S C: CubeMat G (layer 1)= Scalar (Rows 1,N ; Cols 1,N) S1 = ZERO DO 8 L1=1,N DO 8 L2=1,N G(L1,L2,1)=S1 8 CONTINUE C. End SetCubeMat.S IF(KPRINT.GE.0)THEN C Print Start vector data, predescribed precision and max C iteration steps 9 FORMAT('Initial ','data',//) write ( lumon, '(a)' ) ' ' WRITE(LUMON,9) DO 10 J=1,M 11 FORMAT(D13.5,2X) WRITE(LUMON,11)T(J) 12 FORMAT((14X,3(D20.10,1X))) WRITE(LUMON,12)(X(L1),L1=(J-1)*N+1,J*N) 10 CONTINUE 13 FORMAT('N ','=',I2,2X,'M ','=',I2,/, * 'Prescribed ','relative ','precision',D10.2,2X,/, * 'Maximum ','permitted ','number ','of ','iteration ', * 'steps',1X,I3,//,'1') write ( lumon, '(a)' ) ' ' WRITE(LUMON,13)N,M,EPS,ITMAX IF(KPRINT.EQ.0)THEN 14 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,14) 15 FORMAT(4X,'It',4X,'Ny',7X,'Levelf',10X,'Levelx',8X, * 'Rel.Fc.') write ( lumon, '(a)' ) ' ' WRITE(LUMON,15) ENDIF ENDIF C ---------------------------------------------------------- C 1.3 Startup step C ---------------------------------------------------------- C 1.3.1 Computation of the residual vector c CALL BLFCNI(IVPSOL,FCN,BC,N,M,NM,NM1,ITER,KPRINT,HSTART, * FCMIN,T,X,X1,XM,T1,XU,HH,R,TOL,FC,FCOMPT,REDUCT,KFLAG, * KOUNT,INFO,LUMON) C C Main iteration loop C =================== C C: While (expression) 16 IF(INFO.EQ.-1000)THEN C: Begin of Segment Bvpsol.Core C ---------------------------------------------------- C 2 Startup of iteration step IF(.NOT.JACRFR)THEN LEVEL = 0 C -------------------------------------------------- C 2.1 Scaling of variables X(NM) CALL BLSCLE(N,M,NM,NM1,X,XU,XW,XTHR) IF(ITER.NE.0)THEN IRANKA = IRANK C: Begin SetVec.Vec DO 17 L1=1,NM DXQA(L1)=DXQ(L1) 17 CONTINUE C. End SetVec.Vec C: FCNUM = Sum of Formula Elements (for 1,NM) FCNUM = 0.0D0 DO 18 L1=1,NM FCNUM=FCNUM+((DX(L1)/XW(L1))**2) 18 CONTINUE C. End MakeSum.Comp C: FCNMP2 = Sum of Formula Elements (for 1,NM) FCNMP2 = 0.0D0 DO 19 L1=1,NM FCNMP2=FCNMP2+((DXQ(L1)/XW(L1))**2) 19 CONTINUE C. End MakeSum.Comp FCNUMP = FCNUM*FCNMP2 ENDIF IF(ITER.NE.0)THEN IF(IRANK.GE.NB.AND.FC.GT.FCMINH) IRANK = NE IF(IRANK.LT.NB) IRANK = NB TH = FC-ONE C: FCDNM = Sum of Formula Elements (for 1,NM) FCDNM = 0.0D0 DO 20 L1=1,NM FCDNM=FCDNM+(((DXQ(L1)+TH*DX(L1))/XW(L1))**2) 20 CONTINUE C. End MakeSum.Comp FCH = DSQRT(FCNUM/FCDNM)*FC*FC*HALF C ------------------------------------------------ C 2.1.1 Decision criterion for Jacobian updating C technique: C DIFAPP.EQ..TRUE. numerical C differentiation, C DIFAPP.EQ..FALSE. rank1 updating DIFAPP = FC.LT.FCA.AND.NEW.GT.0.OR.FCH.LT.FC* * SIGMA.OR.IRANKA.LT.IRANK.OR.EPH*REDH.GT.EPS FCA = FC IF(NONLIN.GT.0) FC = DMIN1(FCH,ONE) ENDIF C -------------------------------------------------- C 2.2 Difference approximation of jacobian matrix A C ( If Difapp.EQ..TRUE. ) or C Rank-1 update of jacobian matrix A ( If Difapp C .EQ..FALSE. ) CALL BLDERA(BC,N,M,NM,XW,X1,XM,R,T2,A,B,RELDIF) C -------------------------------------------------- C 2.3 Determination of sparse structure of matrices C A and B and determination of internal row C scaling of sensitivity matrix E ISUM = 0 DO 21 I=1,N S = ZERO DO 22 K=1,N TH = DABS(A(I,K))*XW(K) IF(S.LT.TH) S = TH TH = DABS(B(I,K))*XW(K+NM1) IF(S.LT.TH) S = TH 22 CONTINUE IF(S.LT.XTHR) S = XTHR DE(I)=SMALL/S DO 23 K=1,N IF(IA(I,K).LE.0)THEN IF(A(I,K).NE.ZERO)THEN IA(I,K)=1 ISUM = 1 ENDIF ENDIF IF(IB(I,K).LE.0)THEN IF(B(I,K).NE.ZERO)THEN IB(I,K)=1 ISUM = 1 ENDIF ENDIF 23 CONTINUE 21 CONTINUE IF(ISUM.NE.0)THEN C ------------------------------------------------ C 2.3.1 Determination of row and column C permutation vectors DO 24 I=1,N ICOL(I)=I ICOLB(I)=I IROW(I)=I 24 CONTINUE C ------------------------------------------------ C 2.3.2 Search for separable linear boundary C conditions at T(1) NE = N DO 25 I=1,N DO 26 K=1,N IF(IB(I,K).NE.0) GOTO 9996 26 CONTINUE ISUM = 0 DO 27 K=1,N IF(IA(I,K).NE.0)THEN ISUM = ISUM+1 ICA = K ENDIF 27 CONTINUE IF(ISUM.LE.1)THEN DO 28 IS=1,N IH = ICOL(IS) IF(IH.EQ.ICA) ICOL(IS)=ICOL(NE) IH = IROW(IS) IF(IH.EQ.I) IROW(IS)=IROW(NE) 28 CONTINUE ICOL(NE)=ICA IROW(NE)=I NE = NE-1 IF(DABS(R(I)).GT.TEN*EPMACH*DABS(X(ICA))) * THEN INFO = -5 GOTO 9998 ENDIF ENDIF 9996 CONTINUE 25 CONTINUE IF(KPRINT.GE.0.AND.NE.EQ.0)THEN 29 FORMAT('Warning: ','attempt ','to ', * 'solve ','initial ','value ','problem') write ( lumon, '(a)' ) ' ' WRITE(LUMON,29) ENDIF IF(IRANK.GT.NE) IRANK = NE IRANKA = IRANK C ------------------------------------------------ C 2.3.3 Search for separable linear boundary C conditions at T(M) NB = 0 ENDIF IF(ISUM.NE.0.AND.NE.NE.0)THEN DO 30 I=1,NE IR = IROW(I) DO 31 K=1,N IF(IA(IR,K).NE.0) GOTO 9995 31 CONTINUE ISUM = 0 DO 32 K=1,N IF(IB(IR,K).NE.0)THEN ISUM = ISUM+1 ICB = K ENDIF 32 CONTINUE IF(ISUM.LE.1)THEN NB = NB+1 DO 33 IS=1,N IH = ICOLB(IS) IF(IH.EQ.ICB) ICOLB(IS)=ICOLB(NB) 33 CONTINUE ICOLB(NB)=ICB IROW(I)=IROW(NB) IROW(NB)=IR IF(DABS(R(IR)).GT.TEN*EPMACH*DABS(X(ICB+ * NM1)))THEN INFO = -5 GOTO 9998 ENDIF ENDIF 9995 CONTINUE 30 CONTINUE IF(KPRINT.GE.0.AND.NB.EQ.N)THEN 34 FORMAT('Warning: ','attempt ','to ', * 'solve ','initial ','value ','problem') write ( lumon, '(a)' ) ' ' WRITE(LUMON,34) ENDIF C Initial rank strategy for highly nonlinear C problems IF(NB.LT.NE.AND.ITER.EQ.0.AND.NONLIN.GT.2) IRANK * = NB ENDIF ENDIF JACRFR = .FALSE. IF(DIFAPP)THEN NEW = 0 KFLAG = 0 CALL BLDERG(FCN,N,NE,M,M1,NM,NM1,T,X,XU,XW,T2, * TFAIL,G,ICOL,IVPSOL,HSTART,TOL,RELDIF,KFLAG) IF(KFLAG.LT.0)THEN INFO = -3 GOTO 9998 ENDIF IF(M.GT.2) KOUNT = KOUNT+N IF(M.EQ.2) KOUNT = KOUNT+NE ELSE NEW = NEW+1 CALL BLRK1G(N,M,M1,NM,NM1,XW,DX,HH,HHA,T1,G,FCA) ENDIF C ---------------------------------------------------- C 2.3.4 Computation of sensitivity matrix E =-A+B*G(M1 C *..*G(1)) C (projections included) IF(IRANK.NE.0)THEN DO 35 I=1,NE IR = IROW(I) C: Mat E(Row I)= Mat B(Row IR)* Scalar (for 1,N) S1 = DE(IR) DO 36 L1=1,N E(I,L1)=B(IR,L1)*S1 36 CONTINUE C. End SetRow.RowxS 35 CONTINUE DO 37 JJ=1,M1 J = M-JJ DO 38 I=1,NE DO 39 K=1,N S = ZERO DO 40 L=1,N S = S+E(I,L)*G(L,K,J) 40 CONTINUE T1(K)=S 39 CONTINUE DO 41 K=1,N E(I,K)=T1(K) 41 CONTINUE 38 CONTINUE 37 CONTINUE C -------------------------------------------------- C 2.4 Prepare solution of the linear system C -------------------------------------------------- C 2.4.1 internal row and column scaling of matrix A C INTERNAL ROW AND COLUMN SCALING AND PERMUTATION OF C MATRIX E DO 42 K=1,NE KC = ICOL(K) S = XW(KC) DO 43 I=1,NE IR = IROW(I) E(I,K)=-(A(IR,KC)*DE(IR)+E(I,KC))*S 43 CONTINUE 42 CONTINUE C -------------------------------------------------- C 2.4.2 Save matrix E on EH C: Mat EH = Mat E (Rows 1,NE ; Cols 1,NE) DO 44 L1=1,NE DO 44 L2=1,NE EH(L1,L2)=E(L1,L2) 44 CONTINUE C. End SetMat.Mat ENDIF IRANKB = NB C ---------------------------------------------------- C 2.4.3 Monitor for actually applied maximum rank IF(IRKMAX.LT.IRANK) IRKMAX = IRANK C ---------------------------------------------------- C 2.5 Save values of R(N)and HH((M-1)*N) IF(IREPET.EQ.0)THEN C: Begin SetVec.Vec DO 45 L1=1,N RA(L1)=R(L1) 45 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec DO 46 L1=1,NM1 HHA(L1)=HH(L1) 46 CONTINUE C. End SetVec.Vec ENDIF NEXT = .FALSE. C C Pseudo-rank reduction loop C ========================== C C: DO (Until) 47 CONTINUE C -------------------------------------------------- C 3 Main-part of iteration step C -------------------------------------------------- C 3.1 Solution of the linear system C -------------------------------------------------- C 3.1.1 Constrained QR-decomposition of ( ( COMMA NE C NE ) ) - matrix E COND = CONDE IF(IRANK.GT.0) CALL BLDECC(E,N,N,IRANKB,NE,NE, * IRANK,COND,D,PIVOT,IREPET,QE,T1) IF(NONLIN.EQ.0.AND.IRANK.LT.NE)THEN INFO = -8 GOTO 9998 ENDIF C -------------------------------------------------- C 3.1.2 evaluation of subcondition and sensitivity C numbers COND1 = ONE COND2 = ONE SENS1 = ZERO SENS2 = ZERO IF(IRANKB.NE.0)THEN SENS1 = DABS(D(1)) COND1 = SENS1/DABS(D(IRANKB)) ENDIF IF(IRANKB.NE.IRANK)THEN SENS2 = DABS(D(IRANKB+1)) COND2 = SENS2/DABS(D(IRANK)) ENDIF IF(FCA.GE.1.0D0.AND.FC.GE.1.0D0.AND.ITER.NE.0)THEN IF(IRANKB.NE.IRANK.AND.SENS2.LT.(EPS/REDH)*SMALL) * IRANK = IRANKB IF(IRANKB.NE.0.AND.SENS1.LT.(EPS/REDH)*SMALL) * IRANK = 0 ENDIF C -------------------------------------------------- C 3.1.3 (best) (least squares) solution of linear (N, C N)-system CALL BLSOLI(N,M,M1,NM,NM1,LEVEL,NE,NB,IRANK,IRANKB, * IREPET,NYMAX,KPRINT,EPS,REDH,TOLMIN,TOL,RELDIF,EPH, * EPX1H,SIGDEL,SIGDLH,E,EH,HH,DHH,R,A,B,BG,G,QE,U,QU, * DE,DU,T1,T2,US,DX1,D,DDX,DXQ,XW,DR,RF,IROW,ICOL, * ICOLB,PIVOT,NY,INFO,LUMON) IF(INFO.NE.-1000) GOTO 9998 C -------------------------------------------------- C 3.2 Evaluation of scaled natural level function C SUMX C scaled maximum error norm CONV C evaluation of (scaled) standard level function C SUMF ( SUMF only, if KPRINT.GE.0 ) C and computation of ordinary newton corrections C DX(N) CALL BLLVLS(N,M,NM,NM1,XW,DXQ,HH,R,DE,CONV,SUMX, * SUMF,KPRINT) C: Begin SetVec.Vec DO 48 L1=1,NM DX(L1)=DXQ(L1) 48 CONTINUE C. End SetVec.Vec IF(IREPET.EQ.0.AND.IRANK.NE.0)THEN C: Begin SetVec.Vec DO 49 L1=1,IRANK QU(L1)=U(L1) 49 CONTINUE C. End SetVec.Vec ENDIF C: Begin SetVec.Vec DO 50 L1=1,NM XA(L1)=X(L1) 50 CONTINUE C. End SetVec.Vec SUMXA = SUMX CONVA = CONV C -------------------------------------------------- C 3.3 a - priori estimate of relaxation factor FC JRED = 0 REDUCT = .FALSE. IF(ITER.NE.0.AND.NONLIN.NE.0)THEN IF(NEW.LE.0.AND.(IRANK.GE.NE.OR.IRANKA.GE.NE) * .OR.IREPET.NE.0)THEN C ---------------------------------------------- C 3.3.1 Full rank case (independent of preceding C rank) computation of the denominator of C a-priori estimate C: FCDNM = Sum of Formula Elements (for 1,NM) FCDNM = 0.0D0 DO 51 L1=1,NM FCDNM=FCDNM+(((DX(L1)-DXQA(L1))/XW(L1))**2) 51 CONTINUE C. End MakeSum.Comp IF(IRANK.NE.N)THEN C -------------------------------------------- C 3.4 Rank - deficient case ( if previous rank C was full ) C -------------------------------------------- C 3.4.1 Computation of the projected C denominator of a-priori estimate C: Vec T1 = Scalar (for 1,N) S1 = ZERO DO 52 L1=1,N T1(L1)=S1 52 CONTINUE C. End SetVec.S IF(IRANK.NE.0)THEN DO 53 L=1,NE K = ICOL(L) DX1(L)=DXQA(K)/XW(K) 53 CONTINUE C ------------------------------------------ C 3.4.2 Projection for reduced component DX1 C (NE) CALL BLPRJC(N,NE,IRANK,DEL,DX1,D,T2,QE, * PIVOT) DO 54 L=1,NE K = ICOL(L) T1(K)=DX1(L)*XW(K) 54 CONTINUE ENDIF DO 55 J=1,M1 DO 56 I=1,N S = ZERO DO 57 K=1,N S = S+T1(K)*G(I,K,J) 57 CONTINUE T2(I)=S 56 CONTINUE C: Begin SetVec.Vec DO 58 L1=1,N T1(L1)=T2(L1) 58 CONTINUE C. End SetVec.Vec I0 = J*N DO 59 I=1,N ST = ONE/XW(I+I0) S = T1(I) DEL = DEL+S*ST*ST*(S+(DX(I+I0)-DXQA(I+I0)) * *TWO) 59 CONTINUE 55 CONTINUE FCDNM = FCDNM+DEL ENDIF FCDNM = FCDNM*SUMX C ---------------------------------------------- C 3.4.3 New relaxation factor IF(FCDNM.GE.FCNUMP*FCMIN2)THEN MUE = FCA*DSQRT(FCNUMP/FCDNM) FC = DMIN1(MUE,ONE) ELSE FC = ONE ENDIF ENDIF IREPET = 0 REDUCT = FC.LT.FCMIN ENDIF LEVEL = 1 IF(.NOT.REDUCT)THEN C ---------------------------------------------- C 3.5 Save natural level for later computations C of corrector and print iterate FCNUMK = SUMX IF(KPRINT.GE.0)THEN C Print Standard - and natural level IF(KPRINT.GT.0)THEN 60 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,60) 61 FORMAT(4X,'It',4X,'Ny',7X,'Levelf',10X, * 'Levelx',18X,'New',4X,'Rank') write ( lumon, '(a)' ) ' ' WRITE(LUMON,61) ENDIF 62 FORMAT(4X,I2,4X,I2,5X,D10.3,2X,4X,D10.3 * ,2X,13X,I2,6X,I2) write ( lumon, '(a)' ) ' ' WRITE(LUMON,62)ITER,NY,SUMF,SUMXA,NEW,IRANK IF(KPRINT.GT.0)THEN 63 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,63) ENDIF ENDIF C C Relaxation-factor reduction loop C ================================ C C: DO (Until) 64 CONTINUE C -------------------------------------------- C 3.6 Preliminary new iterate C: DO (Until) 65 CONTINUE FCOMPT = .FALSE. C: Vec X = Vec XA + Vec DX * Scalar (for 1,NM) S1 = FC DO 66 L1=1,NM X(L1)=XA(L1)+DX(L1)*S1 66 CONTINUE C. End SetVec.Vec&VecxS IF(ITER.GT.ITMAX)THEN INFO = -2 GOTO 9997 ENDIF C ------------------------------------------ C 3.6.1 Computation of the residual vector CALL BLFCNI(IVPSOL,FCN,BC,N,M,NM,NM1,ITER, * KPRINT,HSTART,FCMIN,T,X,X1,XM,T1,XU,HH,R, * TOL,FC,FCOMPT,REDUCT,KFLAG,KOUNT,INFO,LUMON) IF(INFO.NE.-1000) GOTO 9997 IF(REDUCT) GOTO 9994 IF(.NOT.(FCOMPT)) GOTO 65 C. UNTIL ( expression - negated above) C -------------------------------------------- C 3.6.2 (best) (least squares) solution of C linear (N,N) -system CALL BLSOLI(N,M,M1,NM,NM1,LEVEL,NE,NB,IRANK, * IRANKB,IREPET,NYMAX,KPRINT,EPS,REDH,TOLMIN, * TOL,RELDIF,EPH,EPX1H,SIGDEL,SIGDLH,E,EH,HH, * DHH,R,A,B,BG,G,QE,U,QU,DE,DU,T1,T2,US,DX1,D, * DDX,DXQ,XW,DR,RF,IROW,ICOL,ICOLB,PIVOT,NY, * INFO,LUMON) IF(INFO.NE.-1000) GOTO 9998 C -------------------------------------------- C 3.6.3 Evaluation of scaled natural level C function SUMX C scaled maximum error norm CONV and C evaluation of (scaled) standard level C function SUMF CALL BLLVLS(N,M,NM,NM1,XW,DXQ,HH,R,DE,CONV, * SUMX,SUMF,KPRINT) C -------------------------------------------- C 3.7 Rank independent convergence test IF(CONV.LE.EPS.AND.IRKMAX.EQ.NE)THEN INFO = 0 GOTO 9997 ENDIF C -------------------------------------------- C 3.8 Natural monotonicity test IF(SUMX.GT.SUMXA)THEN C ------------------------------------------ C 3.9 Output of iterate IF(KPRINT.GE.0)THEN C Print Standard - and natural level, and C damping factor IF(KPRINT.GT.0)THEN 67 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,67) 68 FORMAT(4X,'It',4X,'Ny',7X,'Levelf', * 10X,'Levelx',8X,'Rel.Fc.') write ( lumon, '(a)' ) ' ' WRITE(LUMON,68) ENDIF 69 FORMAT(4X,I2,4X,I2,5X,D10.3,2X,4X,D10.3 * ,2X,4X,F5.3) write ( lumon, '(a)' ) ' ' WRITE(LUMON,69)ITER,NY,SUMF,SUMX,FC IF(KPRINT.GT.0)THEN 70 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,70) ENDIF ENDIF JRED = JRED+1 IF(NONLIN.EQ.0)THEN INFO = -4 GOTO 9997 ENDIF C ------------------------------------------ C 3.10 Compute reduced relaxation factor FC TH = FC-ONE C: FCDNM = Sum of Formula Elements (for 1,NM) FCDNM = 0.0D0 DO 71 L1=1,NM FCDNM=FCDNM+(((DXQ(L1)+TH*DX(L1))/XW(L1)) * **2) 71 CONTINUE C. End MakeSum.Comp FC = DSQRT(FCNUMK/FCDNM)*FC*FC*HALF C Rank reduction, if relaxation factor to C small REDUCT = FC.LT.FCMIN.OR.NEW.GT.0.AND.JRED * .GT.1 ELSE NEXT = .TRUE. ENDIF IF(.NOT.(NEXT.OR.REDUCT)) GOTO 64 C. UNTIL ( expression - negated above) C C End of relaxation-factor reduction loop C ======================================= C ENDIF 9994 CONTINUE IF(.NOT.NEXT)THEN C ------------------------------------------------ C 3.11 Restore former values for repeting C iteration step IREPET = 1 C Restore former values LEVEL = 0 C: Begin SetVec.Vec DO 72 L1=1,N R(L1)=RA(L1) 72 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec DO 73 L1=1,N X1(L1)=XA(L1) 73 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec DO 74 L1=1,N XM(L1)=XA(L1+NM1) 74 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec DO 75 L1=1,NM X(L1)=XA(L1) 75 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec&Vec DO 76 L1=1,NM1 XU(L1)=X(L1+N)+HHA(L1) 76 CONTINUE C. End SetVec.Vec&Vec C: Begin SetVec.Vec DO 77 L1=1,NM1 HH(L1)=HHA(L1) 77 CONTINUE C. End SetVec.Vec IF(KPRINT.GE.0)THEN 78 FORMAT(5X,I2,1X,'Not ','accepted ', * 'relaxation ','factor',5X,F5.3,12X,I2) write ( lumon, '(a)' ) ' ' WRITE(LUMON,78)ITER,FC,IRANK ENDIF IF(ITER.EQ.0)THEN FC = FCMIN ENDIF IF(NEW.GT.0)THEN DIFAPP = .TRUE. JACRFR = .TRUE. IRANK = NE IF(IRANK.LT.NB) IRANK = NB GOTO 9998 ENDIF C ------------------------------------------------ C 3.12 Pseudo-rank reduction IREPET = -1 IF(IRANK.EQ.0)THEN INFO = -11 GOTO 9997 ENDIF C: Begin SetVec.Vec DO 79 L1=1,IRANK U(L1)=QU(L1) 79 CONTINUE C. End SetVec.Vec IRANK = IRANK-1 IF(IRANKB.GT.IRANK) IRANKB = IRANK ENDIF IF(.NOT.(NEXT)) GOTO 47 C. UNTIL ( expression - negated above) C C End of pseudo-rank reduction loop C ================================= C C ---------------------------------------------------- C 4 Preparations to start the following iteration step ITER = ITER+1 IRANKA = IRANK C Preliminary pseudo-rank IF(IRANK.GE.NB.AND.FC.GT.FCMINH) IRANK = NE IF(IRANK.LT.NB) IRANK = NB C ---------------------------------------------------- C 4.1 Print values IF(KPRINT.GE.0)THEN C Print Standard - and natural level, and damping C factor IF(KPRINT.GT.0)THEN 80 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,80) 81 FORMAT(4X,'It',4X,'Ny',7X,'Levelf',10X, * 'Levelx',8X,'Rel.Fc.') write ( lumon, '(a)' ) ' ' WRITE(LUMON,81) ENDIF 82 FORMAT(4X,I2,4X,I2,5X,D10.3,2X,4X,D10.3,2X,4X,F5.3) write ( lumon, '(a)' ) ' ' WRITE(LUMON,82)ITER,NY,SUMF,SUMX,FC IF(KPRINT.GT.0)THEN 83 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,83) DO 84 J=1,M 85 FORMAT(D13.5,2X) WRITE(LUMON,85)T(J) 86 FORMAT((14X,3(D20.10,1X))) WRITE(LUMON,86)(X(L1),L1=(J-1)*N+1,J*N) 84 CONTINUE ENDIF ENDIF 9997 CONTINUE C. End of Segment Bvpsol.Core 9998 CONTINUE GOTO 16 ENDIF C. EndWhile C C End of main iteration loop C ========================== C C ---------------------------------------------------------- C 5 Exits C ---------------------------------------------------------- C 5.1 Solution exit IF(INFO.EQ.0)THEN ITER = ITER+1 C: Vec X = Vec X + Vec DXQ (for 1,NM) DO 87 L1=1,NM X(L1)=X(L1)+DXQ(L1) 87 CONTINUE C. End SetVec.&Vec INFO = ITER IF(KPRINT.LT.0)THEN GOTO 9999 ENDIF IF(KPRINT.GT.0)THEN C Print levels, damping factor of last iteration step 88 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,88) 89 FORMAT(4X,'It',4X,'Ny',7X,'Levelf',10X,'Levelx',8X, * 'Rel.Fc.') write ( lumon, '(a)' ) ' ' WRITE(LUMON,89) ENDIF 90 FORMAT(4X,I2,4X,I2,5X,D10.3,2X,4X,D10.3,2X,4X,F5.3) write ( lumon, '(a)' ) ' ' WRITE(LUMON,90)ITER,NY,SUMF,SUMX,FC 91 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,91) 92 FORMAT('1') WRITE(LUMON,92) IF(IRANK.LT.NE)THEN INFO = -1 ELSE C Print solution info 93 FORMAT('Solution ','of',1X,'boundary ','value ', * 'problem',' obtained',/,'BVPSOL',' required',I3,1X, * 'Iteration ','steps ','with',I4,1X,'trajectory', * ' evaluations',//) write ( lumon, '(a)' ) ' ' WRITE(LUMON,93)ITER,KOUNT CALL BLPRCV(LUMON,CONV,EPH) ENDIF ENDIF C ---------------------------------------------------------- C 5.2 Fail exit messages C ---------------------------------------------------------- C 5.2.1 Rank-deficiency : best least squares solution of bvp C obtained IF(INFO.EQ.-1.AND.KPRINT.GE.0)THEN 94 FORMAT('Iteration ','terminates ','at ', * 'stationary ','point',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,94) CALL BLPRCV(LUMON,CONVA,EPH) IF(ITER.NE.0)THEN SKAP = ZERO IF(FCA.EQ.ONE.AND.FC.EQ.ONE.AND.IRANKA.EQ.IRANK) SKAP * = DSQRT(SUMXA/FCNUMK) IF(SKAP.GT.ZERO)THEN 95 FORMAT('Incompatibility ','factor ','kappa',D10.3 * ,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,95)SKAP ENDIF ENDIF ENDIF C ---------------------------------------------------------- C 5.2.2 Termination after more than itmax iterations IF(INFO.EQ.-2.AND.KPRINT.GE.0)THEN 96 FORMAT('Iteration ','terminates ','after ','itmax ', * '=',I3,2X,'iteration ','steps') write ( lumon, '(a)' ) ' ' WRITE(LUMON,96)ITMAX ENDIF C ---------------------------------------------------------- C 5.2.3 Integrator failed in trajectory computation IF(INFO.EQ.-3.AND.KPRINT.GE.0)THEN 97 FORMAT('Integrator ','failed ','in ','trajectory ', * 'computation ',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,97) J1 =-KFLAG 98 FORMAT('BVPSOL ','terminates',/,'Subinterval',I3,1X, * 'possibly ','insert ','new ','node',D20.11,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,98)J1,TFAIL ENDIF C ---------------------------------------------------------- C 5.2.4 Convergence fail of Gauss - Newton method IF(INFO.EQ.-4.AND.KPRINT.GE.0)THEN 99 FORMAT('Gauss ','Newton ','method ','fails ','to ', * 'converge',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,99) ENDIF C ---------------------------------------------------------- C 5.2.5 Inconsistent initial data IF(INFO.EQ.-5.AND.KPRINT.GE.0)THEN 100 FORMAT('Error: ','initial ','data ','and ', * 'boundary ','conditions ','are ','inconsistent',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,100) ENDIF C ---------------------------------------------------------- C 5.2.6 Convergence fail of iterative refinement sweeps IF(INFO.EQ.-6)THEN IF(KPRINT.GE.0)THEN 101 FORMAT('Termination ','since ','iterative ', * 'refinement ','fails ','to ','converge',/,2X,'Insert ', * 'new ','nodes',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,101) ENDIF JN = JN-1 IF(JN.GT.0)THEN 102 FORMAT(8X,'in ','subinterval',2X,I3,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,102)JN ENDIF ENDIF C ---------------------------------------------------------- C 5.2.7 Insufficient error tolerance for integrator IF(INFO.EQ.-7.AND.KPRINT.GE.0)THEN TOLH = EPS/SIGDEL RELDIF = DSQRT(TOLH/SIGDEL) 103 FORMAT('Suggested ','integrator ','accuracy',D10.1 * ,2X,/,'Suggested ','relative ','deviation ', * 'parameter',D10.1,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,103)TOLH,RELDIF 104 FORMAT('Reduce ','relative ','error ','tolerance ', * 'for ','integrator ','to',D10.1,2X,/,2X,'or ','insert ', * 'new ','nodes',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,104)TOLH S = REDH/TOL DO 105 J=1,M1 IF(RF(J).GT.S)THEN 106 FORMAT(2X,'in ','subinterval',I3,/) WRITE(LUMON,106)J ENDIF 105 CONTINUE 107 FORMAT('Reliable ','relative ','accuracy ', * 'greater ','than',1X,D6.1,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,107)1.0D-2 ENDIF C ---------------------------------------------------------- C 5.2.8 ill - conditioned condensed linear system IF(INFO.EQ.-8.AND.KPRINT.GE.0)THEN 108 FORMAT('Gaussian ','block ','elimination ','fails',/, * 2X,'by ','ill ','- ','conditioned ','condensed ', * 'linear ','system',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,108) IF(IRANK.EQ.NE)THEN 109 FORMAT('Relative ','accuracy ','of ','DX1',D10.3 * ,2X) write ( lumon, '(a)' ) ' ' WRITE(LUMON,109)EPX1H ENDIF 110 FORMAT('Possibly set IOPT(3)=1',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,110) ENDIF C ---------------------------------------------------------- C 5.2.9 rank reduced to zero IF(INFO.EQ.-11.AND.KPRINT.GE.0)THEN 11001 FORMAT('rank ','reduction ','failed - ', * 2X,'resulting ','rank ','is ','zero ',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,11001) ENDIF C ---------------------------------------------------------- C 5.3 Common exit IF(KPRINT.GE.0)THEN C J1 = 1 SMALIN = ONE/SMALL IF(IRANKB.NE.0) $ CALL BLPRCD(LUMON,COND1,SENS1,SMALIN,J1,IRANKB) IF(IRANKB.NE.IRANK)THEN J1 = IRANKB+1 CALL BLPRCD(LUMON,COND2,SENS2,SMALIN,J1,IRANK) ENDIF 111 FORMAT('Multiple ','shooting ','condition',D10.3,2X, * /,'1') write ( lumon, '(a)' ) ' ' WRITE(LUMON,111)SIGDLH IF(INFO.GT.0)THEN 112 FORMAT('Solution ','data',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,112) ENDIF IF(INFO.LT.0)THEN 113 FORMAT('Final ','data',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,113) ENDIF DO 114 J=1,M 115 FORMAT(D13.5,2X) WRITE(LUMON,115)T(J) 116 FORMAT((14X,3(D20.10,1X))) WRITE(LUMON,116)(X(L1),L1=(J-1)*N+1,J*N) 114 CONTINUE ENDIF C End of exits 9999 CONTINUE RETURN END SUBROUTINE BVPG(FCN,BC,IVPSOL,N,M,M1,NM,NM1,NMX8,NZ,LICN, *LIRN,LISNQ,NKEEP,T,X,EPS,TOL,RELDIF,NONLIN,ITMAX,INFO,XTHR, *IROW,ICOLA,ICOLB,IA,IB,IW,G,A,B,E,WO,DX,DXQ,DXQA,XA,XW,XU,HH, *DHH,HHA,DE,R,DR,RA,U,DU,X1,XM,T1,T2,RF,IVECT,JVECT,ICN,IRN, *IKEEP,LUMON) c*********************************************************************72 c cc BVPG() carries out the global linear solver approach. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) EXTERNAL FCN,BC,IVPSOL INTEGER N,M,M1,NM,NM1,NMX8,LISNQ,NKEEP,LUMON DOUBLE PRECISION T(M),X(NM) DOUBLE PRECISION EPS DOUBLE PRECISION TOL,RELDIF INTEGER NONLIN,ITMAX INTEGER INFO DOUBLE PRECISION XTHR INTEGER IRN(LIRN),ICN(LICN),IVECT(NZ),JVECT(NZ),IKEEP( *NKEEP) INTEGER IROW(N),ICOLA(N),ICOLB(N) INTEGER IA(N,N),IB(N,N) INTEGER IW(NMX8) DOUBLE PRECISION G(N,N,M1) DOUBLE PRECISION A(N,N),B(N,N) DOUBLE PRECISION E(LICN),WO(NM) DOUBLE PRECISION DX(NM),DXQ(NM),DXQA(NM),XA(NM),XW(NM),XU( *NM1),HH(NM1),DHH(NM1),HHA(NM1),DE(N),R(N),DR(N),RA(N),U(NM), *DU(NM),X1(N),XM(N),T1(N),T2(N),RF(M1) C C Addtional dimensional integer variables: C ======================================== C C M1 M-1 C NM N*M C NM1 N*(M-1) C NMX8 8*N*M C LIRN,LICN,NKEEP,NZ C See driver routine BVPSOL C C Internal real arrays (workspace) : C ================================== C C G(N,N,M1) (N,N) -Wronskian Matrices G(1),...,G(M-1) C . C A(N,N) Wronskian Matrix on left boundary C dBC/dX(X(1,...,N),T(1)). C B(N,N) Wronskian Matrix on right boundary C dBC/dX(X((N-1)*M+1,...,N*M),T(M)). C E(LICN) Holds the values of the Jacobian stored C in sparse mode. C DE(N) Holds row scaling factors for the C boundary conditions part of the Jacobian C matrix. C DHH(NM1) Holds the continuity residuals computed C in BGSOLI . C DR(N) Workspace for subroutine BGSOLI to hold C the boundary residual C BC(DXQ(1,...,N),DXQ((M-1)*N+1,...,M*N))+ C (A*DXQ(1,...,N))+B*DXQ((M-1)*N+1,...,M*N) C . C DU(NM) Used by BGSOLI . Gets the total residual C for the current iterate. C DX(NM) Actual newton correction. C DXQ(NM) Simplified Newton correction J(k-1)*X(k) C with the Jacobian J(k) and the iterate C vector X(k) at the k-th iterate. C DXQA(NM) Previous simplified Newton correction C J(k-2)*X(k-1). C HH(NM1) Elements (J-1)*N+1 to J*N are holding C the values C Y(T(J+1),X((J-1)*N+1,...,J*N))-X(J*N+1, C ...,(J+1)*N) C ( with the trajectory Y in C [T(J),T(J+1)] , J = 1,...,M-1 ). C HHA(NM1) Holds the previous value of HH . C R(N) Value of the boundary condition function C BC for the current iterate. C RA(N) Previous values of R . C RF(M1) Used by BGSOLI . Gets the norms of the C Wronskian matrices. C T1(N) Workspace used for miscellaneous C purposes temporarely. C T2(N) Workspace used for miscellaneous C purposes temporarely. C U(NM) Gets the right hand side of the linear C system to be solved in each iteration C step. Used in BGSOLI . C WO(NM) Workspace needed for sparse solver. Must C not be altered outside the sparse packet C routines. C XA(NM) Previous Newton iterate. C XU(NM1) Elements (J-1)*N+1 to J*N are holding C the values Y(T(J+1),X((J-1)*N+1,...,J*N)) C of the trajectory in the interval C [T(J),T(J+1)] , (for J = 1,...,M-1 ). C XW(NM) Scaling factors for iteration vector. C X1(N) Components of the iteration vector C corresponding to the left boundary C A = T(1). C XM(N) Components of the iteration vector C corresponding to the right boundary C B = T(M). C C Internal integer arrays (workspace) C =================================== C C IROW(N) Row permutations of boundary derivative C matrices A and B . C ICOLA(N) Column permutations of matrix A C (left boundary). C ICOLB(N) Column permutations of matrix B C (right boundary). C IA(N,N) Reflects the sparse structure of matrix C A by values 0, 1. C IB(N,N) Reflects the sparse structure of matrix C B by values 0, 1. C IW(NMX8) Workspace needed for sparse solver C package. C C Internal short integer arrays (workspace) C ========================================= C C IRN(LIRN) Workspace for MA28/MA30 sparse package. C On Input to routine MA28A, it must hold C the row indices of the sparse matrix. C ICN(LICN) Workspace for MA28/MA30 sparse package. C On Input to routine MA28A, it must hold C the column indices of the sparse matrix. C IVECT(NZ) Input to routine MA28B: must hold the C row indices of the sparse matrix. C JVECT(NZ) Input to routine MA28B: must hold the C column indices of the sparse matrix. C IKEEP(NKEEP) Workspace array for MA28 sparse package. C To be preserved across the calls of the C routines MA28A,MA28B,MA28C . C C Internal real variables: C ======================== C C COND Gets the condition of the Jacobian C matrix computed by BGSOLI . C CORR Gets the 1-norm of the residual DU . C Computed by BGSOLI . C CONV Scaled maximum norm of DXQ computed by C subroutine BGLVLS . Used for convergence C test. C CONVA Holds the previous value of CONV . C EPSMIN Smallest reasonable permitted accuracy C EPS that can be prescribed by the user. C FC Actual Gauss Newton iteration damping C factor. C FCA Previous Gauss Newton iteration damping C factor. C FCDNM Used to compute the denominator of the C damping factor FC during computation of C it's predictor, corrector and C aposteriori estimate (in the case of C performing a Rank1 update) . C FCH Temporarely used for storing the new FC C when computing aposteriori estimate. C FCMIN Minimum permitted relaxation factor. If C FC becomes smaller than this value, one C of the following may occur: C a. Recomputation of the sensitivity C matrix by means of difference C approximation (instead of Rank1 C update), if Rank1 - update C previously was used C b. Rank reduction of sensitivity C matrix E , if difference C approximation was used previously C and Rank(E).NE.0 C c. Fail exit otherwise C FCMIN2 FCMIN**2 . Used for FC-predictor C computation. C FCNUM Gets the numerator of the aposteriori C estimate of FC . C FCNUMP Gets the numerator of the predictor C computation of FC . C FCNUMK Gets the numerator of the corrector C computation of FC . C H Actual integrator stepsize. C HMAX Maximum permitted integrator stepsize. C Set to the length of the integration C interval, e.g. the distance of the C effected Shooting points. C HSAVE Stepsize saved across the call of the C integrator. C HSTART Start stepsize for integration used by C subroutines BLFCNI and BLDERG . C MUE Temporary value used during computation C of damping factors predictor. C REDH Multi purpose reduction factor. (???) C RELDIF Relative deviation for numerical C differentation. C SIGMA Decision parameter for Jacobian Rank1 C updates (SIGMA.GT.1) . Rank1 updates are C inhibited, if SIGMA.GT.1/FCMIN is set. C SKAP Used to compute and print out the C incompatibility factor of the nonlinear C boundary value (e.g. least squares) C problem. C SUMF Standard level of the current iterate, C e.g. Norm2(F(X))**2 C with the nonlinear model function F on C which Newton iteration is performed, C arising from the Multiple Shooting C approach. C SUMX Natural level of the current iterate, C e.g. Norm2(DX) C with the Newton correcture DX C (see above). C SUMXA Natural level of the previous iterate. C TFAIL Used to get and print out in case of an C integrator failure the last reached T C value as a proposal for insertion of a C new Shooting point. C TOL Prescribed relative precision for C numerical integration. C TOLH Temporary used for computation of TOL C (may be obmitted|). C TOLMIN Lower bound value for TOL . C XTHR Threshold for scaling. C TJ Used by BLFCNI to hold T(J). C TJ1 Used by BLFCNI to hold T(J+1). C EPH Gets TOL*SIGDEL by BGSOLI . If EPH.GT. C REDH , termination occurs, since C Multiple Shooting condition is too bad. C SIGDEL Used by BGSOLI to compute the required C integrator accuracy from the multiple C shooting condition. C SIGDLH Used by BGSOLI temporary during C determination of SIGDEL . C C Internal integer variables C ========================== C C IAF Indicates, if sparse structure of C Jacobian matrix must be reordered by C sparse solver or not: C 0 Not necessary, can call MA28B to C decompose Jacobian C 1 Must be done, so use MA28A . C IC Permutated index. Used by BGSOLI . C ICA Temporarely used during search for C separable linear boundary conditions at C T(1). C ICB Temporarely used during search for C separable linear boundary conditions at C T(M). C IH Temporarely used during search for C separable linear boundary conditions at C T(1)and T(M). C INZ Count of nonzero Jacobian matrix C elements. C IR Temporary storage for a row permutation C index. C IRANK Rank of Jacobian matrix E of current C iteration step estimated by sparse C solver package (Common block variable). C IS Additional DO loop index. C ISUM Used for determination of sparse C structure of matrices A and B as C nonzeros location counter. C ITER Iteration count. C JJ Used as "reverse DO loop" index: C JJ = IUPB-J in a loop like DO J = 1,IUPB C ... C JRED Damping factor reduction count during an C iterate. C KC Temporary storage for a column C permutation index. C KFLAG Gets the subintervall number of the C failure from subroutine BLDERG , if the C integrator failed. C KOUNT Trajectory evaluations count. C KPRINT Print parameter - copy of input C parameter INFO . C LEVEL Flow control parameter needed by C subroutine BGSOLI : C 0 indicates computation of Newton C correcture, C 1 indicates computation of C simplified Newton correcture C (after computation of the C preliminary new iterate) C NA Number of separable boundary conditions C at T(1): N-NAQ C NAQ Number of not separable boundary C conditions at T(1) C NB Number of separable boundary conditions C at T(M) C NBQ Number of not separable boundary C conditions at T(M): N-NB C NEW Count of subsequent performed Rank1 C (Broyden) updates. C NM Number of rows and columns of the C Jacobian matrix part to be decomposed by C the sparse solver. C NRS N-(NA+NB) C: End Parameter C: EPMACH = relative machine precision DOUBLE PRECISION EPMACH,SMALL DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) DOUBLE PRECISION HALF PARAMETER (HALF=0.5D0) DOUBLE PRECISION REDH PARAMETER (REDH=1.0D-2) DOUBLE PRECISION ONE PARAMETER (ONE=1.0D0) DOUBLE PRECISION TWO PARAMETER (TWO=2.0D0) DOUBLE PRECISION EIGHT PARAMETER (EIGHT=8.0D0) DOUBLE PRECISION TEN PARAMETER (TEN=1.0D1) DOUBLE PRECISION FCMIN PARAMETER (FCMIN=1.0D-2) DOUBLE PRECISION FCMIN2 PARAMETER (FCMIN2=1.0D-4) DOUBLE PRECISION FCNLIN PARAMETER (FCNLIN=1.0D-2) DOUBLE PRECISION SIGMA PARAMETER (SIGMA=2.0D0) INTEGER I,ICA,ICB,ICNCP,IH,IR,IRNCP,IRANK,IS,ISUM,ITER, *J,JRED,J0,J1,K,KFLAG,KOUNT,KPRINT, *L,LEVEL,MINIRN,MINICN,NB,NDIM,NAQ,NEW DOUBLE PRECISION COND,CORR,CONV,CONVA,EPH,EPSMIN, *EPX1H,EPSQ,FC,FCA,FCDNM,FCH,FCNMP2,FCNUM,FCNUMK,FCNUMP, *HSTART,MUE,RESID,RMIN,S,SIGDEL,SIGDLH,SUMF, *SUMX,SUMXA,TFAIL,TH,TOLH,TOLMIN,UQ LOGICAL ABORT1,ABORT2,DIFAPP,FCOMPT,GROW,IBLOCK,IVFAIL, *JACRFR,JACRST,NEXT COMMON /MA28ED/ LP,MP,IBLOCK,GROW COMMON /MA28FD/ EPSQ,RMIN,RESID,IRNCP,ICNCP,MINIRN,MINICN, *IRANK,ABORT1,ABORT2 INTEGER L1,L2 DOUBLE PRECISION S1 C: Begin C: Begin of Segment BVPSOG.Body C ---------------------------------------------------------- C 1 Initialization C ---------------------------------------------------------- C 1.1 Internal parameters C Standard values fixed below C Minimum relative precision of integrator ( to be adapted ) CALL ZIBCONST(EPMACH,SMALL) TOLMIN = EPMACH*TEN*TEN C Maximum permitted number of iterative refinements sweeps C ---------------------------------------------------------- C 1.1.1 Common parameters C Starting value of relaxation factor (FCMIN.LE.FC.LE.1.0) IF(NONLIN.LE.1)THEN C for linear or mildly nonlinear problems FC = ONE ELSE C for highly nonlinear problems FC = FCNLIN ENDIF C Minimum reasonable value for EPS EPSMIN = DSQRT(TEN*EPMACH) IF(EPS.LT.EPSMIN) EPS = EPSMIN C ---------------------------------------------------------- C 1.2 Initial preparations IF(FC.LT.FCMIN) FC = FCMIN IF(FC.GT.ONE) FC = ONE KPRINT = INFO ITER = 0 KOUNT = 0 INFO = -1000 FCA = FC CONV = ZERO JACRFR = .FALSE. JACRST = .FALSE. C: Begin SetVec.Vec DO 6 L1=1,NM XA(L1)=X(L1) 6 CONTINUE C. End SetVec.Vec IF(TOL.LE.ZERO) TOL = EPS/TEN IF(TOL.LT.TOLMIN) TOL = TOLMIN DIFAPP = .TRUE. HSTART =(T(2)-T(1))*REDH LP = 0 MP = 0 IBLOCK = .FALSE. EPSQ = 0.1D0 UQ = 0.1D0 SIGDLH = ZERO SUMF = ONE C: Mat IA = Scalar (Rows 1,N ; Cols 1,N) L1 = 0 DO 7 L2=1,N DO 7 L43=1,N IA(L2,L43)=L1 7 CONTINUE C. End SetIntMat.S C: Mat IB = Scalar (Rows 1,N ; Cols 1,N) L1 = 0 DO 8 L2=1,N DO 8 L43=1,N IB(L2,L43)=L1 8 CONTINUE C. End SetIntMat.S C: CubeMat G (layer 1)= Scalar (Rows 1,N ; Cols 1,N) S1 = ZERO DO 9 L1=1,N DO 9 L2=1,N G(L1,L2,1)=S1 9 CONTINUE C. End SetCubeMat.S IF(KPRINT.GE.0)THEN C Print Start vector data, predescribed precision and max C iteration steps 10 FORMAT('Initial ','data',//) write ( lumon, '(a)' ) ' ' WRITE(LUMON,10) DO 11 J=1,M 12 FORMAT(D13.5,2X) WRITE(LUMON,12)T(J) 13 FORMAT((14X,3(D20.10,1X))) WRITE(LUMON,13)(X(L1),L1=(J-1)*N+1,J*N) 11 CONTINUE 14 FORMAT('N ','=',I2,2X,'M ','=',I2,/, * 'Prescribed ','relative ','precision',D10.2,2X,/, * 'Maximum ','permitted ','number ','of ','iteration ', * 'steps',1X,I3,//,'1') write ( lumon, '(a)' ) ' ' WRITE(LUMON,14)N,M,EPS,ITMAX IF(KPRINT.EQ.0)THEN 15 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,15) 16 FORMAT(4X,'It',7X,'Levelf',10X,'Levelx',8X, * 'Rel.Fc.') write ( lumon, '(a)' ) ' ' WRITE(LUMON,16) ENDIF ENDIF C ---------------------------------------------------------- C 1.3 Startup step C ---------------------------------------------------------- C 1.3.1 Computation of the residual vector CALL BLFCNI(IVPSOL,FCN,BC,N,M,NM,NM1,ITER,KPRINT,HSTART, * FCMIN,T,X,X1,XM,T1,XU,HH,R,TOL,FC,FCOMPT,IVFAIL,KFLAG, * KOUNT,INFO,LUMON) C C Main iteration loop C =================== C C: While (expression) 17 IF(INFO.EQ.-1000)THEN C: Begin of Segment BVPSOG.Core C ---------------------------------------------------- C 2 Startup of iteration step IF(.NOT.(JACRFR.OR.JACRST))THEN LEVEL = 0 C -------------------------------------------------- C 2.1 Scaling of variables X(NM) CALL BLSCLE(N,M,NM,NM1,X,XU,XW,XTHR) IF(ITER.NE.0)THEN C: Begin SetVec.Vec DO 18 L1=1,NM DXQA(L1)=DXQ(L1) 18 CONTINUE C. End SetVec.Vec C: FCNUM = Sum of Formula Elements (for 1,NM) FCNUM = 0.0D0 DO 19 L1=1,NM FCNUM=FCNUM+((DX(L1)/XW(L1))**2) 19 CONTINUE C. End MakeSum.Comp C: FCNMP2 = Sum of Formula Elements (for 1,NM) FCNMP2 = 0.0D0 DO 20 L1=1,NM FCNMP2=FCNMP2+((DXQ(L1)/XW(L1))**2) 20 CONTINUE C. End MakeSum.Comp FCNUMP = FCNUM*FCNMP2 ENDIF IF(ITER.NE.0)THEN TH = FC-ONE C: FCDNM = Sum of Formula Elements (for 1,NM) FCDNM = 0.0D0 DO 21 L1=1,NM FCDNM=FCDNM+(((DXQ(L1)+TH*DX(L1))/XW(L1))**2) 21 CONTINUE C. End MakeSum.Comp FCH = DSQRT(FCNUM/FCDNM)*FC*FC*HALF C ------------------------------------------------ C 2.1.1 Decision criterion for Jacobian updating C technique: C DIFAPP.EQ..TRUE. numerical C differentiation, C DIFAPP.EQ..FALSE. rank1 updating DIFAPP = FC.LT.FCA.AND.NEW.GT.0.OR.FCH.LT.FC* * SIGMA FCA = FC IF(NONLIN.GT.0) FC = DMIN1(FCH,ONE) ENDIF C -------------------------------------------------- C 2.2 Difference approximation of jacobian matrix A C ( If Difapp.EQ..TRUE. ) or C Rank-1 update of jacobian matrix A ( If Difapp C .EQ..FALSE. ) CALL BLDERA(BC,N,M,NM,XW,X1,XM,R,T2,A,B,RELDIF) C -------------------------------------------------- C 2.3 Determination of sparse structure of matrices C A and B IAF = 0 DO 22 I=1,N S = ZERO DO 23 K=1,N TH = DABS(A(I,K))*XW(K) S = S+TH TH = DABS(B(I,K))*XW(K+NM1) S = S+TH 23 CONTINUE IF(S.LT.XTHR) S = XTHR DE(I)=ONE/S DO 24 K=1,N IF(IA(I,K).LE.0)THEN IF(A(I,K).NE.ZERO)THEN IA(I,K)=1 IAF = 1 ENDIF ENDIF IF(IB(I,K).LE.0)THEN IF(B(I,K).NE.ZERO)THEN IB(I,K)=1 IAF = 1 ENDIF ENDIF 24 CONTINUE 22 CONTINUE IF(IAF.NE.0)THEN C ------------------------------------------------ C 2.3.1 Determination of row and column C permutation vectors DO 25 I=1,N ICOLA(I)=I ICOLB(I)=I IROW(I)=I 25 CONTINUE C ------------------------------------------------ C 2.3.2 Search for separable linear boundary C conditions at T(1) NAQ = N DO 26 I=1,N DO 27 K=1,N IF(IB(I,K).NE.0) GOTO 9996 27 CONTINUE ISUM = 0 DO 28 K=1,N IF(IA(I,K).NE.0)THEN ISUM = ISUM+1 ICA = K ENDIF 28 CONTINUE IF(ISUM.LE.1)THEN DO 29 IS=1,N IH = ICOLA(IS) IF(IH.EQ.ICA) ICOLA(IS)=ICOLA(NAQ) IH = IROW(IS) IF(IH.EQ.I) IROW(IS)=IROW(NAQ) 29 CONTINUE ICOLA(NAQ)=ICA IROW(NAQ)=I NAQ = NAQ-1 IF(DABS(R(I)).GT.TEN*EPMACH*DABS(X(ICA))) * THEN INFO = -5 GOTO 9998 ENDIF ENDIF 9996 CONTINUE 26 CONTINUE IF(KPRINT.GE.0.AND.NAQ.EQ.0)THEN 30 FORMAT('Warning: ','attempt ','to ', * 'solve ','initial ','value ','problem') write ( lumon, '(a)' ) ' ' WRITE(LUMON,30) ENDIF NA = N-NAQ C ------------------------------------------------ C 2.3.3 Search for separable linear boundary C conditions at T(M) NB = 0 ENDIF IF(IAF.NE.0.AND.NAQ.NE.0)THEN DO 31 I=1,NAQ IR = IROW(I) DO 32 K=1,N IF(IA(IR,K).NE.0) GOTO 9995 32 CONTINUE ISUM = 0 DO 33 K=1,N IF(IB(IR,K).NE.0)THEN ISUM = ISUM+1 ICB = K ENDIF 33 CONTINUE IF(ISUM.LE.1)THEN NB = NB+1 DO 34 IS=1,N IH = ICOLB(IS) IF(IH.EQ.ICB) ICOLB(IS)=ICOLB(NB) 34 CONTINUE ICOLB(NB)=ICB IROW(I)=IROW(NB) IROW(NB)=IR IF(DABS(R(IR)).GT.TEN*EPMACH*DABS(X(ICB+ * NM1)))THEN INFO = -5 GOTO 9998 ENDIF ENDIF 9995 CONTINUE 31 CONTINUE IF(KPRINT.GE.0.AND.NB.EQ.N)THEN 35 FORMAT('Warning: ','attempt ','to ', * 'solve ','initial ','value ','problem') write ( lumon, '(a)' ) ' ' WRITE(LUMON,35) ENDIF NBQ = N-NB C ------------------------------------------------ C 2.3.4 Count non-zeroes in jacobian and store C their locations NDIM = NM-NA-NB NRS = N-NA-NB INZ = 0 DO 36 I=1,N IF(NAQ.NE.0)THEN DO 37 K=1,NAQ INZ = INZ+1 IVECT(INZ)=I JVECT(INZ)=K 37 CONTINUE ENDIF II = I IF(M.EQ.2)THEN II = 0 IF(NBQ.EQ.0) GOTO 9994 DO 38 L=1,NBQ IF(ICOLB(NB+L).EQ.I) II = L 38 CONTINUE IF(II.EQ.0) GOTO 9994 ENDIF INZ = INZ+1 IVECT(INZ)=I JVECT(INZ)=NAQ+II 9994 CONTINUE 36 CONTINUE IF(M.NE.2)THEN DO 39 J=2,M1 DO 40 I=1,N DO 41 K=1,N INZ = INZ+1 IVECT(INZ)=(J-1)*N+I JVECT(INZ)=(J-1)*N+K-NA 41 CONTINUE II = I IF(J.EQ.M1)THEN II = 0 IF(NBQ.EQ.0) GOTO 9993 DO 42 L=1,NBQ IF(ICOLB(NB+L).EQ.I) II = L 42 CONTINUE IF(II.EQ.0) GOTO 9993 ENDIF INZ = INZ+1 IVECT(INZ)=(J-1)*N+I JVECT(INZ)=J*N+II-NA 9993 CONTINUE 40 CONTINUE 39 CONTINUE ENDIF IF(NRS.NE.0)THEN DO 43 I=1,NRS II = I+NB IF(NAQ.NE.0)THEN DO 44 K=1,NAQ IF(IA(IROW(II),ICOLA(K)).NE.0)THEN INZ = INZ+1 IVECT(INZ)=NM1+I JVECT(INZ)=K ENDIF 44 CONTINUE ENDIF IF(NBQ.NE.0)THEN DO 45 KK=1,NBQ K = KK+NB IF(IB(IROW(II),ICOLB(K)).NE.0)THEN INZ = INZ+1 IVECT(INZ)=NM1+I JVECT(INZ)=NM1+KK-NA ENDIF 45 CONTINUE ENDIF 43 CONTINUE ENDIF ENDIF ENDIF JACRFR = .FALSE. IF(.NOT.JACRST)THEN IF(DIFAPP)THEN NEW = 0 KFLAG = 0 CALL BLDERG(FCN,N,NAQ,M,M1,NM,NM1,T,X,XU,XW,T2, * TFAIL,G,ICOLA,IVPSOL,HSTART,TOL,RELDIF,KFLAG) IF(KFLAG.LT.0)THEN INFO = -3 GOTO 9998 ENDIF IF(M.GT.2) KOUNT = KOUNT+N IF(M.EQ.2) KOUNT = KOUNT+NAQ ELSE NEW = NEW+1 CALL BLRK1G(N,M,M1,NM,NM1,XW,DX,HH,HHA,T1,G,FCA) ENDIF ENDIF JACRST = .FALSE. C ---------------------------------------------------- C 2.3.5 Storing of total sparse jacobian C ( including row and column scaling ) INZ = 0 DO 46 I=1,N IF(NAQ.NE.0)THEN DO 47 K=1,NAQ INZ = INZ+1 E(INZ)=-G(I,ICOLA(K),1)*XW(ICOLA(K))/XW(I+N) 47 CONTINUE ENDIF II = I IF(M.EQ.2)THEN II = 0 IF(NBQ.EQ.0) GOTO 9992 DO 48 L=1,NBQ IF(ICOLB(NB+L).EQ.I) II = L 48 CONTINUE IF(II.EQ.0) GOTO 9992 ENDIF INZ = INZ+1 E(INZ)=ONE 9992 CONTINUE 46 CONTINUE IF(M.NE.2)THEN DO 49 J=2,M1 J0 =(J-1)*N J1 = J0+N DO 50 I=1,N DO 51 K=1,N INZ = INZ+1 E(INZ)=-G(I,K,J)*XW(K+J0)/XW(I+J1) 51 CONTINUE II = I IF(J.EQ.M1)THEN II = 0 IF(NBQ.EQ.0) GOTO 9991 DO 52 L=1,NBQ IF(ICOLB(NB+L).EQ.I) II = L 52 CONTINUE IF(II.EQ.0) GOTO 9991 ENDIF INZ = INZ+1 E(INZ)=ONE 9991 CONTINUE 50 CONTINUE 49 CONTINUE ENDIF IF(NRS.NE.0)THEN DO 53 I=1,NRS II = I+NB IF(NAQ.NE.0)THEN DO 54 K=1,NAQ IF(IA(IROW(II),ICOLA(K)).NE.0)THEN INZ = INZ+1 E(INZ)=-A(IROW(II),ICOLA(K))*DE(IROW(II))* * XW(ICOLA(K)) ENDIF 54 CONTINUE ENDIF IF(NBQ.NE.0)THEN DO 55 KK=1,NBQ K = KK+NB IF(IB(IROW(II),ICOLB(K)).NE.0)THEN INZ = INZ+1 E(INZ)=-B(IROW(II),ICOLB(K))*DE(IROW(II))* * XW(ICOLB(K)+NM1) ENDIF 55 CONTINUE ENDIF 53 CONTINUE ENDIF C ---------------------------------------------------- C 2.4 Save values of R(N)and HH((M-1)*N) C: Begin SetVec.Vec DO 56 L1=1,N RA(L1)=R(L1) 56 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec DO 57 L1=1,NM1 HHA(L1)=HH(L1) 57 CONTINUE C. End SetVec.Vec NEXT = .FALSE. C ---------------------------------------------------- C 3 Main-part of iteration step C ---------------------------------------------------- C 3.1 Solution of the linear system C ---------------------------------------------------- C 3.1.1 Decomposition of (N,N)-matrix A C ---------------------------------------------------- C 3.1.2 LU-decomposition of(NDIM,NDIM)-MATRIX E IF(IAF.EQ.1)THEN DO 58 I=1,INZ IRN(I)=IVECT(I) ICN(I)=JVECT(I) 58 CONTINUE CALL MA28AD(NDIM,INZ,E,LICN,IRN,LIRN,ICN,UQ,IKEEP, * IW,WO,IFLAG) IAF = 0 ELSE CALL MA28BD(NDIM,INZ,E,LICN,IVECT,JVECT,ICN,IKEEP, * IW,WO,IFLAG) IF(RMIN.LT.1.0D-4) IAF = 1 IF(IAF.EQ.1)THEN JACRST = .TRUE. GOTO 9998 ENDIF ENDIF IF(IFLAG.EQ.-1.OR.IFLAG.EQ.-2)THEN INFO = -1 GOTO 9998 ENDIF IF(IFLAG.LE.-3)THEN INFO = -9 GOTO 9998 ENDIF C ---------------------------------------------------- C 3.1.3 Solution of linear (N,N)-system CALL BGSOLI(N,M,M1,NM,NM1,NDIM,LICN,NKEEP,NA,NAQ,NB, * NBQ,NRS,ITER,LEVEL,KPRINT,EPS,REDH,TOLMIN,FC,FCA,TOL, * RELDIF,EPH,EPX1H,SIGDEL,SIGDLH,COND,CORR,HH,DHH,R,A, * B,G,U,DE,DU,T1,DXQ,XW,DR,RF,WO,E,IROW,ICOLA,ICOLB, * ICN,IKEEP,INFO,LUMON) IF(INFO.NE.-1000) GOTO 9998 C ---------------------------------------------------- C 3.2 Evaluation of scaled natural level function SUMX C scaled maximum error norm CONV C evaluation of (scaled) standard level function C SUMF ( SUMF only, if KPRINT.GE.0 ) C and computation of ordinary newton corrections C DX(N) CALL BGLVLS(N,M,NM,NM1,XW,DXQ,HH,R,DE,CONV,SUMX,SUMF, * KPRINT) C: Begin SetVec.Vec DO 59 L1=1,NM DX(L1)=DXQ(L1) 59 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec DO 60 L1=1,NM XA(L1)=X(L1) 60 CONTINUE C. End SetVec.Vec SUMXA = SUMX CONVA = CONV C ---------------------------------------------------- C 3.3 a - priori estimate of relaxation factor FC JRED = 0 IF(ITER.NE.0.AND.NONLIN.NE.0)THEN IF(NEW.EQ.0)THEN C ------------------------------------------------ C 3.3.1 Computation of the denominator of a-priori C estimate C: FCDNM = Sum of Formula Elements (for 1,NM) FCDNM = 0.0D0 DO 61 L1=1,NM FCDNM=FCDNM+(((DX(L1)-DXQA(L1))/XW(L1))**2) 61 CONTINUE C. End MakeSum.Comp C ------------------------------------------------ C 3.3.2 New relaxation factor FCDNM = FCDNM*SUMX IF(FCDNM.GE.FCNUMP*FCMIN2)THEN MUE = FCA*DSQRT(FCNUMP/FCDNM) FC = DMIN1(MUE,ONE) ELSE FC = ONE ENDIF ENDIF IF(FC.LT.FCMIN)THEN INFO = -4 GOTO 9997 ENDIF ENDIF LEVEL = 1 C ---------------------------------------------------- C 3.4 Save natural level for later computations of C corrector and print iterate FCNUMK = SUMX IF(KPRINT.GE.0)THEN C Print Standard - and natural level IF(KPRINT.GT.0)THEN 62 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,62) 63 FORMAT(4X,'It',7X,'Levelf',10X,'Levelx',18X, * 'New') write ( lumon, '(a)' ) ' ' WRITE(LUMON,63) ENDIF 64 FORMAT(4X,I2,5X,D10.3,2X,4X,D10.3,2X,13X,I2) write ( lumon, '(a)' ) ' ' WRITE(LUMON,64)ITER,SUMF,SUMXA,NEW IF(KPRINT.GT.0)THEN 65 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,65) ENDIF ENDIF C C Relaxation-factor reduction loop C ================================ C C: DO (Until) 66 CONTINUE C -------------------------------------------------- C 3.5 Preliminary new iterate C: DO (Until) 67 CONTINUE FCOMPT = .FALSE. C: Vec X = Vec XA + Vec DX * Scalar (for 1,NM) S1 = FC DO 68 L1=1,NM X(L1)=XA(L1)+DX(L1)*S1 68 CONTINUE C. End SetVec.Vec&VecxS IF(ITER.GT.ITMAX)THEN INFO = -2 GOTO 9997 ENDIF C ------------------------------------------------ C 3.5.1 Computation of the residual vector CALL BLFCNI(IVPSOL,FCN,BC,N,M,NM,NM1,ITER,KPRINT, * HSTART,FCMIN,T,X,X1,XM,T1,XU,HH,R,TOL,FC,FCOMPT, * IVFAIL,KFLAG,KOUNT,INFO,LUMON) IF(IVFAIL)THEN INFO = -4 GOTO 9997 ENDIF IF(INFO.NE.-1000) GOTO 9997 IF(.NOT.(FCOMPT)) GOTO 67 C. UNTIL ( expression - negated above) C -------------------------------------------------- C 3.5.2 Solution of linear (N,N)-system CALL BGSOLI(N,M,M1,NM,NM1,NDIM,LICN,NKEEP,NA,NAQ, * NB,NBQ,NRS,ITER,LEVEL,KPRINT,EPS,REDH,TOLMIN,FC, * FCA,TOL,RELDIF,EPH,EPX1H,SIGDEL,SIGDLH,COND,CORR, * HH,DHH,R,A,B,G,U,DE,DU,T1,DXQ,XW,DR,RF,WO,E,IROW, * ICOLA,ICOLB,ICN,IKEEP,INFO,LUMON) IF(INFO.NE.-1000) GOTO 9998 C -------------------------------------------------- C 3.5.3 Evaluation of scaled natural level function C SUMX C scaled maximum error norm CONV and C evaluation of (scaled) standard level C function SUMF CALL BGLVLS(N,M,NM,NM1,XW,DXQ,HH,R,DE,CONV,SUMX, * SUMF,KPRINT) C -------------------------------------------------- C 3.6 Convergence test IF(CONV.LE.EPS)THEN INFO = 0 GOTO 9997 ENDIF C -------------------------------------------------- C 3.7 Natural monotonicity test IF(SUMX.GT.SUMXA)THEN C ------------------------------------------------ C 3.8 Output of iterate IF(KPRINT.GE.0)THEN C Print Standard - and natural level, and C damping factor IF(KPRINT.GT.0)THEN 69 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,69) 70 FORMAT(4X,'It',7X,'Levelf',10X,'Levelx', * 8X,'Rel.Fc.') write ( lumon, '(a)' ) ' ' WRITE(LUMON,70) ENDIF 71 FORMAT(4X,I2,5X,D10.3,2X,4X,D10.3,2X,4X,F5.3) write ( lumon, '(a)' ) ' ' WRITE(LUMON,71)ITER,SUMF,SUMX,FC IF(KPRINT.GT.0)THEN 72 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,72) ENDIF ENDIF JRED = JRED+1 IF(NONLIN.EQ.0)THEN INFO = -4 GOTO 9997 ENDIF C ------------------------------------------------ C 3.9 Compute reduced relaxation factor TH = FC-ONE C: FCDNM = Sum of Formula Elements (for 1,NM) FCDNM = 0.0D0 DO 73 L1=1,NM FCDNM=FCDNM+(((DXQ(L1)+TH*DX(L1))/XW(L1))**2) 73 CONTINUE C. End MakeSum.Comp FC = DSQRT(FCNUMK/FCDNM)*FC*FC*HALF C ------------------------------------------------ C 3.10 Fail exit, if relaxation factor to small JACRFR = FC.LT.FCMIN.OR.NEW.GT.0.AND.JRED.GT.1 IF(JACRFR.AND.NEW.EQ.0)THEN INFO = -4 GOTO 9998 ENDIF ENDIF IF(.NOT.(SUMX.LE.SUMXA.OR.JACRFR)) GOTO 66 C. UNTIL ( expression - negated above) C C End of relaxation-factor reduction loop C ======================================= C IF(JACRFR)THEN C -------------------------------------------------- C 3.11 Restore former values for repeting iteration C step C Restore former values LEVEL = 0 C: Begin SetVec.Vec DO 74 L1=1,N R(L1)=RA(L1) 74 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec DO 75 L1=1,N X1(L1)=XA(L1) 75 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec DO 76 L1=1,N XM(L1)=XA(L1+NM1) 76 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec DO 77 L1=1,NM X(L1)=XA(L1) 77 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec&Vec DO 78 L1=1,NM1 XU(L1)=X(L1+N)+HHA(L1) 78 CONTINUE C. End SetVec.Vec&Vec C: Begin SetVec.Vec DO 79 L1=1,NM1 HH(L1)=HHA(L1) 79 CONTINUE C. End SetVec.Vec IF(KPRINT.GE.0)THEN 80 FORMAT(5X,I2,1X,'Not ','accepted ', * 'relaxation ','factor',5X,F5.3) write ( lumon, '(a)' ) ' ' WRITE(LUMON,80)ITER,FC ENDIF IF(ITER.EQ.0)THEN FC = FCMIN ENDIF DIFAPP = .TRUE. JACRFR = .TRUE. GOTO 9998 ENDIF C ---------------------------------------------------- C 4 Preparations to start the following iteration step ITER = ITER+1 FCA = FC C ---------------------------------------------------- C 4.1 Print values IF(KPRINT.GE.0)THEN C Print Standard - and natural level, and damping C factor IF(KPRINT.GT.0)THEN 81 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,81) 82 FORMAT(4X,'It',7X,'Levelf',10X,'Levelx',8X, * 'Rel.Fc.') write ( lumon, '(a)' ) ' ' WRITE(LUMON,82) ENDIF 83 FORMAT(4X,I2,5X,D10.3,2X,4X,D10.3,2X,4X,F5.3) write ( lumon, '(a)' ) ' ' WRITE(LUMON,83)ITER,SUMF,SUMX,FC IF(KPRINT.GT.0)THEN 84 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,84) DO 85 J=1,M 86 FORMAT(D13.5,2X) WRITE(LUMON,86)T(J) 87 FORMAT((14X,3(D20.10,1X))) WRITE(LUMON,87)(X(L1),L1=(J-1)*N+1,J*N) 85 CONTINUE ENDIF ENDIF 9997 CONTINUE C. End of Segment BVPSOG.Core 9998 CONTINUE GOTO 17 ENDIF C. EndWhile C C End of main iteration loop C ========================== C C ---------------------------------------------------------- C 5 Exits C ---------------------------------------------------------- C 5.1 Solution exit IF(INFO.EQ.0)THEN ITER = ITER+1 C: Vec X = Vec X + Vec DXQ (for 1,NM) DO 88 L1=1,NM X(L1)=X(L1)+DXQ(L1) 88 CONTINUE C. End SetVec.&Vec INFO = ITER IF(KPRINT.LT.0)THEN GOTO 9999 ENDIF IF(KPRINT.GT.0)THEN C Print levels, damping factor of last iteration step C and solution info 89 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,89) 90 FORMAT(4X,'It',7X,'Levelf',10X,'Levelx',8X, * 'Rel.Fc.') write ( lumon, '(a)' ) ' ' WRITE(LUMON,90) ENDIF 91 FORMAT(4X,I2,5X,D10.3,2X,4X,D10.3,2X,4X,F5.3) write ( lumon, '(a)' ) ' ' WRITE(LUMON,91)ITER,SUMF,SUMX,FC 92 FORMAT(1X,66('*')) write ( lumon, '(a)' ) ' ' WRITE(LUMON,92) write ( lumon, '(a)' ) ' ' write ( lumon, '(a)' ) ' ' write ( lumon, '(a)' ) ' ' C Print solution info 94 FORMAT('Solution ','of',1X,'boundary ','value ', * 'problem',' obtained',/,'BVPSOL',' required',I3,1X, * 'Iteration ','steps ','with',I4,1X,'trajectory', * ' evaluations',//) write ( lumon, '(a)' ) ' ' WRITE(LUMON,94)ITER,KOUNT 95 FORMAT('Achieved ','relative ','accuracy',D10.3,2X) write ( lumon, '(a)' ) ' ' WRITE(LUMON,95)CONV IF(EPH.GT.CONV) CONV = EPH 96 FORMAT('Reliable ','relative ','accuracy',D10.3,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,96)CONV S = DLOG(DBLE((M-1)*(2*N+M-1))*EPMACH) DO 97 J=1,M1 S = S+DLOG(RF(J)) 97 CONTINUE IF(S.LT.DLOG(EPS))THEN 98 FORMAT('This ','boundary ','value ','problem ', * 'can ','also ','be ','solved ','with IOPT(3)=0 set',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,98) ENDIF ENDIF C ---------------------------------------------------------- C 5.2 Fail exit messages C ---------------------------------------------------------- C 5.2.1 Gaussian decomposition failed by singular jacobian IF(INFO.EQ.-1.AND.KPRINT.GE.0)THEN 99 FORMAT('Gaussian ','elimination ','failed ','by ', * 'singular ','Jacobian',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,99) ENDIF C ---------------------------------------------------------- C 5.2.2 Termination after more than itmax iterations IF(INFO.EQ.-2.AND.KPRINT.GE.0)THEN 100 FORMAT('Iteration ','terminates ','after ','itmax ', * '=',I3,2X,'iteration ','steps') write ( lumon, '(a)' ) ' ' WRITE(LUMON,100)ITMAX ENDIF C ---------------------------------------------------------- C 5.2.3 Integrator failed in trajectory computation IF(INFO.EQ.-3.AND.KPRINT.GE.0)THEN 101 FORMAT('Integrator ','failed ','in ','trajectory ', * 'computation ',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,101) J1 =-KFLAG 102 FORMAT('BVPSOL ','terminates',/,'Subinterval',I3,1X, * 'possibly ','insert ','new ','node',D20.11,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,102)J1,TFAIL ENDIF C ---------------------------------------------------------- C 5.2.4 Convergence fail of Gauss - Newton method IF(INFO.EQ.-4.AND.KPRINT.GE.0)THEN 103 FORMAT('Gauss ','Newton ','method ','fails ','to ', * 'converge',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,103) ENDIF C ---------------------------------------------------------- C 5.2.5 Inconsistent initial data IF(INFO.EQ.-5.AND.KPRINT.GE.0)THEN 104 FORMAT('Error: ','initial ','data ','and ', * 'boundary ','conditions ','are ','inconsistent',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,104) ENDIF C ---------------------------------------------------------- C 5.2.6 Multiple shooting condition too bad -insert new C nodes IF(INFO.EQ.-6.AND.KPRINT.GE.0)THEN 105 FORMAT('Termination ','since ','Multiple ', * 'Shooting ','condition',/,' ','or ','condition ','of ', * 'Jacobian ','is ','too ','bad',/,'insert ','new ', * 'nodes',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,105) S = REDH/TOL DO 106 J=1,M1 IF(RF(J).GT.S)THEN 107 FORMAT(8X,'in ','subinterval',2X,I3,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,107)J ENDIF 106 CONTINUE ENDIF C ---------------------------------------------------------- C 5.2.7 Insufficient error tolerance for integrator IF(INFO.EQ.-7.AND.KPRINT.GE.0)THEN TOLH = DMIN1(REDH/COND,EPS)/SIGDEL RELDIF = DSQRT(TOLH/SIGDEL) 108 FORMAT('Suggested ','integrator ','accuracy',D10.1 * ,2X,//,'Suggested ','relative ','deviation ', * 'parameter',D10.1,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,108)TOLH,RELDIF 109 FORMAT('Reduce ','relative ','error ','tolerance ', * 'for ','integrator ','to',D10.1,2X,/,2X,'or ','insert ', * 'new ','nodes',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,109)TOLH S = REDH/TOL DO 110 J=1,M1 IF(RF(J).GT.S)THEN 111 FORMAT(8X,'in ','subinterval',2X,I3,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,111)J ENDIF 110 CONTINUE 112 FORMAT('Reliable ','relative ','accuracy ', * 'greater ','than',1X,D6.1,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,112)1.0D-2 ENDIF C ---------------------------------------------------------- C 5.2.8 Too small storage for sparse linear system solver IF(INFO.EQ.-9.AND.KPRINT.GE.0)THEN 113 FORMAT('Too ','small ','storage ','for ','linear ', * 'system ','solver',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,113) ENDIF C ---------------------------------------------------------- C 5.3 Common exit IF(KPRINT.GE.0)THEN C 114 FORMAT('Condition ','of ','Jacobian',D10.3,2X,/, * 'Multiple ','shooting ','condition',D10.3,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,114)COND,SIGDLH write ( lumon, '(a)' ) ' ' IF(INFO.GT.0)THEN 115 FORMAT('Solution ','data',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,115) ENDIF IF(INFO.LT.0)THEN 116 FORMAT('Final ','data',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,116) ENDIF DO 117 J=1,M 118 FORMAT(D13.5,2X) WRITE(LUMON,118)T(J) 119 FORMAT((14X,3(D20.10,1X))) WRITE(LUMON,119)(X(L1),L1=(J-1)*N+1,J*N) 117 CONTINUE ENDIF C End of exits 9999 CONTINUE RETURN END SUBROUTINE BLFCNI(IVPSOL,FCN,BC,N,M,NM,NM1,ITER,KPRINT, *HSTART,FCMIN,T,X,X1,XM,T1,XU,HH,R,TOL,FC,FCOMPT,REDUCT,KFLAG, *KOUNT,INFO,LUMON) c*********************************************************************72 c cc BLFCNI() computes the residual vector. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) EXTERNAL FCN,IVPSOL,BC INTEGER N,M,NM,NM1,ITER,KPRINT DOUBLE PRECISION HSTART,FCMIN INTEGER KOUNT,KFLAG,INFO,LUMON LOGICAL REDUCT DOUBLE PRECISION TOL,FC LOGICAL FCOMPT DOUBLE PRECISION T(M),X(NM) DOUBLE PRECISION XU(NM1),HH(NM1),R(N),X1(N),XM(N),T1(N) C: End Parameter DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) DOUBLE PRECISION HALF PARAMETER (HALF=0.5D0) INTEGER J,J1,KB,KB2 DOUBLE PRECISION HMAX,HSAVE,TJ,TJ1,H INTEGER L1 C: Begin of Segment FcnInt.Body C Computation of the trajectories (solution of M1 initial C value problems) REDUCT = .FALSE. KOUNT = KOUNT+1 HSAVE = HSTART DO 117 J=1,M-1 J1 = J+1 TJ = T(J) TJ1 = T(J1) H = HSAVE HMAX = DABS(TJ1-TJ) KFLAG = 0 KB =(J-1)*N C: Begin SetVec.Vec DO 118 L1=1,N T1(L1)=X(L1+KB) 118 CONTINUE C. End SetVec.Vec CALL IVPSOL(N,FCN,TJ,T1,TJ1,TOL,HMAX,H,KFLAG) HSAVE = H IF(H.EQ.ZERO)THEN C trajectory computation failed IF(ITER.EQ.0)THEN INFO = -3 GOTO 9993 ENDIF IF(KPRINT.GE.0)THEN 119 FORMAT('trajectory ','computation ','failed, ', * 'relaxation ','factor ','or ','pseudo-rank ','reduced',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,119) ENDIF FC = FC*HALF IF(FC.LT.FCMIN)THEN REDUCT = .TRUE. GOTO 9993 ENDIF FCOMPT = .FALSE. GOTO 9993 ENDIF FCOMPT = .TRUE. C continuity conditions C: Begin SetVec.Vec DO 120 L1=1,N XU(L1+KB)=T1(L1) 120 CONTINUE C. End SetVec.Vec KB2 = KB+N C: Begin SetVec.Vec-Vec DO 121 L1=1,N HH(L1+KB)=T1(L1)-X(L1+KB2) 121 CONTINUE C. End SetVec.Vec-Vec 117 CONTINUE C two-point boundary conditions C: Begin SetVec.Vec DO 122 L1=1,N XM(L1)=X(L1+NM1) 122 CONTINUE C. End SetVec.Vec C: Begin SetVec.Vec DO 123 L1=1,N X1(L1)=X(L1) 123 CONTINUE C. End SetVec.Vec CALL BC(X1,XM,R) 9993 CONTINUE C. End of Segment FcnInt.Body RETURN END SUBROUTINE BLSOLI(N,M,M1,NM,NM1,LEVEL,NE,NB,IRANK,IRANKB, *IREPET,NYMAX,KPRINT,EPS,REDH,TOLMIN,TOL,RELDIF,EPH,EPX1H, *SIGDEL,SIGDLH,E,EH,HH,DHH,R,A,B,BG,G,QE,U,QU,DE,DU,T1,T2,US, *DX1,D,DDX,DXQ,XW,DR,RF,IROW,ICOL,ICOLB,PIVOT,NY,INFO,LUMON) c*********************************************************************72 c cc BLSOLI() seeks the least squares solution for the local approach. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) INTEGER N,M,M1,NM,NM1,LEVEL,NE,NB,IRANK,IRANKB,IREPET,NYMAX, *KPRINT DOUBLE PRECISION EPS,REDH,TOLMIN DOUBLE PRECISION TOL,RELDIF,EPH,EPX1H,SIGDEL,SIGDLH INTEGER INFO,NY,LUMON INTEGER IROW(N),ICOL(N),ICOLB(N),PIVOT(N) DOUBLE PRECISION G(N,N,M1) DOUBLE PRECISION A(N,N),B(N,N),BG(N,N),E(N,N),QE(N,N) DOUBLE PRECISION DDX(NM),DXQ(NM),XW(NM),HH(NM1),DHH(NM1),D(N), *DE(N),R(N),DR(N),U(N),DU(N),QU(N),T1(N),T2(N),DX1(N),RF(M), *US(N) DOUBLE PRECISION EH(N,N) C: End Parameter DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) DOUBLE PRECISION HALF PARAMETER (HALF=0.5D0) DOUBLE PRECISION TEN PARAMETER (TEN=1.0D1) INTEGER I,IC,I0,J,JA,JIN,JN,J1,K,K0 DOUBLE PRECISION CORR,EPDX1,EPX1,S,TH,TOLH LOGICAL NOTFND INTEGER L1,L2 DOUBLE PRECISION S1 C: Begin C: Begin of Segment SolvIn.Body IF(IREPET.GE.0)THEN C -------------------------------------------------------- C 1 computation of condensed right - hand side U(NE) IF(IRANK.GT.0) CALL BLRHS1(N,NE,M1,NM1,1,HH,R,B,G,U,DE, * T1,BG,IROW) C -------------------------------------------------------- C 2 saving of right - hand side U IF(IRANK.GE.NE)THEN C: Begin SetVec.Vec DO 124 L1=1,NE US(L1)=U(L1) 124 CONTINUE C. End SetVec.Vec ENDIF ENDIF C ---------------------------------------------------------- C 3 ( best ) constrained least squares solution of linear(NE, C NE) -system C: Vec DX1 = Scalar (for 1,N) S1 = ZERO DO 125 L1=1,N DX1(L1)=S1 125 CONTINUE C. End SetVec.S IF(IRANK.GT.0)THEN CALL BLSOLC(E,N,N,IRANKB,NE,NE,DX1,U,IRANK,D,PIVOT, * IREPET,QE,T1) ENDIF C ---------------------------------------------------------- C 4 iterative refinement of DX1 EPH = EPS IF(IRANK.GE.NE.AND.NE.NE.0)THEN C: Begin SetVec.MatxVec DO 126 L1=1,NE S1=0.0 DO 127 L2=1,NE S1=S1+EH(L1,L2)*DX1(L2) 127 CONTINUE DU(L1)=S1 126 CONTINUE C. End SetVec.MatxVec C: Begin SetVec.Vec-Vec DO 128 L1=1,NE DU(L1)=US(L1)-DU(L1) 128 CONTINUE C. End SetVec.Vec-Vec C Solution of residual equation CALL BLSOLC(E,N,N,IRANKB,NE,NE,T2,DU,IRANK,D,PIVOT, * IREPET,QE,T1) C: EPDX1 = Max.Norm of Vec T2 (for 1,NE) EPDX1 = 0.0 DO 129 L1=1,NE S1 = DABS(T2(L1)) IF(S1.GT.EPDX1) EPDX1=S1 129 CONTINUE C. End MakeMaxNorm.Vec C: Vec DX1 = Vec DX1 + Vec T2 (for 1,NE) DO 130 L1=1,NE DX1(L1)=DX1(L1)+T2(L1) 130 CONTINUE C. End SetVec.&Vec C: EPX1 = Max.Norm of Vec DX1 (for 1,NE) EPX1 = 0.0 DO 131 L1=1,NE S1 = DABS(DX1(L1)) IF(S1.GT.EPX1) EPX1=S1 131 CONTINUE C. End MakeMaxNorm.Vec EPX1H = EPDX1/EPX1 EPH = TEN*EPDX1 IF(EPX1H.GT.HALF)THEN INFO = -8 GOTO 9992 ENDIF ENDIF C ---------------------------------------------------------- C 5 Descaling and back - permutation of solution DX1 C: Permutation ICOL of Vec DXQ = Vec DX1 (for 1,N) DO 132 L1=1,N L2 = ICOL(L1) DXQ(L2) = DX1(L1) 132 CONTINUE C. End SetVecByPermVec.Vec C: Vec DXQ = Vec DXQ * Vec XW (for 1,N) DO 133 L1=1,N DXQ(L1)=DXQ(L1)*XW(L1) 133 CONTINUE C. End SetVec.VecxVec C ---------------------------------------------------------- C 6 Recursive computation of DXQ(N,2), ... , DXQ(N,M) CALL BLRCRS(N,M,M1,NM,NM1,1,HH,G,DXQ,T1,T2) C ---------------------------------------------------------- C 1 Iterative refinement sweeps NY = 1 , ... , NYMAX NY = 0 SIGDEL = TEN*TEN SIGDLH = ZERO IF(EPH.LT.EPS) EPH = EPS IF(NYMAX.NE.0.AND.IRANK.GE.NE.AND.NE.NE.0)THEN IF(KPRINT.GT.0)THEN 134 FORMAT('Iterative ','refinement',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,134) ENDIF C -------------------------------------------------------- C 1.1 Computation of required continuity residuals DHH(N, C M1) JN = 1 JIN = M C -------------------------------------------------------- C 1.2 Computation of boundary residual DR(N) C: DO (Until) 135 CONTINUE C: Begin SetVec.MatxVec DO 136 L1=1,N S1=0.0 DO 137 L2=1,N S1=S1+A(L1,L2)*DXQ(L2) 137 CONTINUE DR(L1)=S1 136 CONTINUE C. End SetVec.MatxVec C: Begin SetVec.MatxVec DO 138 L1=1,N S1=0.0 DO 139 L2=1,N S1=S1+B(L1,L2)*DXQ(L2+NM1) 139 CONTINUE T1(L1)=S1 138 CONTINUE C. End SetVec.MatxVec C: Vec DR = Formula (for 1,N) DO 140 L1=1,N DR(L1)=R(L1)+DR(L1)+T1(L1) 140 CONTINUE C. End SetVec.Comp C ------------------------------------------------------ C 1.3 Computation of condensed residual DU(NE) IF(IRANK.GT.0) CALL BLRHS1(N,NE,M1,NM1,JIN,DHH,DR,B,G, * DU,DE,T1,BG,IROW) C ------------------------------------------------------ C 1.4 Computation of correction DDX(N) C: Vec DX1 = Scalar (for 1,N) S1 = ZERO DO 141 L1=1,N DX1(L1)=S1 141 CONTINUE C. End SetVec.S IF(IRANK.GT.0) CALL BLSOLC(E,N,N,IRANKB,NE,NE,DX1,DU, * IRANK,D,PIVOT,IREPET,QE,T1) C ------------------------------------------------------ C 2 Descaling of DDX(N), refinement of DXQ(N) C: CORR = Max.Norm of Vec DX1 (for 1,N) CORR = 0.0 DO 142 L1=1,N S1 = DABS(DX1(L1)) IF(S1.GT.CORR) CORR=S1 142 CONTINUE C. End MakeMaxNorm.Vec C: Permutation ICOL of Vec T1 = Vec DX1 (for 1,N) DO 143 L1=1,N L2 = ICOL(L1) T1(L2) = DX1(L1) 143 CONTINUE C. End SetVecByPermVec.Vec C: Vec DDX = Vec T1 * Vec XW (for 1,N) DO 144 L1=1,N DDX(L1)=T1(L1)*XW(L1) 144 CONTINUE C. End SetVec.VecxVec C: Vec DXQ = Vec DXQ + Vec DDX (for 1,N) DO 145 L1=1,N DXQ(L1)=DXQ(L1)+DDX(L1) 145 CONTINUE C. End SetVec.&Vec IF(CORR.GE.EPH)THEN EPH = CORR INFO = -8 GOTO 9992 ENDIF RF(1)=CORR C ------------------------------------------------------ C 3 Recursive computation of DDX(N+1), ... , DDX(NM) CALL BLRCRS(N,M,M1,NM,NM1,JIN,DHH,G,DDX,T1,T2) C ------------------------------------------------------ C 3.1 Refinement of DXQ(N+1), ... , DXQ(NM) DO 146 J=2,M I0 =(J-1)*N C: Vec DXQ = Vec DXQ + Vec DDX (for I0+1,I0+N) DO 147 L1=I0+1,I0+N DXQ(L1)=DXQ(L1)+DDX(L1) 147 CONTINUE C. End SetVec.&Vec C: CORR = Max of Formula Elements (for I0+1,I0+N) CORR =-7.23D75 DO 148 L1=I0+1,I0+N S1=DABS(DDX(L1)/XW(L1)) IF(S1.GT.CORR)CORR=S1 148 CONTINUE C. End MakeMax.Comp RF(J)=CORR 146 CONTINUE C ------------------------------------------------------ C 3.2 Determination of sweep index JN JA = JN J = 1 NOTFND = .TRUE. C: While (expression) 149 IF(J.LE.M.AND.NOTFND)THEN IF(RF(J).GT.EPH)THEN NOTFND = .FALSE. ELSE JN = J J = J+1 ENDIF GOTO 149 ENDIF C. EndWhile NY = NY+1 IF(KPRINT.GT.0)THEN 150 FORMAT('Sweep',1X,I3,1X,'starts ','at',1X,I3) write ( lumon, '(a)' ) ' ' WRITE(LUMON,150)NY,JA 151 FORMAT((1X,5(D12.3,1X))) WRITE(LUMON,151)(RF(L1),L1=1,M) ENDIF IF(JN.LE.JA)THEN INFO = -6 GOTO 9992 ENDIF IF(JN.NE.M) JIN = JN C ------------------------------------------------------ C 3.3 Determination and adaptation of parameters TOL AND C RELDIF IF(LEVEL.NE.0.AND.NY.LE.1)THEN DO 152 J=1,M1 S = RF(J+1)/RF(J) IF(SIGDLH.LT.S) SIGDLH = S RF(J)=S 152 CONTINUE IF(KPRINT.GT.0)THEN 153 FORMAT('Norms ','of ','wronskians') write ( lumon, '(a)' ) ' ' WRITE(LUMON,153) 154 FORMAT((1X,5(D12.3,1X))) WRITE(LUMON,154)(RF(L1),L1=1,M1) ENDIF SIGDEL = DMAX1(SIGDLH,SIGDEL) TH = TOL*SIGDEL IF(TH.GT.REDH)THEN INFO = -7 GOTO 9992 ENDIF IF(TH.GT.EPH) EPH = TH TOLH = EPS/SIGDEL IF(TOLH.LT.TOLMIN) TOLH = TOLMIN CWei; TOL = TOLH CWei; RELDIF = DSQRT(TOL/SIGDEL) IF(KPRINT.GE.0)THEN 155 FORMAT('Suggested ','integrator ','accuracy',D10.1 * ,2X,/,'Suggested ','relative ','deviation ', * 'parameter',D10.1,2X,//,'Adapted ','in ', * 'the ','next ','iteration ','step',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,155)TOLH,RELDIF ENDIF ENDIF IF(JN.NE.M)THEN DO 156 J=JN,M1 J1 = J+1 DO 157 I=1,N K0 =(J-1)*N S = HH(I+K0) DO 158 K=1,N S = S+G(I,K,J)*DXQ(K+K0) 158 CONTINUE DHH(I+K0)=S-DXQ(I+K0+N) 157 CONTINUE 156 CONTINUE ENDIF IF(.NOT.(JN.EQ.M)) GOTO 135 C. UNTIL ( expression - negated above) ENDIF C End of iterative refinement sweeps C ---------------------------------------------------------- C 4 Projection of separated linear boundary conditions at T( C M) IF(NB.NE.0)THEN DO 159 K=1,NB IC = ICOLB(K) DXQ(IC+NM1)=ZERO 159 CONTINUE ENDIF 9992 CONTINUE C. End of Segment SolvIn.Body RETURN END SUBROUTINE BGSOLI(N,M,M1,NM,NM1,NDIM,LICN,NKEEP,NA,NAQ,NB, *NBQ,NRS,ITER,LEVEL,KPRINT,EPS,REDH,TOLMIN,FC,FCA,TOL,RELDIF, *EPH,EPX1H,SIGDEL,SIGDLH,COND,CORR,HH,DHH,R,A,B,G,U,DE,DU,T1, *DXQ,XW,DR,RF,WO,E,IROW,ICOLA,ICOLB,ICN,IKEEP,INFO,LUMON) c*********************************************************************72 c cc BGSOLI() seeks the least squares solution for the global approach. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) INTEGER N,M,M1,NM,NM1,NDIM,LICN,NKEEP,NA,NAQ,NB,NBQ,NRS,ITER, *LEVEL,KPRINT,LUMON DOUBLE PRECISION EPS,REDH,TOLMIN,FC,FCA DOUBLE PRECISION TOL,RELDIF,EPH,EPX1H,SIGDEL,SIGDLH,COND, *CORR INTEGER INFO INTEGER IROW(N),ICOLA(N),ICOLB(N) INTEGER ICN(LICN),IKEEP(NKEEP) DOUBLE PRECISION G(N,N,M1) DOUBLE PRECISION A(N,N),B(N,N) DOUBLE PRECISION DXQ(NM),XW(NM),HH(NM1),DHH(NM1),DE(N),R(N), *DR(N),U(NM),DU(NM),T1(N),RF(M),E(LICN),WO(NM) C: End Parameter C: EPMACH = relative machine precision DOUBLE PRECISION EPMACH,SMALL DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) DOUBLE PRECISION HALF PARAMETER (HALF=0.5D0) DOUBLE PRECISION ONE PARAMETER (ONE=1.0D0) DOUBLE PRECISION TEN PARAMETER (TEN=1.0D1) INTEGER I,IC,I0,J,J1,K,MTYPE DOUBLE PRECISION S,TOLH REAL ST INTEGER L1,L2 DOUBLE PRECISION S1 C: Begin C: Begin of Segment SolvIn.Body C ---------------------------------------------------------- C 1 Computation of right - hand side U(NDIM) CALL ZIBCONST(EPMACH,SMALL) DO 127 J=1,M1 J0 =(J-1)*N J1 = J0+N DO 128 I=1,N U(I+J0)=HH(I+J0)/XW(I+J1) 128 CONTINUE 127 CONTINUE IF(NRS.NE.0)THEN DO 129 I=1,NRS I0 = IROW(NB+I) U(NM1+I)=DE(I0)*R(I0) 129 CONTINUE ENDIF C ---------------------------------------------------------- C 2 Solution of linear(NDIM,NDIM) -system MTYPE = 1 CALL MA28CD(NDIM,E,LICN,ICN,IKEEP,U,WO,MTYPE) COND = ZERO IF(NAQ.NE.0)THEN DO 130 I=1,NAQ IC = ICOLA(I) S = U(I) DXQ(IC)=S*XW(IC) COND = COND+DABS(S) 130 CONTINUE ENDIF IF(M.NE.2)THEN DO 131 J=2,M1 J0 =(J-1)*N DO 132 I=1,N S = U(J0-NA+I) DXQ(I+J0)=S*XW(I+J0) COND = COND+DABS(S) 132 CONTINUE 131 CONTINUE ENDIF IF(NBQ.NE.0)THEN DO 133 I=1,NBQ IC = ICOLB(NB+I) S = U(NM1-NA+I) DXQ(IC+NM1)=S*XW(IC+NM1) COND = COND+DABS(S) 133 CONTINUE ENDIF IF(NB.NE.0)THEN DO 134 K=1,NB IC = ICOLB(K) DXQ(IC+NM1)=ZERO 134 CONTINUE ENDIF IF(NA.NE.0)THEN DO 135 K=1,NA IC = ICOLA(NAQ+K) DXQ(IC)=ZERO 135 CONTINUE ENDIF C ---------------------------------------------------------- C 3 Iterative refinement IF(KPRINT.GT.0)THEN 136 FORMAT('Iterative ','refinement',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,136) ENDIF C ---------------------------------------------------------- C 4 Computation of required continuity residuals DHH(NM1) DO 137 J=1,M1 J0 =(J-1)*N J1 = J0+N DO 138 I=1,N S = HH(I+J0) DO 139 K=1,N S = S+G(I,K,J)*DXQ(K+J0) 139 CONTINUE DHH(I+J0)=S-DXQ(I+J1) 138 CONTINUE 137 CONTINUE C ---------------------------------------------------------- C 5 Computation of boundary residual DR(N) C: Begin SetVec.MatxVec DO 140 L1=1,N S1=0.0 DO 141 L2=1,N S1=S1+A(L1,L2)*DXQ(L2) 141 CONTINUE DR(L1)=S1 140 CONTINUE C. End SetVec.MatxVec C: Begin SetVec.MatxVec DO 142 L1=1,N S1=0.0 DO 143 L2=1,N S1=S1+B(L1,L2)*DXQ(L2+NM1) 143 CONTINUE T1(L1)=S1 142 CONTINUE C. End SetVec.MatxVec C: Vec DR = Formula (for 1,N) DO 144 L1=1,N DR(L1)=R(L1)+DR(L1)+T1(L1) 144 CONTINUE C. End SetVec.Comp C ---------------------------------------------------------- C 6 Computation of residual DU(NDIM) DO 145 J=1,M1 J0 =(J-1)*N J1 = J0+N DO 146 I=1,N DU(J0+I)=DHH(I+J0)/XW(I+J1) 146 CONTINUE 145 CONTINUE IF(NRS.NE.0)THEN DO 147 I=1,NRS I0 = IROW(NB+I) DU(NM1+I)=DE(I0)*DR(I0) 147 CONTINUE ENDIF C ---------------------------------------------------------- C 7 Computation of correction MTYPE = 1 CALL MA28CD(NDIM,E,LICN,ICN,IKEEP,DU,WO,MTYPE) C ---------------------------------------------------------- C 8 Refinement of DXQ CORR = ZERO IF(NAQ.NE.0)THEN DO 148 I=1,NAQ IC = ICOLA(I) S = DU(I) DXQ(IC)=DXQ(IC)+S*XW(IC) CORR = CORR+DABS(S) 148 CONTINUE ENDIF IF(M.NE.2)THEN DO 149 J=2,M1 J0 =(J-1)*N DO 150 I=1,N S = DU(J0-NA+I) DXQ(I+J0)=DXQ(I+J0)+S*XW(I+J0) CORR = CORR+DABS(S) 150 CONTINUE 149 CONTINUE ENDIF IF(NBQ.NE.0)THEN DO 151 I=1,NBQ IC = ICOLB(NB+I) S = DU(NM1-NA+I) DXQ(IC+NM1)=DXQ(IC+NM1)+S*XW(IC+NM1) CORR = CORR+DABS(S) 151 CONTINUE ENDIF COND = CORR/(COND*EPMACH) IF(KPRINT.GT.0)THEN 152 FORMAT('Norm ','of ','residual',D12.3,2X) write ( lumon, '(a)' ) ' ' WRITE(LUMON,152)CORR ENDIF C End of iterative refinement C ---------------------------------------------------------- C 9 Determination and adaptation of parameters TOL and C RELDIF IF(LEVEL.NE.0)THEN DO 153 J=1,M1 J0 =(J-1)*N J1 = J0+N SIGDEL = ZERO ST = ZERO DO 154 I=1,N II = N IF(J.EQ.1) II = NAQ S = DABS(DXQ(I+J0))/XW(I+J0) IF(ST.LT.S) ST = S S = ZERO DO 155 K=1,II KK = K IF(J.EQ.1) KK = ICOLA(KK) S = S+G(I,KK,J)*DXQ(KK+J0)/XW(I+J1) 155 CONTINUE S = DABS(S) IF(SIGDEL.LT.S) SIGDEL = S 154 CONTINUE RF(J)=SIGDEL/ST+ONE 153 CONTINUE IF(KPRINT.GT.0)THEN 156 FORMAT('Norms ','of ','wronskians') write ( lumon, '(a)' ) ' ' WRITE(LUMON,156) 157 FORMAT((1X,5(D12.3,1X))) WRITE(LUMON,157)(RF(L1),L1=1,M1) ENDIF SIGDLH = ZERO DO 158 J=1,M1 IF(SIGDLH.LT.RF(J)) SIGDLH = RF(J) 158 CONTINUE SIGDEL = SIGDLH IF(FC.EQ.ONE.AND.FCA.EQ.ONE.AND.ITER.GT.0) SIGDEL = * SIGDLH*COND SIGDEL = DMAX1(SIGDEL,TEN) EPH = TOL*SIGDEL IF(EPH.GT.REDH)THEN INFO = -7 GOTO 9989 ENDIF TOLH = EPS/SIGDEL IF(TOLH.LT.TOLMIN) TOLH = TOLMIN CWEI; TOL = TOLH CWEI; RELDIF = DSQRT(TOL/SIGDEL) IF(KPRINT.GE.0)THEN 159 FORMAT('Suggested ','integrator ','accuracy',D10.1 * ,2X,/,'Suggested ','relative ','deviation ', * 'parameter',D10.1,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,159)TOLH,RELDIF 160 FORMAT('Adapted ','in ','the ','next ', * 'iteration ','step',/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,160) ENDIF ENDIF 9989 CONTINUE C. End of Segment SolvIn.Body RETURN END SUBROUTINE BLPRCD(LUMON,COND,SENS,SMALIN,J,IRANK) c*********************************************************************72 c cc BLPRCD() prints the subcondition and sensitivity information. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) DOUBLE PRECISION COND,SENS,SMALIN INTEGER J,IRANK,LUMON C: End Parameter DOUBLE PRECISION ONE PARAMETER (ONE=1.0D0) DOUBLE PRECISION SENSP IF(SENS.LT.ONE)THEN SENSP = SENS*SMALIN 160 FORMAT('Subcondition (',I2,',',I2,') ',D10.3,2X,/, * 'Sensitivity (',I2,',',I2,') ',D10.3,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,160)J,IRANK,COND,J,IRANK,SENSP ELSE 161 FORMAT('Subcondition ','(',I2,',',I2,') ',D10.3,2X,/, * 'Sensitivity ','(',I2,',',I2,') ',D10.3,2X,' *',D7.0 * ,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,161)J,IRANK,COND,J,IRANK,SENS,SMALIN ENDIF return END SUBROUTINE BLPRCV(LUMON,CONV,EPH) c*********************************************************************72 c cc BLPRCV() prints the relative accuracy information. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) DOUBLE PRECISION CONV,EPH INTEGER LUMON 162 FORMAT('Achieved ','relative ','accuracy',D10.3,2X) write ( lumon, '(a)' ) ' ' WRITE(LUMON,162)CONV IF(EPH.GT.CONV) CONV = EPH 163 FORMAT('Reliable ','relative ','accuracy',D10.3,2X,/) write ( lumon, '(a)' ) ' ' WRITE(LUMON,163)CONV return END SUBROUTINE BLSCLE(N,M,NM,NM1,X,XU,XW,XTHR) c*********************************************************************72 c cc BLSCLE() scales the variables. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) INTEGER N INTEGER M INTEGER NM INTEGER NM1 DOUBLE PRECISION X(NM) DOUBLE PRECISION XW(NM) DOUBLE PRECISION XU(NM1) DOUBLE PRECISION XTHR C: End Parameter C PROVIDES SCALING XW(NM)OF VARIABLES X(NM) C: EPMACH = relative machine precision DOUBLE PRECISION EPMACH,SMALL DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) DOUBLE PRECISION HALF PARAMETER (HALF=0.5D0) DOUBLE PRECISION RED PARAMETER (RED=1.0D-2) INTEGER I,J,J0,J1 DOUBLE PRECISION XMAX INTEGER L1 C: Begin C: Vec XW = Formula (for 1,N) CALL ZIBCONST(EPMACH,SMALL) DO 164 L1=1,N XW(L1)=DABS(X(L1)) 164 CONTINUE C. End SetVec.Comp C ------------------------------------------------------------ C 1 Arithmetic mean for XW(2*N)... XW(M*N) DO 165 J=1,M-1 J0 =(J-1)*N J1 = J0+N DO 166 I=1,N XW(I+J1)=(DABS(X(I+J1))+DABS(XU(I+J0)))*HALF 166 CONTINUE 165 CONTINUE C ------------------------------------------------------------ C 2 Threshold DO 167 I=1,N XMAX = ZERO DO 168 J=0,NM1,N IF(XMAX.LT.XW(I+J)) XMAX = XW(I+J) 168 CONTINUE XMAX = XMAX*RED IF(XMAX.LT.XTHR) XMAX = XTHR DO 169 J=0,NM1,N IF(XW(I+J).LT.XMAX) XW(I+J)=XMAX 169 CONTINUE 167 CONTINUE RETURN END SUBROUTINE BLLVLS(N,M,NM,NM1,XW,DXQ,HH,R,DE,CONV,SUMX,SUMF, *KPRINT) c*********************************************************************72 c cc BLLVLS() evaluates level functions for the local approach. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) C INTEGER N,M,NM,NM1,KPRINT DOUBLE PRECISION XW(NM),DXQ(NM),HH(NM1),R(N),DE(N) DOUBLE PRECISION CONV,SUMX,SUMF C: End Parameter C: SMALL = squareroot of "smallest positive machine number C divided by relative machine precision" DOUBLE PRECISION SMALL PARAMETER (SMALL=4.94D-32) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) INTEGER I,J,J0,J1 DOUBLE PRECISION S INTEGER L1 C: Begin C ------------------------------------------------------------ C 1 Evaluation of scaled natural level function SUMX and C scaled maximum error norm CONV CONV = ZERO SUMX = ZERO DO 170 J=1,NM S = DABS(DXQ(J)/XW(J)) IF(CONV.LT.S) CONV = S SUMX = SUMX+S*S 170 CONTINUE C ------------------------------------------------------------ C 2 Evaluation of (scaled) standard level function sumfs (only C if needed for print) C: SUMF = Sum of Formula Elements (for 1,N) SUMF = 0.0D0 DO 171 L1=1,N SUMF=SUMF+((R(L1)*DE(L1)/SMALL)**2) 171 CONTINUE C. End MakeSum.Comp DO 172 J=1,M-1 J0 =(J-1)*N J1 = J0+N DO 173 I=1,N SUMF = SUMF+(HH(I+J0)/XW(I+J1))**2 173 CONTINUE 172 CONTINUE RETURN END SUBROUTINE BGLVLS(N,M,NM,NM1,XW,DXQ,HH,R,DE,CONV,SUMX,SUMF, *KPRINT) c*********************************************************************72 c cc BGLVLS() evaluates level functions for the global approach. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) C INTEGER N,M,NM,NM1,KPRINT DOUBLE PRECISION XW(NM),DXQ(NM),HH(NM1),R(N),DE(N) DOUBLE PRECISION CONV,SUMX,SUMF C: End Parameter C: SMALL = squareroot of "smallest positive machine number C divided by relative machine precision" DOUBLE PRECISION SMALL PARAMETER (SMALL=4.94D-32) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) INTEGER I,J,J0,J1 DOUBLE PRECISION S INTEGER L1 C: Begin C ------------------------------------------------------------ C 1 Evaluation of scaled natural level function SUMX and C scaled maximum error norm CONV CONV = ZERO SUMX = ZERO DO 167 J=1,NM S = DABS(DXQ(J)/XW(J)) IF(CONV.LT.S) CONV = S SUMX = SUMX+S*S 167 CONTINUE C ------------------------------------------------------------ C 2 Evaluation of (scaled) standard level function sumfs (only C if needed for print) C: SUMF = Sum of Formula Elements (for 1,N) SUMF = 0.0D0 DO 168 L1=1,N SUMF=SUMF+((R(L1)*DE(L1))**2) 168 CONTINUE C. End MakeSum.Comp DO 169 J=1,M-1 J0 =(J-1)*N J1 = J0+N DO 170 I=1,N SUMF = SUMF+(HH(I+J0)/XW(I+J1))**2 170 CONTINUE 169 CONTINUE RETURN END SUBROUTINE BLRHS1(N,NE,M1,NM1,JIN,HH,R,B,G,U,DE,V,BG,IROW) c*********************************************************************72 c cc BLRHS1() computes the condensed right hand side. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) INTEGER N,NE,M1,NM1,JIN DOUBLE PRECISION HH(NM1),R(N) DOUBLE PRECISION B(N,N) DOUBLE PRECISION G(N,N,M1) DOUBLE PRECISION U(N),DE(N),V(N) DOUBLE PRECISION BG(N,N) INTEGER IROW(N) C: End Parameter C Computation of condensed right-hand side U(NE) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) INTEGER I,IR,J,JJ,J1,K,K0,L,M2 DOUBLE PRECISION S,TH C: Begin DO 174 I=1,NE IR = IROW(I) U(I)=DE(IR)*R(IR) 174 CONTINUE IF(JIN.GT.M1)THEN RETURN ENDIF DO 175 I=1,NE IR = IROW(I) S = U(I) K0 = NM1-N DO 176 K=1,N TH = DE(IR)*B(IR,K) BG(I,K)=TH S = S+TH*HH(K+K0) 176 CONTINUE U(I)=S 175 CONTINUE IF(M1.EQ.1.OR.JIN.EQ.M1) RETURN M2 = M1-1 DO 177 JJ=JIN,M2 J = M2+JIN-JJ J1 = J+1 DO 178 I=1,NE DO 179 K=1,N S = ZERO DO 180 L=1,N S = S+BG(I,L)*G(L,K,J1) 180 CONTINUE V(K)=S 179 CONTINUE S = U(I) K0 =(J-1)*N DO 181 K=1,N S = S+V(K)*HH(K+K0) BG(I,K)=V(K) 181 CONTINUE U(I)=S 178 CONTINUE 177 CONTINUE RETURN END SUBROUTINE BLRCRS(N,M,M1,NM,NM1,JIN,HH,G,DX,U,V) c*********************************************************************72 c cc BLRCRS() carries out the recursive solution of a system of equations. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) INTEGER N,M,M1,NM,NM1,JIN DOUBLE PRECISION HH(NM1) DOUBLE PRECISION G(N,N,M1) DOUBLE PRECISION DX(NM),U(N),V(N) C: End Parameter C Recursive solution of m1 linear(N,N)-systems DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) INTEGER I,J,J0,J1,K INTEGER L1 DOUBLE PRECISION S C: Begin C: Begin SetVec.Vec DO 182 L1=1,N U(L1)=DX(L1) 182 CONTINUE C. End SetVec.Vec DO 183 J=1,M1 J0 =(J-1)*N J1 = J0+N DO 184 I=1,N IF(J.GE.JIN)THEN S = HH(I+J0) ELSE S = ZERO ENDIF DO 185 K=1,N S = S+G(I,K,J)*U(K) 185 CONTINUE V(I)=S DX(I+J1)=S 184 CONTINUE C: Begin SetVec.Vec DO 186 L1=1,N U(L1)=V(L1) 186 CONTINUE C. End SetVec.Vec 183 CONTINUE RETURN END SUBROUTINE BLPRJC(N,NE,IRANK,DEL,U,D,V,QE,PIVOT) c*********************************************************************72 c cc BLPRJC() projects the reduced component. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) C INTEGER IRANK,N,NE INTEGER PIVOT(N) DOUBLE PRECISION DEL DOUBLE PRECISION QE(N,N) DOUBLE PRECISION U(N),D(N),V(N) C: End Parameter DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) INTEGER L1,L2 INTEGER I,IRK1 DOUBLE PRECISION S,SH C: Begin DO 187 I=1,NE V(I)=U(PIVOT(I)) 187 CONTINUE IRK1 = IRANK+1 DEL = ZERO DO 188 I=IRK1,NE C: SH = Col I of QE * Vec V (for 1,I-1) SH = 0.0 DO 189 L1=1,I-1 SH = SH+QE(L1,I)*V(L1) 189 CONTINUE C. End MakeSProd.ColxVec S =(V(I)-SH)/D(I) DEL = DEL-S*S V(I)=S 188 CONTINUE DO 190 I=IRK1,NE K = NE+IRK1-I S = V(K) IF(K.NE.NE)THEN C: SH = Row K of QE * Vec V (for K+1,NE) SH = 0.0 DO 191 L1=K+1,NE SH=SH+QE(K,L1)*V(L1) 191 CONTINUE C. End MakeSProd.RowxVec S = S-SH ENDIF S = S/D(K) V(K)=S 190 CONTINUE DO 192 I=1,IRANK C: S = Row I of QE * Vec V (for IRK1,NE) S = 0.0 DO 193 L1=IRK1,NE S=S+QE(I,L1)*V(L1) 193 CONTINUE C. End MakeSProd.RowxVec V(I)=-S 192 CONTINUE C: Permutation PIVOT of Vec U = Vec V (for 1,NE) DO 194 L1=1,NE L2 = PIVOT(L1) U(L2) = V(L1) 194 CONTINUE C. End SetVecByPermVec.Vec RETURN END SUBROUTINE BLDERA(BC,N,M,NM,XW,X1,XM,R,RH,A,B,RELDIF) c*********************************************************************72 c cc BLDERA() approximates boundary derivative matrices using differences. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) EXTERNAL BC INTEGER N,M,NM DOUBLE PRECISION XW(NM),X1(N),XM(N),R(N),RH(N) DOUBLE PRECISION A(N,N),B(N,N) DOUBLE PRECISION RELDIF c C Difference approx. of boundary derivative matrices A(N,N)and C B(N,N) DOUBLE PRECISION ONE PARAMETER (ONE=1.0D0) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) INTEGER I,K,NM1 DOUBLE PRECISION S,XH C: Begin NM1 = N*(M-1) DO 195 K=1,N XH = X1(K) S = RELDIF*XW(K) IF(XH.LT.ZERO) S =-S X1(K)=XH+S CALL BC(X1,XM,RH) X1(K)=XH S = ONE/S DO 196 I=1,N A(I,K)=(RH(I)-R(I))*S 196 CONTINUE XH = XM(K) S = RELDIF*XW(K+NM1) IF(XH.LT.ZERO) S =-S XM(K)=XH+S CALL BC(X1,XM,RH) XM(K)=XH S = ONE/S DO 197 I=1,N B(I,K)=(RH(I)-R(I))*S 197 CONTINUE 195 CONTINUE RETURN END SUBROUTINE BLDERG(FCN,N,NE,M,M1,NM,NM1,T,X,XU,XW,XJ,TJ,G, *ICOL,IVPSOL,HSTART,TOL,RELDIF,KFLAG) c*********************************************************************72 c cc BLDERG() estimates Wronskian matrices using differences. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) INTEGER N INTEGER NE INTEGER M INTEGER M1 INTEGER NM INTEGER NM1 DOUBLE PRECISION T(M) DOUBLE PRECISION X(NM) DOUBLE PRECISION XU(NM1) DOUBLE PRECISION XW(NM) DOUBLE PRECISION XJ(N) DOUBLE PRECISION TJ DOUBLE PRECISION G(N,N,M1) INTEGER ICOL(N) EXTERNAL IVPSOL DOUBLE PRECISION HSTART DOUBLE PRECISION TOL DOUBLE PRECISION RELDIF INTEGER KFLAG C: End Parameter C Difference approximation of Wronskian Matrices G(1),.., G(M1) DOUBLE PRECISION ONE PARAMETER (ONE=1.0D0) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) INTEGER I,IK,J,J0,J1,K DOUBLE PRECISION HMAX,H,HSAVE,S,TJ1,TJA,TH EXTERNAL FCN C: Begin HSAVE = HSTART DO 198 J=1,M-1 J0 =(J-1)*N J1 = J+1 TJA = T(J) TJ1 = T(J1) HMAX = DABS(TJ1-TJA) DO 199 IK=1,N I = ICOL(IK) H = HSAVE IF(J.NE.1.OR.IK.LE.NE)THEN TJ = TJA KFLAG = 0 DO 200 K=1,N XJ(K)=X(K+J0) 200 CONTINUE TH = XJ(I) S = RELDIF*XW(I+J0) IF(TH.LT.ZERO) S =-S XJ(I)=TH+S S = ONE/S CALL IVPSOL(N,FCN,TJ,XJ,TJ1,TOL,HMAX,H,KFLAG) IF(H.EQ.ZERO)THEN KFLAG =-J RETURN ENDIF DO K=1,N G(K,I,J)=S*(XJ(K)-XU(K+J0)) end do ENDIF 199 CONTINUE HSAVE = H 198 CONTINUE KFLAG = 0 RETURN END SUBROUTINE BLRK1G(N,M,M1,NM,NM1,XW,DX,HH,HHA,DXJ,G,FCA) c*********************************************************************72 c cc BLRK1G() performs a rank-1 update of the Wronskian matrices. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c IMPLICIT double precision(S) INTEGER N,M,M1,NM,NM1 DOUBLE PRECISION XW(NM),DX(NM),HH(NM1),HHA(NM1),DXJ(N) DOUBLE PRECISION G(N,N,M1) DOUBLE PRECISION FCA C: End Parameter C RANK-1 UPDATES OF WRONSKIAN MATRICES G(1),..., G(M1) DOUBLE PRECISION ONE PARAMETER (ONE=1.0D0) DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) INTEGER I,J,J0,K DOUBLE PRECISION DNM,FCH,S,T FCH = FCA-ONE DO J=1,M1 J0 =(J-1)*N DNM = ZERO DO I=1,N T = DX(I+J0)/XW(I+J0) DXJ(I)=T/XW(I+J0) DNM = DNM+T*T end do DNM = DNM*FCA IF(DNM.NE.ZERO)THEN DO K=1,N T = DXJ(K)/DNM DO I=1,N S = G(I,K,J) IF(S.NE.ZERO) G(I,K,J)=S+T*(HH(I+J0)+FCH*HHA(I+J0)) end do end do ENDIF end do RETURN END SUBROUTINE BLDFX1 (N,FCN,T,Y,TEND,TOL,HMAX,H,KFLAG) c*********************************************************************72 c cc BLDFX1() - an explicit extrapolation integrator for nonstiff ODE's. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c C* Begin prologue DIFEX1 C C --------------------------------------------------------------------- C C* Title C C Explicit extrapolation integrator for non-stiff systems of C ordinary first-order differential equations. C C* Written by P. Deuflhard, U. Nowak, U. Poehle C Adapted by L. Weimann for use with BVPSOL C* Purpose Solution of systems of initial value problems C* Method Explicit mid-point rule discretization with C h**2-extrapolation C* Category i1a1c1. - System of nonstiff first order C differential equations C* Keywords extrapolation, ODE, explicit mid-point rule, C nonstiff C* Version 1.2 , August 1991 C* Latest Change January 2004 C* Library CodeLib C* Code Fortran 77 C Double Precision C* Environment Standard version for FORTRAN77 environments on C PCs, workstations, and hosts C* Copyright (c) Konrad-Zuse-Zentrum fuer Informationstechnik C Berlin (ZIB) C Takustrasse 7, D-14195 Berlin-Dahlem C phone : + 49/30/84185-0 C fax : + 49/30/84185-125 C* Contact Uwe Poehle C ZIB, Scientific Software Group C phone : + 49/30/84185-241 C fax : + 49/30/84185-107 C e-mail: poehle@zib.de C C --------------------------------------------------------------------- C C* Licence C ------- C C You may use or modify this code for your own non-commercial C purposes for an unlimited time. C In any case you should not deliver this code without a special C permission of ZIB. C In case you intend to use the code commercially, we oblige you C to sign an according licence agreement with ZIB. C C C* Warranty C -------- C C This code has been tested up to a certain level. Defects and C weaknesses, which may be included in the code, do not establish C any warranties by ZIB. ZIB does not take over any liabilities C which may follow from acquisition or application of this code. C C C* Software status C --------------- C C This code is under care of ZIB and belongs to ZIB software C class II. C C C --------------------------------------------------------------------- C C* References: C C /1/ W. B. Gragg: C On Extrapolation Algorithms for Ordinary C Initial Value Problems C SIAM J. Numer. Anal. 2, 384-404 (1965) C C /2/ R. Bulirsch, J. Stoer: C Numerical Treatment of Ordinary Differential Equations C by Extrapolation Methods C Num. Math. 8, 1-13 (1966) C C /3/ P. Deuflhard: C Order and Stepsize Control in Extrapolation Methods C Numer. Math. 41, 399-422 (1983) C C C* External subroutine: (to be supplied by the user) C C FCN EXT Subroutine FCN (N,T,Y,DY) C Right-hand side of first-order C differential equations C N Number of first-order ODEs C T Actual position C Y(N) Values at T C DY(N) Derivatives at T C C C* Parameters: (* marks input/output parameters) C C N I IN Number of first-order ODEs C * T D IN Starting point of integration C OUT Achieved final point of integration C * Y D IN Array of initial values Y(1),...,Y(N) C OUT Array of final values C TEND D IN Prescribed final point of integration C TOL D IN Prescribed relative precision (.GT.0) C HMAX D IN Maximum permitted stepsize C * H D IN Initial stepsize guess C OUT Stepsize proposal for next integration step C (H .EQ. 0. ,if DIFEX1 fails to proceed) C * KFLAG I IN Print parameter C 0 no output C 1 Integration monitor C 2 Intermediate solution points T,Y(I),I=1,N C 3 Integration monitor and solution points C +10 if KFLAG is augmented by 10, a time C monitor is printed additionally (this may C be very expensive in terms of cpu-time) C OUT Error flag C .GE. 0 Successful integration C (KFLAG not altered internally) C -2 More than NSTMAX basic integration steps C per interval have been performed C -3 More than JRMAX stepsize reductions C occurred per basic integration step C -4 Stepsize proposal for next basic C integration step too small C C C* End prologue C ------------ C C C COMMON /STAT/ NFCN, NSTEP, NACCPT, NREJCT, NDEC, NSOL C Internally initialized, for statistical C purposes C NFCN Number of FCN-evaluations C NSTEP Number of integration steps C NACCPT Number of steps accepted C NREJCT Number of steps rejected C NDEC Number of decompositions C NSOL Number of substitutions C C* Type declaration C INTEGER I, J, JK, JL, JM, JMACT, JOPT, JRED, JRMAX, J1, K, KFIN, $ KFLAG, KM, KMACT, KOPT, K1, L, LOUT, M, MAXODE, MDT, M1, N, $ NACCPT, NDEC, NFCN, NJ, NREJCT, NSOL, NSTEP, NSTMAX C DOUBLE PRECISION ALPHA, AWK, D, DMAX1, DMIN1, DT, DUMMY(1), DY, & DYM, $ DZ, BLDFER, EPMACH, EPSAFE, ERR, FC, FCK, FCM, FCO, FMIN, $ FNJ, FN1, FN, H, HALF, HJ, HJ2, HMAX, HMAXU, HN, HREST, $ HRTRN, OMJ, OMJO, ONE, PCT101, PCT90, QUART, RED, RMAX, RO, $ SAFE, SMALL, T, TEND, TN, TOL, TOLH, TOLMIN, U, Y, YK, YM, $ YMAX, YWGT, ZERO C LOGICAL QFIRST, QKONV, QLAST, QPRMON, QPRSOL, QRED C CHARACTER CHGDAT*20, PRODCT*8 C EXTERNAL FCN C C C* Constants problem oriented: (to be supplied by the user) C C MAXODE I K Maximal number of first-order ODEs C PARAMETER ( MAXODE = 1024 ) C C* Constants machine oriented: (to be verified by the user) C C EPMACH D K Relative machine precision C LOUT I K Output is written on logical unit LOUT C SMALL D K Square-root of smallest positive machine number C CRAY (adapted to Cray X-MP) CRAY PARAMETER ( EPMACH = 7.106D-15 , CRAY $ EPSAFE = EPMACH*10.0D0 , CRAY $ LOUT = 6 , CRAY $ SMALL = 6.771D-1234 ) C CIBM (adapted to Siemens 7.865, IBM 370-compatible) CIBM PARAMETER ( EPMACH = 2.22D-16 , CIBM $ EPSAFE = EPMACH*10.0D0 , CIBM $ LOUT = 6 , CIBM $ SMALL = 7.35D-40 ) C CSUN (adapted to sun) C PARAMETER ( EPMACH = 0.1085D-18 , C $ EPSAFE = EPMACH*10.0D0 , C $ LOUT = 6 , C $ SMALL = 0.2223D-161 ) C C* Other Constants: C C HALF D K 1/2 C ONE D K 1 C PCT101 D K 101 Percent C PCT90 D K 90 Percent C QUART D K 1/4 C ZERO D K 0 C PARAMETER ( HALF = 0.5 D0 , $ ONE = 1.0 D0 , $ PCT101 = 1.01 D0 , $ PCT90 = 0.9 D0 , $ QUART = 0.25 D0 , $ ZERO = 0.0 D0 ) C C* Control parameters: (to be supplied by the user) C standard values fixed below C C NSTMAX I K Maximum permitted number of integration steps C per interval = 10000 C JRMAX I K Maximum permitted number of stepsize reductions C KM I K Prescribed maximum column number C JM I K Associated maximum row number C (JM = KM + 1) C MDT I K Associated dimension of DT C BLDFSQ EXT Subroutine BLDFSQ(JM,NJ) C Generate stepsize sequence with respect to /1/ C JM Maximum row number C NJ Array(JM) of stepsize sequence C BLDFSC EXT Subroutine BLDFSC (MODE, Y, N, YOLD, YWGT, C YMAX, THREL, THABS) C Scaling for DIFEX1 C MODE ='INITIAL ' Initial scaling C ='INTERNAL' Scaling during C discretization C ='ACCEPTED' Rescaling if step accepted C Else Error C Y Array of values Y(1),...,Y(N) C N Length of vectors Y, YOLD, YWGT, and YMAX C YOLD Array of old values C YWGT Array of scaled values C YMAX Array of maximum values C THREL Relative threshold value C THABS Absolute threshold value C BLDFER EXT Double Precision function BLDFER(Y, N, YWGT) C Scaled root mean square error C Y Array of values Y(1),...,Y(N) C N Length of vectors Y and YWGT C YWGT Array of scaled values C PARAMETER ( NSTMAX = 10000 , $ JRMAX = 10 , $ KM = 8 , $ JM = KM + 1 , $ MDT = MAXODE*JM ) C C* Internal parameters: (modification not recommended) C C PARAMETER ( FMIN = 1.0 D-3 , $ RMAX = 0.75 D0 , $ RO = QUART , $ SAFE = 0.7 D0 ) C C C* Local variables: (workspace) C C C C QFIRST L V First integration step C QKONV L V Convergence detected C QLAST L V Last integration step C QPRMON L V Print integration monitor C QPRSOL L V Print intermediate solution points C C* Dimensions: C DIMENSION ALPHA(JM,JM), AWK(JM), D(JM,JM), DT(MAXODE,JM), $ DY(MAXODE), DYM(MAXODE), DZ(MAXODE), FCK(KM), NJ(JM), $ Y(MAXODE), YK(MAXODE), YM(MAXODE), YMAX(MAXODE), $ YWGT(MAXODE) C COMMON /DXSTAT/ NFCN, NSTEP, NACCPT, NREJCT, NDEC, NSOL C C******* Revision 1 ******* Latest change: DATA CHGDAT /'August 27, 1991'/ DATA PRODCT /'BLDFX1'/ C*************************** C C C* Modification history C -------------------- C C C 1.0 Feb 9, 1988 First release at ZIB C 1.1 Mar 27, 1991 Vectorize extrapolation loop, C Time monitor C 1.2 Aug 27, 1991 Allow reverse integration direction C C DATA DT/MDT*0.D0/ C C---1. Initial preparations CALL ZIBCONST(EPMACH,SMALL) EPSAFE = EPMACH*10.0D0 LOUT = 6 QPRMON = (KFLAG .EQ. 1 .OR. KFLAG .EQ. 3) QPRSOL = (KFLAG .GE. 2) HRTRN = H C DO 1001 I = 1, N YMAX(I) = ZERO 1001 CONTINUE C HREST = TEND - T H = SIGN (DMIN1 (DABS(H), DABS(HREST)), HREST) HMAX = DABS(HMAX) HMAXU = HMAX FCM = DMAX1 (DABS(H)/HMAX, FMIN) KMACT = KM JMACT = JM CALL BLDFSQ (JM, NJ) FN = DBLE (N) FN1 = DBLE (NJ(1)) TOLH = RO*TOL TOLMIN = EPSAFE*FN IF (TOL .LT. TOLMIN) THEN WRITE (LOUT, 10002) PRODCT, TOL, TOLMIN TOL = TOLMIN ENDIF C C--- Compute amount of work per row of extrapolation tableau AWK(1) = FN1 + ONE DO 101 J=2,JM J1 = J - 1 FNJ = DBLE (NJ(J)) AWK(J) = AWK(J1) + FNJ DO 1011 K=1,J1 D(J,K) = (FNJ / DBLE (NJ(K)))*(FNJ / DBLE (NJ(K))) 1011 CONTINUE C IF (J .NE. 2) THEN DO 1012 K1=2,J1 K = K1 - 1 ALPHA(J1,K) = TOLH**((AWK(K1) - AWK(J)) / $ ((AWK(J) - AWK(1) + ONE)*DBLE(K + K1))) 1012 CONTINUE C ENDIF 101 CONTINUE C C---1.2 Determination of maximum column number in extrapolation C--- tableau (information theoretic concept, ref./3/) KOPT = 1 JOPT = 2 121 CONTINUE C DO WHILE (JOPT .LT. KMACT .AND. C AWK(JOPT+1)*PCT101 .LE. AWK(JOPT)*ALPHA(JOPT,KOPT)) IF (JOPT .GE. KMACT .OR. $ AWK(JOPT+1)*PCT101 .GT. AWK(JOPT)*ALPHA(JOPT,KOPT)) GOTO 122 C Exit 121 KOPT = JOPT JOPT = JOPT + 1 GOTO 121 C ENDDO 122 KMACT = KOPT + 1 JMACT = JOPT IF (QPRMON) WRITE (LOUT, 11221) $ PRODCT, CHGDAT, TOL, KMACT, NJ C IF (QPRSOL) WRITE (LOUT, 11222) NSTEP = 0 QFIRST = .TRUE. QLAST = .FALSE. C NFCN = 0 KFIN = 0 OMJO = ZERO CALL BLDFSC ('INITIAL ', Y, N, DUMMY, YWGT, YMAX, TOL, ONE) C C---2. Basic integration step 2 CONTINUE C DO WHILE (T .NE. TEND) IF (QPRMON) WRITE (LOUT, 12001) NSTEP,NFCN,T,H,KFIN,KOPT IF (QPRSOL) WRITE (LOUT, 12002) NSTEP,NFCN,T,H,(Y(I),I=1,N) JRED = 0 C C--- Explicit euler starting step CALL FCN (N, T, Y, DZ) NFCN = NFCN + 1 C C---3. Basic discretization step 3 CONTINUE C DO WHILE (JRED .LE. JRMAX .AND. .NOT. QKONV) IF (QLAST) THEN TN = TEND ELSE TN = T + H ENDIF IF (TN .EQ. T) THEN C Error 4 IF (QPRMON) WRITE (LOUT, 13001) PRODCT KFLAG = -4 GOTO 9 C Exit to Return ENDIF C C---3.1 Internal discretization DO 31 J=1,JMACT M = NJ(J) M1 = M - 1 KFIN = J - 1 FNJ = DBLE (M) HJ = H / FNJ HJ2 = HJ + HJ DO I=1,N YK(I) = Y(I) YM(I) = Y(I) + HJ*DZ(I) end do C C---3.1.3 Explicit mid-point rule DO 313 K=1,M1 CALL FCN (N, T + HJ*DBLE (K), YM, DY) NFCN = NFCN + 1 DO 3135 I=1,N U = YK(I) + HJ2*DY(I) YK(I) = YM(I) YM(I) = U 3135 CONTINUE 313 CONTINUE C C---3.1.4 Smoothing final step CALL FCN (N, TN, YM, DY) NFCN = NFCN + 1 DO 3141 I = 1,N YM(I) = (YM(I) + YK(I) + HJ*DY(I))*HALF 3141 CONTINUE C C C---3.1.5 Extrapolation DO 3153 I=1,N DY(I) = YM(I) YK(I) = DT(I,1) DT(I,1) = DY(I) 3153 CONTINUE C DO 3158 K=2,J JK = J - K + 1 C DO 3155 I=1,N DYM(I) = (DY(I) - YK(I)) / (D(J,JK) - ONE) DY(I) = D(J,JK)*DYM(I) YK(I) = DT(I,K) DT(I,K) = DYM(I) YM(I) = DYM(I) + YM(I) 3155 CONTINUE C 3158 CONTINUE C IF (J .NE. 1) THEN C C---3.1.6 Convergence monitor CALL BLDFSC ('INTERNAL',YM,N,Y,YWGT,YMAX,TOL,ONE) ERR = BLDFER (DYM, N, YWGT) QKONV = ERR .LE. TOL ERR = ERR / TOLH C C--- Order control K = J - 1 FC = ERR**(ONE / DBLE(K + J)) FCK(K) = FC C C--- Order window IF (J .GE. KOPT .OR. QFIRST .OR. QLAST) THEN IF (QKONV) GOTO 25 C Exit 3 for next basic integration step C C--- Check for possible stepsize reduction RED = ONE / FC QRED = .FALSE. IF (K .EQ. KMACT .OR. K .EQ. JOPT) THEN RED = RED*SAFE QRED = .TRUE. ELSE IF (K .EQ. KOPT) THEN RED = RED*ALPHA(JOPT,KOPT) IF (RED .LT. ONE) THEN RED = ONE / FC QRED = .TRUE. ENDIF ELSE IF (KOPT .EQ. KMACT) THEN RED = RED*ALPHA(KMACT,K) IF (RED .LT. ONE) THEN RED = RED * SAFE QRED = .TRUE. ENDIF ELSE RED = RED*ALPHA(JOPT,K) IF (RED .LT. ONE) THEN RED = ALPHA(KOPT,K) / FC QRED = .TRUE. ENDIF ENDIF ENDIF ENDIF IF (QRED) GOTO 32 C Exit 3.1 to stepsize reduction ENDIF ENDIF 31 CONTINUE C C---3.2 Prepare stepsize reduction 32 CONTINUE C C---3.5 Stepsize reduction RED = DMIN1 (RED, RMAX) H = H*RED IF (NSTEP .GT. 0) QLAST = .FALSE. JRED = JRED + 1 IF (QPRMON) WRITE (LOUT, 13501) JRED,RED, $ KFIN,KOPT,KMACT IF (JRED .GT. JRMAX) THEN C Error 3 IF (QPRMON) WRITE (LOUT, 13502) JRMAX KFLAG = -3 GOTO 9 C Exit to Return ENDIF GOTO 3 C ENDDO C C ************************************************ C---2.5 Preparations for next basic integration step 25 NSTEP = NSTEP + 1 QFIRST = .FALSE. IF (NSTEP .GT. NSTMAX) THEN C Error 2 C Emergency exit, if too many steps taken IF (QPRMON) WRITE (LOUT, 12501) PRODCT, NSTMAX KFLAG = -2 GOTO 9 C Exit to return ENDIF C C--- Restoring DO 251 I=1, N Y(I) = YM(I) 251 CONTINUE C T = TN IF (T .EQ. TEND) GOTO 9 C Exit to return CALL BLDFSC ('ACCEPTED', Y, N, DUMMY, YWGT, YMAX, TOL, ONE) C C---2.7 Order and stepsize selection C C---2.7.1 Stepsize restrictions HMAX = DMIN1(HMAXU, DABS(H)/FMIN) FCM = DABS(H) / HMAX C C---2.7.2 Optimal order determination KOPT = 1 JOPT = 2 FCO = DMAX1 (FCK(1), FCM) OMJO = FCO*AWK(2) IF (KFIN .GE. 2) THEN DO 272 L=2,KFIN JL = L + 1 FC = DMAX1 (FCK(L), FCM) OMJ = FC*AWK(JL) IF (OMJ*PCT101 .LE. OMJO .AND. L .LT. KMACT) THEN KOPT = L JOPT = JL OMJO = OMJ FCO = FC ENDIF 272 CONTINUE ENDIF HREST = TEND - T HN = H / FCO C C---2.7.3 Possible increase of order IF (DABS(HN) .LT. DABS(HREST)) THEN IF ((JRED .EQ. 0 .OR. NSTEP .EQ. 0) .AND. $ KOPT .GE. KFIN .AND. KOPT .NE. KMACT) THEN FC = DMAX1 (FCO/ALPHA(JOPT,KOPT), FCM) JL = JOPT + 1 IF (AWK(JL)*FC*PCT101 .LE. OMJO .AND. $ JOPT .LT. KMACT) THEN FCO = FC HN = H / FCO KOPT = JOPT JOPT = JOPT + 1 ENDIF ENDIF ENDIF C C---2.7.4 Stepsize selection H = HN HRTRN = H IF (DABS(H) .GT. DABS(HREST)*PCT90) THEN H = HREST QLAST = .TRUE. ENDIF GOTO 2 C ENDDO C C---9. Exit 9 HMAX = HMAXU IF (KFLAG .LT. 0) THEN C Fail exit H = ZERO ELSE C Solution exit H = HRTRN IF (QPRMON) WRITE (LOUT, 12001) NSTEP,NFCN,T,H,KFIN,KOPT IF (QPRSOL) WRITE (LOUT, 12002) NSTEP,NFCN,T,H,(Y(I),I=1,N) ENDIF RETURN C C 10001 FORMAT(//,' ',A8,' - Error - ' $ , ' Direction if integration is reverse to convention.') 10002 FORMAT(//,' ',A8,' - Warning -' $ , ' Desired tolerance ', D10.3, ' too small.', /, $ 22X,' tolerance set to ', D10.3, '.') C 11221 FORMAT(1H0,A8,' - ',A20,/, $ 1H0,' Rel.prec. TOL ',D10.3,' max.col. ',I3, /, $ ' sequence ',(1H ,13I4)) 11222 FORMAT(//,5X,'Step',3X,'F-calls',8X,'T',25X,'H',5X,'Y1(T)..',//) 12001 FORMAT(2I9,D20.11,D12.4,I9,I6) 12002 FORMAT(2I9,D20.11,D12.4,4D20.11,/,(50X,4D20.11)) 12501 FORMAT(//,' ',A8,' - Error - ' $ ,' more than NSTMAX=',I3,' integration steps',//) 13001 FORMAT(//,' ',A8,' - Error - ' $ ,' stepsize reduction failed to succeed ',//) 13501 FORMAT(I3,' Stepsize reduction factor ',D10.3, $ ' KFIN',I3,' KOPT',I3,' KMAX',I3) 13502 FORMAT(//,' ',A8,' - Error - ' $ ,' more than JRMAX=',I3,' stepsize reductions per step',/) END SUBROUTINE BLDFSQ(M,NJ) c*********************************************************************72 c cc BLDFSQ() sets the stepsize sequence for DIFEX1. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c INTEGER I, M, NJ DIMENSION NJ(M) NJ(1) = 2 DO I=2,M NJ(I) = NJ(I-1) + 2 end do RETURN END SUBROUTINE BLDFSC (MODE, Y, N, YOLD, YWGT, YMAX, THREL, THABS) c*********************************************************************72 c cc BLDFSC() carries out scaling for DIFEX1. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c C (May be altered for real life applications C by the skillful user) C C C* Parameters: C C MODE C*8 IN ='INITIAL ' Initial scaling C ='INTERNAL' Scaling during discretization C ='ACCEPTED' Rescaling if step accepted C Else Error C Y D IN Array of values Y(1),...,Y(N) C N I IN Length of vectors Y, YOLD, YWGT, and YMAX C YOLD D IN Array of old values C YWGT D OUT Array of scaled values new C YMAX D IN Array of maximum values old C OUT Array of maximum values new C THREL D IN Relative threshold value C THABS D IN Absolute threshold value C C* Local variables: C C YUSER D V User defined array of maximum values C C* Type declaration C INTEGER I, LOUT, MAXODE, N C DOUBLE PRECISION DABS, DMAX1, EPMACH, ONE, THABS, THREL, U, Y, $ YMAX, YOLD, YUSER, YWGT, ZERO, SMALL C CHARACTER MODE*8 C C* Constants: C C EPMACH D K Relative machine precision C LOUT I K Output is written on logical unit LOUT C MAXODE I K Maximal number of first-order ODEs C ONE D K 1.0 C ZERO D K 0.0 C CRAY (adapted to Cray X-MP) CRAY PARAMETER ( EPMACH = 7.106D-15 , C CIBM (adapted to Siemens 7.865, IBM 370-compatible) CIBM PARAMETER ( EPMACH = 2.22D-16 , C CSUN (adapted to sun) PARAMETER ( $ LOUT = 6 , $ MAXODE = 1024 , $ ONE = 1.0 D0 , $ ZERO = 0.0 D0 ) C DIMENSION Y(N), YOLD(N), YWGT(N), YMAX(N), YUSER(MAXODE) SAVE YUSER CALL ZIBCONST(EPMACH,SMALL) IF (MODE .EQ. 'INITIAL ' ) THEN C -------- DO 100 I=1,N YUSER(I) = DABS (YMAX(I)) U = DABS (Y(I)) IF (U .LT. EPMACH) U = ONE YMAX(I) = DMAX1 (U, YUSER(I), THABS) YWGT(I) = YMAX(I) 100 CONTINUE C ELSE IF (MODE .EQ. 'INTERNAL' ) THEN C -------- DO 200 I=1,N YWGT(I) = DMAX1 (YMAX(I)*THREL, DABS(Y(I)), $ DABS(YOLD(I)), YUSER(I), THABS) 200 CONTINUE C ELSE IF (MODE .EQ. 'ACCEPTED' ) THEN C -------- DO I=1,N YMAX(I) = DMAX1 (YMAX(I), DABS(Y(I))) end do ELSE WRITE (LOUT, '(//,A,/)') $ ' D1SCAL - ERROR - Illegal mode' ENDIF RETURN END FUNCTION BLDFER ( Y, N, YWGT ) c*********************************************************************72 c cc BLDFER() evaluates the scaled root mean square error. c c Modified: c c 08 January 2013 c c Reference: c c Josef Stoer, Roland Bulirsch, c Introduction to Numerical Mathematics, c Springer, 2002, c ISBN: 038795452X, c LC: QA297.S8213 c C* Parameters: C C Y D IN Array of values Y(1),...,Y(N) C N I IN Length of Vectors Y and YWGT C YWGT D IN Array of scaled values C C* Type declaration C double precision bldfer INTEGER I, N C DOUBLE PRECISION DBLE, DSQRT, SUM, Y, YWGT, ZERO C C* Constants: C C ZERO D K 0 C PARAMETER ( ZERO = 0.0 D0 ) C DIMENSION Y(N), YWGT(N) C SUM = ZERO DO I=1,N SUM = SUM + (Y(I) / YWGT(I)) * (Y(I) / YWGT(I)) 100 end do BLDFER = DSQRT(SUM / DBLE(N)) RETURN END SUBROUTINE BLDECC (A,NROW,NCOL,MCON,M,N,IRANK,COND,D, 1 PIVOT,KRED,AH,V) c*********************************************************************72 c cc BLDECC() carries out a constrained least squares QR decomposition. c c Modified: c c 08 January 2013 C C* Written by P. Deuflhard C* Purpose Solution of least squares problems, optionally C with equality constraints. C* Method Constrained Least Squares QR-Decomposition C (see references below) C* Category D9b1. - Singular, overdetermined or C underdetermined systems of linear C equations, generalized inverses. C Constrained Least Squares solution C* Keywords Linear Least Square Problems, constrained, C QR-decomposition, pseudo inverse. C* Version 0.9 C* Revision April 1984 C* Latest Change January 1991 C* Library CodeLib C* Code Fortran 77, Double Precision C* Environment Standard Fortran 77 environment on PC's, C workstations and hosts. C* Copyright (c) Konrad-Zuse-Zentrum fuer C Informationstechnik Berlin (ZIB) C Takustrasse 7, D-14195 Berlin-Dahlem C phone : + 49/30/84185-0 C fax : + 49/30/84185-125 C* Contact Lutz Weimann C ZIB, Division Scientific Computing, C Department Numerical Analysis and Modelling C phone : + 49/30/84185-185 C fax : + 49/30/84185-107 C e-mail: weimann@zib.de C C* References: C =========== C C /1/ P.Deuflhard, V.Apostolescu: C An underrelaxed Gauss-Newton method for equality C constrained nonlinear least squares problems. C Lecture Notes Control Inform. Sci. vol. 7, p. C 22-32 (1978) C /2/ P.Deuflhard, W.Sautter: C On rank-deficient pseudoinverses. C J. Lin. Alg. Appl. vol. 29, p. 91-111 (1980) C C* Related Programs: SOLCON (BLSOLC) C C --------------------------------------------------------------- C C* Licence C You may use or modify this code for your own non commercial C purposes for an unlimited time. C In any case you should not deliver this code without a special C permission of ZIB. C In case you intend to use the code commercially, we oblige you C to sign an according licence agreement with ZIB. C C* Warranty C This code has been tested up to a certain level. Defects and C weaknesses, which may be included in the code, do not establish C any warranties by ZIB. ZIB does not take over any liabilities C which may follow from aquisition or application of this code. C C* Software status C This code is under partial care of ZIB and belongs to ZIB software C class 2. C C ------------------------------------------------------------ C C* Summary: C ======== C Constrained QR-decomposition of (M,N)-system with C computation of pseudoinverse in case of rank-defeciency . C First MCON rows belong to equality constraints. C C ------------------------------------------------------------ C C INPUT PARAMETERS (* MARKS INOUT PARAMETERS) C ----------------------------------------------- C C C * A(NROW,NCOL) INPUT MATRIX C A(M,N) CONTAINS ACTUAL INPUT C NROW DECLARED NUMBER OF ROWS OF A AND AH C NCOL DECLARED NUMBER OF COLUMNS OF A AND AH C (*)MCON NUMBER OF EQUALITY CONSTRAINTS (MCON<=N) C INTERNALLY REDUCED IF EQUALITY CONSTRAINTS C ARE LINEARLY DEPENDENT C M TREATED NUMBER OF ROWS OF MATRIX A C N TREATED NUMBER OF COLUMNS OF MATRIX A C (*)IRANK PSEUDO-RANK OF MATRIX A C * COND PERMITTED UPPER BOUND OF DABS(D(1)/D(IRANKC)) C AND OF DABS(D(IRANKC+1))/D(IRANK)) C (SUB-CONDITION NUMBERS OF A) C KRED >=0 HOUSEHOLDER TRIANGULARIZATION C (BUILD UP OF PSEUDO-INVERSE,IF IRANK