C      ALGORITHM 742, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 21, NO. 1, MARCH, 1995, P. 98-110.
C
c*** spsrc.f
C......................................................................
C     THIS SUBROUTINE INSERTS A KNOT, NAMELY THE DATA POINT X(INKNOT),
C      INTO THE SPLINE KNOT SEQUENCE OF THE THE LINEAR SPLINE
C      APPROXIMATION SX(.) DEFINED IN SUBROUTINE L2CXFT.
C
C     CALLS SUBROUTINE PLUS.
C......................................................................
      SUBROUTINE ADDKNT(INKNOT,E,B,D,X,F,AKN,IAKN,IACT,NK,NACT,I1,N)
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E(I1-1:N),B(I1-1:N),D(I1-1:N),X(I1:N),F(I1:N),
     1  AKN(I1-2:N),IAKN(I1-1:N),IACT(I1:N)
C
C.... I N P U T ....
C     INKNOT  INTEGER VARIABLE, SUCH THAT X(INKNOT) IS INSERTED
C             INTO THE KNOT SEQUENCE AKN(.).
C     E,B,D,X,F,AKN,IAKN,IACT,I1,N,NK,NACT  AS THEY ARE
C             DEFINED IN SUBROUTINE L2CXFT.
C
C.... O U T P U T ....
C     E,B,D,AKN,IAKN,IACT  AFTER INSERTING X(INKNOT) INTO AKN(.).
C     NK      INCREASED BY ONE.
C     NACT    DECREASED BY ONE.
C
C.... M E T H O D ....
C     FIND POSITION J IN THE KNOT SEQUENCE AKN(.) SUCH THAT
C      AKN(J-1) < X(INKNOT) < AKN(J) BEFORE INSERTION AND
C      AKN(J)=X(INKNOT) AFTER INSERTION.
C     CALL SUBROUTINE PLUS TO MODIFY THE NORMAL EQUATIONS THAT ARE
C      AFFECTED BY THE KNOT INSERTION. HENCE THE 'HAT' BASIS FUNCTIONS
C      AFFECTED BY THE INSERTION HAVE TO DETERMINED. THE J-TH BASIS
C      FUNCTION IS DEFINED FROM X(IP) TO X(IR), AND IS CENTERED AT
C      X(IQ). ALSO THE (J-1)-TH BASIS FUNCTIONS IS CENTERED AT X(IP) AND
C      IS DEFINED FROM X(IPL) TO X(IQ). THE (J+1)-TH BASIS FUNCTION
C      IS CENTERED AT X(IR) AND IS DEFINED FROM X(IQ) TO X(IRR).
C     REVISE IACT(.), THE INDICES OF THE ACTIVE CONSTRAINTS.
C......................................................................
      I0=I1-1
C
C     GIVEN X(INKNOT), APPLY BISECTION IN [I0,NK) IN ORDER TO FIND IP
C      SUCH THAT  AKN(IP) < X(INKNOT) < AKN(IP+1).
C
      IP=I0
      IR=NK
20    IF (IR-IP.EQ.1) GO TO 30
      IQ=(IP+IR)/2
      IF (INKNOT .LT. IAKN(IQ)) THEN
          IR=IQ
      ELSE
          IP=IQ
      END IF
      GO TO 20
30    CONTINUE
C
      DO 100 I=NK,IP+1,-1
         B(I+1)=B(I)
         E(I+1)=E(I)
         D(I+1)=D(I)
         AKN(I+1)=AKN(I)
100      IAKN(I+1)=IAKN(I)
      J=IP+1
      E(J)=0.
      B(J)=0.
      D(J)=0.
      NK=NK+1
      AKN(J)=X(INKNOT)
      IAKN(J)=INKNOT
C
C     DEFINE SUPPORTS OF NEW BASIS FUNCTIONS THAT FOLLOW THE KNOT
C      INSERTION.
C
      IF (J.EQ.I1) THEN
          IPL=I1
      ELSE
          IPL=IAKN(J-2)
      END IF
C
      IF (J.EQ.NK-1) THEN
          IRR=N
      ELSE
          IRR=IAKN(J+2)
      END IF
C
      IP=IAKN(J-1)
      IQ=IAKN(J)
      IR=IAKN(J+1)
C
C     MODIFY NORMAL EQUATIONS.
C
      CALL PLUS(X,F,AKN,E,B,D,NK,I1,N,IPL,IP,IQ,IR,IRR,J)
C
C     REVISE ACTIVE CONSTRAINTS.
C      NACT>=I1 BECAUSE THE KNOT INSERTION SHOWS THAT THERE EXISTS
C      AT LEAST ONE ACTIVE CONSTRAINT.
C
      IF (NACT.EQ.I1) GO TO 140
C
C     I1<INKNOT<N, K=INKNOT-1=ACTIVE CONSTRAINT TO BE DELETED.
C     GIVEN K, APPLY BISECTION IN [I1,NACT] IN ORDER TO FIND
C      IQ SUCH THAT IACT(IQ).EQ.K .
C
      K=INKNOT-1
      IQ=I1
      IF (K.EQ.IACT(I1)) GO TO 120
      IF (K.EQ.IACT(NACT)) GO TO 140
      IP=I1
      IR=NACT
110   IQ=(IP+IR)/2
      IF (IACT(IQ).EQ.K) GO TO 120
      IF (IACT(IQ).LT.K) THEN
          IP=IQ
      ELSE
          IR=IQ
      END IF
      GO TO 110
120   CONTINUE
      DO 130 I=IQ,NACT-1
130      IACT(I)=IACT(I+1)
140   IACT(NACT)=0
      NACT=NACT-1
C
200   RETURN
      END
C......................................................................
C     THIS SUBROUTINE CALCULATES
C      1) THE COEFFICIENTS COEF(.) OF A LINEAR SPLINE SX(.) OVER
C         KNOTS AKN(.), FROM THE NORMAL EQUATIONS.
C      2) THE SPLINE COMPONENTS SX(.) FROM COEF(.) AND
C      3) THE ASSOCIATED LAGRANGE MULTIPLIERS PAR(.).
C     ALSO THIS SUBROUTINE PROVIDES ITERATIVE REFINEMENT OF
C      1) THE COEFFICIENTS COEF(.) AND
C      2) THE LAGRANGE MULTIPLIERS PAR(.).
C
C     CALLS SUBROUTINES COEFF, SXATXI AND LAGRNG.
C......................................................................
      SUBROUTINE CFLGRF(SX,PAR,COEF,IACT,E,B,D,X,AKN,IAKN,F,DC,EC,
     1  UE,RES,NK,NACT,I1,N,ITREF,ITERTN)
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION SX(I1:N),COEF(I1-1:N),E(I1-1:N),B(I1-1:N),D(I1-1:N),
     1  AKN(I1-2:N),IAKN(I1-1:N),X(I1:N),F(I1:N),DC(I1:N),EC(I1:N),
     2  IACT(I1:N),UE(I1-1:N),PAR(I1:N),RES(I1-1:N)
C
C.... I N P U T ....
C     E,B,D,X,F,AKN,IAKN,IACT,DC,EC,I1,N,NK AND NACT AS THEY ARE
C              DEFINED IN SUBROUTINE L2CXFT.
C     ITERTN  INTEGER VARIABLE THAT IS USED TO CHECK THE FIRST CALL
C              OF THIS SUBROUTINE BY SUBROUTINE L2CXFT.
C     ITREF   INTEGER VARIABLE, IF ITREF=1 THEN ITERATIVE REFINEMENT
C              IS INITIATED.
C
C.... O U T P U T ....
C     COEF(I1-1:NK) REAL ARRAY CONTAINING THE SPLINE COEFFICIENTS
C             AT AKN(I), I=I1-1,I1,...,NK.
C     SX(I1:N) REAL ARRAY CONTAINING THE SPLINE COMPONENTS AT X(I),
C             I=I1,I1+1,...,N.
C     PAR(I1:N-2) REAL ARRAY CONTAINING THE LAGRANGE MULTIPLIERS SO
C             THAT PAR(IACT(I)), I=I1,I1+1,...,NACT ARE THE ACTIVE
C             CONSTRAINT MULTIPLIERS WHILE THE REMAINING MULTIPLIERS
C             EQUAL TO ZERO.
C
C.... W O R K I N G  S P A C E ....
C     UE(I1-1:NK+1) REAL ARRAY THAT PROVIDES WORKING SPACE FOR THE
C             FACTORIZATIONS THAT GIVE COEF(.) AND PAR(.).
C     RES(I1-1:N) REAL ARRAY THAT PROVIDES WORKING SPACE FOR THE
C             RESIDUALS OF THE NORMAL EQUATIONS AT COEF(.) AND FOR
C             THE CORRECTIONS OF COEF(.). IT ALSO PROVIDES WORKING
C             SPACE FOR THE RESIDUALS OF THE LAGRANGE EQUATIONS AT
C             PAR(.) AND FOR THE CORRECTIONS OF PAR(.).
C
C.... M E T H O D ....
C     FOLLOWS THE STEPS DESCRIBED IN THE HEADING COMMENT.
C......................................................................
      I0=I1-1
C
C     IF ITERTN=0 THEN AVOID CALCULATION OF COEF(.), SINCE COEF(.)
C      ARE SUPPLIED BY SUBROUTINE L2ILS.
C
      IF (ITERTN.EQ.0) GO TO 50
C
C.... DERIVE THE SPLINE COEF(.) ....
C
      CALL COEFF(COEF,E,B,D,UE,NK,I1,N)
C
C.... INITIATE ITERATIVE REFINEMENT IF ITREF=1 ....
C
      IF (ITREF.NE.1) GO TO 50
C
C         REFINE COEF(.).
C         SET RESIDUALS OF NORMAL EQUATIONS IN RES(.).
C
          RES(I0)=B(I0)-(E(I0)*COEF(I0)+D(I1)*COEF(I1))
          IF (I1 .EQ. NK) GO TO 30
          DO 20 I=I1,NK-1
20           RES(I)=B(I)-(D(I)*COEF(I-1)+E(I)*COEF(I)+D(I+1)*COEF(I+1))
30        RES(NK)=B(NK)-(D(NK)*COEF(NK-1)+E(NK)*COEF(NK))
C
C         SOLVE FOR RES(.).
C
          CALL COEFF(RES,E,RES,D,UE,NK,I1,N)
C
C         CORRECT COEF(.).
C
          DO 40 I=I0,NK
40           COEF(I)=COEF(I)+RES(I)
C
C.... GENERATE SPLINE SX(.) AT X(.) FROM COEF(.) OVER AKN(.) ....
C
50    CALL SXATXI(SX,X,COEF,AKN,IAKN,I1,N,NK)
C
C.... CALCULATE THE LAGRANGE MULTIPLIERS PAR(.) ....
C
      IF (NACT.LE.I0) GO TO 270
      CALL LAGRNG(PAR,SX,F,IACT,NACT,EC,DC,UE,I1,N)
C
C.... INITIATE ITERATIVE REFINEMENT IF ITREF=1 ....
C
      IF (ITREF.NE.1) GO TO 270
C
C         REFINE PAR(.).
C         SET RESIDUALS OF LAGRANGIAN EQUATIONS IN RES(.).
C
          IF (I1 .EQ. NACT) THEN
              II=IACT(I1)
              RES(II)=SX(II+1)-F(II+1) - EC(II)*PAR(II)
              GO TO 200
          END IF
C
          II=IACT(I1)
          RES(II)=SX(II+1)-F(II+1)-(EC(II)*PAR(II)+DC(II+1)*PAR(II+1))
          IF (I1+1 .EQ. NACT) GO TO 190
          DO 180 I=I1+1,NACT-1
             II=IACT(I)
180          RES(II)=SX(II+1)-F(II+1) -
     1          (DC(II)*PAR(II-1)+EC(II)*PAR(II)+DC(II+1)*PAR(II+1))
190       II=IACT(NACT)
          RES(II)=SX(II+1)-F(II+1)-(DC(II)*PAR(II-1)+EC(II)*PAR(II))
200       CONTINUE
C
C         SOLVE FOR RES(.) AND CORRECT PAR(.).
C          THE FOLLOWING CODE DUPLICATES WITH THAT IN SUBROUTINE LAGRNG.
C
          I=I1
210       K=1
220       IF (I.LT.NACT .AND. IACT(I).EQ.IACT(I+1)-1) THEN
              I=I+1
              K=K+1
              GO TO 220
          END IF
C
230       IF (K.EQ.1) THEN
              II=IACT(I)
              PAR(II)=PAR(II)+RES(II)
          ELSE
              UE(I-K+1)=EC(IACT(I-K+1))
              DO 240 L=I-K+2,I
                 DI=DC(IACT(L))
240              UE(L)=EC(IACT(L))-DI/UE(L-1)*DI
C
              II=IACT(I-K+1)
              RES(II)=RES(II)
              DO 250 L=I-K+2,I
                 II=IACT(L)
250              RES(II)=RES(II)-DC(II)*RES(IACT(L-1))/UE(L-1)
C
              RES(IACT(I))=RES(IACT(I))/UE(I)
              PAR(IACT(I))=PAR(IACT(I))+RES(IACT(I))
              DO 260 L=I-1,I-K+1,-1
                 RES(IACT(L))=RES(IACT(L))/UE(L)-DC(IACT(L+1))*
     1                     RES(IACT(L+1))/UE(L)
260              PAR(IACT(L))=PAR(IACT(L))+RES(IACT(L))

          END IF
          I=I+1
          IF (I-NACT) 210,210,270
C
270   RETURN
      END
C......................................................................
C     THIS SUBROUTINE CALCULATES THE COEFFICIENTS COEF(.) OF
C      A LEAST SQUARES LINEAR SPLINE FROM THE NORMAL EQUATIONS.
C......................................................................
      SUBROUTINE COEFF(COEF,E,B,D,UE,NK,I1,N)
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION COEF(I1-1:N),E(I1-1:N),B(I1-1:N),D(I1-1:N),UE(I1-1:N)
C
C.... I N P U T ....
C     E,B,D,I1,N,NK   AS THEY ARE DEFINED IN SUBROUTINE L2CXFT.
C
C.... O U T P U T ....
C     COEF(I1-1:NK) REAL ARRAY CONTAINING THE SPLINE COEFFICIENTS
C             AT AKN(I), I=I1-1,I1,...,NK.
C
C.... W O R K I N G  S P A C E ....
C     UE(I1-1:NK+1) REAL ARRAY THAT PROVIDES WORKING SPACE FOR THE
C             FACTORIZATION THAT GIVES COEF(.). IT KEEPS THE
C             THE DIAGONAL ELEMENTS OF THE FORWARD SWEEP OF
C             THE FACTORIZATION.
C
C.... M E T H O D (CALCULATE COEFFICIENTS) ....
C     EXPRESS SPLINE SX(.) AS A LINEAR COMBINATION OF B-SPLINES
C             SX(I)=SUM(COEF(J)*N(J,I), J=I1-1,I1,...,NK)
C                                          FOR I=I1,I1+1,...,N
C     WHERE N(J,I)=THE I-TH COMPONENT OF THE J-TH LINEAR B-SPLINE.
C     DENOTE BY <.,.> THE INNER PRODUCT OF TWO VECTORS AND,
C     THE SPLINE COEFFICIENTS COEF(.) ARE DETERMINED BY SOLVING THE
C     NORMAL EQUATIONS
C             D(J)*COEF(J-1)+E(J)*COEF(J)+D(J+1)*COEF(J+1)=B(J)
C                                          FOR J=I1-1,I1,...,NK
C     WHERE
C             D(J)=<N(J,.),N(J-1,.)>, J=I1-1,I1,...,NK
C             E(J)=<N(J,.),N(J,.)>, J=I1-1,I1,...,NK
C             B(J)=<N(J,.),F(.)>, J=I1-1,I1,...,NK.
C......................................................................
      I0=I1-1
C
C.....DERIVE COEF(.) BY FACTORIZATION
C
      UE(I0)=E(I0)
      COEF(I0)=B(I0)
      DO 10 I=I0+1,NK
         UE(I)=E(I)-D(I)/UE(I-1)*D(I)
10       COEF(I)=B(I)-D(I)*COEF(I-1)/UE(I-1)
      COEF(NK)=COEF(NK)/UE(NK)
      DO 20 I=NK-1,I0,-1
20       COEF(I)=(COEF(I)-D(I+1)*COEF(I+1))/UE(I)
C
      RETURN
      END
C.......................................................................
C     THIS SUBROUTINE DELETES THE KNOT AKN(KNOUT) OF THE LINEAR
C      SPLINE APPROXIMATION SX(.) DEFINED IN SUBROUTINE L2CXFT.
C
C     CALLS SUBROUTINE MINUS1.
C.......................................................................
      SUBROUTINE DELKNT(E,B,D,AKN,IAKN,IACT,NK,NACT,I1,N,KNOUT)
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E(I1-1:N),B(I1-1:N),D(I1-1:N),AKN(I1-2:N),IAKN(I1-1:N),
     1  IACT(I1:N)
C
C.... I N P U T ....
C     KNOUT  INTEGER VARIABLE, SUCH THAT AKN(KNOUT) IS THE KNOT TO
C             BE DELETED FROM THE KNOT SEQUENCE AKN(.).
C     E,B,D,AKN,IAKN,IACT,I1,N,NK,NACT  AS THEY ARE DEFINED IN
C             SUBROUTINE L2CXFT.
C
C.... O U T P U T ....
C     E,B,D,AKN,IAKN,IACT  AFTER DELETING THE KNOT AKN(KNOUT).
C     NK      DECREASED BY ONE.
C     NACT    INCREASED BY ONE.
C
C.... M E T H O D ....
C     SINCE AKN(J) IS DELETED THE J-TH NORMAL EQUATION IS DELETED TOO.
C      THEN CALL SUBROUTINE MINUS1 TO MODIFY THE (J-1) AND (J+1) NORMAL
C      EQUATIONS. SHIFT E,B,D,AKN AND IAKN.
C     REVISE IACT(.), THE INDICES OF THE ACTIVE CONSTRAINTS.
C.......................................................................
      I0=I1-1
      J=IAKN(KNOUT)-1
C
C     MODIFY NORMAL EQUATIONS.
C
      CALL MINUS1(E,B,D,AKN,I1,N,KNOUT)
C
C     SHIFT ARRAYS.
C
      DO 5 I=KNOUT,NK-1
         AKN(I)=AKN(I+1)
         IAKN(I)=IAKN(I+1)
         E(I)=E(I+1)
         B(I)=B(I+1)
5        D(I)=D(I+1)
C
C     SET LAST POSITION TO ZERO.
C
      AKN(NK)=0.
      IAKN(NK)=0
      E(NK)=0.
      B(NK)=0.
      D(NK)=0.
      NK=NK-1
C
C     UPDATE IACT(.).
C
C     GIVEN J=IAKN(KNOUT)-1, APPLY BISECTION IN [I1,NACT] IN ORDER
C      TO FIND IR SUCH THAT  IACT(IR-1) < J < IACT(IR) , OR
C      IQ IS OUTSIDE [I1,NACT].
C
      IF (NACT.EQ.I0) THEN
          IR=I1
          GO TO 130
      END IF
      IF (J.LT.IACT(I1)) THEN
          IR=I1
          GO TO 110
      END IF
      IF (J.GT.IACT(NACT)) THEN
          IR=NACT+1
          GO TO 130
      END IF
C
C     APPLY BISECTION.
C
      IP=I1
      IR=NACT
100   IQ=(IP+IR)/2
      IF (IR-IP.EQ.1) GO TO 110
      IF (J.LT.IACT(IQ)) THEN
          IR=IQ
      ELSE
          IP=IQ
      END IF
      GO TO 100
110   CONTINUE
      DO 120 I=NACT,IR,-1
120      IACT(I+1)=IACT(I)
130   IACT(IR)=J
      NACT=NACT+1
C
      RETURN
      END
C......................................................................
C     L2CXFT: A FORTRAN SUBROUTINE FOR LEAST SQUARES DATA FITTING
C             WITH NON-NEGATIVE SECOND DIVIDED DIFFERENCES.
C     BY I.C. DEMETRIOU, DEPARTMENT OF ECONOMICS, UNIVERSITY OF ATHENS,
C             8 PESMAZOGLOU STREET, ATHENS 105 59, GREECE.
C             TEL   : (GR) 01-3242438, 8017732
C             E-MAIL: DEMETRI@AUEB.ARIADNE-T.GR
C......................................................................
C     THE USER GIVES THE NOISY FUNCTION MEASUREMENTS F(I): I=1,2,...N
C      AT THE ABSCISSAE X(1)<X(2)<...<X(N),
C     AND THIS SUBROUTINE CALCULATES THE N-TH COMPONENT VECTOR SX(.)
C
C           THAT MINIMIZES THE SUM OF THE SQUARES OF THE ERRORS
C            SX(I)-F(I): I=1,2,...,N
C           SUBJECT TO NON-NEGATIVITY OF SECOND DIVIDED DIFFERENCES
C            OF THE SMOOTHED VALUES SX(.).
C
C     IN THIS SENSE, SX(.) IS
C                   THE BEST CONVEX FIT TO F(.)
C                 OR THE CONVEX REGRESSION ON F(.)
C
C     THERE ARE EXACTLY N-2 CONSTRAINTS THAT GIVE A HIGHLY STRUCTURED
C      CALCULATION WHICH IS SOLVED BY A SPECIAL DUAL FEASIBLE QUADRATIC
C      PROGRAMMING ALGORITHM INITIALIZED BY A O(N) STARTING PROCEDURE.
C
C     THE CONCAVE FITTING CAN BE TREATED AS CONVEX AFTER AN OVERALL
C      CHANGE OF SIGN.
C
C     REF:  I.C.DEMETRIOU AND M.J.D. POWELL, "THE MINIMUM SUM OF
C           SQUARES CHANGE TO UNIVARIATE DATA THAT GIVES CONVEXITY",
C           IMA J. OF NUMERICAL ANALYSIS, 11(1991), PP. 433-448.
C
C     APPLICATIONS: DATA SMOOTHING, STATISTICAL ANALYSIS, ECONOMICS.
C
C     CALLS SUBROUTINES ADDKNT, CFLGRF, DELKNT, L2ILS AND MESSGE.
C.......................................................................
      SUBROUTINE L2CXFT(X,F,SX,I1,N,AKN,IAKN,COEF,NK,IACT,NACT,PAR,
     1  ITER,SXOLD,E,B,D,Z,U,DC,EC,PAROLD,RES,IPRINT,MODE)
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(I1:N),F(I1:N),SX(I1:N),COEF(I1-1:N),AKN(I1-2:N),
     1  IAKN(I1-1:N),IACT(I1:N),PAR(I1:N),DC(I1:N),EC(I1:N),
     2  E(I1-1:N),B(I1-1:N),D(I1-1:N),Z(I1-1:N),U(I1-1:N),
     3  SXOLD(I1:N),PAROLD(I1:N),RES(I1-1:N),ITER(*)
C
C..... I N P U T .....
C      I1      INTEGER VARIABLE, LOWER DATA INDEX.
C      N       INTEGER VARIABLE, UPPER DATA INDEX.
C      X(I1:N) REAL ARRAY OF ABSCISSAE.
C      F(I1:N) REAL ARRAY OF FUNCTION MEASUREMENTS.
C      IPRINT  INTEGER VARIABLE. IF IPRINT=1 THEN SUBROUTINE MESSGE
C               PRODUCES OUTPUT.
C
C .....O U T P U T .....
C      SX(I1:N) REAL ARRAY CONTAINING THE SPLINE COMPONENTS AT X(I),
C                I=I1,I1+1,...,N IE THE BEST CONVEX FIT TO F(.).
C      NK      INTEGER VARIABLE, UPPER INDEX OF SPLINE KNOTS SEQUENCE.
C      NACT    INTEGER VARIABLE, SUCH THAT EITHER NACT=I1-1 OR
C               (NACT-I1+1) IS EQUAL TO THE NUMBER OF ACTIVE CONSTRAINTS
C               IN THE SOLUTION.
C      AKN(I1-2:NK)  REAL ARRAY CONTAINING THE SPLINE KNOTS AND THE
C               DUMMY COMPONENT AKN(I1-2)=X(I1).
C      IAKN(I1-1:NK)  INTEGER ARRAY CONTAINING THE DATA POINT INDICES
C               OF THE ASSOCIATED KNOTS.
C      COEF(I1-1:NK)  REAL ARRAY CONTAINING THE SPLINE COEFFICIENTS
C               AT AKN(I), I=I1-1,I1,...,NK.
C      IACT(I1:NACT)  INTEGER ARRAY CONTAINING THE INDICES OF THE
C              ACTIVE CONSTRAINTS AT THE SOLUTION.
C      PAR(I1:N-2)  INTEGER ARRAY CONTAINING IN PAR(IACT(.))
C               THE LAGRANGE MULTIPLIERS OF THE ACTIVE CONSTRAINTS.
C      ITER(*) INTEGER ARRAY WHOSE ELEMENTS ARE AS FOLLOWS:
C              ITER(1)  KEEPS THE NUMBER OF ITERATIONS AT THE END
C                        OF THE CALCULATION. IF ITER(1) EXCEEDS
C                        THE NUMBER 3*(N-I1+1) THEN TERMINATION OCCURS
C                        WITH THE MESSAGE
C                         'L2CXFT QUITS ON ITERATION COUNT'
C              ITER(2)  KEEPS THE NUMBER OF ACTIVE CONSTRAINTS AT THE
C                        END OF THE CALL OF SUBROUTINE L2ILS.
C              ITER(3)  KEEPS THE NUMBER OF CALLS OF SUBROUTINE
C                        CFLGRF.
C              ITER(4)  KEEPS THE NUMBER OF ACTIVE CONSTRAINTS WHEN
C                        FIRST FEASIBLE LAGRANGIAN VECTOR IS FOUND
C                        JUST AFTER THE CALL OF SUBROUTINE L2ILS.
C              ITER(5)  KEEPS THE NUMBER OF CONSTRAINT DELETIONS FROM
C                        THE ACTIVE SET.
C              ITER(6)  KEEPS THE NUMBER OF CONSTRAINT ADDITIONS TO
C                        THE ACTIVE SET AFTER FIRST FEASIBLE PAR(.).
C      MODE    INTEGER VARIABLE THAT IS SET AUTOMATICALLY BY THE
C               SUBROUTINE TO INDICATE THE TERMINATION STATUS.
C              1 = SUCCESSFUL TERMINATION
C              2 = TERMINATION DUE TO NON-INCREASE OF OBJECTIVE FUNCTION
C              3 = NUMBER OF 3*(N-I1+1) ITERATIONS HAS EXCEEDED
C              4 = ERROR RETURN BECAUSE (N .LT. I1).
C
C.... M E S S A G E S ....
C     INCLUDED IN SUBROUTINE MESSGE AND ACTIVATED WHEN INPUT PARAMETER
C      IPRINT IS SET BY THE USER TO ONE.
C
C.... W O R K I N G   S P A C E .....
C      E(I1-1:NK)  REAL ARRAY THAT HOLDS THE DIAGONAL OF THE NORMAL
C              EQUATIONS ASSOCIATED WITH SPLINE SX(.).
C      D(I1-1:NK)  REAL ARRAY THAT HOLDS THE SUBDIAGONAL OF THE
C              NORMAL EQUATIONS ASSOCIATED WITH SX(.).
C      B(I1-1:NK)  REAL ARRAY THAT HOLDS THE RIGHT HAND SIDE OF
C              THE NORMAL EQUATIONS.
C      EC(I1:N-2)  REAL ARRAY SUCH THAT EC(I)=(SECOND COMPONENT OF
C              THE I-TH DIVIDED DIFFERENCE AFTER SCALING).
C      DC(I1:N-2)  REAL ARRAY SUCH THAT DC(I+1)=(THIRD COMPONENT OF
C              THE I-TH DIVIDED DIFFERENCE AFTER SCALING).
C      PAROLD(I1:N-2)  REAL ARRAY THAT HOLDS THE LAGRANGE MULTIPLIERS
C              THAT OCCURED AT THE FORMER ITERATION. IT IS
C              PAROLD(IACT(.)).GT.0. .
C      SXOLD(I1:N) REAL ARRAY THAT HOLDS THE SPLINE COMPONENTS
C              THAT OCCURED AT THE FORMER ITERATION.
C      RES(I1-1:NK)  REAL ARRAY THAT IS USED TO KEEP THE RESIDUALS OF
C              THE NORMAL EQUATIONS AT COEF(.) AND THE CORRECTIONS OF
C              COEF(.). IT IS ALSO USED TO KEEP THE RESIDUALS OF THE
C              LANGRANGIAN EQUATIONS AT PAR(.) AND THE CORRECTIONS TO
C              PAR(.).
C      U(I1-1,N) REAL ARRAY PROVIDING SPACE FOR FACTORIZATIONS.
C      Z(I1-1,N) REAL ARRAY PROVIDING SPACE FOR FACTORIZATIONS.
C      ITREF   INTEGER VARIABLE. IF ITREF=1 THEN ITERATIVE REFINEMENT
C               IS INITIATED.
C      ITEST   INTEGER VARIABLE. IF ITEST=1 THEN A RETURN IS CAUSED DUE
C               TO NOT INCREASING OF THE FUNCTION VALUES.
C
C .....M E T H O D (DUAL-FEASIBLE ALGORITHM 3 IN QUOTED REFERENCE).
C      VECTOR SX(.) IS CALCULATED BY THIS SUBROUTINE AS FOLLOWS:
C        STEP 0   CALL SUBROUTINE L2ILS TO DERIVE AN INITIAL
C                  APPROXIMATION TO SX(.).
C        STEP 1   USE THIS APPROXIMATION AS STARTING POINT OF A
C                  SPECIAL QUADRATIC PROGRAMMING (QP) METHOD THAT
C                  COMPLETES THE CALCULATION OF SX(.).
C
C      FEATURES OF THE QP METHOD:
C              A) IT GENERATES INFEASIBLE ESTIMATES OF THE OPTIMUM,
C              B) IT KEEPS A SET OF ACTIVE CONSTRAINTS THAT ARE
C                  ASSOCIATED WITH NON-NEGATIVE LAGRANGE MULTIPLIERS,
C              C) THE OBJECTIVE FUNCTION ACHIEVES A STRICTLY LARGER
C                  VALUE ANY TIME A NEW ITERATION BEGINS.
C              D) POWELL'S SAFEGUARD TERMINATION TEST IS APPLIED
C                  WHEN ACCURACY IS INSUFFICIENT TO MAINTAIN INCREASING
C                  FUNCTION VALUES.
C              E) ITERATIVE REFINEMENT IS CALLED AUTOMATICALLY IN
C                  CONJUNCTION WITH D) IN ORDER TO IMPROVE ACCURACY.
C
C      IMPLEMENTATION OF THE QP METHOD:
C      A) CALCULATION OF AN ESTIMATE SX(.) OF THE OPTIMUM:
C          EXPRESS SPLINE SX(.) AS A LINEAR COMBINATION OF LINEAR
C          B-SPLINES
C              SX(I)=SIGMA(COEF(J)*N(J,I), J=I1-1,I1,...,NK)
C                                           FOR I=I1,I1+1,...,N ,
C          WHERE N(J,I)=THE I-TH COMPONENT OF THE J-TH LINEAR B-SPLINE.
C          DENOTE BY <.,.> THE INNER PRODUCT OF TWO VECTORS AND,
C          THE COEFFICIENTS OF SX(.) ARE DETERMINED BY SOLVING THE
C          THE POSITIVE DEFINITE SYSTEM OF THE NORMAL EQUATIONS
C          EFFICIENTLY AND STABLY BY CHOLESKY FACTORIZATION
C              D(J)*COEF(J-1)+E(J)*COEF(J)+D(J+1)*COEF(J+1)=B(J)
C                                           FOR J=I1-1,I1,...,NK
C          WHERE BY <.,.> IS DENOTED THE SCALAR PRODUCT OF TWO VECTORS
C          AND
C              D(J)=<N(J,.),N(J-1,.)>, J=I1,I1+1,...,NK
C              E(J)=<N(J,.),N(J,.)>, J=I1-1,I1,...,NK
C              B(J)=<N(J,.),F(.)>, J=I1-1,I1,...,NK.
C          (SUBROUTINE CFLGRF).
C
C      B) CALCULATION OF THE LAGRANGE MULTIPLIERS
C          THE LAGRANGE EQUATIONS AT AN ESTIMATE SX(.) ARE
C              2*(SX(II+1)-F(II+1))=SIGMA(PAR(II)*A(II,.))
C                                    FOR II=IACT(I), I=I1,I1+1,...,NACT
C          WHERE A(II,J) IS THE J-TH COMPONENT OF THE NORMAL
C          OF THE II-TH CONSTRAINT AND
C              A(II,J)<>0 ONLY FOR (II.LE.J.LE.II+2).
C          (ARRAY A(.,.) IS NOT ACTUALLY USED IN THIS SUBROUTINE;
C          INSTEAD ARRAYS EC(.) AND DC(.) ARE EMPLOYED.). BECAUSE
C          THIS SYSTEM IS OVERDETERMINED PAR(IACT(.)) IS DERIVED
C          BY CHOOSING FROM THE LAGRANGE EQUATIONS THOSE EQUATIONS THAT
C          CORRESPOND TO THE SECOND NON-ZERO COMPONENT OF A(IACT(.),.);
C          SO A BLOCK DIAGONAL POSITIVE DEFINITE SYSTEM IS OBTAINED THAT
C          IS SOLVED EFFICIENTLY AND STABLY BY CHOLESKY FACTORIZATION.
C          (SUBROUTINE CFLGRF).
C
C      C) PARTICULAR METHODS ARE DEVELOPED FOR INSERTING INTO AND
C          DELETING FROM THE SPLINE APPROXIMATION ONE KNOT ECONOMICALLY.
C          (SUBROUTINES ADDKNT AND DELKNT).
C......................................................................
C     SET CERTAIN INITIAL VALUES.
C
      MODE=1
      I0=I1-1
      ITREF=0
      ITEST=0
      ZERO=0.
      DO 10 I=1,10
10       ITER(I)=0
C
      IF (N.LT.I1) THEN
          MODE=4
          CALL MESSGE(ZERO,0,0,ITER,10)
          GO TO 1000
      END IF
C
C     IF N=I1 OR N=I1+1 THEN THE DATA ITSELF IS THE REQUIRED FIT.
C
      IF (N.EQ.I1) THEN
          SX(I1)=F(I1)
          AKN(I1-1)=X(I1)
          IAKN(I1-1)=I1
          COEF(I1-1)=F(I1)
          CALL MESSGE(ZERO,0,0,ITER,20)
          GO TO 1000
      END IF
C
      IF (N.EQ.I1+1) THEN
          SX(I1)=F(I1)
          SX(N)=F(N)
          AKN(I1-1)=X(I1)
          AKN(I1)=X(N)
          IAKN(I1-1)=I1
          IAKN(I1)=N
          COEF(I1-1)=F(I1)
          COEF(I1)=F(N)
          CALL MESSGE(ZERO,0,0,ITER,20)
          GO TO 1000
      END IF
C
C.... I N I T I A L I Z E .............................................
C
C     THERE ARE N-2 CONSTRAINTS DEFINED BY THE DIVIDED DIFFERENCES.
C      WE ASSOCIATE THE I-TH CONSTRAINT WITH THE I-TH DIVIDED DIFFER-
C      ENCE. SET COEFFICIENTS OF FULL SYSTEM OF EQUATIONS OF LAGRANGE
C      MULTIPLIERS. SINCE GRAD(F)=2*(SX-F), THE COEFFICIENTS ARE
C      SCALED BY 2.
C
      DO 50 I=I1,N-3
         CX2=X(I+2)-X(I+1)
         CX1=X(I+1)-X(I)
         CX3=X(I+2)-X(I)
         CX2=1/CX2
         CX1=1/CX1
         EC(I)=-(CX1+CX2)/2
50       DC(I+1)=CX2/2
      CX2=X(N)-X(N-1)
      CX1=X(N-1)-X(N-2)
      CX3=X(N)-X(N-2)
      CX2=1/CX2
      CX1=1/CX1
      EC(N-2)=-(CX1+CX2)/2
C
      DO 60 I=I1,N-2
60       PAR(I)=0.0
C
C.... CALCULATE STARTING POINT BY L2ILS ...............................
C
      CALL L2ILS(COEF,AKN,IAKN,D,E,B,Z,U,NK,F,X,I1,N)
C
C     INITIALIZE INDICES OF ACTIVE CONSTRAINTS IN IACT(.).
C
      J=I0
      I=I1+1
      IAKN(I0)=I1
      KL=I0
80    IF (IAKN(J).LT.I .AND. I.LT.IAKN(J+1)) THEN
          KL=KL+1
          IACT(KL)=I-1
          I=I+1
          GO TO 80
      END IF
      J=J+1
      I=I+1
      IF (J.LT.NK) GO TO 80
      NACT=KL
      ITER(2)=NACT-I0
C
C.... ITERATE BY ONE, UNTIL FIRST FEASIBLE LAGRANGIAN VECTOR IS FOUND.
C
C     CALCULATE COEF OVER KNOTS, SPLINE OVER ABSCISSAE AND
C      LAGRANGE MULTIPLIERS.
C
90    CALL CFLGRF(SX,PAR,COEF,IACT,E,B,D,X,AKN,IAKN,F,DC,EC,U,
     1  RES,NK,NACT,I1,N,ITREF,ITER(3))
      ITER(3)=ITER(3)+1
C
C     FIND MINIMUM LAGRANGE MULTIPLIER.
C
      IF (NACT .LE. I0) GO TO 150
      AMINLA=PAR(IACT(I1))
      IPAR=IACT(I1)
      IF (I1 .EQ. NACT) GO TO 110
      DO 100 I=I1+1,NACT
         IF (PAR(IACT(I)) .LT. AMINLA) THEN
             AMINLA=PAR(IACT(I))
             IPAR=IACT(I)
         END IF
100   CONTINUE
110   ITER(4)=NACT-I0
      IF (AMINLA .GE. 0.0) GO TO 150
C
C     INSERT A KNOT WHERE THE LEAST NEGATIVE LAGRANGE MULTIPLIER.
C
      INKNOT=IPAR+1
      PAR(IPAR)=0.0
      CALL ADDKNT(INKNOT,E,B,D,X,F,AKN,IAKN,IACT,NK,NACT,I1,N)
      GO TO 90
C
C.... M A I N  I T E R A T I O N ......................................
C
150   CONTINUE
C
C     EXAMINE CONVEXITES AT ALL KNOTS OF SPLINE APPROXIMATION.
C      SPECIAL CASE: IF NK-I0=1 THEN SX(.) IS A STRAIGHT LINE.
C
      IF (NK-I0 .EQ. 1) GO TO 1000
C
C     FIND THE KNOT OF THE MOST VIOLATED CONSTRAINT.
C
      CVMAX=0.
      KL=I0
      K=0
      DO 160 J=I1,NK-1
         JJ=IAKN(J)
         CKNOT=(SX(JJ+1)-SX(JJ))/(X(JJ+1)-X(JJ))-
     1         (SX(JJ)-SX(JJ-1))/(X(JJ)-X(JJ-1))
         SUM2=1/(X(JJ+1)-X(JJ))**2+
     1   ((X(JJ+1)-X(JJ-1))/((X(JJ+1)-X(JJ))*(X(JJ)-X(JJ-1))))**2+
     2   1/(X(JJ)-X(JJ-1))**2
         SUM2=SQRT(SUM2)
         CKNOT=CKNOT/SUM2
         IF (CKNOT .LT. 0.) K=K+1
         IF (CKNOT .LT. CVMAX) THEN
             CVMAX=CKNOT
             KNTOUT=J
             IAKNJ=IAKN(J)
             KL=KL+1
         END IF
160   CONTINUE
C
C     IF KL=I0 THEN SX(.) IS CONVEX. SINCE AT THIS LEVEL PAR(.) ARE
C      ALL NON-NEGATIVE THE KUHN-TUCKER CONDITIONS ARE SATISFIED.
C
C.... T E R M I N A T I O N ...........................................
C
      IF (KL .EQ. I0) GO TO 1000
C
C     MSG: PRINT CVMAX,IAKNJ,K.
C
      IF (IPRINT.NE.0) CALL MESSGE(CVMAX,K,IAKNJ,ITER,160)
C
C     POWELL'S SAFEGUARD TERMINATION TEST:
C      RETURN, IF DUE TO ROUNDING ERRORS, THE CHANGE IN SX(.) MAY
C      NOT INCREASE THE OBJECTIVE FUNCTION.
C
      IF (ITER(1).EQ.0) GO TO 190
          R=0.0
          T=0.0
          DO 170 I=I1,N
             R=R+(SX(I)-SXOLD(I))*(SX(I)+SXOLD(I)-2*F(I))
170          T=T+ABS(SX(I)-SXOLD(I))*(ABS(SX(I)+SXOLD(I))+2*ABS(F(I)))
          IF (T+R.LE.T .OR. T+1.5*R.LE.T+R) THEN
              IF (ITEST.EQ.1) THEN
                  MODE=2
                  IF (IPRINT.NE.0) CALL MESSGE(ZERO,0,0,ITER,170)
C
C                     MSG: L2CXFT RETURNS BECAUSE THE ACCURACY IS
C                          INSUFFICIENT TO MAINTAIN INCREASING
C                          THE FUNCTION VALUES.
C
                  GO TO 1000
              END IF
C
              IF (IPRINT.EQ.0) GO TO 180
              IF (ITREF.EQ.0) CALL MESSGE (ZERO,0,0,ITER,180)
C
C                 MSG: ITERATIVE REFINEMENT IS INITIATED.
C
180           ITREF=1
              ITEST=1
          ELSE
              ITEST=0
          END IF
190   CONTINUE
C
C     KEEP OLD LAGRANGE MULTIPLIERS AND OLD SPLINE APPROXIMATION.
C
      DO 200 I=I1,N-2
200      PAROLD(I)=PAR(I)
      DO 210 I=I1,N
210      SXOLD(I)=SX(I)
C
C     TERMINATION ON ITERATION COUNT.
C
      ITER(1)=ITER(1)+1
      IF (ITER(1).GT.3*(N-I0)) THEN
          MODE=3
          CALL MESSGE(ZERO,0,0,ITER,210)
C
C         MSG: L2CXFT QUITS ON ITERATION COUNT.
C
          GO TO 1000
      END IF
C
C     CALCULATE A SEARCH DIRECTION (DELETE THE KNOT KNTOUT).
C
      CALL DELKNT(E,B,D,AKN,IAKN,IACT,NK,NACT,I1,N,KNTOUT)
      ITER(6)=ITER(6)+1
C
C     SOLVE FOR SX(.) OVER AKN(.).
C
220   CALL CFLGRF(SX,PAR,COEF,IACT,E,B,D,X,AKN,IAKN,F,DC,EC,U,
     1  RES,NK,NACT,I1,N,ITREF,ITER(3))
      ITER(3)=ITER(3)+1
C
      IPAR=IAKNJ-1
      IF (IPRINT.NE.0) CALL MESSGE(PAR(IPAR),0,0,ITER,220)
C
C     MSG: PRINT MULTIPLIER OF INSERTED CONSTRAINT.
C
C     APPLY POWELL'S TECHNIQUE IN ORDER THAT THE INDEX IPAR
C      REMAIN IN ACTIVE SET. SET PAR(IPAR)=MAX(PAR(IPAR),0.0)
C      BECAUSE IN THEORY HOLDS THAT PAR(IPAR)>0.
C
      IF (PAR(IPAR).LT.0.0) PAR(IPAR)=0.0
C
C     EXAMINE SIGN OF LAGRANGE MULTIPLIERS
C
300   IF (NACT .EQ. I0) GO TO 150
      KL=I0
      THETA=2.0
      DO 320 I=I1,NACT
         K=IACT(I)
         IF (PAR(K)) 310,320,320
310      RATIO=PAROLD(K)/(PAROLD(K)-PAR(K))
         IF (THETA .LE. RATIO) GO TO 320
             THETA=RATIO
             IPAR=K
             KL=I1
320      CONTINUE
C
C.... ITERATE IN THE PRIMAL SPACE ......................................
C     IF KL=I0 THEN PAR(.) ARE NON-NEGATIVE. BRANCH WHERE TEST
C      PRIMAL FEASIBILITY.
C
      IF (KL .EQ. I0) GO TO 150
C
C.... RECOVER NON-NEGATIVITY OF LAGRANGE MULTIPLIERS ...................
C
      DO 340 I=I1,NACT
         K=IACT(I)
340      PAROLD(K)=(1-THETA)*PAROLD(K)+THETA*PAR(K)
C
C     INSERT KNOT IPAR+1 INTO THE SPLINE AND SOLVE FOR SX(.).
C
      I1NEW=I1
350   INKNOT=IPAR+1
      PAR(IPAR)=0.
      CALL ADDKNT(INKNOT,E,B,D,X,F,AKN,IAKN,IACT,NK,NACT,I1,N)
      ITER(5)=ITER(5)+1
C
C     TEST THE FEASIBILITY OF PAROLD(.). ALL PAROLD(.) MUST BE
C      NON-NEGATIVE BECAUSE OF THE WAY THETA WAS DETERMINED. OTHERWISE,
C      THIS IS DUE TO ROUND-OFF ERROR. HENCE, ANY NEGATIVITY
C      IS TREATED AS BEING ZERO. SEE LAWSON & HANSON (SOLVING
C      LEAST SQUARES PROBLEMS, PRENTICE-HALL, 1974, p.307).
C
      IF (NACT .EQ. I0) GO TO 370
      IF (I1NEW .GT. NACT) GO TO 370
      DO 360 I=I1NEW,NACT
         IPAR=IACT(I)
         IF (PAROLD(IPAR) .LT. 0.0) THEN
             IF (IPRINT.NE.0) CALL MESSGE(PAROLD(IPAR),IPAR,0,ITER,350)
C
C                MSG: ITERATION IN THE DUAL SPACE DUE TO ROUND OFF.
C                PRINT ACTIVE CONSTRAINT AND LAGRANGE MULTIPLIER.
C
             I1NEW=I
             GO TO 350
         END IF
360   CONTINUE
C
370   CONTINUE
C
      CALL CFLGRF(SX,PAR,COEF,IACT,E,B,D,X,AKN,IAKN,F,DC,EC,U,
     1  RES,NK,NACT,I1,N,ITREF,ITER(3))
      ITER(3)=ITER(3)+1
C
C.... ITERATE IN THE DUAL SPACE ........................................
C
      GO TO 300
1000  CONTINUE
C
      CALL MESSGE(ZERO,NACT-I0,MODE,ITER,1000)
C
C     MSG: PRINT MODE AND
C          INITIAL ACTIVE SET
C          STARTING ACTIVE SET
C          FINAL ACTIVE SET
C          ADDITIONS OF CONSTRAINTS
C          DELETIONS OF CONSTRAINTS
C
C.... END OF SUBROUTINE L2CXFT .........................................
C
      RETURN
      END
C......................................................................
C     GIVEN N NOISY FUNCTION MEASUREMENTS F(I): I=1,2,...N
C      AT THE ABSCISSAE X(1)<X(2)<...<X(N),
C     THIS SUBROUTINE CALCULATES IN ONLY O(N) COMPUTER OPERATIONS AN
C      EXCELLENT STARTING POINT FOR THE FOLLOWING DATA FIT CALCULATION:
C
C      FIND SX(.) THAT MINIMIZES THE SUM OF THE SQUARES OF THE ERRORS
C            SX(I)-F(I): I=1,2,...,N
C      SUBJECT TO NON-NEGATIVITY OF SECOND DIVIDED DIFFERENCES
C            OF THE SMOOTHED VALUES SX(.).
C
C     SX(.) APPEARS DUE A B-SPLINE REPRESENTATION (SEE METHOD BELOW).
C......................................................................
C     THIS CODE HAS SUBMITTED AS REF 1 SHOWS.
C     REF 1: "A O(N) STARTING POINT FORTRAN SUBROUTINE FOR LEAST SQUARES
C             DATA FITTING WITH NON-NEGATIVE SECOND DIVIDED DIFFERENCES"
C             BY I.C. DEMETRIOU, SUBMITTED TO ACM TRANSACTIONS ON
C             MATHEMATICAL SOFTWARE, MARCH 1990.
C     THE ALGORITHM OF THIS CODE CAN BE FOUND AS REF 2 SHOWS.
C     REF 2: "THE MINIMUM SUM OF SQUARES CHANGE TO UNIVARIATE DATA
C             THAT GIVES CONVEXITY"
C             BY I.C.DEMETRIOU AND M.J.D. POWELL, SUBMITTED TO IMA J. OF
C             NUMERICAL ANALYSIS, 1989.
C.......................................................................
        SUBROUTINE L2ILS(COEF,AKN,IAKN,D,E,B,Z,U,NK,F,X,I1,N)
*       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION X(I1:N),F(I1:N),COEF(I1-1:N),AKN(I1-2:N),
     1    IAKN(I1-1:N),D(I1-1:N),E(I1-1:N),B(I1-1:N),Z(I1-1:N),U(I1-1:N)
C
C ..... I N P U T .....
C       I1      INTEGER VARIABLE, LOWER DATA INDEX.
C       N       INTEGER VARIABLE, UPPER DATA INDEX.
C       X(I1:N) REAL ARRAY OF ABSCISSAE.
C       F(I1:N) REAL ARRAY OF MEASUREMENTS.
C
C ..... O U T P U T .....
C       NK      INTEGER VARIABLE, UPPER INDEX OF SPLINE KNOTS SEQUENCE.
C       AKN(I1-2:NK)  REAL ARRAY CONTAINING THE SPLINE KNOTS, WHERE
C               AKN(I1-2)=AKN(I1-1).
C       IAKN(I1-1:NK)  INTEGER ARRAY CONTAINING THE DATA POINT INDICES
C               OF THE ASSOCIATED KNOTS.
C       COEF(I1-1:NK)  REAL ARRAY CONTAINING THE SPLINE COEFFICIENTS
C               AT AKN(I), I=I1-1,I1,...,NK.
C
C ..... W O R K I N G   S P A C E .....
C       E(I1-1:N)  REAL ARRAY THAT HOLDS THE DIAGONAL OF THE NORMAL
C               EQUATIONS ASSOCIATED WITH SPLINE SX(.).
C       D(I1-1:N)  REAL ARRAY THAT HOLDS THE SUBDIAGONAL OF THE
C               NORMAL EQUATIONS ASSOCIATED WITH SX(.).
C       U(I1-1:N)  REAL ARRAY THAT HOLDS THE DIAGONAL OF CHOLESKY
C               FACTORIZATION OF THE NORMAL EQUATIONS.
C       B(I1-1:N)  REAL ARRAY THAT HOLDS THE RIGHT HAND SIDE OF
C               THE NORMAL EQUATIONS.
C       Z(I1-1:N)  REAL ARRAY THAT HOLDS THE RIGHT HAND SIDE OF
C               THE CHOLESKY FACTORIZATION.
C
C ..... M E T H O D (ALGORITHM 1 IN QUOTED REFERENCE) .....
C       EXPRESS SPLINE SX(.) AS A LINEAR COMBINATION OF B-SPLINES
C               SX(T)=SUM(COEF(J)*N(I,T), I=I1-1,I1,...,J)
C                                            FOR T=I1,I1+1,...,N
C       WHERE N(I,T)=THE T-TH COMPONENT OF THE I-TH LINEAR B-SPLINE.
C       DENOTE BY <.,.> THE INNER PRODUCT OF TWO VECTORS AND,
C       THE COEFFICIENTS OF SX(.) ARE DETERMINED BY SOLVING THE
C       NORMAL EQUATIONS
C               D(I)*COEF(I-1)+E(I)*COEF(I)+D(I+1)*COEF(I+1)=B(I)
C                                            FOR I=I1-1,I1,...,J
C       WHERE
C               D(I)=<N(I,.),N(I-1,.)>, I=I1-1,I1,...,J
C               E(I)=<N(I,.),N(I,.)>, I=I1-1,I1,...,J
C               B(I)=<N(I,.),F(.)>, I=I1-1,I1,...,J.
C
C       THEN PROCEED AS IN REFERENCE QUOTED.
C......................................................................
C       NO TEST ON THE CONDITIONS N<I1, N=I1 AND N=I1+1 IS MADE
C        HERE BECAUSE THEY ARE CONSIDERED BY SUBROUTINE L2CXFT.
C
C       SET INITIAL VALUES.
C
        I0=I1-1
100     K=I0
        J=I1
        AKN(I0)=X(I1)
        AKN(I0-1)=AKN(I0)
        AKN(I1)=X(I1+1)
        AKN(I1+1)=AKN(I1)
        IAKN(I0)=I1
        IAKN(I1)=I1+1
C
C       SET NORMAL EQUATIONS.
C
        E(I0)=1.
        E(I1)=1.
        D(I1)=0.
        B(I0)=F(I1)
        B(I1)=F(I1+1)
C
C       SET CHOLESKY FACTORIZATION AND COEFFICIENTS.
C
        U(I0)=1.
        U(I1)=1.
        Z(I0)=B(I0)
        Z(I1)=B(I1)
        COEF(I0)=Z(I0)
        COEF(I1)=Z(I1)
110     IF (K.EQ.N-2) GOTO 170
C
C       INSERT A KNOT AT X(K+2).
C
        J=J+1
        AKN(J)=X(K+3)
        IAKN(J)=K+3
        AKN(J+1)=AKN(J)
C
C       UPDATE NORMAL EQUATIONS AND CHOLESKY FACTORIZATION.
C
        E(J)=1.
        D(J)=0.
        B(J)=F(K+3)
        U(J)=1.
        Z(J)=B(J)
C
C       EXAMINE RIGHTMOST KNOT CONVEXITY.
C
120     IF (J.EQ.I1) THEN
            K=K+1
            GO TO 110
        END IF
C
C       CALCULATE RIGHTMOST THREE COEFFICIENTS.
C
        COEF(J)=Z(J)/U(J)
        DO 130 I=J-1,J-2,-1
130        COEF(I)=(Z(I)-D(I+1)*COEF(I+1))/U(I)
C
C       DETERMINE SIGN OF RIGHTMOST SECOND DIVIDED DIFFERENCE.
C
140     CAKNOT=(COEF(J)-COEF(J-1))/(AKN(J)-AKN(J-1))-
     1  (COEF(J-1)-COEF(J-2))/(AKN(J-1)-AKN(J-2))
        IF (CAKNOT.GT.0.0) THEN
            K=K+1
            GOTO 110
        END IF
C
C       DELETE KNOT AKN(J-1)
C       UPDATE NORMAL EQUATIONS AND CHOLESKY FACTORIZATION.
C
        C1=(AKN(J-1)-AKN(J-2))/(AKN(J)-AKN(J-2))
        C2=(AKN(J)-AKN(J-1))/(AKN(J)-AKN(J-2))
        E(J-2)=E(J-2)+C2*C2*E(J-1)+2*C2*D(J-1)
        D(J-1)=C1*D(J-1)+C2*D(J)+C1*C2*E(J-1)
        E(J-1)=C1*C1*E(J-1)+E(J)+2*C1*D(J)
        B(J-2)=B(J-2)+C2*B(J-1)
        B(J-1)=C1*B(J-1)+B(J)
C
        IF (J.EQ.I0+2) THEN
            U(J-2)=E(J-2)
            Z(J-2)=B(J-2)
        ELSE
            U(J-2)=E(J-2)-D(J-2)*D(J-2)/U(J-3)
            Z(J-2)=B(J-2)-D(J-2)*Z(J-3)/U(J-3)
        END IF
        U(J-1)=E(J-1)-D(J-1)*D(J-1)/U(J-2)
        Z(J-1)=B(J-1)-D(J-1)*Z(J-2)/U(J-2)
        U(J)=0.
        Z(J)=0.
        E(J)=0.
        D(J)=0.
        B(J)=0.
        J=J-1
        AKN(J)=X(K+3)
        IAKN(J)=K+3
        AKN(J+1)=AKN(J)
        GOTO 120
C
C       DERIVE COEFFICIENTS BY BACK SUBSTITUTION.
C
170     NK=J
        COEF(NK)=Z(NK)/U(NK)
        DO 180 I=NK-1,I0,-1
180        COEF(I)=(Z(I)-D(I+1)*COEF(I+1))/U(I)
C
        RETURN
        END
C.......................................................................
C     THIS SUBROUTINE CALCULATES THE LAGRANGE MULTIPLIERS PAR(.)
C      ASSOCIATED WITH THE SPLINE SX(.) IN SUBROUTINE L2CXFT.
C.......................................................................
      SUBROUTINE LAGRNG(PAR,SX,F,IACT,NACT,EC,DC,UE,I1,N)
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION SX(I1:N),F(I1:N),DC(I1:N),EC(I1:N),UE(I1-1:N),
     1  IACT(I1:N),PAR(I1:N)
C
C.... I N P U T ....
C     SX,F,IACT,DC,EC,I1,N,NACT  AS THEY ARE DEFINED IN SUBROUTINE
C             L2CXFT.
C
C.... O U T P U T ....
C     PAR(I1:N-2) REAL ARRAY CONTAINING THE LAGRANGE MULTIPLIERS SO
C             THAT PAR(IACT(I)), I=I1,I1+1,...,NACT ARE THE MULTI-
C             PLIERS OF THE ACTIVE CONSTRAINTS WHILE THE REMAINING
C             MULTIPLIERS EQUAL TO ZERO.
C
C.... W O R K I N G  S P A C E ....
C     UE(I1-1:NK+1) REAL ARRAY THAT PROVIDES WORKING SPACE FOR THE
C             FACTORIZATION THAT GIVES PAR(.). IT KEEPS THE
C             DIAGONAL ELEMENTS OF THE FORWARD SWEEP OF THE
C             FACTORIZATION.
C
C.... M E T H O D (CALCULATE THE LAGRANGE MULTIPLIERS) .....
C     THE FIRST ORDER OPTIMALITY CONDITIONS ARE
C             SUM(PAR(II)*A(II,J))=2*(SX(II+1)-F(II+1)),
C                                   FOR II=IACT(I), I=I1,...,NACT
C      WHERE A(II,J) DENOTES THE J-TH COMPONENT OF THE II-TH SECOND
C      DIVIDED DIFFERENCE. ONLY THE A(II,J) FOR J=II,II+1,II+2 ARE
C      NON-ZERO. THESE CONDITIONS CONSTITUTE AN OVERDETERMINED SYSTEM
C      OF EQUATIONS. CHOOSE THE EQUATIONS THAT CORRESPOND TO THE SECOND
C      NON-ZERO COMPONENT OF EACH A(II,.) AND SOLVE BY GAUSS ELIMINATION
C      WITHOUT PIVOTING THE DERIVED BLOCK DIAGONAL POSITIVE DEFINITE
C      SYSTEM OF EQUATIONS. PROCEED AS IT IS EXPLAINED IN THE TEXTUAL
C      PART THIS ALGORITHM.
C.......................................................................
      IF (NACT.LE.I1-1) GO TO 170
C
C     THE GRADIENT OF THE OBJECTIVE FUNCTION IS DEFINED BY, SAY,
C      GRAD(I)=2*(SX(I+1)-F(I+1)), I=I1,I1+1,...,N-2. HOWEVER, INSTEAD
C      OF AN ARRAY GRAD(.) THE SUBROUTINE USES DIRECTLY THE VALUES OF
C      THE GRADIENT AT X(.) WHEN NECESSARY.
C
C     CALCULATE SIZE K OF A BLOCK OF THE LAGRANGE COEFFICIENT
C      MATRIX.
C
      I=I1
110   K=1
120   IF (I.LT.NACT .AND. IACT(I).EQ.IACT(I+1)-1) THEN
          I=I+1
          K=K+1
          GO TO 120
      END IF
C
130   IF (K.EQ.1) THEN
          II=IACT(I)
          PAR(II)=(SX(II+1)-F(II+1))/EC(II)
      ELSE
C
C         FACTORIZE THE K*K 3/DIAG SYSTEM OF LAGRANGE COEF.
C
          UE(I-K+1)=EC(IACT(I-K+1))
          DO 140 L=I-K+2,I
             DI=DC(IACT(L))
140          UE(L)=EC(IACT(L))-DI/UE(L-1)*DI
C
C         USE THIS FACTORIZATION TO SOLVE FOR PAR(.).
C
          GRAD=SX(IACT(I-K+1)+1)-F(IACT(I-K+1)+1)
          PAR(IACT(I-K+1))=GRAD
          DO 150 L=I-K+2,I
             GRAD=SX(IACT(L)+1)-F(IACT(L)+1)
150          PAR(IACT(L))=GRAD-DC(IACT(L))*PAR(IACT(L-1))/UE(L-1)
          PAR(IACT(I))=PAR(IACT(I))/UE(I)
          DO 160 L=I-1,I-K+1,-1
160          PAR(IACT(L))=PAR(IACT(L))/UE(L)-DC(IACT(L+1))*
     1                                     PAR(IACT(L+1))/UE(L)
      END IF
      I=I+1
      IF (I-NACT) 110,110,170
C
170   RETURN
      END
C.......................................................................
C     MESSAGES ASSOCIATED WITH THE OPERATIONS OF SUBROUTINE L2CXFT.
C.......................................................................
      SUBROUTINE MESSGE(TEMP,ITEMP1,ITEMP2,ITER,LABEL)
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ITER(*)
C
C.... I N P U T ....
C     TEMP,ITEMP1,ITEMP2   VARIABLES THAT ALLOW VALUES TO PASS FROM
C                 SUBROUTINE L2CXFT.
C     ITER(*)    INTEGER ARRAY, HOLDING NUMBER OF ITERATIONS. IT IS
C                 DEFINED IN SUBROUTINE L2CXFT.
C     LABEL      INTEGER VARIABLE, INDICATING THE NEAREST LABEL IN
C                 SUBROUTINE L2CXFT, WHERE THE MESSAGE IS CAUSED.
C
C.... O U T P U T ....
C     VARIOUS MESSAGES THAT ARE ACTIVATED ONLY IF THE ARGUMENT IPRINT
C      OF SUBROUTINE L2CXFT HAS THE VALUE 1.
C......................................................................
C
      IF (LABEL.EQ.10) THEN
          PRINT 1010
1010      FORMAT(//5X,'ERROR RETURN FROM L2CXFT BECAUSE THE ',
     1            /5X,'DATA LOWER BOUND EXCEEDS THE UPPER BOUND.')
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.20) THEN
          PRINT 1020
1020      FORMAT(//5X,'THE DATA ITSELF IS THE REQUIRED APPROXIMATION.')
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.160) THEN
          PRINT 1160,ITER(1),ITEMP1,TEMP,ITEMP2
1160      FORMAT (//5X,'ITERATIONS = ',I5,
     1            /5X,'NUMBER OF VIOLATED CONSTRAINTS    =',I5,
     2            /5X,'MAXIMUM CONSTRAINT  VIOLATION     =',E20.10,
     3            /5X,'                AT DATA POINT     =',I5)
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.170) THEN
          PRINT 1170
1170      FORMAT (//5X,'L2CXFT RETURNS BECAUSE THE ACCURACY IS'
     1      ' INSUFFICIENT',/5X,'TO MAINTAIN INCREASING THE '
     2      'FUNCTION VALUES.')
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.180) THEN
          PRINT 1180
1180      FORMAT (/5X,'ITERATIVE REFINEMENT IS INITIATED.')
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.210) THEN
          PRINT 1210,ITER(1)
1210      FORMAT (//5X,'L2CXFT QUITS ON ITERATION COUNT = ',I5)
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.220) THEN
          PRINT 1220,TEMP
1220      FORMAT (/5X,'MULTIPLIER OF INSERTED CONSTRAINT =',E20.10)
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.350) THEN
          PRINT 1350,ITEMP1,TEMP
1350      FORMAT (/5X,'ITERATION IN THE DUAL SPACE DUE TO ROUND OFF',
     1      /5X,'CONSTRAINT INDEX    = ',I5,
     2      /5X,'LAGRANGE MULTIPLIER = ',E20.10)
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.1000) THEN
          PRINT 1500,ITEMP2,ITER(2),ITER(4),ITEMP1,ITER(6),ITER(5)
1500      FORMAT (//5X,'---------------------------------------',
     1             /5X,'RETURN FROM L2CXFT WITH MODE = ',I5,
     2            //5X,'INITIAL ACTIVE SET  = ',I5,
     3             /5X,'STARTING ACTIVE SET = ',I5,
     4             /5X,'FINAL ACTIVE SET    = ',I5,
     5            //5X,'ADDITIONS OF CONSTRAINTS = ',I5,
     6             /5X,'DELETIONS OF CONSTRAINTS = ',I5,
     7             /5X,'---------------------------------------')
          GO TO 2000
      END IF
C
2000  RETURN
      END
C......................................................................
C     THIS SUBROUTINE REVISES THE (J-1) AND (J+1) NORMAL EQUATIONS
C      OF SX(.) AFTER DELETING THE KNOT AKN(J).
C     THE MODIFIED ELEMENTS REMAIN STORED IN THE SAME POSITIONS.
C......................................................................
      SUBROUTINE MINUS1(E,B,D,AKN,I1,N,J)
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E(I1-1:N),B(I1-1:N),D(I1-1:N),AKN(I1-2:N)
C
C.... I N P U T ....
C     J       INTEGER VARIABLE, SUCH THAT AKN(J) IS THE DELETED KNOT.
C     E,B,D,AKN,I1,N  AS THEY ARE DEFINED IN SUBROUTINE L2CXFT.
C
C.... O U T P U T ....
C     E(J-1),E(J+1),B(J-1),B(J),D(J+1) THE ONLY ELEMENTS OF THE
C             NORMAL EQUATIONS THAT NEED BE MODIFIED BECAUSE OF J-TH
C             KNOT DELETION.
C
C.... M E T H O D ....
C     REF: FORMULAE (2.9) AND (2.10) IN "THE MINIMUM SUM OF
C          SQUARES CHANGE TO UNIVARIATE DATA THAT GIVES CONVEXITY",
C          BY I.C.DEMETRIOU AND M.J.D.POWELL, SUBMITTED TO IMA J. OF
C          NUMERICAL ANALYSIS, 1989.
C......................................................................
      C1=(AKN(J)-AKN(J-1))/(AKN(J+1)-AKN(J-1))
      C2=(AKN(J+1)-AKN(J))/(AKN(J+1)-AKN(J-1))
      E(J-1)=E(J-1)+C2*C2*E(J)+2*C2*D(J)
      E(J+1)=C1*C1*E(J)+E(J+1)+2*C1*D(J+1)
      D(J+1)=C1*D(J)+C2*D(J+1)+C1*C2*E(J)
      B(J-1)=B(J-1)+C2*B(J)
      B(J+1)=C1*B(J)+B(J+1)
C
      RETURN
      END
C.......................................................................
C     THIS SUBROUTINE REVISES THE (J-1), THE J AND THE (J+1) NORMAL
C      EQUATIONS OF SPLINE SX(.), AFTER X(IQ) IS INSERTED IN THE KNOT
C      SEQUENCE AND AKN(J) IS REDEFINED AS AKN(J)=X(IQ).
C.......................................................................
      SUBROUTINE PLUS(X,F,AKN,E,B,D,NK,I1,N,IPL,IP,IQ,IR,IRR,J)
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(I1:N),F(I1:N),AKN(I1-2:N),E(I1-1:N),B(I1-1:N),
     1  D(I1-1:N),SV(24)
C
C.... I N P U T ....
C     E,B,D,AKN,X,F,I1,N,NK  AS THEY ARE DEFINED IN SUBROUTINE L2CXFT.
C     J      INTEGER VARIABLE, AS DEFINED IN SUBROUTINE ADDKNT, INDEX
C            OF NEW INSERTED KNOT.
C     IP     INTEGER VARIABLE, DATA INDEX OF THE KNOT ADJACENT TO
C             THE KNOT AKN(J) AT THE LEFT, IP=IAKN(J-1).
C     IPL    IF J=I1 THEN IPL=I1 ELSE IPL=IAKN(J-1).
C     IQ     INTEGER VARIABLE, DATA INDEX OF INSERTED KNOT, IQ=IAKN(J).
C     IR     INTEGER VARIABLE, DATA INDEX OF THE KNOT ADJACENT TO
C             THE KNOT AKN(J) AT THE RIGHT, IR=IAKN(J+1).
C     IRR    IF J=NK-1 THEN IRR=N ELSE IRR=IAKN(J+1).
C
C.... W O R K I N G  S P A C E ...
C     SV(1:24) REAL ARRAY WHOSE ELEMENTS PROVIDE TEMPORARY STORAGE
C            FOR THE CALCULATIONS OF (E(K),B(K), K=J-1,J,J+1) AND
C            (D(K), K=J,J+1).
C
C.... M E T H O D ....
C     GIVEN THE DATA POINTS X(IPL), X(IP), X(IQ), X(IR) AND X(IRR)
C      AS THE METHOD IN SUBROUTINE ADDKNT DESCRIBES, THEN CALCULATE
C      AB INITIO THE ELEMENTS OF E(K) AND B(K) FOR K=J-1,J,J+1
C      AND D(K) FOR K=J,J+1 (SEE C. DE BOOR, A PRACTICAL GUIDE
C      TO SPLINES, SPRINGER-VERLAG, 1978).
C.......................................................................
C     CALCULATE THE SCALAR PRODUCTS WHICH DEFINE E(.), B(.) AND D(.).
C      SEE COMMENTS IN THE BEGINNING OF SUBROUTINE L2CXFT.
C
      DO 50 I=1,24
50       SV(I)=0.
C
      DO 100 I=IP,IQ-1
         SV(1)=X(I)-AKN(J-1)
         SV(2)=AKN(J)-X(I)
         SV(3)=SV(3)+SV(2)*SV(2)
         SV(4)=SV(4)+SV(1)*SV(1)
         SV(5)=SV(5)+SV(2)*SV(1)
         SV(6)=SV(6)+SV(2)*F(I)
100      SV(7)=SV(7)+SV(1)*F(I)
C
      DO 200 I=IQ,IR-1
         SV(14)=AKN(J+1)-X(I)
         SV(2)=X(I)-AKN(J)
         SV(15)=SV(15)+SV(14)*SV(14)
         SV(16)=SV(16)+SV(2)*SV(2)
         SV(17)=SV(17)+SV(14)*SV(2)
         SV(18)=SV(18)+SV(14)*F(I)
200      SV(19)=SV(19)+SV(2)*F(I)
C
      IF (J.EQ.NK-1) THEN
          SV(11)=(AKN(J+1)-AKN(J))**2
          SV(9)=(AKN(J+1)-AKN(J))*F(N)
      ELSE
          DO 300 I=IR,IRR-1
             SV(24)=AKN(J+2)-X(I)
             SV(9)=SV(9)+SV(24)*F(I)
300          SV(11)=SV(11)+SV(24)*SV(24)
      END IF
C
      IF (I1.LT.J) THEN
          DO 400 I=IPL,IP-1
             SV(23)=X(I)-AKN(J-2)
             SV(8)=SV(8)+SV(23)*F(I)
400          SV(10)=SV(10)+SV(23)*SV(23)
      END IF
C
C     CALCULATE B(J-1),B(J),B(J+1).
C
      SV(2)=AKN(J-1)-AKN(J-2)
      SV(14)=AKN(J)-AKN(J-1)
      SV(21)=AKN(J+1)-AKN(J)
      SV(22)=AKN(J+2)-AKN(J+1)
C
      IF (J.EQ.I1) THEN
          B(J-1)=SV(6)/SV(14)
      ELSE
          B(J-1)=SV(8)/SV(2)+SV(6)/SV(14)
      END IF
      B(J)=SV(7)/SV(14)+SV(18)/SV(21)
C
      IF (J.EQ.NK-1) THEN
          B(J+1)=SV(19)/SV(21)+F(N)
      ELSE
          B(J+1)=SV(19)/SV(21)+SV(9)/SV(22)
      END IF
C
C     CALCULATE E(J-1),E(J),E(J+1),D(J-1),D(J).
C
      SV(2)=SV(2)*SV(2)
      SV(14)=SV(14)*SV(14)
      SV(21)=SV(21)*SV(21)
      SV(22)=SV(22)*SV(22)
C
      IF (J.EQ.I1) THEN
          E(J-1)=SV(3)/SV(14)
      ELSE
          E(J-1)=SV(10)/SV(2)+SV(3)/SV(14)
      END IF
      E(J)=SV(4)/SV(14)+SV(15)/SV(21)
C
      IF (J.EQ.NK-1) THEN
          E(J+1)=SV(16)/SV(21)+1
      ELSE
          E(J+1)=SV(16)/SV(21)+SV(11)/SV(22)
      END IF
      D(J)=SV(5)/SV(14)
      D(J+1)=SV(17)/SV(21)
C
      RETURN
      END
C.......................................................................
C     GENERATE SPLINE SX(.) AT X(.) FROM COEF(.) OVER AKN(.).
C.......................................................................
      SUBROUTINE SXATXI(SX,X,COEF,AKN,IAKN,I1,N,NK)
*     IMPLICIT DOUBLE PRECISION  (A-H,O-Z)
      DIMENSION SX(I1:N),X(I1:N),COEF(I1-1:N),AKN(I1-2:N),IAKN(I1-1:N)
C
C.... I N P U T ....
C     X,COEF,AKN,IAKN,I1,N,NK  AS THEY ARE DEFINED IN SUBROUTINE L2CXFT.
C
C.... O U T P U T ....
C     SX(I1:N) REAL ARRAY CONTAINING THE SPLINE COMPONENTS AT X(I)
C               I=I1,I1+1,...,N.
C
C.... M E T H O D (GENERATE SPLINE COMPONENTS FROM COEFFICIENTS) ....
C     SEE ALGORITHM 5.8 OF L.L. SCHUMAKER, SPLINE FUNCTIONS: BASIC
C     THEORY, J. WILEY AND SONS, 1981, P.194 .
C.......................................................................
      I0=I1-1
      SX(I1)=COEF(I0)
      SX(N)=COEF(NK)
      I=I1+1
      L=I0
50    IF (I.LT.N .AND. IAKN(L).LT.I .AND. I.LE.IAKN(L+1)) THEN
          DENOM=AKN(L+1)-AKN(L)
          A1=(X(I)-AKN(L))/DENOM
          A2=1-A1
          SX(I)=COEF(L)*A2+COEF(L+1)*A1
          I=I+1
          GO TO 50
      END IF
      L=L+1
      IF (L.LT.NK) GO TO 50
C
      RETURN
      END
c*** sptests.f
      PROGRAM PROG1
C
C     SAMPLE PROGRAM THAT DEMONSTRATES SUBROUTINE L2CXFT WHEN THE
C      DATA IS READ FROM, AND THE RESULTS ARE DIRECTED TO, FILES.
C
C     CALLS SUBROUTINE L2CXFT.
C.......................................................................
      PARAMETER (IX1=1,
     1           NX=2000,
     2           IX0=IX1-1,
     3           IFILE=1,
     4           IPRINT=0)
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(IX1:NX),F(IX1:NX),SX(IX1:NX),COEF(IX0:NX),
     2  AKN(IX0-1:NX),IAKN(IX0:NX),IACT(IX1:NX),ITER(10),
     3  PAR(IX1:NX),PARNEW(IX1:NX),RES(IX0:NX),DC(IX0:NX),EC(IX0:NX),
     4  E(IX0:NX),B(IX0:NX),D(IX0:NX),Z(IX0:NX),U(IX0:NX),SXOLD(IX1:NX)
C
C......................................................................
C.... I N P U T ....
C     X,F,SX,COEF,AKN,IAKN,IACT,ITER : ARRAYS THAT ARE EXPLAINED IN
C               SUBROUTINE L2CXFT.
C     I1,N      LOWER,UPPER INDICES OF DATA ARRAYS X(.) AND F(.). THEY
C               ARE EXPLAINED IN SUBROUTINE L2CXFT.
C     'XFDTA'   DATAFILE THAT HOLDS I1 AND N, AND THE DATA PAIRS
C               X(I),F(I) FOR I=I1,I1+1,...,N. 'DATA' IS PROVIDED
C               BY THE USER (OR BY PROG3 FOR TESTING PURPOSES).
C     IFILE     INTEGER PARAMETER WHOSE DEFAULT VALUE IS 0.
C               IF IFILE=1  THEN SOME PRINTOUT IS ACTIVATED AT THE
C                           END OF THE SMOOTHING PROCESS.
C     IPRINT    INTEGER PARAMETER WHOSE DEFAULT VALUE IS O.
C               IF IPRINT=1 THEN SUBROUTINE MESSGE IS CALLED BY SUB-
C                           ROUTINE L2CXFT TO PRINT INTERMEDIATE
C                           RESULTS (USEFUL TO SUBRTN L2CXFT INSTALTN).
C
C.... O U T P U T
C     RESULTS FROM THE CALL OF SUBROUTINE L2CXFT: THE LEAST SQUARES
C      BEST CONVEX FIT TO THE DATA X(.)-F(.).
C
C     FILE 'CXAPPX' KEEPS BEST CONVEX FIT.
C     FILE 'SPLNCF' KEEPS SPLINE KNOTS AND COEFFICIENTS OF BEST FIT.
C......................................................................
C
C      READ (UNFORMATTED) DATA FROM DATAFILE.
C
      OPEN (2,FILE='XFDTA')
            REWIND 2
            READ (2,*) I1,N
            I0=I1-1
            DO 100 I=I1,N
100            READ (2,*) X(I),F(I)
      CLOSE(2)
      PRINT 1110
1110  FORMAT(/5X,'X-F DATA SUPPLY FOR SMOOTHING FROM DATAFILE "XFDTA".')
C
C     CALCULATE BEST CONVEX FIT TO DATA.
C
      CALL L2CXFT(X,F,SX,I1,N,AKN,IAKN,COEF,NK,IACT,NACT,PAR,
     1  ITER,SXOLD,E,B,D,Z,U,DC,EC,PARNEW,RES,IPRINT,MODE)
C
C     SEND RESULTS TO FILES.
C
      OPEN (3,FILE='CXAPPX')
            WRITE (3,120) I1,N
120         FORMAT (2I5)
            DO 130 I=I1,N
130            WRITE (3,140) X(I),SX(I)
140            FORMAT(2E20.10)
      CLOSE(3)
      PRINT 1140
1140  FORMAT(/5X,'THE BEST CONVEX APPROXIMATION TO THE DATA ',
     1               'IS NOW KEPT IN FILE "CXAPPX". ')
C
      OPEN (4,FILE='SPLNCF')
            WRITE (4,150) I1,NK
150         FORMAT (2I5)
            DO 160 I=I1-1,NK
160            WRITE (4,170) IAKN(I),COEF(I)
170            FORMAT(I5,E20.10)
      CLOSE(4)
      PRINT 1170
1170  FORMAT(/5X,'SPLINE KNOTS AND COEFFICIENTS OF THE BEST CONVEX',
     1      /5X,'   APPROXIMATION TO THE DATA ARE NOW KEPT IN FILE ',
     2          '........... "SPLNCF". ')
C
C     PRINTOUT.
C      THE USER: IF N=LARGE, MAY DIRECT THIS PRINTOUT TO A FILE.
C
      IF (IFILE.EQ.1) THEN
          WRITE (*,1200) (J,AKN(J),IAKN(J),COEF(J),J=I0,NK)
1200      FORMAT(////
     1     12X,4HKNOT,6X,10HDATA INDEX,6X,13HB-SPLINE COEF/
     2     12X,5H(AKN),6X,6H(IAKN),10X,6H(COEF)//
     3     (I5,F14.8,I5,10X,E20.10))
C
          DO 1310 I=I1,N
             IF (I.EQ.I1) THEN
                 WRITE (*,1250) I,X(I),F(I),SX(I)
             ELSE IF (I.EQ.N) THEN
                      WRITE (*,1260) I,X(I),F(I),SX(I)
             ELSE
                 WRITE (*,1300) I,X(I),F(I),SX(I),PAR(I-1)
             ENDIF
1250         FORMAT(/////7X,11HDATA POINTS,6X,12HMEASUREMENTS,8X,
     1        13HAPPROXIMATION,7X,13HLAGRANGE MULT/
     2        7X,3H(X),14X,3H(F),17X,4H(SX),16X,5H(PAR)//
     3        (I5,F14.8,E20.10,E20.10))
1260          FORMAT(I5,F14.8,E20.10,E20.10)
1300          FORMAT(I5,F14.8,E20.10,E20.10,E20.10)
1310      CONTINUE
C
         ANORM2=0.0
         DO 1350 I=I1,N
1350        ANORM2=ANORM2+(SX(I)-F(I))*(SX(I)-F(I))
            WRITE(*,1400) ANORM2
1400        FORMAT(/2X,43HVALUE OF OBJECTIVE FUNCTION AT THE OPTIMUM=,
     1       E20.10)
      END IF
C
      STOP
      END
      PROGRAM PROG2
C
C     SAMPLE PROGRAM THAT GENERATES PSEUDORANDOMLY DATA FOR
C      TESTING THE SUBROUTINE L2CXFT.
C
C     CALLS FUNCTION RND.
C.......................................................................
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (I1=1,
     1           N=75,
     2           XA=0.0,
     3           XB=1.0,
     4           SIZE=50.0/N**2)
      DIMENSION X(I1:N),F(I1:N)
C
C.... I N P U T  AND  P A R A M E T E R S  (SET BY THE USER) ....
C     I1        LOWER DATA INDEX.
C     N         UPPER DATA INDEX.
C     XA        LEFT LIMIT OF X(.), X(I1)=XA.
C     XB        RIGHT LIMIT OF X(.), X(N)=XB.
C     SIZE      RELATIVE MAGNITUDE OF NOISE.
C                RECOMMENTED VALUES FOR SIZE THAT COVER A WIDE RANGE
C                OF DATA BEHAVIOUR ARE (WHEN I1=1)
C                = 0
C                = 4/N**2
C                = 10/N**2
C                = 50/N**2
C                = 250/N**2
C                = 1000/N**2
C
C.... O U T P U T  (DIRECTED TO THE FILE "XFDTA") ....
C     X(I1:N)   DATA POINTS (ABSCISSAE).
C     F(I1:N)   FUNCTION MEASUREMENTS (INCLUDING NOISE).
C
C.... M E T H O D  (USER INTERFACE) ....
C     THIS PROGRAM GENERATES RANDOMLY DATA AS FOLLOWS:
C      A CONTINOUS (+CONVEX) FUNCTION F(X) DEFINED ON [XA,XB]
C      IS SUPPLIED BY THE USER. THEN F(X) IS EVALUATED ON THE
C      GRID XA=X(I1)<X(I1+1)<...<X(N)=XB AND THEN WE ADD
C      TO THE FUNCTION VALUES RANDOM NUMBERS FROM THE UNIFORM
C      DISTRIBUTION OVER THE INTERVAL [-SIZE,SIZE].
C      THE TYPE OF THE FUNCTION, THE LIMITS OF THE DATA AND THE SIZE
C      OF THE MAGNITUDE OF THE RANDOM NUMBERS ARE UPON THE USER'S
C      DECISION . THE GRID NEED NOT BE EQUALLY SPACED.
C
C     THE METHOD THAT PRODUCES THE PSEUDO-RANDOM NUMBERS MAY BE
C      FOUND IN FUNCTION RND(.).
C
C     THE DATA ARE GENERATED AND THEN DIRECTED TO FILE "XFDTA". THIS
C      FILE IS INPUT TO PROGRAM PROG1.
C......................................................................
C
C     SET DATA POINTS.
C
      STEP=(XB-XA)/(N-I1)
      X(I1)=XA
      A=X(I1)
      DO 20 I=I1+1,N-1
         X(I)=A+STEP
20       A=X(I)
      X(N)=XB
C
C     GENERATE DATA BY ADDING NOISE TO FUNCTION MEASUREMENTS OF THE
C      RUNGE FUNCTION FROM 0.0 TO 1.0.
C     ISEED IS SET BY THE USER. IT CAN BE ANY INTEGER SUCH THAT
C       1 .LE. ISEED .LE. 65535 .
C
      ISEED=1
      DO 50 I=I1,N
50       F(I)=1/(1+10*X(I)*X(I))+SIZE*RND(ISEED)
C
C     SEND (X,F) DATA PAIRS TO A DATAFILE.
C
      OPEN (2,FILE='XFDTA')
            WRITE (2,120) I1,N
120         FORMAT (2I5)
            DO 130 I=I1,N
130            WRITE (2,140) X(I),F(I)
140            FORMAT (2E20.10)
      CLOSE(2)
      PRINT 1140
1140  FORMAT(/5X,'DATA HAS BEEN GENERATED AND KEPT IN FILE "XFDTA".')
C
      STOP
      END
      PROGRAM PROGMD
C......................................................................
C     THIS PROGRAM IS A DRIVER OF SUBROUTINE L2CXFT.
C      IT USES REPEATEDLY THE SIMPLE DRIVER PROGRAM PROG2 BY ADDING
C      RUNDOM NUMBERS (NOISE) TO GIVEN FUNCTION MEASUREMENTS SO THAT
C      SUBSTANTIALLY DIFFERENT DATA SETS ARE PRODUCED AND TESTED.
C      THE NUMBER OF DATA AND THE RELATIVE MAGNITUDE OF THE NOISE ARE
C      VARIABLES IN NESTED LOOPS.
C
C     CALLS SUBROUTINE L2CXFT.
C.......................................................................
      PARAMETER (IX1=1,
     1           NX=2000,
     2           IX0=IX1-1,
     3           IPRINT=0)
C
*     IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(IX1:NX),F(IX1:NX),SX(IX1:NX),COEF(IX0:NX),
     1  AKN(IX0-1:NX),IAKN(IX0:NX),IACT(IX1:NX),ITER(10),
     2  PAR(IX1:NX),PAROLD(IX1:NX),RES(IX0:NX),DC(IX0:NX),EC(IX0:NX),
     3  E(IX0:NX),B(IX0:NX),D(IX0:NX),Z(IX0:NX),U(IX0:NX),SXOLD(IX1:NX)
C
C.... U S E R  I N T E R F A C E  (INPUT AND PARAMETERS) ....
C
C     THE USER HAS ONLY TO SPECIFY THE PARAMETERS IX1 AND NX, A FUNCTION
C      TO PROVIDE MEASUREMENTS AND LET THIS PROGRAM RUN.
C
C     IX1        LOWER DATA INDEX (SET BY THE USER).
C     NX         UPPER DATA INDEX (SET BY THE USER).
C
C     ALL ARRAYS ARE EXPLAINED IN SUBROUTINE L2CXFT.
C     IPRINT    INTEGER PARAMETER WHOSE DEFAULT VALUE IS O.
C               IF IPRINT=1 THEN SUBROUTINE MESSGE IS CALLED BY SUB-
C                           ROUTINE L2CXFT TO PRINT INTERMEDIATE
C                           RESULTS (USEFUL TO SUBRTN L2CXFT INSTALTN).
C
C.... O U T P U T ....
C     RESULTS FROM THE CALL OF SUBROUTINE L2CXFT CONCERNING ACTIVE
C      SET CHANGES. SEE SUBROUTINE MESSGE.
C
C.... M E T H O D  FOR PRODUCING THE RANDOM NUMBERS ....
C     SEE PROG2.
C......................................................................
      I1=IX1
      I0=I1-1
C
C     SIZE      RELATIVE MAGNITUDE OF NOISE.
C                RECOMMENTED VALUES FOR SIZE THAT COVER A WIDE RANGE
C                OF DATA BEHAVIOUR ARE (WHEN I1=1)
C                = 0
C                = 4/N**2
C                = 10/N**2
C                = 50/N**2
C                = 250/N**2
C                = 1000/N**2
C
      DO 160 JJ=0,6
         IF (JJ.EQ.0) ISIZE=0
         IF (JJ.EQ.1) ISIZE=4
         IF (JJ.EQ.2) ISIZE=10
         IF (JJ.EQ.3) ISIZE=50
         IF (JJ.EQ.4) ISIZE=100
         IF (JJ.EQ.5) ISIZE=250
         IF (JJ.EQ.6) ISIZE=1000
C
      DO 150 NN=1,5
         IF (NN.EQ.1)  N=20
         IF (NN.EQ.2)  N=100
         IF (NN.EQ.3)  N=500
         IF (NN.EQ.4)  N=1000
         IF (NN.EQ.5)  N=2000
C
C     SET DATA POINTS.
C
      XA=0.
      XB=1.
      STEP=(XB-XA)/(N-I1)
      X(I1)=XA
      A=X(I1)
      DO 20 I=I1+1,N-1
         X(I)=A+STEP
20       A=X(I)
      X(N)=XB
C
      PRINT *, '--------------------------------------------'
      PRINT *, '   DBL PRECISION'
      PRINT *, '   RUNGE FUNCTION + NOISE '
      PRINT *, '--------------------------------------------'
      PRINT *, ' NUMBER OF DATA = ',N,' ERROR MULTIPLIER= ',ISIZE
      SIZE=FLOAT(ISIZE)/(N*N)
      PRINT *, ' ERROR SIZE     = ',SIZE
C
C     GENERATE MEASUREMENTS WITH NOISE.
C
      ISEED=1
      DO 30 I=I1,N
C        ADD ERROR TO MEASUREMENTS OF THE RUNGE FUNCTION.
         F(I)=1/(1+10.*X(I)*X(I))+SIZE*RND(ISEED)
C
C        ADD ERROR TO MEASUREMENTS OF THE BERNSTEIN FUNCTION.
*         F(I)=X(I)**4-3./4.*X(I)**3+3.*X(I)*X(I)-5./2.*X(I)+
*     1                                      SIZE*RND(ISEED)
30    CONTINUE
C
C.... CALCULATE THE CONVEX FIT ........................................
C
      CALL L2CXFT(X,F,SX,I1,N,AKN,IAKN,COEF,NK,IACT,NACT,PAR,
     1  ITER,SXOLD,E,B,D,Z,U,DC,EC,PAROLD,RES,IPRINT,MODE)
C
C     PRINTOUT
C
      ANORM2=0.0
      DO 50 I=I1,N
50       ANORM2=ANORM2+(SX(I)-F(I))*(SX(I)-F(I))
         PRINT 1050, ANORM2
1050     FORMAT(/2X,43HVALUE OF OBJECTIVE FUNCTION AT THE OPTIMUM=,
     1     E20.10)
C     END OF N LOOP.
150   CONTINUE
C     END OF ISIZE LOOP.
160   CONTINUE
C
      STOP
      END
C.....................................................................
C     THIS FUNCTION GENERATES RANDOM NUMBERS SUCH THAT
C               -1 < RND < 1
C      WHICH ARE ADDED TO GIVEN FUNCTION MEASUREMENTS, IN ORDER
C       TO CONSTRUCT DATA SETS FOR TESTING THE SUBROUTINE L2CXFT.
C.....................................................................
      FUNCTION RND(ISEED)
*     DOUBLE PRECISION RND
C
C     FUNCTION RND CAN GENERATE 65536 "RANDOM" NUMBERS BEFORE
C      REPEATING ITSELF. IT WILL WORK ON MOST COMPUTERS FOR
C      WHICH 2**31-1.LT.MAXINT.
C
C.... I N P U T ....
C     ISEED     INTEGER VARIABLE BETWEEN 1 AND 65535 .
C
C.... O U T P U T ....
C     RND       THE RETURN OF THE FUNCTION EVALUATION,
C               BETWEEN -1 AND 1.
C
C.... M E T H O D (LINEAR CONGRUENTIAL) ....
C      EACH NUMBER IN THE RANDOM SEQUENCE, R(K) SAY,
C      IS CALCULATED FROM ITS PREDECESSOR, R(K-1), USING THE FORMULA:
C               R(K)=(MULTIPLIER * R(K-1) + INCREMENT) MOD 2**16
C      WHERE MULTIPLIER=25173 AND INCREMENT=13849.
C
C..... REF: P.GROGONO, PROGRAMMING IN PASCAL, ADDISON-WESLEY PUB.CO.,
C          1980, P.118.
C.....................................................................
      RND=FLOAT(ISEED)/65535.
      RND=2*RND-1
      ISEED=25173*ISEED+13849
      ISEED=MOD(ISEED,65536)
C
      RETURN
      END
c*** dpsrc.f
C......................................................................
C     THIS SUBROUTINE INSERTS A KNOT, NAMELY THE DATA POINT X(INKNOT),
C      INTO THE SPLINE KNOT SEQUENCE OF THE THE LINEAR SPLINE
C      APPROXIMATION SX(.) DEFINED IN SUBROUTINE L2CXFT.
C
C     CALLS SUBROUTINE PLUS.
C......................................................................
      SUBROUTINE ADDKNT(INKNOT,E,B,D,X,F,AKN,IAKN,IACT,NK,NACT,I1,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E(I1-1:N),B(I1-1:N),D(I1-1:N),X(I1:N),F(I1:N),
     1  AKN(I1-2:N),IAKN(I1-1:N),IACT(I1:N)
C
C.... I N P U T ....
C     INKNOT  INTEGER VARIABLE, SUCH THAT X(INKNOT) IS INSERTED
C             INTO THE KNOT SEQUENCE AKN(.).
C     E,B,D,X,F,AKN,IAKN,IACT,I1,N,NK,NACT  AS THEY ARE
C             DEFINED IN SUBROUTINE L2CXFT.
C
C.... O U T P U T ....
C     E,B,D,AKN,IAKN,IACT  AFTER INSERTING X(INKNOT) INTO AKN(.).
C     NK      INCREASED BY ONE.
C     NACT    DECREASED BY ONE.
C
C.... M E T H O D ....
C     FIND POSITION J IN THE KNOT SEQUENCE AKN(.) SUCH THAT
C      AKN(J-1) < X(INKNOT) < AKN(J) BEFORE INSERTION AND
C      AKN(J)=X(INKNOT) AFTER INSERTION.
C     CALL SUBROUTINE PLUS TO MODIFY THE NORMAL EQUATIONS THAT ARE
C      AFFECTED BY THE KNOT INSERTION. HENCE THE 'HAT' BASIS FUNCTIONS
C      AFFECTED BY THE INSERTION HAVE TO DETERMINED. THE J-TH BASIS
C      FUNCTION IS DEFINED FROM X(IP) TO X(IR), AND IS CENTERED AT
C      X(IQ). ALSO THE (J-1)-TH BASIS FUNCTIONS IS CENTERED AT X(IP) AND
C      IS DEFINED FROM X(IPL) TO X(IQ). THE (J+1)-TH BASIS FUNCTION
C      IS CENTERED AT X(IR) AND IS DEFINED FROM X(IQ) TO X(IRR).
C     REVISE IACT(.), THE INDICES OF THE ACTIVE CONSTRAINTS.
C......................................................................
      I0=I1-1
C
C     GIVEN X(INKNOT), APPLY BISECTION IN [I0,NK) IN ORDER TO FIND IP
C      SUCH THAT  AKN(IP) < X(INKNOT) < AKN(IP+1).
C
      IP=I0
      IR=NK
20    IF (IR-IP.EQ.1) GO TO 30
      IQ=(IP+IR)/2
      IF (INKNOT .LT. IAKN(IQ)) THEN
          IR=IQ
      ELSE
          IP=IQ
      END IF
      GO TO 20
30    CONTINUE
C
      DO 100 I=NK,IP+1,-1
         B(I+1)=B(I)
         E(I+1)=E(I)
         D(I+1)=D(I)
         AKN(I+1)=AKN(I)
100      IAKN(I+1)=IAKN(I)
      J=IP+1
      E(J)=0.
      B(J)=0.
      D(J)=0.
      NK=NK+1
      AKN(J)=X(INKNOT)
      IAKN(J)=INKNOT
C
C     DEFINE SUPPORTS OF NEW BASIS FUNCTIONS THAT FOLLOW THE KNOT
C      INSERTION.
C
      IF (J.EQ.I1) THEN
          IPL=I1
      ELSE
          IPL=IAKN(J-2)
      END IF
C
      IF (J.EQ.NK-1) THEN
          IRR=N
      ELSE
          IRR=IAKN(J+2)
      END IF
C
      IP=IAKN(J-1)
      IQ=IAKN(J)
      IR=IAKN(J+1)
C
C     MODIFY NORMAL EQUATIONS.
C
      CALL PLUS(X,F,AKN,E,B,D,NK,I1,N,IPL,IP,IQ,IR,IRR,J)
C
C     REVISE ACTIVE CONSTRAINTS.
C      NACT>=I1 BECAUSE THE KNOT INSERTION SHOWS THAT THERE EXISTS
C      AT LEAST ONE ACTIVE CONSTRAINT.
C
      IF (NACT.EQ.I1) GO TO 140
C
C     I1<INKNOT<N, K=INKNOT-1=ACTIVE CONSTRAINT TO BE DELETED.
C     GIVEN K, APPLY BISECTION IN [I1,NACT] IN ORDER TO FIND
C      IQ SUCH THAT IACT(IQ).EQ.K .
C
      K=INKNOT-1
      IQ=I1
      IF (K.EQ.IACT(I1)) GO TO 120
      IF (K.EQ.IACT(NACT)) GO TO 140
      IP=I1
      IR=NACT
110   IQ=(IP+IR)/2
      IF (IACT(IQ).EQ.K) GO TO 120
      IF (IACT(IQ).LT.K) THEN
          IP=IQ
      ELSE
          IR=IQ
      END IF
      GO TO 110
120   CONTINUE
      DO 130 I=IQ,NACT-1
130      IACT(I)=IACT(I+1)
140   IACT(NACT)=0
      NACT=NACT-1
C
200   RETURN
      END
C......................................................................
C     THIS SUBROUTINE CALCULATES
C      1) THE COEFFICIENTS COEF(.) OF A LINEAR SPLINE SX(.) OVER
C         KNOTS AKN(.), FROM THE NORMAL EQUATIONS.
C      2) THE SPLINE COMPONENTS SX(.) FROM COEF(.) AND
C      3) THE ASSOCIATED LAGRANGE MULTIPLIERS PAR(.).
C     ALSO THIS SUBROUTINE PROVIDES ITERATIVE REFINEMENT OF
C      1) THE COEFFICIENTS COEF(.) AND
C      2) THE LAGRANGE MULTIPLIERS PAR(.).
C
C     CALLS SUBROUTINES COEFF, SXATXI AND LAGRNG.
C......................................................................
      SUBROUTINE CFLGRF(SX,PAR,COEF,IACT,E,B,D,X,AKN,IAKN,F,DC,EC,
     1  UE,RES,NK,NACT,I1,N,ITREF,ITERTN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION SX(I1:N),COEF(I1-1:N),E(I1-1:N),B(I1-1:N),D(I1-1:N),
     1  AKN(I1-2:N),IAKN(I1-1:N),X(I1:N),F(I1:N),DC(I1:N),EC(I1:N),
     2  IACT(I1:N),UE(I1-1:N),PAR(I1:N),RES(I1-1:N)
C
C.... I N P U T ....
C     E,B,D,X,F,AKN,IAKN,IACT,DC,EC,I1,N,NK AND NACT AS THEY ARE
C              DEFINED IN SUBROUTINE L2CXFT.
C     ITERTN  INTEGER VARIABLE THAT IS USED TO CHECK THE FIRST CALL
C              OF THIS SUBROUTINE BY SUBROUTINE L2CXFT.
C     ITREF   INTEGER VARIABLE, IF ITREF=1 THEN ITERATIVE REFINEMENT
C              IS INITIATED.
C
C.... O U T P U T ....
C     COEF(I1-1:NK) REAL ARRAY CONTAINING THE SPLINE COEFFICIENTS
C             AT AKN(I), I=I1-1,I1,...,NK.
C     SX(I1:N) REAL ARRAY CONTAINING THE SPLINE COMPONENTS AT X(I),
C             I=I1,I1+1,...,N.
C     PAR(I1:N-2) REAL ARRAY CONTAINING THE LAGRANGE MULTIPLIERS SO
C             THAT PAR(IACT(I)), I=I1,I1+1,...,NACT ARE THE ACTIVE
C             CONSTRAINT MULTIPLIERS WHILE THE REMAINING MULTIPLIERS
C             EQUAL TO ZERO.
C
C.... W O R K I N G  S P A C E ....
C     UE(I1-1:NK+1) REAL ARRAY THAT PROVIDES WORKING SPACE FOR THE
C             FACTORIZATIONS THAT GIVE COEF(.) AND PAR(.).
C     RES(I1-1:N) REAL ARRAY THAT PROVIDES WORKING SPACE FOR THE
C             RESIDUALS OF THE NORMAL EQUATIONS AT COEF(.) AND FOR
C             THE CORRECTIONS OF COEF(.). IT ALSO PROVIDES WORKING
C             SPACE FOR THE RESIDUALS OF THE LAGRANGE EQUATIONS AT
C             PAR(.) AND FOR THE CORRECTIONS OF PAR(.).
C
C.... M E T H O D ....
C     FOLLOWS THE STEPS DESCRIBED IN THE HEADING COMMENT.
C......................................................................
      I0=I1-1
C
C     IF ITERTN=0 THEN AVOID CALCULATION OF COEF(.), SINCE COEF(.)
C      ARE SUPPLIED BY SUBROUTINE L2ILS.
C
      IF (ITERTN.EQ.0) GO TO 50
C
C.... DERIVE THE SPLINE COEF(.) ....
C
      CALL COEFF(COEF,E,B,D,UE,NK,I1,N)
C
C.... INITIATE ITERATIVE REFINEMENT IF ITREF=1 ....
C
      IF (ITREF.NE.1) GO TO 50
C
C         REFINE COEF(.).
C         SET RESIDUALS OF NORMAL EQUATIONS IN RES(.).
C
          RES(I0)=B(I0)-(E(I0)*COEF(I0)+D(I1)*COEF(I1))
          IF (I1 .EQ. NK) GO TO 30
          DO 20 I=I1,NK-1
20           RES(I)=B(I)-(D(I)*COEF(I-1)+E(I)*COEF(I)+D(I+1)*COEF(I+1))
30        RES(NK)=B(NK)-(D(NK)*COEF(NK-1)+E(NK)*COEF(NK))
C
C         SOLVE FOR RES(.).
C
          CALL COEFF(RES,E,RES,D,UE,NK,I1,N)
C
C         CORRECT COEF(.).
C
          DO 40 I=I0,NK
40           COEF(I)=COEF(I)+RES(I)
C
C.... GENERATE SPLINE SX(.) AT X(.) FROM COEF(.) OVER AKN(.) ....
C
50    CALL SXATXI(SX,X,COEF,AKN,IAKN,I1,N,NK)
C
C.... CALCULATE THE LAGRANGE MULTIPLIERS PAR(.) ....
C
      IF (NACT.LE.I0) GO TO 270
      CALL LAGRNG(PAR,SX,F,IACT,NACT,EC,DC,UE,I1,N)
C
C.... INITIATE ITERATIVE REFINEMENT IF ITREF=1 ....
C
      IF (ITREF.NE.1) GO TO 270
C
C         REFINE PAR(.).
C         SET RESIDUALS OF LAGRANGIAN EQUATIONS IN RES(.).
C
          IF (I1 .EQ. NACT) THEN
              II=IACT(I1)
              RES(II)=SX(II+1)-F(II+1) - EC(II)*PAR(II)
              GO TO 200
          END IF
C
          II=IACT(I1)
          RES(II)=SX(II+1)-F(II+1)-(EC(II)*PAR(II)+DC(II+1)*PAR(II+1))
          IF (I1+1 .EQ. NACT) GO TO 190
          DO 180 I=I1+1,NACT-1
             II=IACT(I)
180          RES(II)=SX(II+1)-F(II+1) -
     1          (DC(II)*PAR(II-1)+EC(II)*PAR(II)+DC(II+1)*PAR(II+1))
190       II=IACT(NACT)
          RES(II)=SX(II+1)-F(II+1)-(DC(II)*PAR(II-1)+EC(II)*PAR(II))
200       CONTINUE
C
C         SOLVE FOR RES(.) AND CORRECT PAR(.).
C          THE FOLLOWING CODE DUPLICATES WITH THAT IN SUBROUTINE LAGRNG.
C
          I=I1
210       K=1
220       IF (I.LT.NACT .AND. IACT(I).EQ.IACT(I+1)-1) THEN
              I=I+1
              K=K+1
              GO TO 220
          END IF
C
230       IF (K.EQ.1) THEN
              II=IACT(I)
              PAR(II)=PAR(II)+RES(II)
          ELSE
              UE(I-K+1)=EC(IACT(I-K+1))
              DO 240 L=I-K+2,I
                 DI=DC(IACT(L))
240              UE(L)=EC(IACT(L))-DI/UE(L-1)*DI
C
              II=IACT(I-K+1)
              RES(II)=RES(II)
              DO 250 L=I-K+2,I
                 II=IACT(L)
250              RES(II)=RES(II)-DC(II)*RES(IACT(L-1))/UE(L-1)
C
              RES(IACT(I))=RES(IACT(I))/UE(I)
              PAR(IACT(I))=PAR(IACT(I))+RES(IACT(I))
              DO 260 L=I-1,I-K+1,-1
                 RES(IACT(L))=RES(IACT(L))/UE(L)-DC(IACT(L+1))*
     1                     RES(IACT(L+1))/UE(L)
260              PAR(IACT(L))=PAR(IACT(L))+RES(IACT(L))

          END IF
          I=I+1
          IF (I-NACT) 210,210,270
C
270   RETURN
      END
C......................................................................
C     THIS SUBROUTINE CALCULATES THE COEFFICIENTS COEF(.) OF
C      A LEAST SQUARES LINEAR SPLINE FROM THE NORMAL EQUATIONS.
C......................................................................
      SUBROUTINE COEFF(COEF,E,B,D,UE,NK,I1,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION COEF(I1-1:N),E(I1-1:N),B(I1-1:N),D(I1-1:N),UE(I1-1:N)
C
C.... I N P U T ....
C     E,B,D,I1,N,NK   AS THEY ARE DEFINED IN SUBROUTINE L2CXFT.
C
C.... O U T P U T ....
C     COEF(I1-1:NK) REAL ARRAY CONTAINING THE SPLINE COEFFICIENTS
C             AT AKN(I), I=I1-1,I1,...,NK.
C
C.... W O R K I N G  S P A C E ....
C     UE(I1-1:NK+1) REAL ARRAY THAT PROVIDES WORKING SPACE FOR THE
C             FACTORIZATION THAT GIVES COEF(.). IT KEEPS THE
C             THE DIAGONAL ELEMENTS OF THE FORWARD SWEEP OF
C             THE FACTORIZATION.
C
C.... M E T H O D (CALCULATE COEFFICIENTS) ....
C     EXPRESS SPLINE SX(.) AS A LINEAR COMBINATION OF B-SPLINES
C             SX(I)=SUM(COEF(J)*N(J,I), J=I1-1,I1,...,NK)
C                                          FOR I=I1,I1+1,...,N
C     WHERE N(J,I)=THE I-TH COMPONENT OF THE J-TH LINEAR B-SPLINE.
C     DENOTE BY <.,.> THE INNER PRODUCT OF TWO VECTORS AND,
C     THE SPLINE COEFFICIENTS COEF(.) ARE DETERMINED BY SOLVING THE
C     NORMAL EQUATIONS
C             D(J)*COEF(J-1)+E(J)*COEF(J)+D(J+1)*COEF(J+1)=B(J)
C                                          FOR J=I1-1,I1,...,NK
C     WHERE
C             D(J)=<N(J,.),N(J-1,.)>, J=I1-1,I1,...,NK
C             E(J)=<N(J,.),N(J,.)>, J=I1-1,I1,...,NK
C             B(J)=<N(J,.),F(.)>, J=I1-1,I1,...,NK.
C......................................................................
      I0=I1-1
C
C.....DERIVE COEF(.) BY FACTORIZATION
C
      UE(I0)=E(I0)
      COEF(I0)=B(I0)
      DO 10 I=I0+1,NK
         UE(I)=E(I)-D(I)/UE(I-1)*D(I)
10       COEF(I)=B(I)-D(I)*COEF(I-1)/UE(I-1)
      COEF(NK)=COEF(NK)/UE(NK)
      DO 20 I=NK-1,I0,-1
20       COEF(I)=(COEF(I)-D(I+1)*COEF(I+1))/UE(I)
C
      RETURN
      END
C.......................................................................
C     THIS SUBROUTINE DELETES THE KNOT AKN(KNOUT) OF THE LINEAR
C      SPLINE APPROXIMATION SX(.) DEFINED IN SUBROUTINE L2CXFT.
C
C     CALLS SUBROUTINE MINUS1.
C.......................................................................
      SUBROUTINE DELKNT(E,B,D,AKN,IAKN,IACT,NK,NACT,I1,N,KNOUT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E(I1-1:N),B(I1-1:N),D(I1-1:N),AKN(I1-2:N),IAKN(I1-1:N),
     1  IACT(I1:N)
C
C.... I N P U T ....
C     KNOUT  INTEGER VARIABLE, SUCH THAT AKN(KNOUT) IS THE KNOT TO
C             BE DELETED FROM THE KNOT SEQUENCE AKN(.).
C     E,B,D,AKN,IAKN,IACT,I1,N,NK,NACT  AS THEY ARE DEFINED IN
C             SUBROUTINE L2CXFT.
C
C.... O U T P U T ....
C     E,B,D,AKN,IAKN,IACT  AFTER DELETING THE KNOT AKN(KNOUT).
C     NK      DECREASED BY ONE.
C     NACT    INCREASED BY ONE.
C
C.... M E T H O D ....
C     SINCE AKN(J) IS DELETED THE J-TH NORMAL EQUATION IS DELETED TOO.
C      THEN CALL SUBROUTINE MINUS1 TO MODIFY THE (J-1) AND (J+1) NORMAL
C      EQUATIONS. SHIFT E,B,D,AKN AND IAKN.
C     REVISE IACT(.), THE INDICES OF THE ACTIVE CONSTRAINTS.
C.......................................................................
      I0=I1-1
      J=IAKN(KNOUT)-1
C
C     MODIFY NORMAL EQUATIONS.
C
      CALL MINUS1(E,B,D,AKN,I1,N,KNOUT)
C
C     SHIFT ARRAYS.
C
      DO 5 I=KNOUT,NK-1
         AKN(I)=AKN(I+1)
         IAKN(I)=IAKN(I+1)
         E(I)=E(I+1)
         B(I)=B(I+1)
5        D(I)=D(I+1)
C
C     SET LAST POSITION TO ZERO.
C
      AKN(NK)=0.
      IAKN(NK)=0
      E(NK)=0.
      B(NK)=0.
      D(NK)=0.
      NK=NK-1
C
C     UPDATE IACT(.).
C
C     GIVEN J=IAKN(KNOUT)-1, APPLY BISECTION IN [I1,NACT] IN ORDER
C      TO FIND IR SUCH THAT  IACT(IR-1) < J < IACT(IR) , OR
C      IQ IS OUTSIDE [I1,NACT].
C
      IF (NACT.EQ.I0) THEN
          IR=I1
          GO TO 130
      END IF
      IF (J.LT.IACT(I1)) THEN
          IR=I1
          GO TO 110
      END IF
      IF (J.GT.IACT(NACT)) THEN
          IR=NACT+1
          GO TO 130
      END IF
C
C     APPLY BISECTION.
C
      IP=I1
      IR=NACT
100   IQ=(IP+IR)/2
      IF (IR-IP.EQ.1) GO TO 110
      IF (J.LT.IACT(IQ)) THEN
          IR=IQ
      ELSE
          IP=IQ
      END IF
      GO TO 100
110   CONTINUE
      DO 120 I=NACT,IR,-1
120      IACT(I+1)=IACT(I)
130   IACT(IR)=J
      NACT=NACT+1
C
      RETURN
      END
C......................................................................
C     L2CXFT: A FORTRAN SUBROUTINE FOR LEAST SQUARES DATA FITTING
C             WITH NON-NEGATIVE SECOND DIVIDED DIFFERENCES.
C     BY I.C. DEMETRIOU, DEPARTMENT OF ECONOMICS, UNIVERSITY OF ATHENS,
C             8 PESMAZOGLOU STREET, ATHENS 105 59, GREECE.
C             TEL   : (GR) 01-3242438, 8017732
C             E-MAIL: DEMETRI@AUEB.ARIADNE-T.GR
C......................................................................
C     THE USER GIVES THE NOISY FUNCTION MEASUREMENTS F(I): I=1,2,...N
C      AT THE ABSCISSAE X(1)<X(2)<...<X(N),
C     AND THIS SUBROUTINE CALCULATES THE N-TH COMPONENT VECTOR SX(.)
C
C           THAT MINIMIZES THE SUM OF THE SQUARES OF THE ERRORS
C            SX(I)-F(I): I=1,2,...,N
C           SUBJECT TO NON-NEGATIVITY OF SECOND DIVIDED DIFFERENCES
C            OF THE SMOOTHED VALUES SX(.).
C
C     IN THIS SENSE, SX(.) IS
C                   THE BEST CONVEX FIT TO F(.)
C                 OR THE CONVEX REGRESSION ON F(.)
C
C     THERE ARE EXACTLY N-2 CONSTRAINTS THAT GIVE A HIGHLY STRUCTURED
C      CALCULATION WHICH IS SOLVED BY A SPECIAL DUAL FEASIBLE QUADRATIC
C      PROGRAMMING ALGORITHM INITIALIZED BY A O(N) STARTING PROCEDURE.
C
C     THE CONCAVE FITTING CAN BE TREATED AS CONVEX AFTER AN OVERALL
C      CHANGE OF SIGN.
C
C     REF:  I.C.DEMETRIOU AND M.J.D. POWELL, "THE MINIMUM SUM OF
C           SQUARES CHANGE TO UNIVARIATE DATA THAT GIVES CONVEXITY",
C           IMA J. OF NUMERICAL ANALYSIS, 11(1991), PP. 433-448.
C
C     APPLICATIONS: DATA SMOOTHING, STATISTICAL ANALYSIS, ECONOMICS.
C
C     CALLS SUBROUTINES ADDKNT, CFLGRF, DELKNT, L2ILS AND MESSGE.
C.......................................................................
      SUBROUTINE L2CXFT(X,F,SX,I1,N,AKN,IAKN,COEF,NK,IACT,NACT,PAR,
     1  ITER,SXOLD,E,B,D,Z,U,DC,EC,PAROLD,RES,IPRINT,MODE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(I1:N),F(I1:N),SX(I1:N),COEF(I1-1:N),AKN(I1-2:N),
     1  IAKN(I1-1:N),IACT(I1:N),PAR(I1:N),DC(I1:N),EC(I1:N),
     2  E(I1-1:N),B(I1-1:N),D(I1-1:N),Z(I1-1:N),U(I1-1:N),
     3  SXOLD(I1:N),PAROLD(I1:N),RES(I1-1:N),ITER(*)
C
C..... I N P U T .....
C      I1      INTEGER VARIABLE, LOWER DATA INDEX.
C      N       INTEGER VARIABLE, UPPER DATA INDEX.
C      X(I1:N) REAL ARRAY OF ABSCISSAE.
C      F(I1:N) REAL ARRAY OF FUNCTION MEASUREMENTS.
C      IPRINT  INTEGER VARIABLE. IF IPRINT=1 THEN SUBROUTINE MESSGE
C               PRODUCES OUTPUT.
C
C .....O U T P U T .....
C      SX(I1:N) REAL ARRAY CONTAINING THE SPLINE COMPONENTS AT X(I),
C                I=I1,I1+1,...,N IE THE BEST CONVEX FIT TO F(.).
C      NK      INTEGER VARIABLE, UPPER INDEX OF SPLINE KNOTS SEQUENCE.
C      NACT    INTEGER VARIABLE, SUCH THAT EITHER NACT=I1-1 OR
C               (NACT-I1+1) IS EQUAL TO THE NUMBER OF ACTIVE CONSTRAINTS
C               IN THE SOLUTION.
C      AKN(I1-2:NK)  REAL ARRAY CONTAINING THE SPLINE KNOTS AND THE
C               DUMMY COMPONENT AKN(I1-2)=X(I1).
C      IAKN(I1-1:NK)  INTEGER ARRAY CONTAINING THE DATA POINT INDICES
C               OF THE ASSOCIATED KNOTS.
C      COEF(I1-1:NK)  REAL ARRAY CONTAINING THE SPLINE COEFFICIENTS
C               AT AKN(I), I=I1-1,I1,...,NK.
C      IACT(I1:NACT)  INTEGER ARRAY CONTAINING THE INDICES OF THE
C              ACTIVE CONSTRAINTS AT THE SOLUTION.
C      PAR(I1:N-2)  INTEGER ARRAY CONTAINING IN PAR(IACT(.))
C               THE LAGRANGE MULTIPLIERS OF THE ACTIVE CONSTRAINTS.
C      ITER(*) INTEGER ARRAY WHOSE ELEMENTS ARE AS FOLLOWS:
C              ITER(1)  KEEPS THE NUMBER OF ITERATIONS AT THE END
C                        OF THE CALCULATION. IF ITER(1) EXCEEDS
C                        THE NUMBER 3*(N-I1+1) THEN TERMINATION OCCURS
C                        WITH THE MESSAGE
C                         'L2CXFT QUITS ON ITERATION COUNT'
C              ITER(2)  KEEPS THE NUMBER OF ACTIVE CONSTRAINTS AT THE
C                        END OF THE CALL OF SUBROUTINE L2ILS.
C              ITER(3)  KEEPS THE NUMBER OF CALLS OF SUBROUTINE
C                        CFLGRF.
C              ITER(4)  KEEPS THE NUMBER OF ACTIVE CONSTRAINTS WHEN
C                        FIRST FEASIBLE LAGRANGIAN VECTOR IS FOUND
C                        JUST AFTER THE CALL OF SUBROUTINE L2ILS.
C              ITER(5)  KEEPS THE NUMBER OF CONSTRAINT DELETIONS FROM
C                        THE ACTIVE SET.
C              ITER(6)  KEEPS THE NUMBER OF CONSTRAINT ADDITIONS TO
C                        THE ACTIVE SET AFTER FIRST FEASIBLE PAR(.).
C      MODE    INTEGER VARIABLE THAT IS SET AUTOMATICALLY BY THE
C               SUBROUTINE TO INDICATE THE TERMINATION STATUS.
C              1 = SUCCESSFUL TERMINATION
C              2 = TERMINATION DUE TO NON-INCREASE OF OBJECTIVE FUNCTION
C              3 = NUMBER OF 3*(N-I1+1) ITERATIONS HAS EXCEEDED
C              4 = ERROR RETURN BECAUSE (N .LT. I1).
C
C.... M E S S A G E S ....
C     INCLUDED IN SUBROUTINE MESSGE AND ACTIVATED WHEN INPUT PARAMETER
C      IPRINT IS SET BY THE USER TO ONE.
C
C.... W O R K I N G   S P A C E .....
C      E(I1-1:NK)  REAL ARRAY THAT HOLDS THE DIAGONAL OF THE NORMAL
C              EQUATIONS ASSOCIATED WITH SPLINE SX(.).
C      D(I1-1:NK)  REAL ARRAY THAT HOLDS THE SUBDIAGONAL OF THE
C              NORMAL EQUATIONS ASSOCIATED WITH SX(.).
C      B(I1-1:NK)  REAL ARRAY THAT HOLDS THE RIGHT HAND SIDE OF
C              THE NORMAL EQUATIONS.
C      EC(I1:N-2)  REAL ARRAY SUCH THAT EC(I)=(SECOND COMPONENT OF
C              THE I-TH DIVIDED DIFFERENCE AFTER SCALING).
C      DC(I1:N-2)  REAL ARRAY SUCH THAT DC(I+1)=(THIRD COMPONENT OF
C              THE I-TH DIVIDED DIFFERENCE AFTER SCALING).
C      PAROLD(I1:N-2)  REAL ARRAY THAT HOLDS THE LAGRANGE MULTIPLIERS
C              THAT OCCURED AT THE FORMER ITERATION. IT IS
C              PAROLD(IACT(.)).GT.0. .
C      SXOLD(I1:N) REAL ARRAY THAT HOLDS THE SPLINE COMPONENTS
C              THAT OCCURED AT THE FORMER ITERATION.
C      RES(I1-1:NK)  REAL ARRAY THAT IS USED TO KEEP THE RESIDUALS OF
C              THE NORMAL EQUATIONS AT COEF(.) AND THE CORRECTIONS OF
C              COEF(.). IT IS ALSO USED TO KEEP THE RESIDUALS OF THE
C              LANGRANGIAN EQUATIONS AT PAR(.) AND THE CORRECTIONS TO
C              PAR(.).
C      U(I1-1,N) REAL ARRAY PROVIDING SPACE FOR FACTORIZATIONS.
C      Z(I1-1,N) REAL ARRAY PROVIDING SPACE FOR FACTORIZATIONS.
C      ITREF   INTEGER VARIABLE. IF ITREF=1 THEN ITERATIVE REFINEMENT
C               IS INITIATED.
C      ITEST   INTEGER VARIABLE. IF ITEST=1 THEN A RETURN IS CAUSED DUE
C               TO NOT INCREASING OF THE FUNCTION VALUES.
C
C .....M E T H O D (DUAL-FEASIBLE ALGORITHM 3 IN QUOTED REFERENCE).
C      VECTOR SX(.) IS CALCULATED BY THIS SUBROUTINE AS FOLLOWS:
C        STEP 0   CALL SUBROUTINE L2ILS TO DERIVE AN INITIAL
C                  APPROXIMATION TO SX(.).
C        STEP 1   USE THIS APPROXIMATION AS STARTING POINT OF A
C                  SPECIAL QUADRATIC PROGRAMMING (QP) METHOD THAT
C                  COMPLETES THE CALCULATION OF SX(.).
C
C      FEATURES OF THE QP METHOD:
C              A) IT GENERATES INFEASIBLE ESTIMATES OF THE OPTIMUM,
C              B) IT KEEPS A SET OF ACTIVE CONSTRAINTS THAT ARE
C                  ASSOCIATED WITH NON-NEGATIVE LAGRANGE MULTIPLIERS,
C              C) THE OBJECTIVE FUNCTION ACHIEVES A STRICTLY LARGER
C                  VALUE ANY TIME A NEW ITERATION BEGINS.
C              D) POWELL'S SAFEGUARD TERMINATION TEST IS APPLIED
C                  WHEN ACCURACY IS INSUFFICIENT TO MAINTAIN INCREASING
C                  FUNCTION VALUES.
C              E) ITERATIVE REFINEMENT IS CALLED AUTOMATICALLY IN
C                  CONJUNCTION WITH D) IN ORDER TO IMPROVE ACCURACY.
C
C      IMPLEMENTATION OF THE QP METHOD:
C      A) CALCULATION OF AN ESTIMATE SX(.) OF THE OPTIMUM:
C          EXPRESS SPLINE SX(.) AS A LINEAR COMBINATION OF LINEAR
C          B-SPLINES
C              SX(I)=SIGMA(COEF(J)*N(J,I), J=I1-1,I1,...,NK)
C                                           FOR I=I1,I1+1,...,N ,
C          WHERE N(J,I)=THE I-TH COMPONENT OF THE J-TH LINEAR B-SPLINE.
C          DENOTE BY <.,.> THE INNER PRODUCT OF TWO VECTORS AND,
C          THE COEFFICIENTS OF SX(.) ARE DETERMINED BY SOLVING THE
C          THE POSITIVE DEFINITE SYSTEM OF THE NORMAL EQUATIONS
C          EFFICIENTLY AND STABLY BY CHOLESKY FACTORIZATION
C              D(J)*COEF(J-1)+E(J)*COEF(J)+D(J+1)*COEF(J+1)=B(J)
C                                           FOR J=I1-1,I1,...,NK
C          WHERE BY <.,.> IS DENOTED THE SCALAR PRODUCT OF TWO VECTORS
C          AND
C              D(J)=<N(J,.),N(J-1,.)>, J=I1,I1+1,...,NK
C              E(J)=<N(J,.),N(J,.)>, J=I1-1,I1,...,NK
C              B(J)=<N(J,.),F(.)>, J=I1-1,I1,...,NK.
C          (SUBROUTINE CFLGRF).
C
C      B) CALCULATION OF THE LAGRANGE MULTIPLIERS
C          THE LAGRANGE EQUATIONS AT AN ESTIMATE SX(.) ARE
C              2*(SX(II+1)-F(II+1))=SIGMA(PAR(II)*A(II,.))
C                                    FOR II=IACT(I), I=I1,I1+1,...,NACT
C          WHERE A(II,J) IS THE J-TH COMPONENT OF THE NORMAL
C          OF THE II-TH CONSTRAINT AND
C              A(II,J)<>0 ONLY FOR (II.LE.J.LE.II+2).
C          (ARRAY A(.,.) IS NOT ACTUALLY USED IN THIS SUBROUTINE;
C          INSTEAD ARRAYS EC(.) AND DC(.) ARE EMPLOYED.). BECAUSE
C          THIS SYSTEM IS OVERDETERMINED PAR(IACT(.)) IS DERIVED
C          BY CHOOSING FROM THE LAGRANGE EQUATIONS THOSE EQUATIONS THAT
C          CORRESPOND TO THE SECOND NON-ZERO COMPONENT OF A(IACT(.),.);
C          SO A BLOCK DIAGONAL POSITIVE DEFINITE SYSTEM IS OBTAINED THAT
C          IS SOLVED EFFICIENTLY AND STABLY BY CHOLESKY FACTORIZATION.
C          (SUBROUTINE CFLGRF).
C
C      C) PARTICULAR METHODS ARE DEVELOPED FOR INSERTING INTO AND
C          DELETING FROM THE SPLINE APPROXIMATION ONE KNOT ECONOMICALLY.
C          (SUBROUTINES ADDKNT AND DELKNT).
C......................................................................
C     SET CERTAIN INITIAL VALUES.
C
      MODE=1
      I0=I1-1
      ITREF=0
      ITEST=0
      ZERO=0.
      DO 10 I=1,10
10       ITER(I)=0
C
      IF (N.LT.I1) THEN
          MODE=4
          CALL MESSGE(ZERO,0,0,ITER,10)
          GO TO 1000
      END IF
C
C     IF N=I1 OR N=I1+1 THEN THE DATA ITSELF IS THE REQUIRED FIT.
C
      IF (N.EQ.I1) THEN
          SX(I1)=F(I1)
          AKN(I1-1)=X(I1)
          IAKN(I1-1)=I1
          COEF(I1-1)=F(I1)
          CALL MESSGE(ZERO,0,0,ITER,20)
          GO TO 1000
      END IF
C
      IF (N.EQ.I1+1) THEN
          SX(I1)=F(I1)
          SX(N)=F(N)
          AKN(I1-1)=X(I1)
          AKN(I1)=X(N)
          IAKN(I1-1)=I1
          IAKN(I1)=N
          COEF(I1-1)=F(I1)
          COEF(I1)=F(N)
          CALL MESSGE(ZERO,0,0,ITER,20)
          GO TO 1000
      END IF
C
C.... I N I T I A L I Z E .............................................
C
C     THERE ARE N-2 CONSTRAINTS DEFINED BY THE DIVIDED DIFFERENCES.
C      WE ASSOCIATE THE I-TH CONSTRAINT WITH THE I-TH DIVIDED DIFFER-
C      ENCE. SET COEFFICIENTS OF FULL SYSTEM OF EQUATIONS OF LAGRANGE
C      MULTIPLIERS. SINCE GRAD(F)=2*(SX-F), THE COEFFICIENTS ARE
C      SCALED BY 2.
C
      DO 50 I=I1,N-3
         CX2=X(I+2)-X(I+1)
         CX1=X(I+1)-X(I)
         CX3=X(I+2)-X(I)
         CX2=1/CX2
         CX1=1/CX1
         EC(I)=-(CX1+CX2)/2
50       DC(I+1)=CX2/2
      CX2=X(N)-X(N-1)
      CX1=X(N-1)-X(N-2)
      CX3=X(N)-X(N-2)
      CX2=1/CX2
      CX1=1/CX1
      EC(N-2)=-(CX1+CX2)/2
C
      DO 60 I=I1,N-2
60       PAR(I)=0.0
C
C.... CALCULATE STARTING POINT BY L2ILS ...............................
C
      CALL L2ILS(COEF,AKN,IAKN,D,E,B,Z,U,NK,F,X,I1,N)
C
C     INITIALIZE INDICES OF ACTIVE CONSTRAINTS IN IACT(.).
C
      J=I0
      I=I1+1
      IAKN(I0)=I1
      KL=I0
80    IF (IAKN(J).LT.I .AND. I.LT.IAKN(J+1)) THEN
          KL=KL+1
          IACT(KL)=I-1
          I=I+1
          GO TO 80
      END IF
      J=J+1
      I=I+1
      IF (J.LT.NK) GO TO 80
      NACT=KL
      ITER(2)=NACT-I0
C
C.... ITERATE BY ONE, UNTIL FIRST FEASIBLE LAGRANGIAN VECTOR IS FOUND.
C
C     CALCULATE COEF OVER KNOTS, SPLINE OVER ABSCISSAE AND
C      LAGRANGE MULTIPLIERS.
C
90    CALL CFLGRF(SX,PAR,COEF,IACT,E,B,D,X,AKN,IAKN,F,DC,EC,U,
     1  RES,NK,NACT,I1,N,ITREF,ITER(3))
      ITER(3)=ITER(3)+1
C
C     FIND MINIMUM LAGRANGE MULTIPLIER.
C
      IF (NACT .LE. I0) GO TO 150
      AMINLA=PAR(IACT(I1))
      IPAR=IACT(I1)
      IF (I1 .EQ. NACT) GO TO 110
      DO 100 I=I1+1,NACT
         IF (PAR(IACT(I)) .LT. AMINLA) THEN
             AMINLA=PAR(IACT(I))
             IPAR=IACT(I)
         END IF
100   CONTINUE
110   ITER(4)=NACT-I0
      IF (AMINLA .GE. 0.0) GO TO 150
C
C     INSERT A KNOT WHERE THE LEAST NEGATIVE LAGRANGE MULTIPLIER.
C
      INKNOT=IPAR+1
      PAR(IPAR)=0.0
      CALL ADDKNT(INKNOT,E,B,D,X,F,AKN,IAKN,IACT,NK,NACT,I1,N)
      GO TO 90
C
C.... M A I N  I T E R A T I O N ......................................
C
150   CONTINUE
C
C     EXAMINE CONVEXITES AT ALL KNOTS OF SPLINE APPROXIMATION.
C      SPECIAL CASE: IF NK-I0=1 THEN SX(.) IS A STRAIGHT LINE.
C
      IF (NK-I0 .EQ. 1) GO TO 1000
C
C     FIND THE KNOT OF THE MOST VIOLATED CONSTRAINT.
C
      CVMAX=0.
      KL=I0
      K=0
      DO 160 J=I1,NK-1
         JJ=IAKN(J)
         CKNOT=(SX(JJ+1)-SX(JJ))/(X(JJ+1)-X(JJ))-
     1         (SX(JJ)-SX(JJ-1))/(X(JJ)-X(JJ-1))
         SUM2=1/(X(JJ+1)-X(JJ))**2+
     1   ((X(JJ+1)-X(JJ-1))/((X(JJ+1)-X(JJ))*(X(JJ)-X(JJ-1))))**2+
     2   1/(X(JJ)-X(JJ-1))**2
         SUM2=SQRT(SUM2)
         CKNOT=CKNOT/SUM2
         IF (CKNOT .LT. 0.) K=K+1
         IF (CKNOT .LT. CVMAX) THEN
             CVMAX=CKNOT
             KNTOUT=J
             IAKNJ=IAKN(J)
             KL=KL+1
         END IF
160   CONTINUE
C
C     IF KL=I0 THEN SX(.) IS CONVEX. SINCE AT THIS LEVEL PAR(.) ARE
C      ALL NON-NEGATIVE THE KUHN-TUCKER CONDITIONS ARE SATISFIED.
C
C.... T E R M I N A T I O N ...........................................
C
      IF (KL .EQ. I0) GO TO 1000
C
C     MSG: PRINT CVMAX,IAKNJ,K.
C
      IF (IPRINT.NE.0) CALL MESSGE(CVMAX,K,IAKNJ,ITER,160)
C
C     POWELL'S SAFEGUARD TERMINATION TEST:
C      RETURN, IF DUE TO ROUNDING ERRORS, THE CHANGE IN SX(.) MAY
C      NOT INCREASE THE OBJECTIVE FUNCTION.
C
      IF (ITER(1).EQ.0) GO TO 190
          R=0.0
          T=0.0
          DO 170 I=I1,N
             R=R+(SX(I)-SXOLD(I))*(SX(I)+SXOLD(I)-2*F(I))
170          T=T+ABS(SX(I)-SXOLD(I))*(ABS(SX(I)+SXOLD(I))+2*ABS(F(I)))
          IF (T+R.LE.T .OR. T+1.5*R.LE.T+R) THEN
              IF (ITEST.EQ.1) THEN
                  MODE=2
                  IF (IPRINT.NE.0) CALL MESSGE(ZERO,0,0,ITER,170)
C
C                     MSG: L2CXFT RETURNS BECAUSE THE ACCURACY IS
C                          INSUFFICIENT TO MAINTAIN INCREASING
C                          THE FUNCTION VALUES.
C
                  GO TO 1000
              END IF
C
              IF (IPRINT.EQ.0) GO TO 180
              IF (ITREF.EQ.0) CALL MESSGE (ZERO,0,0,ITER,180)
C
C                 MSG: ITERATIVE REFINEMENT IS INITIATED.
C
180           ITREF=1
              ITEST=1
          ELSE
              ITEST=0
          END IF
190   CONTINUE
C
C     KEEP OLD LAGRANGE MULTIPLIERS AND OLD SPLINE APPROXIMATION.
C
      DO 200 I=I1,N-2
200      PAROLD(I)=PAR(I)
      DO 210 I=I1,N
210      SXOLD(I)=SX(I)
C
C     TERMINATION ON ITERATION COUNT.
C
      ITER(1)=ITER(1)+1
      IF (ITER(1).GT.3*(N-I0)) THEN
          MODE=3
          CALL MESSGE(ZERO,0,0,ITER,210)
C
C         MSG: L2CXFT QUITS ON ITERATION COUNT.
C
          GO TO 1000
      END IF
C
C     CALCULATE A SEARCH DIRECTION (DELETE THE KNOT KNTOUT).
C
      CALL DELKNT(E,B,D,AKN,IAKN,IACT,NK,NACT,I1,N,KNTOUT)
      ITER(6)=ITER(6)+1
C
C     SOLVE FOR SX(.) OVER AKN(.).
C
220   CALL CFLGRF(SX,PAR,COEF,IACT,E,B,D,X,AKN,IAKN,F,DC,EC,U,
     1  RES,NK,NACT,I1,N,ITREF,ITER(3))
      ITER(3)=ITER(3)+1
C
      IPAR=IAKNJ-1
      IF (IPRINT.NE.0) CALL MESSGE(PAR(IPAR),0,0,ITER,220)
C
C     MSG: PRINT MULTIPLIER OF INSERTED CONSTRAINT.
C
C     APPLY POWELL'S TECHNIQUE IN ORDER THAT THE INDEX IPAR
C      REMAIN IN ACTIVE SET. SET PAR(IPAR)=MAX(PAR(IPAR),0.0)
C      BECAUSE IN THEORY HOLDS THAT PAR(IPAR)>0.
C
      IF (PAR(IPAR).LT.0.0) PAR(IPAR)=0.0
C
C     EXAMINE SIGN OF LAGRANGE MULTIPLIERS
C
300   IF (NACT .EQ. I0) GO TO 150
      KL=I0
      THETA=2.0
      DO 320 I=I1,NACT
         K=IACT(I)
         IF (PAR(K)) 310,320,320
310      RATIO=PAROLD(K)/(PAROLD(K)-PAR(K))
         IF (THETA .LE. RATIO) GO TO 320
             THETA=RATIO
             IPAR=K
             KL=I1
320      CONTINUE
C
C.... ITERATE IN THE PRIMAL SPACE ......................................
C     IF KL=I0 THEN PAR(.) ARE NON-NEGATIVE. BRANCH WHERE TEST
C      PRIMAL FEASIBILITY.
C
      IF (KL .EQ. I0) GO TO 150
C
C.... RECOVER NON-NEGATIVITY OF LAGRANGE MULTIPLIERS ...................
C
      DO 340 I=I1,NACT
         K=IACT(I)
340      PAROLD(K)=(1-THETA)*PAROLD(K)+THETA*PAR(K)
C
C     INSERT KNOT IPAR+1 INTO THE SPLINE AND SOLVE FOR SX(.).
C
      I1NEW=I1
350   INKNOT=IPAR+1
      PAR(IPAR)=0.
      CALL ADDKNT(INKNOT,E,B,D,X,F,AKN,IAKN,IACT,NK,NACT,I1,N)
      ITER(5)=ITER(5)+1
C
C     TEST THE FEASIBILITY OF PAROLD(.). ALL PAROLD(.) MUST BE
C      NON-NEGATIVE BECAUSE OF THE WAY THETA WAS DETERMINED. OTHERWISE,
C      THIS IS DUE TO ROUND-OFF ERROR. HENCE, ANY NEGATIVITY
C      IS TREATED AS BEING ZERO. SEE LAWSON & HANSON (SOLVING
C      LEAST SQUARES PROBLEMS, PRENTICE-HALL, 1974, p.307).
C
      IF (NACT .EQ. I0) GO TO 370
      IF (I1NEW .GT. NACT) GO TO 370
      DO 360 I=I1NEW,NACT
         IPAR=IACT(I)
         IF (PAROLD(IPAR) .LT. 0.0) THEN
             IF (IPRINT.NE.0) CALL MESSGE(PAROLD(IPAR),IPAR,0,ITER,350)
C
C                MSG: ITERATION IN THE DUAL SPACE DUE TO ROUND OFF.
C                PRINT ACTIVE CONSTRAINT AND LAGRANGE MULTIPLIER.
C
             I1NEW=I
             GO TO 350
         END IF
360   CONTINUE
C
370   CONTINUE
C
      CALL CFLGRF(SX,PAR,COEF,IACT,E,B,D,X,AKN,IAKN,F,DC,EC,U,
     1  RES,NK,NACT,I1,N,ITREF,ITER(3))
      ITER(3)=ITER(3)+1
C
C.... ITERATE IN THE DUAL SPACE ........................................
C
      GO TO 300
1000  CONTINUE
C
      CALL MESSGE(ZERO,NACT-I0,MODE,ITER,1000)
C
C     MSG: PRINT MODE AND
C          INITIAL ACTIVE SET
C          STARTING ACTIVE SET
C          FINAL ACTIVE SET
C          ADDITIONS OF CONSTRAINTS
C          DELETIONS OF CONSTRAINTS
C
C.... END OF SUBROUTINE L2CXFT .........................................
C
      RETURN
      END
C......................................................................
C     GIVEN N NOISY FUNCTION MEASUREMENTS F(I): I=1,2,...N
C      AT THE ABSCISSAE X(1)<X(2)<...<X(N),
C     THIS SUBROUTINE CALCULATES IN ONLY O(N) COMPUTER OPERATIONS AN
C      EXCELLENT STARTING POINT FOR THE FOLLOWING DATA FIT CALCULATION:
C
C      FIND SX(.) THAT MINIMIZES THE SUM OF THE SQUARES OF THE ERRORS
C            SX(I)-F(I): I=1,2,...,N
C      SUBJECT TO NON-NEGATIVITY OF SECOND DIVIDED DIFFERENCES
C            OF THE SMOOTHED VALUES SX(.).
C
C     SX(.) APPEARS DUE A B-SPLINE REPRESENTATION (SEE METHOD BELOW).
C......................................................................
C     THIS CODE HAS SUBMITTED AS REF 1 SHOWS.
C     REF 1: "A O(N) STARTING POINT FORTRAN SUBROUTINE FOR LEAST SQUARES
C             DATA FITTING WITH NON-NEGATIVE SECOND DIVIDED DIFFERENCES"
C             BY I.C. DEMETRIOU, SUBMITTED TO ACM TRANSACTIONS ON
C             MATHEMATICAL SOFTWARE, MARCH 1990.
C     THE ALGORITHM OF THIS CODE CAN BE FOUND AS REF 2 SHOWS.
C     REF 2: "THE MINIMUM SUM OF SQUARES CHANGE TO UNIVARIATE DATA
C             THAT GIVES CONVEXITY"
C             BY I.C.DEMETRIOU AND M.J.D. POWELL, SUBMITTED TO IMA J. OF
C             NUMERICAL ANALYSIS, 1989.
C.......................................................................
        SUBROUTINE L2ILS(COEF,AKN,IAKN,D,E,B,Z,U,NK,F,X,I1,N)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION X(I1:N),F(I1:N),COEF(I1-1:N),AKN(I1-2:N),
     1    IAKN(I1-1:N),D(I1-1:N),E(I1-1:N),B(I1-1:N),Z(I1-1:N),U(I1-1:N)
C
C ..... I N P U T .....
C       I1      INTEGER VARIABLE, LOWER DATA INDEX.
C       N       INTEGER VARIABLE, UPPER DATA INDEX.
C       X(I1:N) REAL ARRAY OF ABSCISSAE.
C       F(I1:N) REAL ARRAY OF MEASUREMENTS.
C
C ..... O U T P U T .....
C       NK      INTEGER VARIABLE, UPPER INDEX OF SPLINE KNOTS SEQUENCE.
C       AKN(I1-2:NK)  REAL ARRAY CONTAINING THE SPLINE KNOTS, WHERE
C               AKN(I1-2)=AKN(I1-1).
C       IAKN(I1-1:NK)  INTEGER ARRAY CONTAINING THE DATA POINT INDICES
C               OF THE ASSOCIATED KNOTS.
C       COEF(I1-1:NK)  REAL ARRAY CONTAINING THE SPLINE COEFFICIENTS
C               AT AKN(I), I=I1-1,I1,...,NK.
C
C ..... W O R K I N G   S P A C E .....
C       E(I1-1:N)  REAL ARRAY THAT HOLDS THE DIAGONAL OF THE NORMAL
C               EQUATIONS ASSOCIATED WITH SPLINE SX(.).
C       D(I1-1:N)  REAL ARRAY THAT HOLDS THE SUBDIAGONAL OF THE
C               NORMAL EQUATIONS ASSOCIATED WITH SX(.).
C       U(I1-1:N)  REAL ARRAY THAT HOLDS THE DIAGONAL OF CHOLESKY
C               FACTORIZATION OF THE NORMAL EQUATIONS.
C       B(I1-1:N)  REAL ARRAY THAT HOLDS THE RIGHT HAND SIDE OF
C               THE NORMAL EQUATIONS.
C       Z(I1-1:N)  REAL ARRAY THAT HOLDS THE RIGHT HAND SIDE OF
C               THE CHOLESKY FACTORIZATION.
C
C ..... M E T H O D (ALGORITHM 1 IN QUOTED REFERENCE) .....
C       EXPRESS SPLINE SX(.) AS A LINEAR COMBINATION OF B-SPLINES
C               SX(T)=SUM(COEF(J)*N(I,T), I=I1-1,I1,...,J)
C                                            FOR T=I1,I1+1,...,N
C       WHERE N(I,T)=THE T-TH COMPONENT OF THE I-TH LINEAR B-SPLINE.
C       DENOTE BY <.,.> THE INNER PRODUCT OF TWO VECTORS AND,
C       THE COEFFICIENTS OF SX(.) ARE DETERMINED BY SOLVING THE
C       NORMAL EQUATIONS
C               D(I)*COEF(I-1)+E(I)*COEF(I)+D(I+1)*COEF(I+1)=B(I)
C                                            FOR I=I1-1,I1,...,J
C       WHERE
C               D(I)=<N(I,.),N(I-1,.)>, I=I1-1,I1,...,J
C               E(I)=<N(I,.),N(I,.)>, I=I1-1,I1,...,J
C               B(I)=<N(I,.),F(.)>, I=I1-1,I1,...,J.
C
C       THEN PROCEED AS IN REFERENCE QUOTED.
C......................................................................
C       NO TEST ON THE CONDITIONS N<I1, N=I1 AND N=I1+1 IS MADE
C        HERE BECAUSE THEY ARE CONSIDERED BY SUBROUTINE L2CXFT.
C
C       SET INITIAL VALUES.
C
        I0=I1-1
100     K=I0
        J=I1
        AKN(I0)=X(I1)
        AKN(I0-1)=AKN(I0)
        AKN(I1)=X(I1+1)
        AKN(I1+1)=AKN(I1)
        IAKN(I0)=I1
        IAKN(I1)=I1+1
C
C       SET NORMAL EQUATIONS.
C
        E(I0)=1.
        E(I1)=1.
        D(I1)=0.
        B(I0)=F(I1)
        B(I1)=F(I1+1)
C
C       SET CHOLESKY FACTORIZATION AND COEFFICIENTS.
C
        U(I0)=1.
        U(I1)=1.
        Z(I0)=B(I0)
        Z(I1)=B(I1)
        COEF(I0)=Z(I0)
        COEF(I1)=Z(I1)
110     IF (K.EQ.N-2) GOTO 170
C
C       INSERT A KNOT AT X(K+2).
C
        J=J+1
        AKN(J)=X(K+3)
        IAKN(J)=K+3
        AKN(J+1)=AKN(J)
C
C       UPDATE NORMAL EQUATIONS AND CHOLESKY FACTORIZATION.
C
        E(J)=1.
        D(J)=0.
        B(J)=F(K+3)
        U(J)=1.
        Z(J)=B(J)
C
C       EXAMINE RIGHTMOST KNOT CONVEXITY.
C
120     IF (J.EQ.I1) THEN
            K=K+1
            GO TO 110
        END IF
C
C       CALCULATE RIGHTMOST THREE COEFFICIENTS.
C
        COEF(J)=Z(J)/U(J)
        DO 130 I=J-1,J-2,-1
130        COEF(I)=(Z(I)-D(I+1)*COEF(I+1))/U(I)
C
C       DETERMINE SIGN OF RIGHTMOST SECOND DIVIDED DIFFERENCE.
C
140     CAKNOT=(COEF(J)-COEF(J-1))/(AKN(J)-AKN(J-1))-
     1  (COEF(J-1)-COEF(J-2))/(AKN(J-1)-AKN(J-2))
        IF (CAKNOT.GT.0.0) THEN
            K=K+1
            GOTO 110
        END IF
C
C       DELETE KNOT AKN(J-1)
C       UPDATE NORMAL EQUATIONS AND CHOLESKY FACTORIZATION.
C
        C1=(AKN(J-1)-AKN(J-2))/(AKN(J)-AKN(J-2))
        C2=(AKN(J)-AKN(J-1))/(AKN(J)-AKN(J-2))
        E(J-2)=E(J-2)+C2*C2*E(J-1)+2*C2*D(J-1)
        D(J-1)=C1*D(J-1)+C2*D(J)+C1*C2*E(J-1)
        E(J-1)=C1*C1*E(J-1)+E(J)+2*C1*D(J)
        B(J-2)=B(J-2)+C2*B(J-1)
        B(J-1)=C1*B(J-1)+B(J)
C
        IF (J.EQ.I0+2) THEN
            U(J-2)=E(J-2)
            Z(J-2)=B(J-2)
        ELSE
            U(J-2)=E(J-2)-D(J-2)*D(J-2)/U(J-3)
            Z(J-2)=B(J-2)-D(J-2)*Z(J-3)/U(J-3)
        END IF
        U(J-1)=E(J-1)-D(J-1)*D(J-1)/U(J-2)
        Z(J-1)=B(J-1)-D(J-1)*Z(J-2)/U(J-2)
        U(J)=0.
        Z(J)=0.
        E(J)=0.
        D(J)=0.
        B(J)=0.
        J=J-1
        AKN(J)=X(K+3)
        IAKN(J)=K+3
        AKN(J+1)=AKN(J)
        GOTO 120
C
C       DERIVE COEFFICIENTS BY BACK SUBSTITUTION.
C
170     NK=J
        COEF(NK)=Z(NK)/U(NK)
        DO 180 I=NK-1,I0,-1
180        COEF(I)=(Z(I)-D(I+1)*COEF(I+1))/U(I)
C
        RETURN
        END
C.......................................................................
C     THIS SUBROUTINE CALCULATES THE LAGRANGE MULTIPLIERS PAR(.)
C      ASSOCIATED WITH THE SPLINE SX(.) IN SUBROUTINE L2CXFT.
C.......................................................................
      SUBROUTINE LAGRNG(PAR,SX,F,IACT,NACT,EC,DC,UE,I1,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION SX(I1:N),F(I1:N),DC(I1:N),EC(I1:N),UE(I1-1:N),
     1  IACT(I1:N),PAR(I1:N)
C
C.... I N P U T ....
C     SX,F,IACT,DC,EC,I1,N,NACT  AS THEY ARE DEFINED IN SUBROUTINE
C             L2CXFT.
C
C.... O U T P U T ....
C     PAR(I1:N-2) REAL ARRAY CONTAINING THE LAGRANGE MULTIPLIERS SO
C             THAT PAR(IACT(I)), I=I1,I1+1,...,NACT ARE THE MULTI-
C             PLIERS OF THE ACTIVE CONSTRAINTS WHILE THE REMAINING
C             MULTIPLIERS EQUAL TO ZERO.
C
C.... W O R K I N G  S P A C E ....
C     UE(I1-1:NK+1) REAL ARRAY THAT PROVIDES WORKING SPACE FOR THE
C             FACTORIZATION THAT GIVES PAR(.). IT KEEPS THE
C             DIAGONAL ELEMENTS OF THE FORWARD SWEEP OF THE
C             FACTORIZATION.
C
C.... M E T H O D (CALCULATE THE LAGRANGE MULTIPLIERS) .....
C     THE FIRST ORDER OPTIMALITY CONDITIONS ARE
C             SUM(PAR(II)*A(II,J))=2*(SX(II+1)-F(II+1)),
C                                   FOR II=IACT(I), I=I1,...,NACT
C      WHERE A(II,J) DENOTES THE J-TH COMPONENT OF THE II-TH SECOND
C      DIVIDED DIFFERENCE. ONLY THE A(II,J) FOR J=II,II+1,II+2 ARE
C      NON-ZERO. THESE CONDITIONS CONSTITUTE AN OVERDETERMINED SYSTEM
C      OF EQUATIONS. CHOOSE THE EQUATIONS THAT CORRESPOND TO THE SECOND
C      NON-ZERO COMPONENT OF EACH A(II,.) AND SOLVE BY GAUSS ELIMINATION
C      WITHOUT PIVOTING THE DERIVED BLOCK DIAGONAL POSITIVE DEFINITE
C      SYSTEM OF EQUATIONS. PROCEED AS IT IS EXPLAINED IN THE TEXTUAL
C      PART THIS ALGORITHM.
C.......................................................................
      IF (NACT.LE.I1-1) GO TO 170
C
C     THE GRADIENT OF THE OBJECTIVE FUNCTION IS DEFINED BY, SAY,
C      GRAD(I)=2*(SX(I+1)-F(I+1)), I=I1,I1+1,...,N-2. HOWEVER, INSTEAD
C      OF AN ARRAY GRAD(.) THE SUBROUTINE USES DIRECTLY THE VALUES OF
C      THE GRADIENT AT X(.) WHEN NECESSARY.
C
C     CALCULATE SIZE K OF A BLOCK OF THE LAGRANGE COEFFICIENT
C      MATRIX.
C
      I=I1
110   K=1
120   IF (I.LT.NACT .AND. IACT(I).EQ.IACT(I+1)-1) THEN
          I=I+1
          K=K+1
          GO TO 120
      END IF
C
130   IF (K.EQ.1) THEN
          II=IACT(I)
          PAR(II)=(SX(II+1)-F(II+1))/EC(II)
      ELSE
C
C         FACTORIZE THE K*K 3/DIAG SYSTEM OF LAGRANGE COEF.
C
          UE(I-K+1)=EC(IACT(I-K+1))
          DO 140 L=I-K+2,I
             DI=DC(IACT(L))
140          UE(L)=EC(IACT(L))-DI/UE(L-1)*DI
C
C         USE THIS FACTORIZATION TO SOLVE FOR PAR(.).
C
          GRAD=SX(IACT(I-K+1)+1)-F(IACT(I-K+1)+1)
          PAR(IACT(I-K+1))=GRAD
          DO 150 L=I-K+2,I
             GRAD=SX(IACT(L)+1)-F(IACT(L)+1)
150          PAR(IACT(L))=GRAD-DC(IACT(L))*PAR(IACT(L-1))/UE(L-1)
          PAR(IACT(I))=PAR(IACT(I))/UE(I)
          DO 160 L=I-1,I-K+1,-1
160          PAR(IACT(L))=PAR(IACT(L))/UE(L)-DC(IACT(L+1))*
     1                                     PAR(IACT(L+1))/UE(L)
      END IF
      I=I+1
      IF (I-NACT) 110,110,170
C
170   RETURN
      END
C.......................................................................
C     MESSAGES ASSOCIATED WITH THE OPERATIONS OF SUBROUTINE L2CXFT.
C.......................................................................
      SUBROUTINE MESSGE(TEMP,ITEMP1,ITEMP2,ITER,LABEL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ITER(*)
C
C.... I N P U T ....
C     TEMP,ITEMP1,ITEMP2   VARIABLES THAT ALLOW VALUES TO PASS FROM
C                 SUBROUTINE L2CXFT.
C     ITER(*)    INTEGER ARRAY, HOLDING NUMBER OF ITERATIONS. IT IS
C                 DEFINED IN SUBROUTINE L2CXFT.
C     LABEL      INTEGER VARIABLE, INDICATING THE NEAREST LABEL IN
C                 SUBROUTINE L2CXFT, WHERE THE MESSAGE IS CAUSED.
C
C.... O U T P U T ....
C     VARIOUS MESSAGES THAT ARE ACTIVATED ONLY IF THE ARGUMENT IPRINT
C      OF SUBROUTINE L2CXFT HAS THE VALUE 1.
C......................................................................
C
      IF (LABEL.EQ.10) THEN
          PRINT 1010
1010      FORMAT(//5X,'ERROR RETURN FROM L2CXFT BECAUSE THE ',
     1            /5X,'DATA LOWER BOUND EXCEEDS THE UPPER BOUND.')
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.20) THEN
          PRINT 1020
1020      FORMAT(//5X,'THE DATA ITSELF IS THE REQUIRED APPROXIMATION.')
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.160) THEN
          PRINT 1160,ITER(1),ITEMP1,TEMP,ITEMP2
1160      FORMAT (//5X,'ITERATIONS = ',I5,
     1            /5X,'NUMBER OF VIOLATED CONSTRAINTS    =',I5,
     2            /5X,'MAXIMUM CONSTRAINT  VIOLATION     =',E20.10,
     3            /5X,'                AT DATA POINT     =',I5)
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.170) THEN
          PRINT 1170
1170      FORMAT (//5X,'L2CXFT RETURNS BECAUSE THE ACCURACY IS'
     1      ' INSUFFICIENT',/5X,'TO MAINTAIN INCREASING THE '
     2      'FUNCTION VALUES.')
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.180) THEN
          PRINT 1180
1180      FORMAT (/5X,'ITERATIVE REFINEMENT IS INITIATED.')
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.210) THEN
          PRINT 1210,ITER(1)
1210      FORMAT (//5X,'L2CXFT QUITS ON ITERATION COUNT = ',I5)
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.220) THEN
          PRINT 1220,TEMP
1220      FORMAT (/5X,'MULTIPLIER OF INSERTED CONSTRAINT =',E20.10)
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.350) THEN
          PRINT 1350,ITEMP1,TEMP
1350      FORMAT (/5X,'ITERATION IN THE DUAL SPACE DUE TO ROUND OFF',
     1      /5X,'CONSTRAINT INDEX    = ',I5,
     2      /5X,'LAGRANGE MULTIPLIER = ',E20.10)
          GO TO 2000
      END IF
C
      IF (LABEL.EQ.1000) THEN
          PRINT 1500,ITEMP2,ITER(2),ITER(4),ITEMP1,ITER(6),ITER(5)
1500      FORMAT (//5X,'---------------------------------------',
     1             /5X,'RETURN FROM L2CXFT WITH MODE = ',I5,
     2            //5X,'INITIAL ACTIVE SET  = ',I5,
     3             /5X,'STARTING ACTIVE SET = ',I5,
     4             /5X,'FINAL ACTIVE SET    = ',I5,
     5            //5X,'ADDITIONS OF CONSTRAINTS = ',I5,
     6             /5X,'DELETIONS OF CONSTRAINTS = ',I5,
     7             /5X,'---------------------------------------')
          GO TO 2000
      END IF
C
2000  RETURN
      END
C......................................................................
C     THIS SUBROUTINE REVISES THE (J-1) AND (J+1) NORMAL EQUATIONS
C      OF SX(.) AFTER DELETING THE KNOT AKN(J).
C     THE MODIFIED ELEMENTS REMAIN STORED IN THE SAME POSITIONS.
C......................................................................
      SUBROUTINE MINUS1(E,B,D,AKN,I1,N,J)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION E(I1-1:N),B(I1-1:N),D(I1-1:N),AKN(I1-2:N)
C
C.... I N P U T ....
C     J       INTEGER VARIABLE, SUCH THAT AKN(J) IS THE DELETED KNOT.
C     E,B,D,AKN,I1,N  AS THEY ARE DEFINED IN SUBROUTINE L2CXFT.
C
C.... O U T P U T ....
C     E(J-1),E(J+1),B(J-1),B(J),D(J+1) THE ONLY ELEMENTS OF THE
C             NORMAL EQUATIONS THAT NEED BE MODIFIED BECAUSE OF J-TH
C             KNOT DELETION.
C
C.... M E T H O D ....
C     REF: FORMULAE (2.9) AND (2.10) IN "THE MINIMUM SUM OF
C          SQUARES CHANGE TO UNIVARIATE DATA THAT GIVES CONVEXITY",
C          BY I.C.DEMETRIOU AND M.J.D.POWELL, SUBMITTED TO IMA J. OF
C          NUMERICAL ANALYSIS, 1989.
C......................................................................
      C1=(AKN(J)-AKN(J-1))/(AKN(J+1)-AKN(J-1))
      C2=(AKN(J+1)-AKN(J))/(AKN(J+1)-AKN(J-1))
      E(J-1)=E(J-1)+C2*C2*E(J)+2*C2*D(J)
      E(J+1)=C1*C1*E(J)+E(J+1)+2*C1*D(J+1)
      D(J+1)=C1*D(J)+C2*D(J+1)+C1*C2*E(J)
      B(J-1)=B(J-1)+C2*B(J)
      B(J+1)=C1*B(J)+B(J+1)
C
      RETURN
      END
C.......................................................................
C     THIS SUBROUTINE REVISES THE (J-1), THE J AND THE (J+1) NORMAL
C      EQUATIONS OF SPLINE SX(.), AFTER X(IQ) IS INSERTED IN THE KNOT
C      SEQUENCE AND AKN(J) IS REDEFINED AS AKN(J)=X(IQ).
C.......................................................................
      SUBROUTINE PLUS(X,F,AKN,E,B,D,NK,I1,N,IPL,IP,IQ,IR,IRR,J)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(I1:N),F(I1:N),AKN(I1-2:N),E(I1-1:N),B(I1-1:N),
     1  D(I1-1:N),SV(24)
C
C.... I N P U T ....
C     E,B,D,AKN,X,F,I1,N,NK  AS THEY ARE DEFINED IN SUBROUTINE L2CXFT.
C     J      INTEGER VARIABLE, AS DEFINED IN SUBROUTINE ADDKNT, INDEX
C            OF NEW INSERTED KNOT.
C     IP     INTEGER VARIABLE, DATA INDEX OF THE KNOT ADJACENT TO
C             THE KNOT AKN(J) AT THE LEFT, IP=IAKN(J-1).
C     IPL    IF J=I1 THEN IPL=I1 ELSE IPL=IAKN(J-1).
C     IQ     INTEGER VARIABLE, DATA INDEX OF INSERTED KNOT, IQ=IAKN(J).
C     IR     INTEGER VARIABLE, DATA INDEX OF THE KNOT ADJACENT TO
C             THE KNOT AKN(J) AT THE RIGHT, IR=IAKN(J+1).
C     IRR    IF J=NK-1 THEN IRR=N ELSE IRR=IAKN(J+1).
C
C.... W O R K I N G  S P A C E ...
C     SV(1:24) REAL ARRAY WHOSE ELEMENTS PROVIDE TEMPORARY STORAGE
C            FOR THE CALCULATIONS OF (E(K),B(K), K=J-1,J,J+1) AND
C            (D(K), K=J,J+1).
C
C.... M E T H O D ....
C     GIVEN THE DATA POINTS X(IPL), X(IP), X(IQ), X(IR) AND X(IRR)
C      AS THE METHOD IN SUBROUTINE ADDKNT DESCRIBES, THEN CALCULATE
C      AB INITIO THE ELEMENTS OF E(K) AND B(K) FOR K=J-1,J,J+1
C      AND D(K) FOR K=J,J+1 (SEE C. DE BOOR, A PRACTICAL GUIDE
C      TO SPLINES, SPRINGER-VERLAG, 1978).
C.......................................................................
C     CALCULATE THE SCALAR PRODUCTS WHICH DEFINE E(.), B(.) AND D(.).
C      SEE COMMENTS IN THE BEGINNING OF SUBROUTINE L2CXFT.
C
      DO 50 I=1,24
50       SV(I)=0.
C
      DO 100 I=IP,IQ-1
         SV(1)=X(I)-AKN(J-1)
         SV(2)=AKN(J)-X(I)
         SV(3)=SV(3)+SV(2)*SV(2)
         SV(4)=SV(4)+SV(1)*SV(1)
         SV(5)=SV(5)+SV(2)*SV(1)
         SV(6)=SV(6)+SV(2)*F(I)
100      SV(7)=SV(7)+SV(1)*F(I)
C
      DO 200 I=IQ,IR-1
         SV(14)=AKN(J+1)-X(I)
         SV(2)=X(I)-AKN(J)
         SV(15)=SV(15)+SV(14)*SV(14)
         SV(16)=SV(16)+SV(2)*SV(2)
         SV(17)=SV(17)+SV(14)*SV(2)
         SV(18)=SV(18)+SV(14)*F(I)
200      SV(19)=SV(19)+SV(2)*F(I)
C
      IF (J.EQ.NK-1) THEN
          SV(11)=(AKN(J+1)-AKN(J))**2
          SV(9)=(AKN(J+1)-AKN(J))*F(N)
      ELSE
          DO 300 I=IR,IRR-1
             SV(24)=AKN(J+2)-X(I)
             SV(9)=SV(9)+SV(24)*F(I)
300          SV(11)=SV(11)+SV(24)*SV(24)
      END IF
C
      IF (I1.LT.J) THEN
          DO 400 I=IPL,IP-1
             SV(23)=X(I)-AKN(J-2)
             SV(8)=SV(8)+SV(23)*F(I)
400          SV(10)=SV(10)+SV(23)*SV(23)
      END IF
C
C     CALCULATE B(J-1),B(J),B(J+1).
C
      SV(2)=AKN(J-1)-AKN(J-2)
      SV(14)=AKN(J)-AKN(J-1)
      SV(21)=AKN(J+1)-AKN(J)
      SV(22)=AKN(J+2)-AKN(J+1)
C
      IF (J.EQ.I1) THEN
          B(J-1)=SV(6)/SV(14)
      ELSE
          B(J-1)=SV(8)/SV(2)+SV(6)/SV(14)
      END IF
      B(J)=SV(7)/SV(14)+SV(18)/SV(21)
C
      IF (J.EQ.NK-1) THEN
          B(J+1)=SV(19)/SV(21)+F(N)
      ELSE
          B(J+1)=SV(19)/SV(21)+SV(9)/SV(22)
      END IF
C
C     CALCULATE E(J-1),E(J),E(J+1),D(J-1),D(J).
C
      SV(2)=SV(2)*SV(2)
      SV(14)=SV(14)*SV(14)
      SV(21)=SV(21)*SV(21)
      SV(22)=SV(22)*SV(22)
C
      IF (J.EQ.I1) THEN
          E(J-1)=SV(3)/SV(14)
      ELSE
          E(J-1)=SV(10)/SV(2)+SV(3)/SV(14)
      END IF
      E(J)=SV(4)/SV(14)+SV(15)/SV(21)
C
      IF (J.EQ.NK-1) THEN
          E(J+1)=SV(16)/SV(21)+1
      ELSE
          E(J+1)=SV(16)/SV(21)+SV(11)/SV(22)
      END IF
      D(J)=SV(5)/SV(14)
      D(J+1)=SV(17)/SV(21)
C
      RETURN
      END
C.......................................................................
C     GENERATE SPLINE SX(.) AT X(.) FROM COEF(.) OVER AKN(.).
C.......................................................................
      SUBROUTINE SXATXI(SX,X,COEF,AKN,IAKN,I1,N,NK)
      IMPLICIT DOUBLE PRECISION  (A-H,O-Z)
      DIMENSION SX(I1:N),X(I1:N),COEF(I1-1:N),AKN(I1-2:N),IAKN(I1-1:N)
C
C.... I N P U T ....
C     X,COEF,AKN,IAKN,I1,N,NK  AS THEY ARE DEFINED IN SUBROUTINE L2CXFT.
C
C.... O U T P U T ....
C     SX(I1:N) REAL ARRAY CONTAINING THE SPLINE COMPONENTS AT X(I)
C               I=I1,I1+1,...,N.
C
C.... M E T H O D (GENERATE SPLINE COMPONENTS FROM COEFFICIENTS) ....
C     SEE ALGORITHM 5.8 OF L.L. SCHUMAKER, SPLINE FUNCTIONS: BASIC
C     THEORY, J. WILEY AND SONS, 1981, P.194 .
C.......................................................................
      I0=I1-1
      SX(I1)=COEF(I0)
      SX(N)=COEF(NK)
      I=I1+1
      L=I0
50    IF (I.LT.N .AND. IAKN(L).LT.I .AND. I.LE.IAKN(L+1)) THEN
          DENOM=AKN(L+1)-AKN(L)
          A1=(X(I)-AKN(L))/DENOM
          A2=1-A1
          SX(I)=COEF(L)*A2+COEF(L+1)*A1
          I=I+1
          GO TO 50
      END IF
      L=L+1
      IF (L.LT.NK) GO TO 50
C
      RETURN
      END
c*** dptests.f
      PROGRAM PROG1
C
C     SAMPLE PROGRAM THAT DEMONSTRATES SUBROUTINE L2CXFT WHEN THE
C      DATA IS READ FROM, AND THE RESULTS ARE DIRECTED TO, FILES.
C
C     CALLS SUBROUTINE L2CXFT.
C.......................................................................
      PARAMETER (IX1=1,
     1           NX=2000,
     2           IX0=IX1-1,
     3           IFILE=1,
     4           IPRINT=0)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(IX1:NX),F(IX1:NX),SX(IX1:NX),COEF(IX0:NX),
     2  AKN(IX0-1:NX),IAKN(IX0:NX),IACT(IX1:NX),ITER(10),
     3  PAR(IX1:NX),PARNEW(IX1:NX),RES(IX0:NX),DC(IX0:NX),EC(IX0:NX),
     4  E(IX0:NX),B(IX0:NX),D(IX0:NX),Z(IX0:NX),U(IX0:NX),SXOLD(IX1:NX)
C
C......................................................................
C.... I N P U T ....
C     X,F,SX,COEF,AKN,IAKN,IACT,ITER : ARRAYS THAT ARE EXPLAINED IN
C               SUBROUTINE L2CXFT.
C     I1,N      LOWER,UPPER INDICES OF DATA ARRAYS X(.) AND F(.). THEY
C               ARE EXPLAINED IN SUBROUTINE L2CXFT.
C     'XFDTA'   DATAFILE THAT HOLDS I1 AND N, AND THE DATA PAIRS
C               X(I),F(I) FOR I=I1,I1+1,...,N. 'DATA' IS PROVIDED
C               BY THE USER (OR BY PROG3 FOR TESTING PURPOSES).
C     IFILE     INTEGER PARAMETER WHOSE DEFAULT VALUE IS 0.
C               IF IFILE=1  THEN SOME PRINTOUT IS ACTIVATED AT THE
C                           END OF THE SMOOTHING PROCESS.
C     IPRINT    INTEGER PARAMETER WHOSE DEFAULT VALUE IS O.
C               IF IPRINT=1 THEN SUBROUTINE MESSGE IS CALLED BY SUB-
C                           ROUTINE L2CXFT TO PRINT INTERMEDIATE
C                           RESULTS (USEFUL TO SUBRTN L2CXFT INSTALTN).
C
C.... O U T P U T
C     RESULTS FROM THE CALL OF SUBROUTINE L2CXFT: THE LEAST SQUARES
C      BEST CONVEX FIT TO THE DATA X(.)-F(.).
C
C     FILE 'CXAPPX' KEEPS BEST CONVEX FIT.
C     FILE 'SPLNCF' KEEPS SPLINE KNOTS AND COEFFICIENTS OF BEST FIT.
C......................................................................
C
C      READ (UNFORMATTED) DATA FROM DATAFILE.
C
      OPEN (2,FILE='XFDTA')
            REWIND 2
            READ (2,*) I1,N
            I0=I1-1
            DO 100 I=I1,N
100            READ (2,*) X(I),F(I)
      CLOSE(2)
      PRINT 1110
1110  FORMAT(/5X,'X-F DATA SUPPLY FOR SMOOTHING FROM DATAFILE "XFDTA".')
C
C     CALCULATE BEST CONVEX FIT TO DATA.
C
      CALL L2CXFT(X,F,SX,I1,N,AKN,IAKN,COEF,NK,IACT,NACT,PAR,
     1  ITER,SXOLD,E,B,D,Z,U,DC,EC,PARNEW,RES,IPRINT,MODE)
C
C     SEND RESULTS TO FILES.
C
      OPEN (3,FILE='CXAPPX')
            WRITE (3,120) I1,N
120         FORMAT (2I5)
            DO 130 I=I1,N
130            WRITE (3,140) X(I),SX(I)
140            FORMAT(2E20.10)
      CLOSE(3)
      PRINT 1140
1140  FORMAT(/5X,'THE BEST CONVEX APPROXIMATION TO THE DATA ',
     1               'IS NOW KEPT IN FILE "CXAPPX". ')
C
      OPEN (4,FILE='SPLNCF')
            WRITE (4,150) I1,NK
150         FORMAT (2I5)
            DO 160 I=I1-1,NK
160            WRITE (4,170) IAKN(I),COEF(I)
170            FORMAT(I5,E20.10)
      CLOSE(4)
      PRINT 1170
1170  FORMAT(/5X,'SPLINE KNOTS AND COEFFICIENTS OF THE BEST CONVEX',
     1      /5X,'   APPROXIMATION TO THE DATA ARE NOW KEPT IN FILE ',
     2          '........... "SPLNCF". ')
C
C     PRINTOUT.
C      THE USER: IF N=LARGE, MAY DIRECT THIS PRINTOUT TO A FILE.
C
      IF (IFILE.EQ.1) THEN
          WRITE (*,1200) (J,AKN(J),IAKN(J),COEF(J),J=I0,NK)
1200      FORMAT(////
     1     12X,4HKNOT,6X,10HDATA INDEX,6X,13HB-SPLINE COEF/
     2     12X,5H(AKN),6X,6H(IAKN),10X,6H(COEF)//
     3     (I5,F14.8,I5,10X,E20.10))
C
          DO 1310 I=I1,N
             IF (I.EQ.I1) THEN
                 WRITE (*,1250) I,X(I),F(I),SX(I)
             ELSE IF (I.EQ.N) THEN
                      WRITE (*,1260) I,X(I),F(I),SX(I)
             ELSE
                 WRITE (*,1300) I,X(I),F(I),SX(I),PAR(I-1)
             ENDIF
1250         FORMAT(/////7X,11HDATA POINTS,6X,12HMEASUREMENTS,8X,
     1        13HAPPROXIMATION,7X,13HLAGRANGE MULT/
     2        7X,3H(X),14X,3H(F),17X,4H(SX),16X,5H(PAR)//
     3        (I5,F14.8,E20.10,E20.10))
1260          FORMAT(I5,F14.8,E20.10,E20.10)
1300          FORMAT(I5,F14.8,E20.10,E20.10,E20.10)
1310      CONTINUE
C
         ANORM2=0.0
         DO 1350 I=I1,N
1350        ANORM2=ANORM2+(SX(I)-F(I))*(SX(I)-F(I))
            WRITE(*,1400) ANORM2
1400        FORMAT(/2X,43HVALUE OF OBJECTIVE FUNCTION AT THE OPTIMUM=,
     1       E20.10)
      END IF
C
      STOP
      END
      PROGRAM PROG2
C
C     SAMPLE PROGRAM THAT GENERATES PSEUDORANDOMLY DATA FOR
C      TESTING THE SUBROUTINE L2CXFT.
C
C     CALLS FUNCTION RND.
C.......................................................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (I1=1,
     1           N=75,
     2           XA=0.0,
     3           XB=1.0,
     4           SIZE=50.0/N**2)
      DIMENSION X(I1:N),F(I1:N)
C
C.... I N P U T  AND  P A R A M E T E R S  (SET BY THE USER) ....
C     I1        LOWER DATA INDEX.
C     N         UPPER DATA INDEX.
C     XA        LEFT LIMIT OF X(.), X(I1)=XA.
C     XB        RIGHT LIMIT OF X(.), X(N)=XB.
C     SIZE      RELATIVE MAGNITUDE OF NOISE.
C                RECOMMENTED VALUES FOR SIZE THAT COVER A WIDE RANGE
C                OF DATA BEHAVIOUR ARE (WHEN I1=1)
C                = 0
C                = 4/N**2
C                = 10/N**2
C                = 50/N**2
C                = 250/N**2
C                = 1000/N**2
C
C.... O U T P U T  (DIRECTED TO THE FILE "XFDTA") ....
C     X(I1:N)   DATA POINTS (ABSCISSAE).
C     F(I1:N)   FUNCTION MEASUREMENTS (INCLUDING NOISE).
C
C.... M E T H O D  (USER INTERFACE) ....
C     THIS PROGRAM GENERATES RANDOMLY DATA AS FOLLOWS:
C      A CONTINOUS (+CONVEX) FUNCTION F(X) DEFINED ON [XA,XB]
C      IS SUPPLIED BY THE USER. THEN F(X) IS EVALUATED ON THE
C      GRID XA=X(I1)<X(I1+1)<...<X(N)=XB AND THEN WE ADD
C      TO THE FUNCTION VALUES RANDOM NUMBERS FROM THE UNIFORM
C      DISTRIBUTION OVER THE INTERVAL [-SIZE,SIZE].
C      THE TYPE OF THE FUNCTION, THE LIMITS OF THE DATA AND THE SIZE
C      OF THE MAGNITUDE OF THE RANDOM NUMBERS ARE UPON THE USER'S
C      DECISION . THE GRID NEED NOT BE EQUALLY SPACED.
C
C     THE METHOD THAT PRODUCES THE PSEUDO-RANDOM NUMBERS MAY BE
C      FOUND IN FUNCTION RND(.).
C
C     THE DATA ARE GENERATED AND THEN DIRECTED TO FILE "XFDTA". THIS
C      FILE IS INPUT TO PROGRAM PROG1.
C......................................................................
C
C     SET DATA POINTS.
C
      STEP=(XB-XA)/(N-I1)
      X(I1)=XA
      A=X(I1)
      DO 20 I=I1+1,N-1
         X(I)=A+STEP
20       A=X(I)
      X(N)=XB
C
C     GENERATE DATA BY ADDING NOISE TO FUNCTION MEASUREMENTS OF THE
C      RUNGE FUNCTION FROM 0.0 TO 1.0.
C     ISEED IS SET BY THE USER. IT CAN BE ANY INTEGER SUCH THAT
C       1 .LE. ISEED .LE. 65535 .
C
      ISEED=1
      DO 50 I=I1,N
50       F(I)=1/(1+10*X(I)*X(I))+SIZE*RND(ISEED)
C
C     SEND (X,F) DATA PAIRS TO A DATAFILE.
C
      OPEN (2,FILE='XFDTA')
            WRITE (2,120) I1,N
120         FORMAT (2I5)
            DO 130 I=I1,N
130            WRITE (2,140) X(I),F(I)
140            FORMAT (2E20.10)
      CLOSE(2)
      PRINT 1140
1140  FORMAT(/5X,'DATA HAS BEEN GENERATED AND KEPT IN FILE "XFDTA".')
C
      STOP
      END
      PROGRAM PROGMD
C......................................................................
C     THIS PROGRAM IS A DRIVER OF SUBROUTINE L2CXFT.
C      IT USES REPEATEDLY THE SIMPLE DRIVER PROGRAM PROG2 BY ADDING
C      RUNDOM NUMBERS (NOISE) TO GIVEN FUNCTION MEASUREMENTS SO THAT
C      SUBSTANTIALLY DIFFERENT DATA SETS ARE PRODUCED AND TESTED.
C      THE NUMBER OF DATA AND THE RELATIVE MAGNITUDE OF THE NOISE ARE
C      VARIABLES IN NESTED LOOPS.
C
C     CALLS SUBROUTINE L2CXFT.
C.......................................................................
      PARAMETER (IX1=1,
     1           NX=2000,
     2           IX0=IX1-1,
     3           IPRINT=0)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(IX1:NX),F(IX1:NX),SX(IX1:NX),COEF(IX0:NX),
     1  AKN(IX0-1:NX),IAKN(IX0:NX),IACT(IX1:NX),ITER(10),
     2  PAR(IX1:NX),PAROLD(IX1:NX),RES(IX0:NX),DC(IX0:NX),EC(IX0:NX),
     3  E(IX0:NX),B(IX0:NX),D(IX0:NX),Z(IX0:NX),U(IX0:NX),SXOLD(IX1:NX)
C
C.... U S E R  I N T E R F A C E  (INPUT AND PARAMETERS) ....
C
C     THE USER HAS ONLY TO SPECIFY THE PARAMETERS IX1 AND NX, A FUNCTION
C      TO PROVIDE MEASUREMENTS AND LET THIS PROGRAM RUN.
C
C     IX1        LOWER DATA INDEX (SET BY THE USER).
C     NX         UPPER DATA INDEX (SET BY THE USER).
C
C     ALL ARRAYS ARE EXPLAINED IN SUBROUTINE L2CXFT.
C     IPRINT    INTEGER PARAMETER WHOSE DEFAULT VALUE IS O.
C               IF IPRINT=1 THEN SUBROUTINE MESSGE IS CALLED BY SUB-
C                           ROUTINE L2CXFT TO PRINT INTERMEDIATE
C                           RESULTS (USEFUL TO SUBRTN L2CXFT INSTALTN).
C
C.... O U T P U T ....
C     RESULTS FROM THE CALL OF SUBROUTINE L2CXFT CONCERNING ACTIVE
C      SET CHANGES. SEE SUBROUTINE MESSGE.
C
C.... M E T H O D  FOR PRODUCING THE RANDOM NUMBERS ....
C     SEE PROG2.
C......................................................................
      I1=IX1
      I0=I1-1
C
C     SIZE      RELATIVE MAGNITUDE OF NOISE.
C                RECOMMENTED VALUES FOR SIZE THAT COVER A WIDE RANGE
C                OF DATA BEHAVIOUR ARE (WHEN I1=1)
C                = 0
C                = 4/N**2
C                = 10/N**2
C                = 50/N**2
C                = 250/N**2
C                = 1000/N**2
C
      DO 160 JJ=0,6
         IF (JJ.EQ.0) ISIZE=0
         IF (JJ.EQ.1) ISIZE=4
         IF (JJ.EQ.2) ISIZE=10
         IF (JJ.EQ.3) ISIZE=50
         IF (JJ.EQ.4) ISIZE=100
         IF (JJ.EQ.5) ISIZE=250
         IF (JJ.EQ.6) ISIZE=1000
C
      DO 150 NN=1,5
         IF (NN.EQ.1)  N=20
         IF (NN.EQ.2)  N=100
         IF (NN.EQ.3)  N=500
         IF (NN.EQ.4)  N=1000
         IF (NN.EQ.5)  N=2000
C
C     SET DATA POINTS.
C
      XA=0.
      XB=1.
      STEP=(XB-XA)/(N-I1)
      X(I1)=XA
      A=X(I1)
      DO 20 I=I1+1,N-1
         X(I)=A+STEP
20       A=X(I)
      X(N)=XB
C
      PRINT *, '--------------------------------------------'
      PRINT *, '   DBL PRECISION'
      PRINT *, '   RUNGE FUNCTION + NOISE '
      PRINT *, '--------------------------------------------'
      PRINT *, ' NUMBER OF DATA = ',N,' ERROR MULTIPLIER= ',ISIZE
      SIZE=FLOAT(ISIZE)/(N*N)
      PRINT *, ' ERROR SIZE     = ',SIZE
C
C     GENERATE MEASUREMENTS WITH NOISE.
C
      ISEED=1
      DO 30 I=I1,N
C        ADD ERROR TO MEASUREMENTS OF THE RUNGE FUNCTION.
         F(I)=1/(1+10.*X(I)*X(I))+SIZE*RND(ISEED)
C
C        ADD ERROR TO MEASUREMENTS OF THE BERNSTEIN FUNCTION.
*         F(I)=X(I)**4-3./4.*X(I)**3+3.*X(I)*X(I)-5./2.*X(I)+
*     1                                      SIZE*RND(ISEED)
30    CONTINUE
C
C.... CALCULATE THE CONVEX FIT ........................................
C
      CALL L2CXFT(X,F,SX,I1,N,AKN,IAKN,COEF,NK,IACT,NACT,PAR,
     1  ITER,SXOLD,E,B,D,Z,U,DC,EC,PAROLD,RES,IPRINT,MODE)
C
C     PRINTOUT
C
      ANORM2=0.0
      DO 50 I=I1,N
50       ANORM2=ANORM2+(SX(I)-F(I))*(SX(I)-F(I))
         PRINT 1050, ANORM2
1050     FORMAT(/2X,43HVALUE OF OBJECTIVE FUNCTION AT THE OPTIMUM=,
     1     E20.10)
C     END OF N LOOP.
150   CONTINUE
C     END OF ISIZE LOOP.
160   CONTINUE
C
      STOP
      END
C.....................................................................
C     THIS FUNCTION GENERATES RANDOM NUMBERS SUCH THAT
C               -1 < RND < 1
C      WHICH ARE ADDED TO GIVEN FUNCTION MEASUREMENTS, IN ORDER
C       TO CONSTRUCT DATA SETS FOR TESTING THE SUBROUTINE L2CXFT.
C.....................................................................
      FUNCTION RND(ISEED)
      DOUBLE PRECISION RND
C
C     FUNCTION RND CAN GENERATE 65536 "RANDOM" NUMBERS BEFORE
C      REPEATING ITSELF. IT WILL WORK ON MOST COMPUTERS FOR
C      WHICH 2**31-1.LT.MAXINT.
C
C.... I N P U T ....
C     ISEED     INTEGER VARIABLE BETWEEN 1 AND 65535 .
C
C.... O U T P U T ....
C     RND       THE RETURN OF THE FUNCTION EVALUATION,
C               BETWEEN -1 AND 1.
C
C.... M E T H O D (LINEAR CONGRUENTIAL) ....
C      EACH NUMBER IN THE RANDOM SEQUENCE, R(K) SAY,
C      IS CALCULATED FROM ITS PREDECESSOR, R(K-1), USING THE FORMULA:
C               R(K)=(MULTIPLIER * R(K-1) + INCREMENT) MOD 2**16
C      WHERE MULTIPLIER=25173 AND INCREMENT=13849.
C
C..... REF: P.GROGONO, PROGRAMMING IN PASCAL, ADDISON-WESLEY PUB.CO.,
C          1980, P.118.
C.....................................................................
      RND=FLOAT(ISEED)/65535.
      RND=2*RND-1
      ISEED=25173*ISEED+13849
      ISEED=MOD(ISEED,65536)
C
      RETURN
      END
c*** l2cxft.txt
         ============================================================
                 DESCRIPTION AND USE OF THE FORTRAN CODES FOR

         "L2CXFT: A FORTRAN SUBROUTINE FOR LEAST SQUARES DATA FITTING
                 WITH NON-NEGATIVE SECOND DIVIDED DIFFERENCES"

                                I.C. DEMETRIOU

              An algorithm of ACM TOMS, Vol.21, No.1, March 1995
         ============================================================

                             University of Athens
                            Department of Economics
                      Unit of Mathematics and Informatics
                  8 Pesmazoglou Street, 105 59 Athens, Greece
                       E-mail: demetri@aueb.ariadne-t.gr


This file provides information about:

      DOCUMENTATION, LISTINGS, DRIVER PROGRAMS, INSTALLATION, TESTING AND
                 MAINTENANCE OF THE FORTRAN SUBROUTINE L2CXFT

SUMMARY
This file is complementary to the ACM  TOMS  submission  entitled  "L2CXFT:  A
Fortran subroutine for least squares data  fitting  with  non-negative  second
divided differences" by I.C. Demetriou [1]. It consists of three sections  and
four appendices organized as follows. Section 1 states the problem, Section  2
presents the purpose of each  subroutine  and  Section  3  outlines  the  user
interface. Appendix A keeps the Fortran listings of the  subroutines  together
with symbol explanations. Appendix B presents driver programs  and  output  of
test examples in order to help the usage of L2CXFT. Appendix C gives technical
details about installation, compilation, linking and running  of  the  Fortran
codes. Appendix D considers some coding details of each subroutine.
The entire set code occupies 1469 Fortran 77 lines including comments.  Single
and double precision versions have  been  developed  for  practical  use.  All
subroutines begin with comments that explain the input and  output  arguments,
the working space and the method followed.
Details about the method of calculation, the user interface and the subroutine
purpose  are  presented  in  [1].  The  underlying  method  is  mathematically
established by [2].
APPLICATIONS
Since the underlying method may be attended as a convex regression  technique,
it may be particularly useful to  statistical  analysis  and  economics  where
marginal utility, productivity, supply and demand functions are assumed to  be
convex or concave. The concave smoothing can be treated  as  convex  after  an
overall change of sign.
REFERENCES
[1] Demetriou, I.C., L2CXFT: A  Fortran  subroutine  for  least  squares  data
    fitting with nonnegative second divided differences. ACM  Transactions  on
    Mathematical Software, Vol.21, No.1, March 1995.
[2] Demetriou, I.C. and M.J.D.Powell, The minimum sum  of  squares  change  to
    univariate data that  gives  convexity.  IMA  J.  of  Numerical  Analysis,
    11(1991), pp. 433-448.
[3] Powell, M.J.D., "Approximation theory and methods".  Cambridge  University
    Press, Cambridge, 1981.

                              THIS FILE CONTA