share/colnew/fortran/consts.f

   1       SUBROUTINE CONSTS (K, RHO, COEF)
   2 C
   3 C**********************************************************************
   4 C
   5 C   purpose
   6 C            assign (once) values to various array constants.
   7 C
   8 C   arrays assigned during compilation:
   9 C     cnsts1 - weights for extrapolation error estimate
  10 C     cnsts2 - weights for mesh selection
  11 C              (the above weights come from the theoretical form for
  12 C              the collocation error -- see [3])
  13 C
  14 C   arrays assigned during execution:
  15 C     wgterr - the particular values of cnsts1 used for current run
  16 C              (depending on k, m)
  17 C     wgtmsh - gotten from the values of cnsts2 which in turn are
  18 C              the constants in the theoretical expression for the
  19 C              errors. the quantities in wgtmsh are 10x the values
  20 C              in cnsts2 so that the mesh selection algorithm
  21 C              is aiming for errors .1x as large as the user
  22 C              requested tolerances.
  23 C     jtol   - components of differential system to which tolerances
  24 C              refer (viz, if ltol(i) refers to a derivative of u(j),
  25 C              then jtol(i)=j)
  26 C     root   - reciprocals of expected rates of convergence of compo-
  27 C              nents of z(j) for which tolerances are specified
  28 C     rho    - the k collocation points on (0,1)
  29 C     coef   -
  30 C     acol  -  the runge-kutta coefficients values at collocation
  31 C              points
  32 C
  33 C**********************************************************************
  34 C
  35       IMPLICIT REAL*8 (A-H,O-Z)
  36       DIMENSION RHO(7), COEF(K,1), CNSTS1(28), CNSTS2(28), DUMMY(1)
  37 C
  38       COMMON /COLORD/ KDUM, NCOMP, MSTAR, KD, MMAX, M(20)
  39       COMMON /COLBAS/ B(28), ACOL(28,7), ASAVE(28,4)
  40       COMMON /COLEST/ TOL(40), WGTMSH(40), WGTERR(40), TOLIN(40),
  41      1                ROOT(40), JTOL(40), LTOL(40), NTOL
  42 C
  43       DATA CNSTS1 /    .25D0,     .625D-1,  7.2169D-2, 1.8342D-2,
  44      1     1.9065D-2, 5.8190D-2, 5.4658D-3, 5.3370D-3, 1.8890D-2,
  45      2     2.7792D-2, 1.6095D-3, 1.4964D-3, 7.5938D-3, 5.7573D-3,
  46      3     1.8342D-2, 4.673D-3,  4.150D-4,  1.919D-3,  1.468D-3,
  47      4     6.371D-3,  4.610D-3,  1.342D-4,  1.138D-4,  4.889D-4,
  48      5     4.177D-4,  1.374D-3,  1.654D-3,  2.863D-3  /
  49       DATA CNSTS2 /   1.25D-1,   2.604D-3,  8.019D-3,  2.170D-5,
  50      1     7.453D-5,  5.208D-4,  9.689D-8,  3.689D-7,  3.100D-6,
  51      2     2.451D-5,  2.691D-10, 1.120D-9,  1.076D-8,  9.405D-8,
  52      3     1.033D-6,  5.097D-13, 2.290D-12, 2.446D-11, 2.331D-10,
  53      4     2.936D-9,  3.593D-8,  7.001D-16, 3.363D-15, 3.921D-14,
  54      5     4.028D-13, 5.646D-12, 7.531D-11, 1.129D-9  /
  55 C
  56 C...  assign weights for error estimate
  57 C
  58       KOFF = K * ( K + 1 ) / 2
  59       IZ = 1
  60       DO 10 J = 1, NCOMP
  61            MJ = M(J)
  62            DO 10 L = 1, MJ
  63              WGTERR(IZ) = CNSTS1(KOFF - MJ + L)
  64              IZ = IZ + 1
  65    10 CONTINUE
  66 C
  67 C...  assign array values for mesh selection: wgtmsh, jtol, and root
  68 C
  69       JCOMP = 1
  70       MTOT = M(1)
  71       DO 40 I = 1, NTOL
  72            LTOLI = LTOL(I)
  73    20      CONTINUE
  74            IF ( LTOLI .LE. MTOT )                   GO TO 30
  75            JCOMP = JCOMP + 1
  76            MTOT = MTOT + M(JCOMP)
  77            GO TO 20
  78    30      CONTINUE
  79            JTOL(I) = JCOMP
  80            WGTMSH(I) = 1.D1 * CNSTS2(KOFF+LTOLI-MTOT) / TOLIN(I)
  81            ROOT(I) = 1.D0 / DFLOAT(K+MTOT-LTOLI+1)
  82    40 CONTINUE
  83 C
  84 C...  specify collocation points
  85 C
  86       GO TO (50,60,70,80,90,100,110), K
  87    50 RHO(1) = 0.D0
  88       GO TO 120
  89    60 RHO(2) = .57735026918962576451D0
  90       RHO(1) = - RHO(2)
  91       GO TO 120
  92    70 RHO(3) = .77459666924148337704D0
  93       RHO(2) = .0D0
  94       RHO(1) = - RHO(3)
  95       GO TO 120
  96    80 RHO(4) = .86113631159405257523D0
  97       RHO(3) = .33998104358485626480D0
  98       RHO(2) = - RHO(3)
  99       RHO(1) = - RHO(4)
 100       GO TO 120
 101    90 RHO(5) = .90617984593866399280D0
 102       RHO(4) = .53846931010568309104D0
 103       RHO(3) = .0D0
 104       RHO(2) = - RHO(4)
 105       RHO(1) = - RHO(5)
 106       GO TO 120
 107   100 RHO(6) = .93246951420315202781D0
 108       RHO(5) = .66120938646626451366D0
 109       RHO(4) = .23861918608319690863D0
 110       RHO(3) = -RHO(4)
 111       RHO(2) = -RHO(5)
 112       RHO(1) = -RHO(6)
 113       GO TO 120
 114   110 RHO(7) = .949107991234275852452D0
 115       RHO(6) = .74153118559939443986D0
 116       RHO(5) = .40584515137739716690D0
 117       RHO(4) = 0.D0
 118       RHO(3) = -RHO(5)
 119       RHO(2) = -RHO(6)
 120       RHO(1) = -RHO(7)
 121   120 CONTINUE
 122 C
 123 C...  map (-1,1) to (0,1) by  t = .5 * (1. + x)
 124 C
 125       DO 130 J = 1, K
 126          RHO(J) = .5D0 * (1.D0 + RHO(J))
 127   130 CONTINUE
 128 C
 129 C...  now find runge-kutta coeffitients b, acol and asave
 130 C...  the values of asave are to be used in  newmsh  and errchk .
 131 C
 132       DO 140 J = 1, K
 133          DO 135 I = 1, K
 134   135      COEF(I,J) = 0.D0
 135          COEF(J,J) = 1.D0
 136          CALL VMONDE (RHO, COEF(1,J), K)
 137   140 CONTINUE
 138       CALL RKBAS ( 1.D0, COEF, K, MMAX, B, DUMMY, 0)
 139       DO 150 I = 1, K
 140          CALL RKBAS ( RHO(I), COEF, K, MMAX, ACOL(1,I), DUMMY, 0)
 141   150 CONTINUE
 142       CALL RKBAS ( 1.D0/6.D0, COEF, K, MMAX, ASAVE(1,1), DUMMY, 0)
 143       CALL RKBAS ( 1.D0/3.D0, COEF, K, MMAX, ASAVE(1,2), DUMMY, 0)
 144       CALL RKBAS ( 2.D0/3.D0, COEF, K, MMAX, ASAVE(1,3), DUMMY, 0)
 145       CALL RKBAS ( 5.D0/6.D0, COEF, K, MMAX, ASAVE(1,4), DUMMY, 0)
 146       RETURN
 147       END