← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/optim/fcube.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/optim/fcube.f
 
1
      subroutine fcube(t,f,fp,ta,fa,fpa,tlower,tupper)
 
2
c     Copyright INRIA
 
3
      implicit double precision (a-h,o-z)
 
4
c
 
5
c           Using f and fp at t and ta, computes new t by cubic formula
 
6
c           safeguarded inside [tlower,tupper].
 
7
c
 
8
      z1=fp+fpa-3.d0*(fa-f)/(ta-t)
 
9
      b=z1+fp
 
10
c
 
11
c              first compute the discriminant (without overflow)
 
12
c
 
13
      if (dabs(z1).le.1.d0) then
 
14
          discri=z1*z1-fp*fpa
 
15
        else
 
16
          discri=fp/z1
 
17
          discri=discri*fpa
 
18
          discri=z1-discri
 
19
          if (z1.ge.0.d0 .and. discri.ge.0.d0) then
 
20
              discri=dsqrt(z1)*dsqrt(discri)
 
21
              go to 200
 
22
          endif
 
23
          if (z1.le.0.d0 .and. discri.le.0.d0) then
 
24
              discri=dsqrt(-z1)*dsqrt(-discri)
 
25
              go to 200
 
26
          endif
 
27
          discri=-1.d0
 
28
      endif
 
29
      if (discri.lt.0.d0) then
 
30
          if (fp.lt.0.d0) t=tupper
 
31
          if (fp.ge.0.d0) t=tlower
 
32
          go to 990
 
33
      endif
 
34
c
 
35
c       discriminant nonnegative, stable solution formula
 
36
c
 
37
      discri=dsqrt(discri)
 
38
  200 if (t-ta.lt.0.d0) discri=-discri
 
39
      sign=(t-ta)/dabs(t-ta)
 
40
      if (b*sign.gt.0.) then
 
41
          anum=(ta-t)*fp
 
42
          den=b+discri
 
43
        else
 
44
          den=z1+b+fpa
 
45
          anum=(ta-t)*(b-discri)
 
46
      endif
 
47
c
 
48
c               now compute the ratio (without overflow)
 
49
c
 
50
      if (dabs(den).ge.1.d0) then
 
51
          t=t+anum/den
 
52
        else
 
53
          if (dabs(anum).lt.(tupper-tlower)*dabs(den)) then
 
54
              t=t+anum/den
 
55
            else
 
56
              if (fp.lt.0.d0) t=tupper
 
57
              if (fp.ge.0.d0) t=tlower
 
58
          endif
 
59
      endif
 
60
c
 
61
c                       finally, safeguard
 
62
c
 
63
      t=dmax1(t,tlower)
 
64
      t=dmin1(t,tupper)
 
65
  990 return
 
66
      end