source/unres/src_CSA/rmdd.f

   1 c     algorithm 611, collected algorithms from acm.
   2 c     algorithm appeared in acm-trans. math. software, vol.9, no. 4,
   3 c     dec., 1983, p. 503-524.
   4       integer function imdcon(k)
   5 c
   6       integer k
   7 c
   8 c  ***  return integer machine-dependent constants  ***
   9 c
  10 c     ***  k = 1 means return standard output unit number.   ***
  11 c     ***  k = 2 means return alternate output unit number.  ***
  12 c     ***  k = 3 means return  input unit number.            ***
  13 c          (note -- k = 2, 3 are used only by test programs.)
  14 c
  15 c  +++  port version follows...
  16 c     external i1mach
  17 c     integer i1mach
  18 c     integer mdperm(3)
  19 c     data mdperm(1)/2/, mdperm(2)/4/, mdperm(3)/1/
  20 c     imdcon = i1mach(mdperm(k))
  21 c  +++  end of port version  +++
  22 c
  23 c  +++  non-port version follows...
  24       integer mdcon(3)
  25       data mdcon(1)/6/, mdcon(2)/8/, mdcon(3)/5/
  26       imdcon = mdcon(k)
  27 c  +++  end of non-port version  +++
  28 c
  29  999  return
  30 c  ***  last card of imdcon follows  ***
  31       end
  32       double precision function rmdcon(k)
  33 c
  34 c  ***  return machine dependent constants used by nl2sol  ***
  35 c
  36 c +++  comments below contain data statements for various machines.  +++
  37 c +++  to convert to another machine, place a c in column 1 of the   +++
  38 c +++  data statement line(s) that correspond to the current machine +++
  39 c +++  and remove the c from column 1 of the data statement line(s)  +++
  40 c +++  that correspond to the new machine.                           +++
  41 c
  42       integer k
  43 c
  44 c  ***  the constant returned depends on k...
  45 c
  46 c  ***        k = 1... smallest pos. eta such that -eta exists.
  47 c  ***        k = 2... square root of eta.
  48 c  ***        k = 3... unit roundoff = smallest pos. no. machep such
  49 c  ***                 that 1 + machep .gt. 1 .and. 1 - machep .lt. 1.
  50 c  ***        k = 4... square root of machep.
  51 c  ***        k = 5... square root of big (see k = 6).
  52 c  ***        k = 6... largest machine no. big such that -big exists.
  53 c
  54       double precision big, eta, machep
  55       integer bigi(4), etai(4), machei(4)
  56 c/+
  57       double precision dsqrt
  58 c/
  59       equivalence (big,bigi(1)), (eta,etai(1)), (machep,machei(1))
  60 c
  61 c  +++  ibm 360, ibm 370, or xerox  +++
  62 c
  63 c     data big/z7fffffffffffffff/, eta/z0010000000000000/,
  64 c    1     machep/z3410000000000000/
  65 c
  66 c  +++  data general  +++
  67 c
  68 c     data big/0.7237005577d+76/, eta/0.5397605347d-78/,
  69 c    1     machep/2.22044605d-16/
  70 c
  71 c  +++  dec 11  +++
  72 c
  73 c     data big/1.7d+38/, eta/2.938735878d-39/, machep/2.775557562d-17/
  74 c
  75 c  +++  hp3000  +++
  76 c
  77 c     data big/1.157920892d+77/, eta/8.636168556d-78/,
  78 c    1     machep/5.551115124d-17/
  79 c
  80 c  +++  honeywell  +++
  81 c
  82 c     data big/1.69d+38/, eta/5.9d-39/, machep/2.1680435d-19/
  83 c
  84 c  +++  dec10  +++
  85 c
  86 c     data big/"377777100000000000000000/,
  87 c    1     eta/"002400400000000000000000/,
  88 c    2     machep/"104400000000000000000000/
  89 c
  90 c  +++  burroughs  +++
  91 c
  92 c     data big/o0777777777777777,o7777777777777777/,
  93 c    1     eta/o1771000000000000,o7770000000000000/,
  94 c    2     machep/o1451000000000000,o0000000000000000/
  95 c
  96 c  +++  control data  +++
  97 c
  98 c     data big/37767777777777777777b,37167777777777777777b/,
  99 c    1     eta/00014000000000000000b,00000000000000000000b/,
 100 c    2     machep/15614000000000000000b,15010000000000000000b/
 101 c
 102 c  +++  prime  +++
 103 c
 104 c     data big/1.0d+9786/, eta/1.0d-9860/, machep/1.4210855d-14/
 105 c
 106 c  +++  univac  +++
 107 c
 108 c     data big/8.988d+307/, eta/1.2d-308/, machep/1.734723476d-18/
 109 c
 110 c  +++  vax  +++
 111 c
 112       data big/1.7d+38/, eta/2.939d-39/, machep/1.3877788d-17/
 113 c
 114 c  +++  cray 1  +++
 115 c
 116 c     data bigi(1)/577767777777777777777b/,
 117 c    1     bigi(2)/000007777777777777776b/,
 118 c    2     etai(1)/200004000000000000000b/,
 119 c    3     etai(2)/000000000000000000000b/,
 120 c    4     machei(1)/377224000000000000000b/,
 121 c    5     machei(2)/000000000000000000000b/
 122 c
 123 c  +++  port library -- requires more than just a data statement... +++
 124 c
 125 c     external d1mach
 126 c     double precision d1mach, zero
 127 c     data big/0.d+0/, eta/0.d+0/, machep/0.d+0/, zero/0.d+0/
 128 c     if (big .gt. zero) go to 1
 129 c        big = d1mach(2)
 130 c        eta = d1mach(1)
 131 c        machep = d1mach(4)
 132 c1    continue
 133 c
 134 c  +++ end of port +++
 135 c
 136 c-------------------------------  body  --------------------------------
 137 c
 138       go to (10, 20, 30, 40, 50, 60), k
 139 c
 140  10   rmdcon = eta
 141       go to 999
 142 c
 143  20   rmdcon = dsqrt(256.d+0*eta)/16.d+0
 144       go to 999
 145 c
 146  30   rmdcon = machep
 147       go to 999
 148 c
 149  40   rmdcon = dsqrt(machep)
 150       go to 999
 151 c
 152  50   rmdcon = dsqrt(big/256.d+0)*16.d+0
 153       go to 999
 154 c
 155  60   rmdcon = big
 156 c
 157  999  return
 158 c  ***  last card of rmdcon follows  ***
 159       end