Line data    Source code

       1           92 :       subroutine aramax(x,a,ia,nn,xext,aext,d2aext)
       2              : c
       3              : c  subroutine armax finds the maximum aext of the function 
       4              : c  a = abs(a(x)),
       5              : c  and the point xext where the maximum occurs. x(n) contains x and
       6              : c  a(n) contains a, for n = 1,...,nn. ia is the first dimension of a.
       7              : c  the extremum is found by fitting a parabola through
       8              : c  the functional values at the three points (x(n),a(n)) closest to
       9              : c  the extremum, if this is not at n = 1 or n = nn.
      10              : c  if the extremum is in the interior of the interval considered
      11              : c  d2aext is set to the second derivative of a wrt x at the
      12              : c  extremum. otherwise d2aext is set to 0.
      13              : c
      14              : c  note that x, a, xext, aext and d2aext are real.
      15              : c
      16              : c  this is a version made to overcome problems with entry in
      17              : c  s/r armax.
      18              : c
      19              : c  calls leq.
      20              : c  **********
      21              : c
      22              : c  ......................................................................
      23              : c
      24              : c  Double precision version.
      25              : c  ++++++++++++++++++++++++
      26              : c
      27              : c  Dated 10/3/90.
      28              : c
      29              : c  Note: this double precision version of the routine has same name
      30              : c  as single precision version, but is distinguished by its file name
      31              : c
      32              :       implicit double precision(a-h,o-z)
      33              :       logical max,abm
      34              :       dimension x(nn),a(ia,nn),wrk(3,3),b(3)
      35              : c
      36              : c  common defining standard input and output
      37              : c
      38              :       common/cstdio/ istdin, istdou, istdpr, istder
      39              : c
      40           92 :       max=.true.
      41           92 :       abm=.true.
      42              : c
      43              : c  find extremal value of a
      44           92 :    12 axt=a(1,1)
      45           92 :       if(abm) axt=abs(axt)
      46           92 :       nxt=1
      47       135056 :       do 20 n=2,nn
      48       134964 :       aa=a(1,n)
      49       134964 :       if(abm) aa=abs(aa)
      50       134964 :       if(max) go to 15
      51            0 :       if(aa-axt) 17,17,20
      52       134964 :    15 if(axt-aa) 17,17,20
      53        41896 :    17 nxt=n
      54       134964 :       axt=aa
      55           92 :    20 continue
      56           92 :       if(nxt.ne.1.and.nxt.ne.nn) go to 30
      57           46 :       xext=x(nxt)
      58           46 :       aext=axt
      59           46 :       d2aext=0
      60           46 :       return
      61              : c  fit parabola around x(nxt)
      62           46 :    30 xtr=x(nxt)
      63           46 :       scl=x(nxt+1)-x(nxt-1)
      64           46 :       if(scl) 35,90,35
      65          184 :    35 do 40 i=1,3
      66          138 :       t=1
      67          138 :       wrk(i,1)=t
      68          138 :       i1=nxt-2+i
      69          138 :       xx=(x(i1)-xtr)/scl
      70          414 :       do 38 k=2,3
      71          276 :       t=t*xx
      72          414 :    38 wrk(i,k)=t
      73          184 :    40 b(i)=a(1,i1)
      74              : c
      75           46 :       call leq(wrk,b,3,1,3,1,det)
      76              : c
      77              : c  find extremum
      78           46 :       if(det) 42,90,42
      79           46 :    42 if(b(3)) 45,95,45
      80           46 :    45 xext=-b(2)/(2*b(3))
      81           46 :       aext=b(1)+xext*(b(2)+xext*b(3))
      82           46 :       if(abm) aext=abs(aext)
      83           46 :       xext=xtr+scl*xext
      84           46 :       d2aext=2*b(3)/(scl*scl)
      85           46 :       return
      86            0 :    90 if(istdpr.gt.0) write(istdpr,100) x(nxt-1),x(nxt),x(nxt+1)
      87            0 :       d2aext=0
      88            0 :       return
      89            0 :    95 if(istdpr.gt.0) write(istdpr,110) b
      90            0 :       d2aext=0
      91            0 :       return
      92              :   100 format(///1x,10('*'),' points coincide in armax'/
      93              :      *  ' x(nxt-1),x(nxt),x(nxt+1):',1p3e15.5)
      94              :   110 format(///1x,10('*'),' coefficient to x**2 is zero in armax'/
      95              :      *  ' b =',1p3e15.5)
      96              :       end