NAG.f

Go to the documentation of this file.

      REAL*8 FUNCTION s21bbf(X, Y, Z, IFAIL)

C     MARK 8 RELEASE. NAG COPYRIGHT 1979.

C

C     ******************************************************************

C

C      CALCULATES AN APPROXIMATE VALUE FOR THE ELLIPTIC INTEGRAL OF

C      THE FIRST KIND, RF(X,Y,Z), AS DEFINED BY B.C.CARLSON.

C      -SPECIAL FUNCTIONS OF APPLIED MATHS-ACADEMIC PRESS(1977).

C      THE ALGORITHM IS ALSO DUE TO CARLSON BUT HAS BEEN RECODED BY

C      J.L.SCHONFELDER TO AVOID INTERMEDIATE UNDER AND OVERFLOWS

C

C      RF(X,Y,Z)=INTEGRAL(ZERO TO INFINITY) OF

C          0.5/SQRT((T+X)*(T+Y)*(T+Z)) DT

C

C      FOR X,Y,Z ALL .GE. ZERO AND AT MOST ONE EQUAL TO ZERO

C

C     **************************************************************** *

C

C     ******************************************************************

C     PRECISION-DEPENDENT CONSTANT

C

C        ACC  CONTROLS FINAL TRUNCATION ACCURACY. AN ACCURACY

C                 CLOSE TO FULL MACHINE PRECISION CAN BE OBTAINED

C                 FOR ALL LEGAL ARGUMENTS IF ACC IS CHOSEN SO THAT

C           0.25*ACC**6/(1.0-ACC).LE.X02AAF/X02BAF

C

C     TO SELECT THE CORRECT VALUE FOR A PARTICULAR MACHINE-RANGE,

C     ACTIVATE THE STATEMENT CONTAINED IN A COMMENT BEGINNING  CDD ,

C     WHERE  DD  IS THE APPROXIMATE NUMBER OF SIGNIFICANT DECIMAL

C     DIGITS REPRESENTED BY THE MACHINE.

C     ******************************************************************

C

C     .. SCALAR ARGUMENTS ..

*      REAL X, Y, Z

      REAL*8 x, y, z

      INTEGER ifail

C     ..

C     .. LOCAL SCALARS ..

*+SELF,IF=IBM.

      CHARACTER*8 srname

*+SELF.

*      REAL ACC, C1, C2, C3, CXN, CYN, CZN, E2, E3, LAMDA, LOLIM,

      REAL*8 acc, c1, c2, c3, cxn, cyn, czn, e2, e3, lamda, lolim,

     * mu, rscale, rtx, rty, rtz, uplim, xn, yn, zn

      INTEGER ind

C     .. FUNCTION REFERENCES ..

*      REAL SQRT, X02ABE, X02ACE

      REAL*8 sqrt, x02abf, x02acf

      INTEGER p01aae

*

      REAL*8 xarg

      sqrt(xarg)=dsqrt(xarg)

*

C     ..

*      DATA SRNAME /8H S21BBE /

      DATA srname /8h s21bbf /

*+SELF,IF=E09.

*      DATA ACC/5.7E-2/

*+SELF,IF=E10.

*      DATA ACC/3.9E-2/

*+SELF,IF=E13.

*      DATA ACC/1.2E-2/

*+SELF,IF=E16.

      DATA acc/3.9e-3/

*+SELF,IF=E18.

*      DATA ACC/1.8E-3/

*+SELF,IF=E20.

*      DATA ACC/8.5E-4/

*+SELF.

C

C     ORDER X,Y,Z INTO XN,YN,ZN  ST. XN.LE.YN.LE.ZN

      IF (x.LE.y) go to 20

      xn = y

      yn = x

      go to 40

   20 xn = x

      yn = y

   40 zn = z

      IF (yn.LE.zn) go to 60

      zn = yn

      yn = z

      IF (xn.LE.yn) go to 60

      yn = xn

      xn = z

C

C     TEST FOR VALID ARGUMENTS

   60 ind = 1

      IF (xn.LT.0.0) go to 180

      ind = 2

      IF (yn.LE.0.0) go to 180

C

C     VALID CALL

      ifail = 0

      rscale = 1.0

*      LOLIM = 16.0*X02ABE(0.0)

      lolim = 16.0*x02abf(0.0d0)

*      UPLIM = 0.0625*X02ACE(0.0)

      uplim = 0.0625*x02acf(0.0d0)

C

C     FOR EXTREME ARGUMENTS SCALE TO AVOID UNDER AND OVERFLOWS

      IF (zn.LE.uplim) go to 120

      rscale = 0.25

      zn = zn*0.0625

      IF (yn.LE.lolim) go to 80

      yn = yn*0.0625

      IF (xn.LE.lolim) go to 100

      xn = xn*0.0625

      go to 140

   80 lamda = (sqrt(xn)+sqrt(yn))*(sqrt(zn)*0.25)

      xn = lamda*0.25

      yn = xn

      zn = (zn+lamda)*0.25

      go to 140

  100 rtz = sqrt(zn)

      rty = sqrt(yn)

      lamda = rtz*rty + ((rtz+rty)*0.25)*sqrt(xn)

      xn = lamda*0.25

      yn = (yn+lamda)*0.25

      zn = (zn+lamda)*0.25

      go to 140

  120 IF (zn.GT.lolim) go to 140

      rscale = 4.0

      xn = xn*16.0

      yn = yn*16.0

      zn = zn*16.0

C

C     MAIN RECURSION

  140 mu = (xn+yn+zn)/3.0

      czn = 2.0 - (zn+mu)/mu

      cxn = 2.0 - (xn+mu)/mu

      IF (dmax1(cxn,-czn).LE.acc) go to 160

      rtx = sqrt(xn)

      rty = sqrt(yn)

      rtz = sqrt(zn)

      lamda = rtz*(rtx+rty) + rtx*rty

      xn = (xn+lamda)*0.25

      yn = (yn+lamda)*0.25

      zn = (zn+lamda)*0.25

      go to 140

C

C     POWER SERIES EXPANSION

  160 c1 = 1.0/24.0

      c2 = 3.0/44.0

      c3 = 1.0/14.0

      cyn = -cxn - czn

      e2 = cxn*cyn - czn*czn

      e3 = cxn*czn*cyn

*      S21BBE = RSCALE*(1.0+(C1*E2-0.1-C2*E3)*E2+C3*E3)/SQRT(MU)

      s21bbf = rscale*(1.0+(c1*e2-0.1-c2*e3)*e2+c3*e3)/sqrt(mu)

      go to 200

C

C     FAILURE EXITS

  180 ifail = p01aae(ifail,ind,srname)

*      S21BBE = 0.0

      s21bbf = 0.0

  200 RETURN

      END



      REAL*8 FUNCTION s21bcf(X, Y, Z, IFAIL)

C     MARK 8 RELEASE. NAG COPYRIGHT 1979.

C

C     ******************************************************************

C

C      CALCULATES AN APPROXIMATE VALUE FOR THE ELLIPTIC INTEGRAL OF

C      THE SECOND KIND, RD(X,Y,Z), AS DEFINED BY B.C.CARLSON.

C      -SPECIAL FUNCTIONS OF APPLIED MATHS-ACADEMIC PRESS(1977).

C      THE ALGORITHM IS ALSO DUE TO CARLSON.

C

C      RD(X,Y,Z)=INTEGRAL(ZERO TO INFINITY) OF

C          1.5/SQRT((T+X)*(T+Y)*(T+Z)**3) DT

C

C      FOR X.GE.0.0 , Y.GE.0.0, Z.GT.0.0 AND AT MOST ONE OF X

C      AND Y EQUAL TO 0.0.

C

C     ******************************************************************

C

C     ******************************************************************

C     PRECISION-DEPENDENT CONSTANT

C

C        ACC  CONTROLS FINAL TRUNCATION ACCURACY. AN ACCURACY

C                 CLOSE TO FULL MACHINE PRECISION CAN BE OBTAINED

C                 FOR ALL LEGAL ARGUMENTS IF ACC IS CHOSEN SO THAT

C           3.0*ACC**6/(1.0-ACC)**1.5.LE.X02AAF/X02BAF

C

C     TO SELECT THE CORRECT VALUE FOR A PARTICULAR MACHINE-RANGE,

C     ACTIVATE THE STATEMENT CONTAINED IN A COMMENT BEGINNING  CDD ,

C     WHERE  DD  IS THE APPROXIMATE NUMBER OF SIGNIFICANT DECIMAL

C     DIGITS REPRESENTED BY THE MACHINE.

C     ******************************************************************

C

C     .. SCALAR ARGUMENTS ..

*      REAL X, Y, Z

      REAL*8 x, y, z

      INTEGER ifail

C     ..

C     .. LOCAL SCALARS ..

*+SELF,IF=IBM.

      CHARACTER*8 srname

*+SELF.

*      REAL ACC, C1, C2, C3, C4, C5, C6, CXN, CYN, CZN, EA, EB, EC,

      REAL*8 acc, c1, c2, c3, c4, c5, c6, cxn, cyn, czn, ea, eb, ec,

     * ed, ef, lamda, lolim, mu, pow4, rtx, rty, rtz, sigma, uplim,

     * xn, yn, zn

      INTEGER ind

C     .. FUNCTION REFERENCES ..

*      REAL SQRT, X02ABE, X02ACE

      REAL*8 sqrt, x02abf, x02acf

      INTEGER p01aae

*

      REAL*8 xarg

      sqrt(xarg)=dsqrt(xarg)

*

C     ..

*      DATA SRNAME /8H S21BCE /

      DATA srname /8h s21bcf /

*+SELF,IF=E09.

*      DATA ACC/3.8E-2/

*+SELF,IF=E10.

*      DATA ACC/2.6E-2/

*+SELF,IF=E13.

*      DATA ACC/8.3E-3/

*+SELF,IF=E16.

      DATA acc/2.6e-3/

*+SELF,IF=E18.

*      DATA ACC/1.2E-3/

*+SELF,IF=E20.

*      DATA ACC/5.6E-4/

*+SELF.

C

C     TEST FOR VALID ARGUMENTS

      ind = 1

      IF (x.LT.0.0 .OR. y.LT.0.0 .OR. x+y.EQ.0.0) go to 60

      ind = 2

      IF (z.LE.0.0) go to 60

      ind = 3

*      LOLIM = 2.0/X02ACE(0.0)**0.66667

      lolim = 2.0/x02acf(0.0d0)**0.66667

      IF (dmin1(x+y,z).LT.lolim) go to 60

      ind = 4

*      UPLIM = (0.1*ACC/X02ABE(0.0))**0.66667

      uplim = (0.1*acc/x02abf(0.0d0))**0.66667

      IF (dmax1(x,y,z).GE.uplim) go to 60

C

C     VALID CALL

      ifail = 0

      xn = x

      yn = y

      zn = z

      sigma = 0.0

      pow4 = 1.0

C

C     MAIN RECURSION

   20 mu = (xn+yn+3.0*zn)*0.2

      cxn = (mu-xn)/mu

      cyn = (mu-yn)/mu

      czn = (mu-zn)/mu

      IF (dmax1(abs(cxn),abs(cyn),abs(czn)).LT.acc) go to 40

      rtx = sqrt(xn)

      rty = sqrt(yn)

      rtz = sqrt(zn)

      lamda = rtx*rty + rty*rtz + rtz*rtx

      zn = zn + lamda

      sigma = sigma + pow4/(rtz*zn)

      pow4 = pow4*0.25

      zn = zn*0.25

      xn = (xn+lamda)*0.25

      yn = (yn+lamda)*0.25

      go to 20

C

C     POWER SERIES EVALUATION

   40 c1 = 3.0/14.0

      c2 = 9.0/88.0

      c3 = 9.0/52.0

      c4 = 1.0/6.0

      c5 = 9.0/22.0

      c6 = 3.0/26.0

      ea = cxn*cyn

      eb = czn*czn

      ec = 3.0*ea - 8.0*eb

      ed = ea - eb

      ef = ea - 6.0*eb

*      S21BCE = 3.0*SIGMA + POW4*(1.0+EF*(C2*EF-C1-C3*EC*CZN)+CZN*

      s21bcf = 3.0*sigma + pow4*(1.0+ef*(c2*ef-c1-c3*ec*czn)+czn*

     * (c4*ec+czn*(czn*c6*ea-c5*ed)))/(mu*sqrt(mu))

      go to 80

C

C     FAILURE EXITS

   60 ifail = p01aae(ifail,ind,srname)

*      S21BCE = 0.0

      s21bcf = 0.0

   80 RETURN

      END

      INTEGER FUNCTION p01aae(IFAIL, ERROR, SRNAME)

C     MARK 1 RELEASE.  NAG COPYRIGHT 1971

C     MARK 3 REVISED

C     MARK 4A REVISED, IER-45

C     MARK 4.5 REVISED

C     MARK 7 REVISED (DEC 1978)

C     RETURNS THE VALUE OF ERROR OR TERMINATES THE PROGRAM.

      INTEGER error, ifail, nout

C+SELF,IF=IBM.

      CHARACTER*8 srname

C+SELF.

C     TEST IF NO ERROR DETECTED

      IF (error.EQ.0) go to 20

C     DETERMINE OUTPUT UNIT FOR MESSAGE

      CALL x04aae(0,nout)

C     TEST FOR SOFT FAILURE

      IF (mod(ifail,10).EQ.1) go to 10

C     HARD FAILURE

      WRITE (nout,99999) srname, error

C     STOPPING MECHANISM MAY ALSO DIFFER

      stop

C     SOFT FAIL

C     TEST IF ERROR MESSAGES SUPPRESSED

   10 IF (mod(ifail/10,10).EQ.0) go to 20

      WRITE (nout,99999) srname, error

   20 p01aae = error

      RETURN

99999 FORMAT (1h0, 38herror detected by nag library routine , a8,

     * 11h - ifail = , i5//)

      END



      REAL*8 FUNCTION x02abf(X)

C     NAG COPYRIGHT 1975

C     MARK 4.5 RELEASE

*      REAL X, Z

      REAL*8 x, z

*+SELF,IF=R01.

*     DATA Z/Z00100000/

*+SELF,IF=R04.

*      DATA Z/0001 4000 0000 0000 0000B/

      DATA z/1.d-15/

*+SELF.

C     * RMIN *

C     RETURNS THE VALUE OF THE SMALLEST POSITIVE REAL FLOATING-

C     POINT NUMBER EXACTLY REPRESENTABLE ON THE COMPUTER

C     THE X PARAMETER IS NOT USED

*      X02ABE = Z

      x02abf = z

      RETURN

      END

*+DECK,X02ACE.

*      REAL FUNCTION X02ACE(X)

      REAL*8 FUNCTION x02acf(X)

C     NAG COPYRIGHT 1975

C     MARK 4.5 RELEASE

*      REAL X, Z

      REAL*8 x, z

*+SELF,IF=R01.

*      DATA Z/Z7FFFFFFF/

*+SELF,IF=R04.

*    DATA Z/3776 7777 7777 7777 7777B/

      DATA z/1.d+15/

*+SELF.

C     * RMAX *

C     RETURNS THE VALUE OF THE LARGEST POSITIVE REAL  FLOATING-

C     POINT NUMBER REPRESENTABLE ON THE COMPUTER

*      X02ACE = Z

      x02acf = z

      RETURN

      END



      SUBROUTINE x04aae(I,NERR)

C     MARK 7 RELEASE. NAG COPYRIGHT 1978

C     MARK 7C REVISED IER-190 (MAY 1979)

C     IF I = 0, SETS NERR TO CURRENT ERROR MESSAGE UNIT NUMBER

C     (STORED IN NERR1).

C     IF I = 1, CHANGES CURRENT ERROR MESSAGE UNIT NUMBER TO

C     VALUE SPECIFIED BY NERR.

C

C     *** NOTE ***

C     THIS ROUTINE ASSUMES THAT THE VALUE OF NERR1 IS SAVED

C     BETWEEN CALLS.  IN SOME IMPLEMENTATIONS IT MAY BE

C     NECESSARY TO STORE NERR1 IN A LABELLED COMMON

C     BLOCK /AX04AA/ TO ACHIEVE THIS.

C

C     .. SCALAR ARGUMENTS ..

      INTEGER i, nerr

C     ..

C     .. LOCAL SCALARS ..

      INTEGER nerr1

C     ..

C     DATA NERR1 /5/

      DATA nerr1 /6/

      IF (i.EQ.0) nerr = nerr1

      IF (i.EQ.1) nerr1 = nerr

      RETURN

      END

*     SUBROUTINE F02AFE(A, IA, N, RR, RI, INTGER, LFAIL)

      SUBROUTINE f02aff(A, IA, N, RR, RI, INTGER, LFAIL)

C     MARK 2 RELEASE. NAG COPYRIGHT 1972

C     MARK 3 REVISED.

C     MARK 4.5 REVISED

C

C     EIGENVALUES OF A REAL UNSYMMETRIC MATRIX

C     1ST AUGUST 1971

C

      INTEGER p01aae, isave, lfail, n, ia, intger(n)

      CHARACTER*8 srname

      REAL*8 acc, xxxx, a(ia,n), rr(n), ri(n), x02aae

      DATA srname /'  F02AFE'/

      isave = lfail

      lfail = 1

      acc = x02aae(xxxx)

      CALL f01akf(n, 1, n, a, ia, intger)

      CALL f02apf(n, acc, a, ia, rr, ri, intger, lfail)

      IF (lfail.NE.0) lfail = p01aae(isave,lfail,srname)

      RETURN

      END

*     SUBROUTINE F02APE(NN, ACC, H, IH, WR, WI, ICNT, LFAIL)

      SUBROUTINE f02apf(NN, ACC, H, IH, WR, WI, ICNT, LFAIL)

C     MARK 2 RELEASE. NAG COPYRIGHT 1972

C     MARK 3 REVISED.

C     MARK 4 REVISED.

C     MARK 4.5 REVISED

C     MARK 7D REVISED IER-197 (JUL 1979)

C     MARK 7E REVISED IER-203 (JUL 1979)

C     MARK 8 REVISED. IER-234 (APR 1980).

C     MARK 9 REVISED. IER-326 (SEP 1981).

C

C     HQR

C     FINDS THE EIGENVALUES OF A REAL UPPER HESSENBERG MATRIX, H,

C     STORED IN THE ARRAY H(N,N), AND STORES THE REAL PARTS IN

C     THE ARRAY WR(N) AND THE IMAGINARY PARTS IN THE ARRAY

C     WI(N). ACC IS THE RELATIVE MACHINE PRECISION. THE

C     SUBROUTINE FAILS IF ALL EIGENVALUES TAKE MORE THAN 30*N

C     ITERATIONS.

C     1ST DECEMBER 1971

C

      INTEGER p01aae, isave, lfail, n, nn, its, na, l, ll, i, m,

     * n2, mm, m2, m3, k, j, ih, icnt(nn), nhs, itn

      CHARACTER*8 srname

      REAL*8 t, acc, x, y, w, s, z, r, p, q, h(ih,nn), wr(nn), wi(nn)

      REAL*8 norm, x02ade

      LOGICAL notlst

      DATA srname /'  F02APE'/

      isave = lfail

      lfail = 0

      t = 0.0

      n = nn

      itn = 30*n

C     COMPUTE MATRIX NORM

      norm = 0.0

      k = 1

      DO 30 i = 1,n

         DO 25 j = k,n

            norm = norm + abs(h(i,j))

   25    CONTINUE

         k = i

   30 CONTINUE

      nhs = n*(n+1)/2 + n - 1

   20 IF (n.EQ.0) go to 600

      its = 0

      na = n - 1

C     LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT

   40 l = n + 1

      IF (n.LT.2) go to 80

      DO 60 ll=2,n

         l = l - 1

         s = abs(h(l-1,l-1)) + abs(h(l,l))

         IF (s.LT.x02ade(0.0d0)) s = norm/float(nhs)

         IF (abs(h(l,l-1)).LE.acc*s) go to 100

   60 CONTINUE

   80 l = 1

  100 x = h(n,n)

      IF (l.EQ.n) go to 520

      y = h(na,na)

      w = h(n,na)*h(na,n)

      IF (l.EQ.na) go to 540

      IF (itn.GT.0) go to 120

      lfail = p01aae(isave,1,srname)

      RETURN

  120 IF ((its.EQ.10) .OR. (its.EQ.20)) go to 140

      go to 200

C     FORM EXCEPTIONAL SHIFT

  140 t = t + x

      IF (na.LT.1) go to 180

      DO 160 i=1,na

         h(i,i) = h(i,i) - x

  160 CONTINUE

  180 h(n,n) = 0.0

      s = abs(h(n,na)) + abs(h(na,n-2))

      x = 0.75*s

      y = 0.75*s

      w = -0.4375*s**2

  200 its = its + 1

      itn = itn - 1

C     LOOK FOR TWO SMALL CONSECUTIVE SUB-DIAGONAL ELEMENTS

      IF (l.GT.(n-2)) go to 240

      m = n - 1

      n2 = n - 2

      DO 220 mm=l,n2

         m = m - 1

         z = h(m,m)

         r = x - z

         s = y - z

         p = (r*s-w)/h(m+1,m) + h(m,m+1)

         q = h(m+1,m+1) - z - r - s

         r = h(m+2,m+1)

         s = abs(p) + abs(q) + abs(r)

         p = p/s

         q = q/s

         r = r/s

         IF (m.EQ.l) go to 240

         IF ((abs(h(m,m-1))*(abs(q)+abs(r))).LE.(acc*abs(p)*

     *    (abs(h(m-1,m-1))+abs(z)+abs(h(m+1,m+1))))) go to 240

  220 CONTINUE

  240 m2 = m + 2

      IF (m2.GT.n) go to 280

      DO 260 i=m2,n

         h(i,i-2) = 0.0

  260 CONTINUE

  280 m3 = m + 3

      IF (m3.GT.n) go to 320

      DO 300 i=m3,n

         h(i,i-3) = 0.0

  300 CONTINUE

  320 IF (m.GT.na) go to 40

C     DOUBLE QR STEP INVOLVING ROWS L TO N AND COLUMNS M TO N

      DO 500 k=m,na

         notlst = .true.

         IF (k.EQ.na) notlst = .false.

         IF (k.EQ.m) go to 340

         p = h(k,k-1)

         q = h(k+1,k-1)

         r = 0.0

         IF (notlst) r = h(k+2,k-1)

         x = abs(p) + abs(q) + abs(r)

         IF (x.EQ.0.0) go to 500

         p = p/x

         q = q/x

         r = r/x

  340    s = sqrt(p**2+q**2+r**2)

         IF (p.LT.0.0) s = -s

         IF (k.NE.m) go to 360

         IF (l.NE.m) h(k,k-1) = -h(k,k-1)

         go to 380

  360    h(k,k-1) = -s*x

  380    p = p + s

         x = p/s

         y = q/s

         z = r/s

         q = q/p

         r = r/p

         IF (k.GT.n) go to 440

C     ROW MODIFICATION

         DO 420 j=k,n

            p = h(k,j) + q*h(k+1,j)

            IF (.NOT.notlst) go to 400

            p = p + r*h(k+2,j)

            h(k+2,j) = h(k+2,j) - p*z

  400       h(k+1,j) = h(k+1,j) - p*y

            h(k,j) = h(k,j) - p*x

  420    CONTINUE

  440    j = n

         IF ((k+3).LT.n) j = k + 3

         IF (l.GT.j) go to 500

C     COLUMN MODIFICATION

         DO 480 i=l,j

            p = x*h(i,k) + y*h(i,k+1)

            IF (.NOT.notlst) go to 460

            p = p + z*h(i,k+2)

            h(i,k+2) = h(i,k+2) - p*r

  460       h(i,k+1) = h(i,k+1) - p*q

            h(i,k) = h(i,k) - p

  480    CONTINUE

  500 CONTINUE

      go to 40

C     ONE ROOT FOUND

  520 wr(n) = x + t

      wi(n) = 0.0

      icnt(n) = its

      n = na

      go to 20

C     TWO ROOTS FOUND

  540 p = (y-x)/2.0

      q = p**2 + w

      y = sqrt(abs(q))

      x = x + t

      icnt(n) = -its

      icnt(na) = its

      IF (q.GT.0.0) go to 560

C     COMPLEX PAIR

      wr(na) = x + p

      wr(n) = x + p

      wi(na) = y

      wi(n) = -y

      go to 580

C     REAL PAIR

  560 IF (p.LT.0.0) y = -y

      y = p + y

      wr(na) = x + y

      wr(n) = x - w/y

      wi(n) = 0.0

      wi(na) = 0.0

  580 n = n - 2

      go to 20

  600 RETURN

      END

*     SUBROUTINE F01AKE(N, K, L, A, IA, INTGER)

      SUBROUTINE f01akf(N, K, L, A, IA, INTGER)

C     MARK 2 RELEASE. NAG COPYRIGHT 1972

C     MARK 4 REVISED.

C     MARK 4.5 REVISED

C     MARK 8 REVISED. IER-248 (JUN 1980).

C

C     DIRHES

C     AUGUST 1ST, 1971 .

C     GIVEN THE UNSYMMETRIC MATRIX, A, STORED IN THE ARRAY A(N,N),

C     THIS SUBROUTINE REDUCES THE SUB-MATRIX OF ORDER L - K + 1,

C     WHICH STARTS AT THE ELEMENT A(K,K) AND FINISHES AT THE

C     ELEMENT A(L,L), TO HESSENBERG FORM, H, BY THE DIRECT

C     METHOD(AN = NH). THE MATRIX H IS OVERWRITTEN ON A WITH

C     DETAILS OF THE TRANSFORMATIONS (N) STORED IN THE REMAINING

C     TRIANGLE UNDER H AND IN ELEMENTS K TO L OF THE ARRAY

C     INTGER(N).

C     1ST AUGUST 1971

C

      INTEGER ifail1, k1, k, l, j, m, i, n, j1, ia, ij, intger(n)

      REAL*8 x, y, d2, a(ia,n)

      ifail1 = 0

      k1 = k + 1

      IF (k1.GT.n) RETURN

      DO 200 j=k1,n

         m = j

         x = 0.0

         IF (j.GT.l) go to 120

         DO 20 i=j,l

            IF (abs(a(i,j-1)).LE.abs(x)) go to 20

            x = a(i,j-1)

            m = i

   20    CONTINUE

         intger(j) = m

         IF (m.EQ.j) go to 80

C     INTERCHANGE ROWS AND COLUMNS OF A.

         DO 40 i=k,n

            y = a(m,i)

            a(m,i) = a(j,i)

            a(j,i) = y

   40    CONTINUE

         DO 60 i=1,l

            y = a(i,m)

            a(i,m) = a(i,j)

            a(i,j) = y

   60    CONTINUE

   80    IF (x.EQ.0.0) go to 120

         IF (j.EQ.l) go to 120

         j1 = j + 1

         DO 100 i=j1,l

            a(i,j-1) = a(i,j-1)/x

  100    CONTINUE

  120    DO 180 i=1,l

            y = a(i,j)

            d2 = 0.0

            IF (x.EQ.0.0) go to 140

            IF (j.GE.l) go to 140

            j1 = j + 1

            CALL x03aaf(a(i,j1), (n-j)*ia-i+1, a(j1,j-1),

     *       (n-j+2)*ia-j, l-j, ia, 1, y, d2, y, d2, .true., ifail1)

  140       ij = j

            IF (i.LE.j) ij = i - 1

            y = -y

            d2 = -d2

            IF ((k+1).GT.ij) go to 160

            CALL x03aaf(a(i,k), (n-k+1)*ia-i+1, a(k1,j),

     *       (n-j+1)*ia-k, ij-k, ia, 1, y, d2, y, d2, .true.,

     *       ifail1)

  160       a(i,j) = -y

  180    CONTINUE

  200 CONTINUE

      RETURN

      END

*+PATCH,TEXTX03.

*+DECK,X03AAE.

*     SUBROUTINE X03AAE(A, ISIZEA, B, ISIZEB, N, ISTEPA, ISTEPB,

      SUBROUTINE x03aaf(A, ISIZEA, B, ISIZEB, N, ISTEPA, ISTEPB,

     * c1, c2, d1, d2, sw, ifail)

C     NAG COPYRIGHT 1975

C     MARK 4.5 RELEASE

C     MARK 6 REVISED

C

C     CALCULATES THE VALUE OF A SCALAR PRODUCT USING BASIC

C     OR ADDITIONAL PRECISION AND ADDS IT TO A BASIC OR ADDITIONAL

C     PRECISION INITIAL VALUE.

C

      INTEGER p01aae, isave, isizea, isizeb, istepa, istepb, ifail,

     * is, it, n, i

      DOUBLE PRECISION sum

*+SELF,IF=IBM.

      CHARACTER*8 srname

*+SELF.

      REAL*8 a(isizea), b(isizeb), c1, c2, d1, d2, x

      LOGICAL sw

      DATA srname /8h x03aae /

      isave = ifail

      ifail = 0

      IF (istepa.GT.0 .AND. istepb.GT.0) go to 20

      ifail = p01aae(isave,1,srname)

      RETURN

   20 is = 1 - istepa

      it = 1 - istepb

      IF (sw) go to 80

      x = 0.0

      IF (n.LT.1) go to 60

      DO 40 i=1,n

         is = is + istepa

         it = it + istepb

         x = x + a(is)*b(it)

   40 CONTINUE

   60 d1 = x + (c1+c2)

      d2 = 0.0

      RETURN

   80 sum = 0.0d0

      IF (n.LT.1) go to 120

      DO 100 i=1,n

         is = is + istepa

         it = it + istepb

         sum = sum + dble(a(is))*b(it)

  100 CONTINUE

  120 sum = sum + (dble(c1)+c2)

      d1 = sum + sum - dble(sngl(sum))

C     THE LAST STATEMENT ASSUMES THAT THE MACHINE SIMPLY

C     TRUNCATES WHEN ASSIGNING A DOUBLE PRECISION QUANTITY

C     TO A SINGLE PRECISION VARIABLE. IF INSTEAD THE MACHINE

C     ROUNDS, REPLACE THE LAST STATEMENT BY

C     D1 = SUM

      d2 = sum - d1

      RETURN

      END

*PATCH,TEXTX02.

*DECK,X02AAE.

      REAL*8 FUNCTION x02aae(X)

C     NAG COPYRIGHT 1975

C     MARK 4.5 RELEASE

      REAL*8 x, z

*SELF,IF=R01.

C      DATA Z/Z3C100000/

*SELF,IF=R04.

C      DATA Z/1641 4000 0000 0000 0000B/

      DATA z/1.d-15/

*SELF.

C     * EPS *

C     RETURNS THE VALUE EPS WHERE EPS IS THE SMALLEST

C     POSITIVE

C     NUMBER SUCH THAT 1.0 + EPS > 1.0

C     THE X PARAMETER IS NOT USED

      x02aae = z

      RETURN

      END

*+DECK,X02ADE.

      REAL*8 FUNCTION x02ade(X)

C     NAG COPYRIGHT 1975

C     MARK 4.5 RELEASE

      REAL*8 x, z

*+SELF,IF=R01.

*       DATA Z/Z05100000/

*+SELF,IF=R04.

*      DATA Z/0060 4000 0000 0000 0000B/

*+SELF.

      DATA z/1.d-15/

C     * TOL *

C     RETURNS THE RATIO OF THE SMALLEST POSITIVE REAL FLOATING-

C     POINT NUMBER REPRESENTABLE ON THE COMPUTER TO EPS

      x02ade = z

      RETURN

      END

*+DECK,X02ADE.

      SUBROUTINE e04abe(FUN, EPS, T, A, B, MAXCAL, X, F, IFAIL)

C

C     MARK 6 RELEASE NAG COPYRIGHT 1977

C     MARK 8 REVISED. IER-231 (MAR 1980).

C     MARK 8D REVISED. IER-272 (DEC 1980).

C

C     **************************************************************

C

C     E04ABE ATTEMPTS TO FIND A MINIMUM IN AN INTERVAL A .LE. X .LE.

C     B OF A FUNCTION F(X) OF THE SCALAR X, USING FUNCTION VALUES

C     ONLY.

C

C     IT IS BASED ON THE SUBROUTINE UNIFUN IN THE NPL ALGORITHMS

C     LIBRARY (REF. NO. E4/13/F). THE FUNCTION F(X) IS DEFINED BY

C     THE USER-SUPPLIED SUBROUTINE FUN. T AND EPS DEFINE A TOLERANCE

C     TOL = EPS * ABS(X) + T, AND FUN IS NEVER EVALUATED AT TWO

C     POINTS CLOSER THAN TOL. IF FUN IS DELTA-UNIMODAL, FOR SOME

C     DELTA LESS THAN TOL, THEN X APPROXIMATES THE GLOBAL MINIMUM OF

C     FUN WITH AN ERROR LESS THAN 3*TOL. IF FUN IS NOT DELTA-

C     UNIMODAL ON (A, B), THEN X MAY APPROXIMATE A LOCAL, BUT NON

C     GLOBAL, MINIMUM. EPS SHOULD BE NO SMALLER THAN 2*EPSMCH, AND

C     PREFERABLY NOT MUCH LESS THAN SQRT(EPSMCH), WHERE EPSMCH IS

C     THE RELATIVE MACHINE PRECISION. T SHOULD BE POSITIVE. NOTE

C     THAT, FOR CONSISTENCY WITH OTHER E04 DOCUMENTATION, THE NAME

C     FUNCT IS USED INSTEAD OF FUN IN THE WRITE-UP.

C

C     PHILIP E. GILL, WALTER MURRAY, SUSAN M. PICKEN, HAZEL M.

C     BARBER AND MARGARET H. WRIGHT, D.N.A.C., NATIONAL PHYSICAL

C     LABORATORY, ENGLAND

C

C     **************************************************************

C

C     .. SCALAR ARGUMENTS ..

      REAL*8 a, b, eps, f, t, x

      INTEGER ifail, maxcal

C     .. SUBROUTINE ARGUMENTS ..

C     FUN

C     ..

C     .. LOCAL SCALARS ..

      CHARACTER*8 srname

      REAL*8 b1, d, e, epsmch, f1, f2, fa, fu, fv, fw, gtest1,

     * gtest2, gu, oldf, pt2, pt4, pt6, rr, rteps, scxbd, ss, tol,

     * u, x1, x2, xlamda, xv, xw, sftbnd

      INTEGER iflag, iloc, isave, numf

C     .. FUNCTION REFERENCES ..

      REAL*8 sqrt, xarg, x02aae

      INTEGER p01aae

        sqrt(xarg)=dsqrt(xarg)



C     .. SUBROUTINE REFERENCES ..

C     ABZE04

C     ..

      DATA srname /'  E04ABE'/

      isave = ifail

C

C     A MACHINE-DEPENDENT CONSTANT IS SET HERE. EPSMCH IS THE

C     SMALLEST POSITIVE REAL NUMBER SUCH THAT 1.0 + EPSMCH .GT. 1.0

C

      epsmch = x02aae(u)

      rteps = sqrt(epsmch)

      IF (eps.LT.epsmch) eps = rteps

      IF (t.LT.epsmch) t = rteps

C

C     ERROR IN INPUT PARAMETERS

C

      ifail = 1

      IF (a+t.GE.b .OR. maxcal.LT.3 .OR. ifail.LT.0 .OR.

     * ifail.GT.1) go to 140

      sftbnd = a

      pt2 = (b-a)*2.0e-1

      pt4 = pt2 + pt2

      pt6 = pt2 + pt4

      x1 = a + pt4

      CALL fun(x1, f1)

      x2 = b - pt4

      CALL fun(x2, f2)

      xlamda = b

      IF (f1.GT.f2) go to 20

      x = x1

      a = -pt4

      b = pt2

      xw = pt2

      b1 = pt2

      rr = 1.0e+0

      d = -pt2

      fw = f2

      fv = f1

      f = f1

C

C     SET STEP TO NEW POINT

C

      u = -pt2

      go to 40

   20 x = x2

      a = -pt2

      b = pt4 + eps*abs(xlamda) + t

      xw = -pt2

      b1 = b

      rr = -1.0e+0

      d = pt2

      fw = f1

      fv = f2

      f = f2

C

C     SET STEP TO NEW POINT

C

      u = pt2

   40 xv = 0.0e+0

      scxbd = pt4

      e = pt6

      ss = 0.0e+0

      fa = fw + t

      oldf = fa

      gtest1 = 0.0e+0

      gtest2 = 0.0e+0

      tol = eps*abs(x) + t

      CALL fun(x+u, fu)

      gu = 0.0e+0

      numf = 3

C

C     SET ILOC TO 3 SO THAT THE MAIN SECTION OF ABZE04 IS EXECUTED

C     AS THE INITIAL 3 POINTS HAVE ALREADY BEEN SET UP

C

      iloc = 3

   60 CALL abze04(eps, t, 0.0d+0, sftbnd, xlamda, u, fu, gu, x, f,

     * xw, fw, xv, fv, a, fa, b, oldf, b1, scxbd, e, d, rr, ss,

     * gtest1, gtest2, tol, iloc, iflag)

      IF (iflag.NE.1) go to 100

      IF (numf.GE.maxcal) go to 80

      CALL fun(x+u, fu)

      numf = numf + 1

      go to 60

   80 ifail = 2

      go to 120

  100 ifail = 0

  120 maxcal = numf

      a = a + x

      b = b + x

  140 CONTINUE

      IF (ifail.EQ.0) RETURN

      ifail = p01aae(isave,ifail,srname)

      RETURN

C

C     END OF E04ABE (UNIFUN)

C

      END

      SUBROUTINE abze04(EPS, T, ETA, SFTBND, XLAMDA, U, FU, GU,

     * xmin, fmin, xw, fw, xv, fv, a, fa, b, oldf, b1, scxbd, e, d,

     * rr, ss, gtest1, gtest2, tol, iloc, itest)

C

C     MARK 6 RELEASE NAG COPYRIGHT 1977

C     MARK 7 REISSUE

C     MARK 8 REVISED. IER-239 (APR 1980).

C     MARK 8 REVISED. IER-244 (MAY 1980).

C     MARK 9 REVISED. IER-317 (SEP 1981).

C

C     **************************************************************

C

C     ABZE04 (NEWPTQ), AN ALGORITHM FOR FINDING A STEPLENGTH, CALLED

C     REPEATEDLY BY NPL OPTIMIZATION ROUTINES WHICH REQUIRE A STEP-

C     LENGTH TO BE COMPUTED USING QUADRATIC INTERPOLATION.

C     THE PARAMETERS SET UP BEFORE THE CALL OF ABZE04 CONTAIN

C     INFORMATION ABOUT THE INTERVAL IN WHICH A LOWER POINT IS TO BE

C     FOUND AND FROM THIS ABZE04 PRODUCES A POINT AT WHICH THE

C     FUNCTION CAN BE EVALUATED OUTSIDE THIS SUBROUTINE.

C     THE VALUE OF THE INTEGER PARAMETER ILOC DETERMINES THE PATH

C     TAKEN THROUGH THE CODE. FOR A FURTHER DESCRIPTION OF ILOC

C     AND THE OTHER PARAMETERS SEE NPL ALGORITHMS LIBRARY REF. NO.

C     E4/15/F.

C

C     PHILIP E. GILL, WALTER MURRAY, SUSAN M. PICKEN,

C     MARGARET H. WRIGHT AND ENID M. LONG

C     D.N.A.C. NATIONAL PHYSICAL LABORATORY, ENGLAND

C

C     **************************************************************

C

C     .. SCALAR ARGUMENTS ..

      REAL*8 a, b1, b, d, e, eps, eta, fa, fmin, fu, fv, fw, gtest1,

     * gtest2, gu, oldf, rr, scxbd, sftbnd, ss, t, tol, u, xlamda,

     * xmin, xv, xw

      INTEGER iloc, itest

C     ..

C     .. LOCAL SCALARS ..

      REAL*8 a1, d1, d2, q, r, s, t2, xm

C     .. FUNCTION REFERENCES ..

      REAL*8 sqrt, xarg

        sqrt(xarg)=dsqrt(xarg)

C     ..

C

C     BRANCH TO APPROPRIATE SECTION OF CODE DEPENDING ON THE

C     VALUE OF THE FLAG ILOC

C     THE SIGNIFICANCE OF THE FLAGS ILOC AND ITEST ARE DESCRIBED IN

C     NPL ALGORITHMS LIBRARY DOCUMENT REF. NO. E4/15/F.

C

      go to(20, 40, 40, 460, 440), iloc

C

C     ILOC = 1

C

C     CHECK INPUT PARAMETERS

C

   20 itest = 2

      tol = t

      t2 = tol + tol

      IF (u.LE.0.0e+0 .OR. xlamda.LE.t2 .OR. gu.GT.0.0e+0) RETURN

      itest = 1

C

C     A AND B DEFINE THE INTERVAL OF UNCERTAINTY. XMIN DENOTES

C     THE LOWEST POINT OBTAINED SO FAR, XW THE LAST VALUE

C     OF XMIN AND XV THE SCALED VALUE OF ALPHA CORRESPONDING TO

C     THE HIGHEST FUNCTION VALUE OF THE THREE POINTS THROUGH

C     WHICH A PARABOLA MAY BE FITTED. INITIALIZE A, XV, XW, XMIN

C     AT ORIGIN AND CORRESPONDING FUNCTION VALUES AT LATEST

C     ESTIMATE OF MINIMUM.

C

      xmin = 0.0e+0

      xw = 0.0e+0

      xv = 0.0e+0

      a = 0.0e+0

      oldf = fu

      fmin = fu

      fw = fu

      fv = fu

      fa = fu

      d = u

C

C     THE PARAMETER RR HAS TWO USES DURING THE EXECUTION OF THIS

C     SUBROUTINE.  INITIALLY THE SIGN OF RR INDICATES WHETHER OR NOT

C     THE MINIMUM HAS BEEN BRACKETED. LATER, WHEN A POINT SATISFYING

C     THE GTEST2 CRITERION HAS BEEN FOUND, RR IS USED TO COMPUTE A

C     STEPLENGTH WHICH SATISFIES THE SECOND CRITERION INVOLVING

C     GTEST1.

C

      rr = -1.0e+0

C

C     SET UP XBND AS A BOUND ON THE STEP TO BE TAKEN. (XBND IS NOT

C     COMPUTED EXPLICITLY BUT SCXBD IS ITS SCALED VALUE.) SET THE

C     UPPER BOUND ON THE INTERVAL OF UNCERTAINTY INITIALLY TO

C     XLAMDA + TOL(XLAMDA).

C

      scxbd = xlamda

      b = scxbd + eps*abs(scxbd) + t

      e = 2.0e+0*b

      b1 = b

C

C     COMPUTE THE CONSTANTS REQUIRED FOR THE TWO CONVERGENCE

C     CRITERIA.

C

      gtest1 = -1.0e-4*gu

      gtest2 = -eta*gu

C

C     SET ILOC TO INDICATE THAT ONLY TWO POINTS ARE AVAILABLE

C

      iloc = 2

      go to 380

C

C     ILOC = 2 OR 3

C

C     UPDATE A, B, XV, XW, AND XMIN.

C

   40 IF (fu.GT.fmin) go to 100

C

C     IF FUNCTION VALUE NOT INCREASED, NEW POINT BECOMES

C     NEXT ORIGIN AND OTHER POINTS ARE SCALED ACCORDINGLY.

C

      IF (u.LT.0.0e+0) go to 60

      a = 0.0e+0

      fa = fmin

      IF (xw.EQ.xv .AND. fmin.EQ.fu) rr = 1.0e+0

      go to 80

   60 b = 0.0e+0

      rr = 1.0e+0

   80 xv = xw

      fv = fw

      fw = fmin

      fmin = fu

      xmin = xmin + u

      a = a - u

      b = b - u

      xv = xv - u

      xw = 0.0e+0 - u

C

C     THIS MAY BE CHANGED TO XW = - U IF THE COMPUTER IS

C     SUCH THAT  - U AND 0.0 - U ARE IDENTICAL.

C

      scxbd = scxbd - u

      tol = eps*abs(xmin) + t

      go to 180

C

C     IF FUNCTION VALUE INCREASED, ORIGIN REMAINS UNCHANGED

C     BUT OTHER POINTS MAY BE INTERCHANGED.

C

  100 IF (u.GE.0.0e+0) go to 120

      a = u

      fa = fu

      go to 140

  120 b = u

      rr = 1.0e+0

  140 IF (fu.GT.fw .AND. xw.NE.0.0e+0) go to 160

      xv = xw

      fv = fw

      xw = u

      fw = fu

      go to 180

  160 xv = u

      fv = fu

  180 t2 = tol + tol

      xm = 5.0e-1*(a+b)

C

C     CHECK TERMINATION CRITERIA.

C

      IF (abs(xm).LE.t2-5.0e-1*(b-a) .OR. xmin+b.LE.sftbnd .OR.

     * fa-fmin.LE.abs(a)*gtest2 .AND. fmin.LT.oldf .AND.

     * (abs(xmin-xlamda).GT.tol .OR. rr.LT.0.0e+0)) go to 420

      r = 0.0e+0

      q = 0.0e+0

      s = 0.0e+0

      IF (abs(e).LE.tol) go to 240

C

C     FIT PARABOLA THROUGH XMIN, XV, XW.

C

      IF (iloc.NE.2) go to 200

C

C     SPECIAL CASE. ONLY TWO POINTS ARE AVAILABLE FOR

C     QUADRATIC INTERPOLATION

C

      q = 2.0e+0*(fw-fmin-xw*gu)

      s = gu*xw*xw

      IF (xmin.NE.0.0e+0) s = (2.0e+0*(fmin-fw)+xw*gu)*xw

      go to 220

  200 r = xw*(fv-fmin)

      q = xv*(fw-fmin)

      s = r*xw - q*xv

      q = 2.0e+0*(q-r)

  220 IF (q.GT.0.0e+0) s = -s

      IF (q.LE.0.0e+0) q = -q

      r = e

C

C     IF THE LAST STEP EXPANDED THE INTERVAL OR THE MINIMUM HAS

C     ALREADY BEEN BRACKETED SET E AS THE LAST STEP TAKEN

C

      IF (d.NE.b1 .OR. rr.GT.0.0e+0) e = d

C

C     CONSTRUCT AN ARTIFICIAL BOUND ON THE ESTIMATED STEPLENGTH.

C

  240 a1 = a

      b1 = b

      IF (xmin.NE.a) go to 260

      d = xm

      go to 330

  260 IF (rr.GT.0.0e+0) go to 280

      d = -4.0e+0*a

      IF (d.GE.scxbd) d = scxbd

      go to 320

C

C     DETERMINE INTERVAL OF LENGTH D2 IN WHICH TO SET

C     ARTIFICIAL BOUND.

C

  280 d1 = a

      d2 = b

      IF (abs(d2).GT.tol .AND. (xw.LE.0.0e+0 .OR. abs(d1).LE.tol))

     * go to 300

      u = d1

      d1 = d2

      d2 = u

  300 u = -d1/d2

      IF (u.GE.1.0e+0) d = 5.0e+0*d2*(1.0e-1+1.0e+0/u)/1.1e+1

      IF (u.LT.1.0e+0) d = 5.0e-1*d2*sqrt(u)

C

C     IF THE MINIMUM IS BRACKETED BY XV AND XW THE STEP MUST LIE

C     WITHIN (A, B).

C

  320 IF (xw.LT.0.0e+0 .AND. xv.GT.0.0e+0 .OR. xw.GT.0.0e+0 .AND.

     * xv.LT.0.0e+0) go to 330

C

C     IF THE MINIMUM IS NOT BRACKETED BY XV AND XW THE STEP MUST LIE

C     WITHIN (A1, B1).

C

      IF (d.LE.0.0e+0) a1 = d

      IF (d.GT.0.0e+0) b1 = d

C

C     REJECT THE STEP OBTAINED BY INTERPOLATION IF IT LIES OUTSIDE

C     THE REQUIRED INTERVAL OR IT IS GREATER THAN HALF THAT

C     OBTAINED DURING THE LAST-BUT-ONE ITERATION.

C

  330 IF (abs(s).GE.abs(5.0e-1*q*r) .OR. s.LE.q*a1 .OR. s.GE.q*b1)

     * go to 340

C

C     A PARABOLIC INTERPOLATION STEP.

C

      d = s/q

C

C     F MUST NOT BE EVALUATED TOO CLOSE TO A OR B.

C

      IF (d-a.GE.t2 .AND. b-d.GE.t2) go to 360

      d = tol

      IF (xm.LE.0.0e+0) d = -tol

      go to 360

C

C     A NON-INTERPOLATION STEP.

C

  340 e = b

      IF (xm.LE.0.0e+0) e = a

  360 iloc = 3

C

C     CHECK THAT THE NEW STEP LENGTH WILL NOT BE GREATER THAN

C     XLAMDA.

C

  380 IF (d.LT.scxbd) go to 400

C

C     REPLACE THE STEP LENGTH BY THE SCALED BOUND (SO AS TO COMPUTE

C     THE NEW POINT ON THE BOUNDARY.

C

      d = scxbd

C

C     MOVE SCXBD TO THE LEFT SO THAT NEWXBND + TOL(NEWXBND) = XBND.

C

      scxbd = (scxbd - tol)/(1.0e+0 + eps)

  400 u = d

      IF (abs(d).LT.tol .AND. d.LE.0.0e+0) u = -tol

      IF (abs(d).LT.tol .AND. d.GT.0.0e+0) u = tol

      itest = 1

      RETURN

C

C     THE FIRST CONVERGENCE CRITERION HAS BEEN SATISFIED. NOW CHECK

C     THAT THE FUNCTION VALUE HAS BEEN REDUCED SUFFICIENTLY. THE

C     VARIABLE RR IS NOW USED TO REDUCE THE STEP LENGTH.

C

  420 d = rr

      rr = xmin

      ss = 5.0e-1

      fu = fmin

      IF (xmin.EQ.0.0e+0) xmin = t

C

C     IF XMIN LIES WITHIN TOL OF THE BOUNDARY AND THE MINIMUM HAS

C     BEEN BRACKETED, THEN RECOMPUTE THE POINT ON THE BOUNDARY.

C

  440 IF (abs(xmin-xlamda).GE.tol .OR. xmin.EQ.t) go to 460

      IF (scxbd.LT.0.0e+0 .AND. xw.LT.0.0e+0 .AND. xv.LT.0.0e+0)

     * xmin = xlamda

      IF (d.LT.0.0e+0) go to 460

      u = 0.0e+0

      iloc = 4

      itest = 1

      RETURN

C

C     CHECK THAT THE NEW POINT SATISFIES SAFEGUARD CONDITIONS.

C     IF NECESSARY ATTEMPT TO FIND A SUFFICIENTLY LOWER POINT

C     BY SUCCESSIVELY DECREASING THE STEPLENGTH.

C

  460 IF (xmin+b.GT.sftbnd) go to 480

      itest = 4

C

C     ITEST = 4 IMPLIES THAT THE REQUIRED STEP LENGTH IS SMALLER

C     THAN SFTBND.

C

      RETURN

  480 IF (oldf-fu.LE.gtest1*xmin) go to 500

      fmin = fu

      itest = 0

C

C     THE ALGORITHM HAS SUCCESSFULLY FOUND A SUFFICIENTLY LOWER

C     POINT.

C

      RETURN

  500 IF (xmin.NE.t) go to 520

      itest = 3

C

C     DESPITE REPEATED REDUCTIONS IN THE STEP SIZE, A LOWER POINT

C     COULD NOT BE FOUND.

C

      RETURN

C

C     A SUFFICIENT REDUCTION IN THE FUNCTION VALUE HAS NOT YET BEEN

C     FOUND, TRY A FURTHER REDUCTION IN THE STEP LENGTH.

C

  520 xmin = rr*ss

      ss = ss*ss

      IF (xmin.LT.t) xmin = t

      itest = 1

      u = 0.0e+0

      iloc = 5

      RETURN

C

C     END OF ABZE04 (NEWPTQ)

C

      END

C

      SUBROUTINE e01baf(M, X, Y, K, C, LCK, WRK, LWRK, IFAIL)

C     MARK 8 RELEASE. NAG COPYRIGHT 1979.

C

C     ******************************************************

C

C     NPL ALGORITHMS LIBRARY ROUTINE SP3INT

C

C     CREATED 16/5/79.                        RELEASE 00/00

C

C     AUTHORS ... GERALD T. ANTHONY, MAURICE G.COX

C                 J.GEOFFREY HAYES AND MICHAEL A. SINGER.

C     NATIONAL PHYSICAL LABORATORY, TEDDINGTON,

C     MIDDLESEX TW11 OLW, ENGLAND

C

C     ******************************************************

C

C     E01BAE.  AN ALGORITHM, WITH CHECKS, TO DETERMINE THE

C     COEFFICIENTS IN THE B-SPLINE REPRESENTATION OF A CUBIC

C     SPLINE WHICH INTERPOLATES (PASSES EXACTLY THROUGH) A

C     GIVEN SET OF POINTS.

C

C     INPUT PARAMETERS

C        M        THE NUMBER OF DISTINCT POINTS WHICH THE

C                    SPLINE IS TO INTERPOLATE.

C                    (M MUST BE AT LEAST 4.)

C        X        ARRAY CONTAINING THE DISTINCT VALUES OF THE

C                    INDEPENDENT VARIABLE. NB X(I) MUST BE

C                    STRICTLY GREATER THAN X(J) WHENEVER I IS

C                    STRICTLY GREATER THAN J.

C        Y        ARRAY CONTAINING THE VALUES OF THE DEPENDENT

C                    VARIABLE.

C        LCK      THE SMALLER OF THE ACTUALLY DECLARED DIMENSIONS

C                    OF K AND C. MUST BE AT LEAST M + 4.

C

C     OUTPUT PARAMETERS

C        K        ON SUCCESSFUL EXIT, K CONTAINS THE KNOTS

C                    SET UP BY THE ROUTINE. IF THE SPLINE IS

C                    TO BE EVALUATED (BY NPL ROUTINE E02BEF,

C                    FOR EXAMPLE) THE ARRAY K MUST NOT BE

C                    ALTERED BEFORE CALLING THAT ROUTINE.

C        C        ON SUCCESSFUL EXIT, C CONTAINS THE B-SPLINE

C                    COEFFICIENTS OF THE INTERPOLATING SPLINE.

C                    THESE ARE ALSO REQUIRED BY THE EVALUATING

C                    ROUTINE E02BEF.

C        IFAIL    FAILURE INDICATOR

C                    0 - SUCCESSFUL TERMINATION.

C                    1 - ONE OF THE FOLLOWING CONDITIONS HAS

C                        BEEN VIOLATED -

C                        M AT LEAST 4

C                        LK AT LEAST M + 4

C                        LWORK AT LEAST 6 * M + 16

C                    2 - THE VALUES OF THE INDEPENDENT VARIABLE

C                        ARE DISORDERED. IN OTHER WORDS, THE

C                        CONDITION MENTIONED UNDER X IS NOT

C                        SATISFIED.

C

C     WORKSPACE (AND ASSOCIATED DIMENSION) PARAMETERS

C        WRK     WORKSPACE ARRAY, OF LENGTH LWRK.

C        LWRK    ACTUAL DECLARED DIMENSION OF WRK.

C                    MUST BE AT LEAST 6 * M + 16.

C

C     .. SCALAR ARGUMENTS ..

      INTEGER ifail, lck, lwrk, m

C     .. ARRAY ARGUMENTS ..

      REAL*8 c(lck), k(lck), wrk(lwrk), x(m), y(m)

C     ..

C     .. LOCAL SCALARS ..



c      DOUBLE PRECISION SRNAME

      CHARACTER*8 srname



      REAL*8 one, ss

      INTEGER i, ierror, m1, m2

C     .. FUNCTION REFERENCES ..

      INTEGER p01aae

C     .. SUBROUTINE REFERENCES ..

C     E02BAE

C     ..

      DATA one /1.0e+0/

      DATA srname /'  E01BAE'/

      ierror = 1

C

C     TESTS FOR ADEQUACY OF ARRAY LENGTHS AND THAT M IS GREATER

C     THAN 4.

C

      IF (lwrk.LT.6*m+16 .OR. m.LT.4) go to 80

      IF (lck.LT.m+4) go to 80

C

C     TESTS FOR THE CORRECT ORDERING OF THE X(I)

C

      ierror = 2

      DO 20 i=2,m

         IF (x(i).LE.x(i-1)) go to 80

   20 CONTINUE

C

C     INITIALISE THE ARRAY OF KNOTS AND THE ARRAY OF WEIGHTS

C

      wrk(1) = one

      wrk(2) = one

      wrk(3) = one

      wrk(4) = one

      IF (m.EQ.4) go to 60

      DO 40 i=5,m

         k(i) = x(i-2)

         wrk(i) = one

   40 CONTINUE

   60 m1 = m + 1

      m2 = m1 + m

C

C     CALL THE SPLINE FITTING ROUTINE

C

      ierror = 0

      CALL e02baf(m, m+4, x, y, wrk, k, wrk(m1), wrk(m2), c, ss,

     * ierror)

C

C     ALL THE TESTS PERFORMED BY E02BAE ARE REDUNDANT

C     BECAUSE OF THE ABOVE TESTS AND ASSIGNMENTS, AND SO

C     IERROR = 0 AFTER THIS CALL.

C

   80 ifail = p01aae(ifail,ierror,srname)

      RETURN

C

C     END OF E01BAE.

C

      END

      SUBROUTINE e02bcf(NCAP7, K, C, X, LEFT, S, IFAIL)

C     MARK 7 RELEASE. NAG COPYRIGHT 1978.

C

C     **************************************************

C     *                                                *

C     *         NAG LIBRARY SUBROUTINE  E02BCE         *

C     *                                                *

C     *  EVALUATION OF CUBIC SPLINE AND ITS            *

C     *  DERIVATIVES FROM ITS B-SPLINE REPRESENTATION  *

C     *                                                *

C     *  ROUTINE CREATED ... 17 NOV 1977               *

C     *  LATEST UPDATE ....  24 APR 1978               *

C     *  RELEASE NUMBER ...  01                        *

C     *  AUTHORS ... MAURICE G. COX AND                *

C     *              J. GEOFFREY HAYES, N.P.L.         *

C     *                                                *

C     **************************************************

C

C     .. SCALAR ARGUMENTS ..

      REAL*8 x

      INTEGER ifail, left, ncap7

C     .. ARRAY ARGUMENTS ..

      REAL*8 c(ncap7), k(ncap7), s(4)

C     ..

C     .. LOCAL SCALARS ..

c      DOUBLE PRECISION SRNAME

      CHARACTER*8 srname

      REAL*8 c1, c2, c3, c4, d1n41, d1n42, d1n43, d1n44, d2n41,

     * d2n42, d2n43, d2n44, d3n41, d3n42, d3n43, d3n44, e2, e3, e4,

     * e5, half, k1, k2, k3, k4, k5, k6, m11, m21, m22, m32, n41,

     * n42, n43, n44, p4, p5, p6, six, three

      INTEGER ierror, j1, j, l

C     .. FUNCTION REFERENCES ..

      INTEGER p01aae

C     ..

      DATA srname /'  E02BCE'/

      half = 0.5e+00

      three = 3.0e+00

      six = 6.0e+00

C

C     *********  DATA VALIDATION  *********

C

C     CHECK WHETHER AT LEAST ONE INTERVAL HAS BEEN SPECIFIED -

C

      ierror = 1

      IF (ncap7.LT.8) go to 80

C

C     CHECK WHETHER THE RANGE OF DEFINITION OF THE SPLINE IS

C     STRICTLY POSITIVE IN LENGTH -

C

      ierror = 2

      IF (k(4).GE.k(ncap7-3)) go to 80

C

C     CHECK WHETHER  X  IS A VALID ARGUMENT -

C

      IF (x.LT.k(4) .OR. x.GT.k(ncap7-3)) go to 80

C

C     *********  COMPUTATION  *********

C

C     BINARY SEARCH FOR INTERVAL CONTAINING  X  -

C

C     IF RIGHT-HAND DERIVATIVES ARE REQUIRED  (LEFT .NE. 1)

C     SEARCH FOR  J  SATISFYING

C     K(J + 3) .LE. X .LT. K(J + 4)

C     (SETTING  J = NCAP  IN THE EXCEPTIONAL CASE  X = K(NCAP +

C     4)).

C

C     IF LEFT-HAND DERIVATIVES ARE REQUIRED  (LEFT .EQ. 1)

C     SEARCH FOR  J  SATISFYING

C     K(J + 3) .LT. X .LE. K(J + 4)

C     (SETTING  J = 1  IN THE EXCEPTIONAL CASE  X = K(4)).

C

      ierror = 0

      j1 = 4

      j = ncap7 - 3

   20 l = (j1+j)/2

      IF (j-j1.LE.1) go to 60

      IF (left.NE.1 .AND. x.GE.k(l)) go to 40

      IF (left.EQ.1 .AND. x.GT.k(l)) go to 40

      j = l

      go to 20

   40 j1 = l

      go to 20

   60 j = j - 4

C

C     FORM CERTAIN CONSTANTS -

C

      k1 = k(j+1)

      k2 = k(j+2)

      k3 = k(j+3)

      k4 = k(j+4)

      k5 = k(j+5)

      k6 = k(j+6)

      e2 = x - k2

      e3 = x - k3

      e4 = k4 - x

      e5 = k5 - x

      p4 = k4 - k1

      p5 = k5 - k2

      p6 = k6 - k3

C

C     FORM BASIS FUNCTIONS AND THEIR DERIVATIVES -

C

C     THE VALUES OF THE NON-ZERO UN-NORMALIZED B-SPLINES OF ORDER R

C     ARE DENOTED BY  MR1, MR2,...  .  THE CORRESPONDING NORMALIZED

C     B-SPLINES ARE DENOTED BY  NR1, NR2,...  AND THEIR DERIVATIVES

C     OF ORDER  L  BY  DLNR1, DLNR2,...  .

C

      m11 = six/(k4-k3)

      m21 = -m11/(k4-k2)

      m22 = m11/(k5-k3)

      d3n41 = m21/p4

      m32 = (m21-m22)/p5

      d3n44 = m22/p6

      d3n42 = -d3n41 - m32

      d3n43 = m32 - d3n44

      m21 = -e4*m21

      m22 = e3*m22

      d2n41 = m21/p4

      m32 = (m21-m22)/p5

      d2n44 = m22/p6

      d2n42 = -d2n41 - m32

      d2n43 = m32 - d2n44

      m21 = half*m21

      m22 = half*m22

      d1n41 = -e4*m21/p4

      m32 = (e2*m21+e5*m22)/p5

      d1n44 = e3*m22/p6

      d1n42 = -d1n41 - m32

      d1n43 = m32 - d1n44

      n41 = -e4*d1n41/three

      n42 = (-(x-k1)*d1n41+e5*m32)/three

      n43 = (e2*m32+(k6-x)*d1n44)/three

      n44 = e3*d1n44/three

C

C     FORM THE VALUES OF THE CUBIC SPLINE AND ITS DERIVATIVES -

C

      c1 = c(j)

      c2 = c(j+1)

      c3 = c(j+2)

      c4 = c(j+3)

      s(1) = c1*n41 + c2*n42 + c3*n43 + c4*n44

      s(2) = c1*d1n41 + c2*d1n42 + c3*d1n43 + c4*d1n44

      s(3) = c1*d2n41 + c2*d2n42 + c3*d2n43 + c4*d2n44

      s(4) = c1*d3n41 + c2*d3n42 + c3*d3n43 + c4*d3n44

C

C     *********  ERROR DIAGNOSTICS  *********

C

   80 ifail = p01aae(ifail,ierror,srname)

      RETURN

C

C     END OF SUBROUTINE  E02BCF

C

      END

      SUBROUTINE e02baf(M, NCAP7, X, Y, W, K, WORK1, WORK2, C, SS,

     * ifail)

C     NAG COPYRIGHT 1975

C     MARK 5 RELEASE

C     MARK 6 REVISED  IER-84

C     MARK 8 RE-ISSUE. IER-224 (APR 1980).

C     MARK 9A REVISED. IER-356 (NOV 1981)

C

C     NAG LIBRARY SUBROUTINE  E02BAE

C

C     E02BAE  COMPUTES A WEIGHTED LEAST-SQUARES APPROXIMATION

C     TO AN ARBITRARY SET OF DATA POINTS BY A CUBIC SPLINE

C     WITH KNOTS PRESCRIBED BY THE USER.  CUBIC SPLINE

C     INTERPOLATION CAN ALSO BE CARRIED OUT.

C

C     COX-DE BOOR METHOD FOR EVALUATING B-SPLINES WITH

C     ADAPTATION OF GENTLEMAN*S PLANE ROTATION SCHEME FOR

C     SOLVING OVER-DETERMINED LINEAR SYSTEMS.

C

C     USES NAG LIBRARY ROUTINE  P01AAE.

C

C     STARTED - 1973.

C     COMPLETED - 1976.

C     AUTHOR - MGC AND JGH.

C

C     REDESIGNED TO USE CLASSICAL GIVENS ROTATIONS IN

C     ORDER TO AVOID THE OCCASIONAL UNDERFLOW (AND HENCE

C     OVERFLOW) PROBLEMS EXPERIENCED BY GENTLEMAN*S 3-

C     MULTIPLICATION PLANE ROTATION SCHEME

C

C     WORK1  AND  WORK2  ARE WORKSPACE AREAS.

C     WORK1(R)  CONTAINS THE VALUE OF THE  R TH  DISTINCT DATA

C     ABSCISSA AND, SUBSEQUENTLY, FOR  R = 1, 2, 3, 4,  THE

C     VALUES OF THE NON-ZERO B-SPLINES FOR EACH SUCCESSIVE

C     ABSCISSA VALUE.

C     WORK2(L, J)  CONTAINS, FOR  L = 1, 2, 3, 4,  THE VALUE OF

C     THE  J TH  ELEMENT IN THE  L TH  DIAGONAL OF THE

C     UPPER TRIANGULAR MATRIX OF BANDWIDTH  4  IN THE

C     TRIANGULAR SYSTEM DEFINING THE B-SPLINE COEFFICIENTS.

C

C     .. SCALAR ARGUMENTS ..

      REAL*8 ss

      INTEGER ifail, m, ncap7

C     .. ARRAY ARGUMENTS ..

      REAL*8 c(ncap7), k(ncap7), w(m), work1(m), work2(4,ncap7),

     * x(m), y(m)

C     ..

C     .. LOCAL SCALARS ..

c      DOUBLE PRECISION SRNAME

      CHARACTER*8 srname

      REAL*8 acol, arow, ccol, cosine, crow, d4, d5, d6, d7, d8, d9,

     * d, dprime, e2, e3, e4, e5, k0, k1, k2, k3, k4, k5, k6, n1,

     * n2, n3, relemt, s, sigma, sine, wi, xi

      INTEGER i, ierror, iplusj, iu, j, jold, jplusl, jrev, l4,

     * l, lplus1, lplusu, ncap3, ncap, ncapm1, r

C     .. FUNCTION REFERENCES ..

      REAL*8 sqrt

      INTEGER p01aae

C     ..

      DATA srname /'  E02BAE'/

      ierror = 4

C     CHECK THAT THE VALUES OF  M  AND  NCAP7  ARE REASONABLE

      IF (ncap7.LT.8 .OR. m.LT.ncap7-4) go to 420

      ncap = ncap7 - 7

      ncapm1 = ncap - 1

      ncap3 = ncap + 3

C

C     IN ORDER TO DEFINE THE FULL B-SPLINE BASIS, AUGMENT THE

C     PRESCRIBED INTERIOR KNOTS BY KNOTS OF MULTIPLICITY FOUR

C     AT EACH END OF THE DATA RANGE.

C

      DO 20 j=1,4

         i = ncap3 + j

         k(j) = x(1)

         k(i) = x(m)

   20 CONTINUE

C

C     TEST THE VALIDITY OF THE DATA.

C

C     CHECK THAT THE KNOTS ARE ORDERED AND ARE INTERIOR

C     TO THE DATA INTERVAL.

C

      ierror = 1

      IF (k(5).LE.x(1) .OR. k(ncap3).GE.x(m)) go to 420

      DO 40 j=4,ncap3

         IF (k(j).GT.k(j+1)) go to 420

   40 CONTINUE

C

C     CHECK THAT THE WEIGHTS ARE STRICTLY POSITIVE.

C

      ierror = 2

      DO 60 i=1,m

         IF (w(i).LE.0.0) go to 420

   60 CONTINUE

C

C     CHECK THAT THE DATA ABSCISSAE ARE ORDERED, THEN FORM THE

C     ARRAY  WORK1  FROM THE ARRAY  X.  THE ARRAY  WORK1  CONTAINS

C     THE

C     SET OF DISTINCT DATA ABSCISSAE.

C

      ierror = 3

      work1(1) = x(1)

      j = 2

      DO 80 i=2,m

         IF (x(i).LT.work1(j-1)) go to 420

         IF (x(i).EQ.work1(j-1)) go to 80

         work1(j) = x(i)

         j = j + 1

   80 CONTINUE

      r = j - 1

C

C     CHECK THAT THERE ARE SUFFICIENT DISTINCT DATA ABSCISSAE FOR

C     THE PRESCRIBED NUMBER OF KNOTS.

C

      ierror = 4

      IF (r.LT.ncap3) go to 420

C

C     CHECK THE FIRST  S  AND THE LAST  S  SCHOENBERG-WHITNEY

C     CONDITIONS ( S = MIN(NCAP - 1, 4) ).

C

      ierror = 5

      DO 100 j=1,4

         IF (j.GE.ncap) go to 160

         i = ncap3 - j + 1

         l = r - j + 1

         IF (work1(j).GE.k(j+4) .OR. k(i).GE.work1(l)) go to 420

  100 CONTINUE

C

C     CHECK ALL THE REMAINING SCHOENBERG-WHITNEY CONDITIONS.

C

      IF (ncap.LE.5) go to 160

      r = r - 4

      i = 4

      DO 140 j=5,ncapm1

         k0 = k(j+4)

         k4 = k(j)

  120    i = i + 1

         IF (work1(i).LE.k4) go to 120

         IF (i.GT.r .OR. work1(i).GE.k0) go to 420

  140 CONTINUE

C

C     INITIALISE A BAND TRIANGULAR SYSTEM (I.E. A

C     MATRIX AND A RIGHT HAND SIDE) TO ZERO. THE

C     PROCESSING OF EACH DATA POINT IN TURN RESULTS

C     IN AN UPDATING OF THIS SYSTEM. THE SUBSEQUENT

C     SOLUTION OF THE RESULTING BAND TRIANGULAR SYSTEM

C     YIELDS THE COEFFICIENTS OF THE B-SPLINES.

C

  160 DO 200 i=1,ncap3

         DO 180 l=1,4

            work2(l,i) = 0.0

  180    CONTINUE

         c(i) = 0.0

  200 CONTINUE

      sigma = 0.0

      j = 0

      jold = 0

      DO 340 i=1,m

C

C     FOR THE DATA POINT  (X(I), Y(I))  DETERMINE AN INTERVAL

C     K(J + 3) .LE. X .LT. K(J + 4)  CONTAINING  X(I).  (IN THE

C     CASE  J + 4 .EQ. NCAP  THE SECOND EQUALITY IS RELAXED TO

C     INCLUDE

C     EQUALITY).

C

         wi = w(i)

         xi = x(i)

  220    IF (xi.LT.k(j+4) .OR. j.GT.ncapm1) go to 240

         j = j + 1

         go to 220

  240    IF (j.EQ.jold) go to 260

C

C     SET CERTAIN CONSTANTS RELATING TO THE INTERVAL

C     K(J + 3) .LE. X .LE. K(J + 4).

C

         k1 = k(j+1)

         k2 = k(j+2)

         k3 = k(j+3)

         k4 = k(j+4)

         k5 = k(j+5)

         k6 = k(j+6)

         d4 = 1.0/(k4-k1)

         d5 = 1.0/(k5-k2)

         d6 = 1.0/(k6-k3)

         d7 = 1.0/(k4-k2)

         d8 = 1.0/(k5-k3)

         d9 = 1.0/(k4-k3)

         jold = j

C

C     COMPUTE AND STORE IN  WORK1(L) (L = 1, 2, 3, 4)  THE VALUES

C     OF

C     THE FOUR NORMALIZED CUBIC B-SPLINES WHICH ARE NON-ZERO AT

C     X=X(I).

C

  260    e5 = k5 - xi

         e4 = k4 - xi

         e3 = xi - k3

         e2 = xi - k2

         n1 = wi*d9

         n2 = e3*n1*d8

         n1 = e4*n1*d7

         n3 = e3*n2*d6

         n2 = (e2*n1+e5*n2)*d5

         n1 = e4*n1*d4

         work1(4) = e3*n3

         work1(3) = e2*n2 + (k6-xi)*n3

         work1(2) = (xi-k1)*n1 + e5*n2

         work1(1) = e4*n1

         crow = y(i)*wi

C

C     ROTATE THIS ROW INTO THE BAND TRIANGULAR SYSTEM USING PLANE

C     ROTATIONS.

C

         DO 320 lplus1=1,4

            l = lplus1 - 1

            relemt = work1(lplus1)

            IF (relemt.EQ.0.0) go to 320

            jplusl = j + l

            l4 = 4 - l

            d = work2(1,jplusl)

            IF (abs(relemt).GE.d) dprime =

     *       abs(relemt)*sqrt(1.0+(d/relemt)**2)

            IF (abs(relemt).LT.d) dprime = d*sqrt(1.0+(relemt/d)**2)

            work2(1,jplusl) = dprime

            cosine = d/dprime

            sine = relemt/dprime

            IF (l4.LT.2) go to 300

            DO 280 iu=2,l4

               lplusu = l + iu

               acol = work2(iu,jplusl)

               arow = work1(lplusu)

               work2(iu,jplusl) = cosine*acol + sine*arow

               work1(lplusu) = cosine*arow - sine*acol

  280       CONTINUE

  300       ccol = c(jplusl)

            c(jplusl) = cosine*ccol + sine*crow

            crow = cosine*crow - sine*ccol

  320    CONTINUE

         sigma = sigma + crow**2

  340 CONTINUE

      ss = sigma

C

C     SOLVE THE BAND TRIANGULAR SYSTEM FOR THE B-SPLINE

C     COEFFICIENTS. IF A DIAGONAL ELEMENT IS ZERO, AND HENCE

C     THE TRIANGULAR SYSTEM IS SINGULAR, THE IMPLICATION IS

C     THAT THE SCHOENBERG-WHITNEY CONDITIONS ARE ONLY JUST

C     SATISFIED. THUS IT IS APPROPRIATE TO EXIT IN THIS

C     CASE WITH THE SAME VALUE  (IFAIL=5)  OF THE ERROR

C     INDICATOR.

C

      l = -1

      DO 400 jrev=1,ncap3

         j = ncap3 - jrev + 1

         d = work2(1,j)

         IF (d.EQ.0.0) go to 420

         IF (l.LT.3) l = l + 1

         s = c(j)

         IF (l.EQ.0) go to 380

         DO 360 i=1,l

            iplusj = i + j

            s = s - work2(i+1,j)*c(iplusj)

  360    CONTINUE

  380    c(j) = s/d

  400 CONTINUE

      ierror = 0

  420 IF (ierror) 440, 460, 440

  440 ifail = p01aae(ifail,ierror,srname)

      RETURN

  460 ifail = 0

      RETURN

      END