### Initial revision of COLNEW and additional support routines from

`Linpack and BLAS.`
parent 272dc22e
 SUBROUTINE APPROX (I, X, ZVAL, A, COEF, XI, N, Z, DMZ, K, 1 NCOMP, MMAX, M, MSTAR, MODE, DMVAL, MODM ) C C********************************************************************** C C purpose C (1) (m1-1) (mncomp-1) C evaluate z(u(x))=(u (x),u (x),...,u (x),...,u (x) ) C 1 1 1 mncomp C at one point x. C C variables C a - array of mesh independent rk-basis coefficients C basm - array of mesh dependent monomial coefficients C xi - the current mesh (having n subintervals) C z - the current solution vector C dmz - the array of mj-th derivatives of the current solution C mode - determines the amount of initialization needed C = 4 forms z(u(x)) using z, dmz and ha C = 3 as in =4, but computes local rk-basis C = 2 as in =3, but determines i such that C xi(i) .le. x .lt. xi(i+1) (unless x=xi(n+1)) C = 1 retrieve z=z(u(x(i))) directly C C********************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION ZVAL(1), DMVAL(1), XI(1), M(1), A(7,1), DM(7) DIMENSION Z(1), DMZ(1), BM(4), COEF(1) C COMMON /COLOUT/ PRECIS, IOUT, IPRINT C GO TO (10, 30, 80, 90), MODE C C... mode = 1 , retrieve z( u(x) ) directly for x = xi(i). C 10 X = XI(I) IZ = (I-1) * MSTAR DO 20 J = 1, MSTAR IZ = IZ + 1 ZVAL(J) = Z(IZ) 20 CONTINUE RETURN C C... mode = 2 , locate i so xi(i) .le. x .lt. xi(i+1) C 30 CONTINUE IF ( X .GE. XI(1)-PRECIS .AND. X .LE. XI(N+1)+PRECIS ) 1 GO TO 40 IF (IPRINT .LT. 1) WRITE(IOUT,900) X, XI(1), XI(N+1) IF ( X .LT. XI(1) ) X = XI(1) IF ( X .GT. XI(N+1) ) X = XI(N+1) 40 IF ( I .GT. N .OR. I .LT. 1 ) I = (N+1) / 2 ILEFT = I IF ( X .LT. XI(ILEFT) ) GO TO 60 DO 50 L = ILEFT, N I = L IF ( X .LT. XI(L+1) ) GO TO 80 50 CONTINUE GO TO 80 60 IRIGHT = ILEFT - 1 DO 70 L = 1, IRIGHT I = IRIGHT + 1 - L IF ( X .GE. XI(I) ) GO TO 80 70 CONTINUE C C... mode = 2 or 3 , compute mesh independent rk-basis. C 80 CONTINUE S = (X - XI(I)) / (XI(I+1) - XI(I)) CALL RKBAS ( S, COEF, K, MMAX, A, DM, MODM ) C C... mode = 2, 3, or 4 , compute mesh dependent rk-basis. C 90 CONTINUE BM(1) = X - XI(I) DO 95 L = 2, MMAX BM(L) = BM(1) / DFLOAT(L) 95 CONTINUE C C... evaluate z( u(x) ). C 100 IR = 1 IZ = (I-1) * MSTAR + 1 IDMZ = (I-1) * K * NCOMP DO 140 JCOMP = 1, NCOMP MJ = M(JCOMP) IR = IR + MJ IZ = IZ + MJ DO 130 L = 1, MJ IND = IDMZ + JCOMP ZSUM = 0.D0 DO 110 J = 1, K ZSUM = ZSUM + A(J,L) * DMZ(IND) 110 IND = IND + NCOMP DO 120 LL = 1, L LB = L + 1 - LL 120 ZSUM = ZSUM * BM(LB) + Z(IZ-LL) 130 ZVAL(IR-L) = ZSUM 140 CONTINUE IF ( MODM .EQ. 0 ) RETURN C C... for modm = 1 evaluate dmval(j) = mj-th derivative of uj. C DO 150 JCOMP = 1, NCOMP 150 DMVAL(JCOMP) = 0.D0 IDMZ = IDMZ + 1 DO 170 J = 1, K FACT = DM(J) DO 160 JCOMP = 1, NCOMP DMVAL(JCOMP) = DMVAL(JCOMP) + FACT * DMZ(IDMZ) IDMZ = IDMZ + 1 160 CONTINUE 170 CONTINUE RETURN C-------------------------------------------------------------------- 900 FORMAT(37H ****** DOMAIN ERROR IN APPROX ****** 1 /4H X =,D20.10, 10H ALEFT =,D20.10, 2 11H ARIGHT =,D20.10) END
 C C---------------------------------------------------------------------- C p a r t 4 C polynomial and service routines C---------------------------------------------------------------------- C SUBROUTINE APPSLN (X, Z, FSPACE, ISPACE) C C***************************************************************** C C purpose C C set up a standard call to approx to evaluate the C approximate solution z = z( u(x) ) at a point x C (it has been computed by a call to colnew ). C the parameters needed for approx are retrieved C from the work arrays ispace and fspace . C C***************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Z(1), FSPACE(1), ISPACE(1), A(28), DUMMY(1) IS6 = ISPACE(6) IS5 = ISPACE(1) + 2 IS4 = IS5 + ISPACE(4) * (ISPACE(1) + 1) I = 1 CALL APPROX (I, X, Z, A, FSPACE(IS6), FSPACE(1), ISPACE(1), 1 FSPACE(IS5), FSPACE(IS4), ISPACE(2), ISPACE(3), 2 ISPACE(5), ISPACE(8), ISPACE(4), 2, DUMMY, 0) RETURN END
This diff is collapsed.
 SUBROUTINE CONSTS (K, RHO, COEF) C C********************************************************************** C C purpose C assign (once) values to various array constants. C C arrays assigned during compilation: C cnsts1 - weights for extrapolation error estimate C cnsts2 - weights for mesh selection C (the above weights come from the theoretical form for C the collocation error -- see ) C C arrays assigned during execution: C wgterr - the particular values of cnsts1 used for current run C (depending on k, m) C wgtmsh - gotten from the values of cnsts2 which in turn are C the constants in the theoretical expression for the C errors. the quantities in wgtmsh are 10x the values C in cnsts2 so that the mesh selection algorithm C is aiming for errors .1x as large as the user C requested tolerances. C jtol - components of differential system to which tolerances C refer (viz, if ltol(i) refers to a derivative of u(j), C then jtol(i)=j) C root - reciprocals of expected rates of convergence of compo- C nents of z(j) for which tolerances are specified C rho - the k collocation points on (0,1) C coef - C acol - the runge-kutta coefficients values at collocation C points C C********************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION RHO(7), COEF(K,1), CNSTS1(28), CNSTS2(28), DUMMY(1) C COMMON /COLORD/ KDUM, NCOMP, MSTAR, KD, MMAX, M(20) COMMON /COLBAS/ B(28), ACOL(28,7), ASAVE(28,4) COMMON /COLEST/ TOL(40), WGTMSH(40), WGTERR(40), TOLIN(40), 1 ROOT(40), JTOL(40), LTOL(40), NTOL C DATA CNSTS1 / .25D0, .625D-1, 7.2169D-2, 1.8342D-2, 1 1.9065D-2, 5.8190D-2, 5.4658D-3, 5.3370D-3, 1.8890D-2, 2 2.7792D-2, 1.6095D-3, 1.4964D-3, 7.5938D-3, 5.7573D-3, 3 1.8342D-2, 4.673D-3, 4.150D-4, 1.919D-3, 1.468D-3, 4 6.371D-3, 4.610D-3, 1.342D-4, 1.138D-4, 4.889D-4, 5 4.177D-4, 1.374D-3, 1.654D-3, 2.863D-3 / DATA CNSTS2 / 1.25D-1, 2.604D-3, 8.019D-3, 2.170D-5, 1 7.453D-5, 5.208D-4, 9.689D-8, 3.689D-7, 3.100D-6, 2 2.451D-5, 2.691D-10, 1.120D-9, 1.076D-8, 9.405D-8, 3 1.033D-6, 5.097D-13, 2.290D-12, 2.446D-11, 2.331D-10, 4 2.936D-9, 3.593D-8, 7.001D-16, 3.363D-15, 3.921D-14, 5 4.028D-13, 5.646D-12, 7.531D-11, 1.129D-9 / C C... assign weights for error estimate C KOFF = K * ( K + 1 ) / 2 IZ = 1 DO 10 J = 1, NCOMP MJ = M(J) DO 10 L = 1, MJ WGTERR(IZ) = CNSTS1(KOFF - MJ + L) IZ = IZ + 1 10 CONTINUE C C... assign array values for mesh selection: wgtmsh, jtol, and root C JCOMP = 1 MTOT = M(1) DO 40 I = 1, NTOL LTOLI = LTOL(I) 20 CONTINUE IF ( LTOLI .LE. MTOT ) GO TO 30 JCOMP = JCOMP + 1 MTOT = MTOT + M(JCOMP) GO TO 20 30 CONTINUE JTOL(I) = JCOMP WGTMSH(I) = 1.D1 * CNSTS2(KOFF+LTOLI-MTOT) / TOLIN(I) ROOT(I) = 1.D0 / DFLOAT(K+MTOT-LTOLI+1) 40 CONTINUE C C... specify collocation points C GO TO (50,60,70,80,90,100,110), K 50 RHO(1) = 0.D0 GO TO 120 60 RHO(2) = .57735026918962576451D0 RHO(1) = - RHO(2) GO TO 120 70 RHO(3) = .77459666924148337704D0 RHO(2) = .0D0 RHO(1) = - RHO(3) GO TO 120 80 RHO(4) = .86113631159405257523D0 RHO(3) = .33998104358485626480D0 RHO(2) = - RHO(3) RHO(1) = - RHO(4) GO TO 120 90 RHO(5) = .90617984593866399280D0 RHO(4) = .53846931010568309104D0 RHO(3) = .0D0 RHO(2) = - RHO(4) RHO(1) = - RHO(5) GO TO 120 100 RHO(6) = .93246951420315202781D0 RHO(5) = .66120938646626451366D0 RHO(4) = .23861918608319690863D0 RHO(3) = -RHO(4) RHO(2) = -RHO(5) RHO(1) = -RHO(6) GO TO 120 110 RHO(7) = .949107991234275852452D0 RHO(6) = .74153118559939443986D0 RHO(5) = .40584515137739716690D0 RHO(4) = 0.D0 RHO(3) = -RHO(5) RHO(2) = -RHO(6) RHO(1) = -RHO(7) 120 CONTINUE C C... map (-1,1) to (0,1) by t = .5 * (1. + x) C DO 130 J = 1, K RHO(J) = .5D0 * (1.D0 + RHO(J)) 130 CONTINUE C C... now find runge-kutta coeffitients b, acol and asave C... the values of asave are to be used in newmsh and errchk . C DO 140 J = 1, K DO 135 I = 1, K 135 COEF(I,J) = 0.D0 COEF(J,J) = 1.D0 CALL VMONDE (RHO, COEF(1,J), K) 140 CONTINUE CALL RKBAS ( 1.D0, COEF, K, MMAX, B, DUMMY, 0) DO 150 I = 1, K CALL RKBAS ( RHO(I), COEF, K, MMAX, ACOL(1,I), DUMMY, 0) 150 CONTINUE CALL RKBAS ( 1.D0/6.D0, COEF, K, MMAX, ASAVE(1,1), DUMMY, 0) CALL RKBAS ( 1.D0/3.D0, COEF, K, MMAX, ASAVE(1,2), DUMMY, 0) CALL RKBAS ( 2.D0/3.D0, COEF, K, MMAX, ASAVE(1,3), DUMMY, 0) CALL RKBAS ( 5.D0/6.D0, COEF, K, MMAX, ASAVE(1,4), DUMMY, 0) RETURN END
 SUBROUTINE CONTRL (XI, XIOLD, Z, DMZ, RHS, DELZ, DELDMZ, 1 DQZ, DQDMZ, G, W, V, VALSTR, SLOPE, SCALE, DSCALE, 2 ACCUM, IPVTG, INTEGS, IPVTW, NFXPNT, FIXPNT, IFLAG, 3 FSUB, DFSUB, GSUB, DGSUB, GUESS ) C C********************************************************************** C C purpose C this subroutine is the actual driver. the nonlinear iteration C strategy is controlled here ( see  ). upon convergence, errchk C is called to test for satisfaction of the requested tolerances. C C variables C C check - maximum tolerance value, used as part of criteria for C checking for nonlinear iteration convergence C relax - the relaxation factor for damped newton iteration C relmin - minimum allowable value for relax (otherwise the C jacobian is considered singular). C rlxold - previous relax C rstart - initial value for relax when problem is sensitive C ifrz - number of fixed jacobian iterations C lmtfrz - maximum value for ifrz before performing a reinversion C iter - number of iterations (counted only when jacobian C reinversions are performed). C xi - current mesh C xiold - previous mesh C ipred = 0 if relax is determined by a correction C = 1 if relax is determined by a prediction C ifreez = 0 if the jacobian is to be updated C = 1 if the jacobian is currently fixed (frozen) C iconv = 0 if no previous convergence has been obtained C = 1 if convergence on a previous mesh has been obtained C icare =-1 no convergence occurred (used for regular problems) C = 0 a regular problem C = 1 a sensitive problem C = 2 used for continuation (see description of ipar(10) C in colnew). C rnorm - norm of rhs (right hand side) for current iteration C rnold - norm of rhs for previous iteration C anscl - scaled norm of newton correction C anfix - scaled norm of newton correction at next step C anorm - scaled norm of a correction obtained with jacobian fixed C nz - number of components of z (see subroutine approx) C ndmz - number of components of dmz (see subroutine approx) C imesh - a control variable for subroutines newmsh and errchk C = 1 the current mesh resulted from mesh selection C or is the initial mesh. C = 2 the current mesh resulted from doubling the C previous mesh C C********************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION XI(1), XIOLD(1), Z(1), DMZ(1), RHS(1) DIMENSION G(1), W(1), V(1), VALSTR(1), SLOPE(1), ACCUM(1) DIMENSION DELZ(1), DELDMZ(1), DQZ(1), DQDMZ(1) , FIXPNT(1) DIMENSION DUMMY(1), SCALE(1), DSCALE(1) DIMENSION INTEGS(1), IPVTG(1), IPVTW(1) C COMMON /COLOUT/ PRECIS, IOUT, IPRINT COMMON /COLORD/ K, NCOMP, MSTAR, KD, MMAX, M(20) COMMON /COLAPR/ N, NOLD, NMAX, NZ, NDMZ COMMON /COLMSH/ MSHFLG, MSHNUM, MSHLMT, MSHALT COMMON /COLSID/ ZETA(40), ALEFT, ARIGHT, IZETA, IDUM COMMON /COLNLN/ NONLIN, ITER, LIMIT, ICARE, IGUESS COMMON /COLEST/ TOL(40), WGTMSH(40), WGTERR(40), TOLIN(40), 1 ROOT(40), JTOL(40), LTOL(40), NTOL C EXTERNAL FSUB, DFSUB, GSUB, DGSUB, GUESS C C... constants for control of nonlinear iteration C RELMIN = 1.D-3 RSTART = 1.D-2 LMTFRZ = 4 C C... compute the maximum tolerance C CHECK = 0.D0 DO 10 I = 1, NTOL 10 CHECK = DMAX1 ( TOLIN(I), CHECK ) IMESH = 1 ICONV = 0 IF ( NONLIN .EQ. 0 ) ICONV = 1 ICOR = 0 NOCONV = 0 MSING = 0 C C... the main iteration begins here . C... loop 20 is executed until error tolerances are satisfied or C... the code fails (due to a singular matrix or storage limitations) C 20 CONTINUE C C... initialization for a new mesh C ITER = 0 IF ( NONLIN .GT. 0 ) GO TO 50 C C... the linear case. C... set up and solve equations C CALL LSYSLV (MSING, XI, XIOLD, DUMMY, DUMMY, Z, DMZ, G, 1 W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 0, 2 FSUB, DFSUB, GSUB, DGSUB, GUESS ) C C... check for a singular matrix C IF ( MSING .EQ. 0 ) GO TO 400 30 IF ( MSING .LT. 0 ) GO TO 40 IF ( IPRINT .LT. 1 ) WRITE (IOUT,495) GO TO 460 40 IF ( IPRINT .LT. 1 ) WRITE (IOUT,490) IFLAG = 0 RETURN C C... iteration loop for nonlinear case C... define the initial relaxation parameter (= relax) C 50 RELAX = 1.D0 C C... check for previous convergence and problem sensitivity C IF ( ICARE .EQ. 1 .OR. ICARE .EQ. (-1) ) RELAX = RSTART IF ( ICONV .EQ. 0 ) GO TO 160 C C... convergence on a previous mesh has been obtained. thus C... we have a very good initial approximation for the newton C... process. proceed with one full newton and then iterate C... with a fixed jacobian. C IFREEZ = 0 C C... evaluate right hand side and its norm and C... find the first newton correction C CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G, 1 W, V, RHS, DQDMZ, INTEGS, IPVTG, IPVTW, RNOLD, 1, 2 FSUB, DFSUB, GSUB, DGSUB, GUESS ) C IF ( IPRINT .LT. 0 ) WRITE(IOUT,530) IF ( IPRINT .LT. 0 ) WRITE (IOUT,510) ITER, RNOLD GO TO 70 C C... solve for the next iterate . C... the value of ifreez determines whether this is a full C... newton step (=0) or a fixed jacobian iteration (=1). C 60 IF ( IPRINT .LT. 0 ) WRITE (IOUT,510) ITER, RNORM RNOLD = RNORM CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G, 1 W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 2 3+IFREEZ, FSUB, DFSUB, GSUB, DGSUB, GUESS ) C C... check for a singular matrix C 70 IF ( MSING .NE. 0 ) GO TO 30 IF ( IFREEZ .EQ. 1 ) GO TO 80 C C... a full newton step C ITER = ITER + 1 IFRZ = 0 80 CONTINUE C C... update z and dmz , compute new rhs and its norm C DO 90 I = 1, NZ Z(I) = Z(I) + DELZ(I) 90 CONTINUE DO 100 I = 1, NDMZ DMZ(I) = DMZ(I) + DELDMZ(I) 100 CONTINUE CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G, 1 W, V, RHS, DUMMY, INTEGS, IPVTG, IPVTW, RNORM, 2, 2 FSUB, DFSUB, GSUB, DGSUB, GUESS ) C C... check monotonicity. if the norm of rhs gets smaller, C... proceed with a fixed jacobian; else proceed cautiously, C... as if convergence has not been obtained before (iconv=0). C IF ( RNORM .LT. PRECIS ) GO TO 390 IF ( RNORM .GT. RNOLD ) GO TO 130 IF ( IFREEZ .EQ. 1 ) GO TO 110 IFREEZ = 1 GO TO 60 C C... verify that the linear convergence with fixed jacobian C... is fast enough. C 110 IFRZ = IFRZ + 1 IF ( IFRZ .GE. LMTFRZ ) IFREEZ = 0 IF ( RNOLD .LT. 4.D0*RNORM ) IFREEZ = 0 C C... check convergence (iconv = 1). C DO 120 IT = 1, NTOL INZ = LTOL(IT) DO 120 IZ = INZ, NZ, MSTAR IF ( DABS(DELZ(IZ)) .GT. 1 TOLIN(IT) * (DABS(Z(IZ)) + 1.D0)) GO TO 60 120 CONTINUE C C... convergence obtained C IF ( IPRINT .LT. 1 ) WRITE (IOUT,560) ITER GO TO 400 C C... convergence of fixed jacobian iteration failed. C 130 IF ( IPRINT .LT. 0 ) WRITE (IOUT,510) ITER, RNORM IF ( IPRINT .LT. 0 ) WRITE (IOUT,540) ICONV = 0 RELAX = RSTART DO 140 I = 1, NZ Z(I) = Z(I) - DELZ(I) 140 CONTINUE DO 150 I = 1, NDMZ DMZ(I) = DMZ(I) - DELDMZ(I) 150 CONTINUE C C... update old mesh C NP1 = N + 1 DO 155 I = 1, NP1 155 XIOLD(I) = XI(I) NOLD = N C ITER = 0 C C... no previous convergence has been obtained. proceed C... with the damped newton method. C... evaluate rhs and find the first newton correction. C 160 IF(IPRINT .LT. 0) WRITE (IOUT,500) CALL LSYSLV (MSING, XI, XIOLD, Z, DMZ, DELZ, DELDMZ, G, 1 W, V, RHS, DQDMZ, INTEGS, IPVTG, IPVTW, RNOLD, 1, 2 FSUB, DFSUB, GSUB, DGSUB, GUESS ) C C... check for a singular matrix C