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 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)=, J=I1-1,I1,...,NK C E(J)=, J=I1-1,I1,...,NK C B(J)=, 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) 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)=, J=I1,I1+1,...,NK C E(J)=, J=I1-1,I1,...,NK C B(J)=, 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) 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)=, I=I1-1,I1,...,J C E(I)=, I=I1-1,I1,...,J C B(I)=, I=I1-1,I1,...,J. C C THEN PROCEED AS IN REFERENCE QUOTED. C...................................................................... C NO TEST ON THE CONDITIONS N=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 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)=, J=I1-1,I1,...,NK C E(J)=, J=I1-1,I1,...,NK C B(J)=, 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) 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)=, J=I1,I1+1,...,NK C E(J)=, J=I1-1,I1,...,NK C B(J)=, 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) 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)=, I=I1-1,I1,...,J C E(I)=, I=I1-1,I1,...,J C B(I)=, I=I1-1,I1,...,J. C C THEN PROCEED AS IN REFERENCE QUOTED. C...................................................................... C NO TEST ON THE CONDITIONS N THE INNER PRODUCT OF TWO VECTORS AND, 154 C THE COEFFICIENTS OF SX(.) ARE DETERMINED BY SOLVING THE 155 C THE POSITIVE DEFINITE SYSTEM OF THE NORMAL EQUATIONS 156 C EFFICIENTLY AND STABLY BY CHOLESKY FACTORIZATION 157 C D(J)*COEF(J-1)+E(J)*COEF(J)+D(J+1)*COEF(J+1)=B(J) 158 C FOR J=I1-1,I1,...,NK 159 C WHERE BY <.,.> IS DENOTED THE SCALAR PRODUCT OF TWO VECTORS 160 C AND 161 C D(J)=, J=I1,I1+1,...,NK 162 C E(J)=, J=I1-1,I1,...,NK 163 C B(J)=, J=I1-1,I1,...,NK. 164 C (SUBROUTINE CFLGRF). 165 C 166 C B) CALCULATION OF THE LAGRANGE MULTIPLIERS 167 C THE LAGRANGE EQUATIONS AT AN ESTIMATE SX(.) ARE 168 C 2*(SX(II+1)-F(II+1))=SIGMA(PAR(II)*A(II,.)) 169 C FOR II=IACT(I), I=I1,I1+1,...,NACT 170 C WHERE A(II,J) IS THE J-TH COMPONENT OF THE NORMAL 171 C OF THE II-TH CONSTRAINT AND 172 C A(II,J)<>0 ONLY FOR (II.LE.J.LE.II+2). 173 C (ARRAY A(.,.) IS NOT ACTUALLY USED IN THIS SUBROUTINE; 174 C INSTEAD ARRAYS EC(.) AND DC(.) ARE EMPLOYED.). BECAUSE 175 C THIS SYSTEM IS OVERDETERMINED PAR(IACT(.)) IS DERIVED 176 C BY CHOOSING FROM THE LAGRANGE EQUATIONS THOSE EQUATIONS THAT 177 C CORRESPOND TO THE SECOND NON-ZERO COMPONENT OF A(IACT(.),.); 178 C SO A BLOCK DIAGONAL POSITIVE DEFINITE SYSTEM IS OBTAINED THA 179 C IS SOLVED EFFICIENTLY AND STABLY BY CHOLESKY FACTORIZATION. 180 C (SUBROUTINE CFLGRF). 181 C 182 C C) PARTICULAR METHODS ARE DEVELOPED FOR INSERTING INTO AND 183 C DELETING FROM THE SPLINE APPROXIMATION ONE KNOT ECONOMICALLY 184 C (SUBROUTINES ADDKNT AND DELKNT). 185 C...................................................................... 186 C SET CERTAIN INITIAL VALUES. 187 C 188 MODE=1 189 I0=I1-1 190 ITREF=0 191 ITEST=0 192 ZERO=0. 193 DO 10 I=1,10 194 10 ITER(I)=0 195 C 196 IF (N.LT.I1) THEN 197 MODE=4 198 CALL MESSGE(ZERO,0,0,ITER,10) 199 GO TO 1000 200 END IF 201 C 202 C IF N=I1 OR N=I1+1 THEN THE DATA ITSELF IS THE REQUIRED FIT. 203 C 204 IF (N.EQ.I1) THEN 205 SX(I1)=F(I1) 206 AKN(I1-1)=X(I1) 207 IAKN(I1-1)=I1 208 COEF(I1-1)=F(I1) 209 CALL MESSGE(ZERO,0,0,ITER,20) 210 GO TO 1000 211 END IF 212 C 213 IF (N.EQ.I1+1) THEN 214 SX(I1)=F(I1) 215 SX(N)=F(N) 216 AKN(I1-1)=X(I1) 217 AKN(I1)=X(N) 218 IAKN(I1-1)=I1 219 IAKN(I1)=N 220 COEF(I1-1)=F(I1) 221 COEF(I1)=F(N) 222 CALL MESSGE(ZERO,0,0,ITER,20) 223 GO TO 1000 224 END IF 225 C 226 C.... I N I T I A L I Z E ............................................. 227 C 228 C THERE ARE N-2 CONSTRAINTS DEFINED BY THE DIVIDED DIFFERENCES. 229 C WE ASSOCIATE THE I-TH CONSTRAINT WITH THE I-TH DIVIDED DIFFER- 230 C ENCE. SET COEFFICIENTS OF FULL SYSTEM OF EQUATIONS OF LAGRANGE 231 C MULTIPLIERS. SINCE GRAD(F)=2*(SX-F), THE COEFFICIENTS ARE 232 C SCALED BY 2. 233 C 234 DO 50 I=I1,N-3 235 CX2=X(I+2)-X(I+1) 236 CX1=X(I+1)-X(I) 237 CX3=X(I+2)-X(I) 238 CX2=1/CX2 239 CX1=1/CX1 240 EC(I)=-(CX1+CX2)/2 241 50 DC(I+1)=CX2/2 242 CX2=X(N)-X(N-1) 243 CX1=X(N-1)-X(N-2) 244 CX3=X(N)-X(N-2) 245 CX2=1/CX2 246 CX1=1/CX1 247 EC(N-2)=-(CX1+CX2)/2 248 C 249 DO 60 I=I1,N-2 250 60 PAR(I)=0.0 251 C 252 C.... CALCULATE STARTING POINT BY L2ILS ............................... 253 C 254 CALL L2ILS(COEF,AKN,IAKN,D,E,B,Z,U,NK,F,X,I1,N) 255 C 256 C INITIALIZE INDICES OF ACTIVE CONSTRAINTS IN IACT(.). 257 C 258 J=I0 259 I=I1+1 260 IAKN(I0)=I1 261 KL=I0 262 80 IF (IAKN(J).LT.I .AND. I.LT.IAKN(J+1)) THEN 263 KL=KL+1 264 IACT(KL)=I-1 265 I=I+1 266 GO TO 80 267 END IF 268 J=J+1 269 I=I+1 270 IF (J.LT.NK) GO TO 80 271 NACT=KL 272 ITER(2)=NACT-I0 273 C 274 C.... ITERATE BY ONE, UNTIL FIRST FEASIBLE LAGRANGIAN VECTOR IS FOUND. 275 C 276 C CALCULATE COEF OVER KNOTS, SPLINE OVER ABSCISSAE AND 277 C LAGRANGE MULTIPLIERS. 278 C 279 90 CALL CFLGRF(SX,PAR,COEF,IACT,E,B,D,X,AKN,IAKN,F,DC,EC,U, 280 1 RES,NK,NACT,I1,N,ITREF,ITER(3)) 281 ITER(3)=ITER(3)+1 282 C 283 C FIND MINIMUM LAGRANGE MULTIPLIER. 284 C 285 IF (NACT .LE. I0) GO TO 150 286 AMINLA=PAR(IACT(I1)) 287 IPAR=IACT(I1) 288 IF (I1 .EQ. NACT) GO TO 110 289 DO 100 I=I1+1,NACT 290 IF (PAR(IACT(I)) .LT. AMINLA) THEN 291 AMINLA=PAR(IACT(I)) 292 IPAR=IACT(I) 293 END IF 294 100 CONTINUE 295 110 ITER(4)=NACT-I0 296 IF (AMINLA .GE. 0.0) GO TO 150 297 C 298 C INSERT A KNOT WHERE THE LEAST NEGATIVE LAGRANGE MULTIPLIER. 299 C 300 INKNOT=IPAR+1 301 PAR(IPAR)=0.0 302 CALL ADDKNT(INKNOT,E,B,D,X,F,AKN,IAKN,IACT,NK,NACT,I1,N) 303 GO TO 90 304 C 305 C.... M A I N I T E R A T I O N ...................................... 306 C 307 150 CONTINUE 308 C 309 C EXAMINE CONVEXITES AT ALL KNOTS OF SPLINE APPROXIMATION. 310 C SPECIAL CASE: IF NK-I0=1 THEN SX(.) IS A STRAIGHT LINE. 311 C 312 IF (NK-I0 .EQ. 1) GO TO 1000 313 C 314 C FIND THE KNOT OF THE MOST VIOLATED CONSTRAINT. 315 C 316 CVMAX=0. 317 KL=I0 318 K=0 319 DO 160 J=I1,NK-1 320 JJ=IAKN(J) 321 CKNOT=(SX(JJ+1)-SX(JJ))/(X(JJ+1)-X(JJ))- 322 1 (SX(JJ)-SX(JJ-1))/(X(JJ)-X(JJ-1)) 323 SUM2=1/(X(JJ+1)-X(JJ))**2+ 324 1 ((X(JJ+1)-X(JJ-1))/((X(JJ+1)-X(JJ))*(X(JJ)-X(JJ-1))))**2+ 325 2 1/(X(JJ)-X(JJ-1))**2 326 SUM2=SQRT(SUM2) 327 CKNOT=CKNOT/SUM2 328 IF (CKNOT .LT. 0.) K=K+1 329 IF (CKNOT .LT. CVMAX) THEN 330 CVMAX=CKNOT 331 KNTOUT=J 332 IAKNJ=IAKN(J) 333 KL=KL+1 334 END IF 335 160 CONTINUE 336 C 337 C IF KL=I0 THEN SX(.) IS CONVEX. SINCE AT THIS LEVEL PAR(.) ARE 338 C ALL NON-NEGATIVE THE KUHN-TUCKER CONDITIONS ARE SATISFIED. 339 C 340 C.... T E R M I N A T I O N ........................................... 341 C 342 IF (KL .EQ. I0) GO TO 1000 343 C 344 C MSG: PRINT CVMAX,IAKNJ,K. 345 C 346 IF (IPRINT.NE.0) CALL MESSGE(CVMAX,K,IAKNJ,ITER,160) 347 C 348 C POWELL'S SAFEGUARD TERMINATION TEST: 349 C RETURN, IF DUE TO ROUNDING ERRORS, THE CHANGE IN SX(.) MAY 350 C NOT INCREASE THE OBJECTIVE FUNCTION. 351 C 352 IF (ITER(1).EQ.0) GO TO 190 353 R=0.0 354 T=0.0 355 DO 170 I=I1,N 356 R=R+(SX(I)-SXOLD(I))*(SX(I)+SXOLD(I)-2*F(I)) 357 170 T=T+ABS(SX(I)-SXOLD(I))*(ABS(SX(I)+SXOLD(I))+2*ABS(F(I))) 358 IF (T+R.LE.T .OR. T+1.5*R.LE.T+R) THEN 359 IF (ITEST.EQ.1) THEN 360 MODE=2 361 IF (IPRINT.NE.0) CALL MESSGE(ZERO,0,0,ITER,170) 362 C 363 C MSG: L2CXFT RETURNS BECAUSE THE ACCURACY IS 364 C INSUFFICIENT TO MAINTAIN INCREASING 365 C THE FUNCTION VALUES. 366 C 367 GO TO 1000 368 END IF 369 C 370 IF (IPRINT.EQ.0) GO TO 180 371 IF (ITREF.EQ.0) CALL MESSGE (ZERO,0,0,ITER,180) 372 C 373 C MSG: ITERATIVE REFINEMENT IS INITIATED. 374 C 375 180 ITREF=1 376 ITEST=1 377 ELSE 378 ITEST=0 379 END IF 380 190 CONTINUE 381 C 382 C KEEP OLD LAGRANGE MULTIPLIERS AND OLD SPLINE APPROXIMATION. 383 C 384 DO 200 I=I1,N-2 385 200 PAROLD(I)=PAR(I) 386 DO 210 I=I1,N 387 210 SXOLD(I)=SX(I) 388 C 389 C TERMINATION ON ITERATION COUNT. 390 C 391 ITER(1)=ITER(1)+1 392 IF (ITER(1).GT.3*(N-I0)) THEN 393 MODE=3 394 CALL MESSGE(ZERO,0,0,ITER,210) 395 C 396 C MSG: L2CXFT QUITS ON ITERATION COUNT. 397 C 398 GO TO 1000 399 END IF 400 C 401 C CALCULATE A SEARCH DIRECTION (DELETE THE KNOT KNTOUT). 402 C 403 CALL DELKNT(E,B,D,AKN,IAKN,IACT,NK,NACT,I1,N,KNTOUT) 404 ITER(6)=ITER(6)+1 405 C 406 C SOLVE FOR SX(.) OVER AKN(.). 407 C 408 220 CALL CFLGRF(SX,PAR,COEF,IACT,E,B,D,X,AKN,IAKN,F,DC,EC,U, 409 1 RES,NK,NACT,I1,N,ITREF,ITER(3)) 410 ITER(3)=ITER(3)+1 411 C 412 IPAR=IAKNJ-1 413 IF (IPRINT.NE.0) CALL MESSGE(PAR(IPAR),0,0,ITER,220) 414 C 415 C MSG: PRINT MULTIPLIER OF INSERTED CONSTRAINT. 416 C 417 C APPLY POWELL'S TECHNIQUE IN ORDER THAT THE INDEX IPAR 418 C REMAIN IN ACTIVE SET. SET PAR(IPAR)=MAX(PAR(IPAR),0.0) 419 C BECAUSE IN THEORY HOLDS THAT PAR(IPAR)>0. 420 C 421 IF (PAR(IPAR).LT.0.0) PAR(IPAR)=0.0 422 C 423 C EXAMINE SIGN OF LAGRANGE MULTIPLIERS 424 C 425 300 IF (NACT .EQ. I0) GO TO 150 426 KL=I0 427 THETA=2.0 428 DO 320 I=I1,NACT 429 K=IACT(I) 430 IF (PAR(K)) 310,320,320 431 310 RATIO=PAROLD(K)/(PAROLD(K)-PAR(K)) 432 IF (THETA .LE. RATIO) GO TO 320 433 THETA=RATIO 434 IPAR=K 435 KL=I1 436 320 CONTINUE 437 C 438 C.... ITERATE IN THE PRIMAL SPACE ..................................... 439 C IF KL=I0 THEN PAR(.) ARE NON-NEGATIVE. BRANCH WHERE TEST 440 C PRIMAL FEASIBILITY. 441 C 442 IF (KL .EQ. I0) GO TO 150 443 C 444 C.... RECOVER NON-NEGATIVITY OF LAGRANGE MULTIPLIERS .................. 445 C 446 DO 340 I=I1,NACT 447 K=IACT(I) 448 340 PAROLD(K)=(1-THETA)*PAROLD(K)+THETA*PAR(K) 449 C 450 C INSERT KNOT IPAR+1 INTO THE SPLINE AND SOLVE FOR SX(.). 451 C 452 I1NEW=I1 453 350 INKNOT=IPAR+1 454 PAR(IPAR)=0. 455 CALL ADDKNT(INKNOT,E,B,D,X,F,AKN,IAKN,IACT,NK,NACT,I1,N) 456 ITER(5)=ITER(5)+1 457 C 458 C TEST THE FEASIBILITY OF PAROLD(.). ALL PAROLD(.) MUST BE 459 C NON-NEGATIVE BECAUSE OF THE WAY THETA WAS DETERMINED. OTHERWISE, 460 C THIS IS DUE TO ROUND-OFF ERROR. HENCE, ANY NEGATIVITY 461 C IS TREATED AS BEING ZERO. SEE LAWSON & HANSON (SOLVING 462 C LEAST SQUARES PROBLEMS, PRENTICE-HALL, 1974, p.307). 463 C 464 IF (NACT .EQ. I0) GO TO 370 465 IF (I1NEW .GT. NACT) GO TO 370 466 DO 360 I=I1NEW,NACT 467 IPAR=IACT(I) 468 IF (PAROLD(IPAR) .LT. 0.0) THEN 469 IF (IPRINT.NE.0) CALL MESSGE(PAROLD(IPAR),IPAR,0,ITER,350) 470 C 471 C MSG: ITERATION IN THE DUAL SPACE DUE TO ROUND OFF. 472 C PRINT ACTIVE CONSTRAINT AND LAGRANGE MULTIPLIER. 473 C 474 I1NEW=I 475 GO TO 350 476 END IF 477 360 CONTINUE 478 C 479 370 CONTINUE 480 C 481 CALL CFLGRF(SX,PAR,COEF,IACT,E,B,D,X,AKN,IAKN,F,DC,EC,U, 482 1 RES,NK,NACT,I1,N,ITREF,ITER(3)) 483 ITER(3)=ITER(3)+1 484 C 485 C.... ITERATE IN THE DUAL SPACE ....................................... 486 C 487 GO TO 300 488 1000 CONTINUE 489 C 490 CALL MESSGE(ZERO,NACT-I0,MODE,ITER,1000) 491 C 492 C MSG: PRINT MODE AND 493 C INITIAL ACTIVE SET 494 C STARTING ACTIVE SET 495 C FINAL ACTIVE SET 496 C ADDITIONS OF CONSTRAINTS 497 C DELETIONS OF CONSTRAINTS 498 C 499 C.... END OF SUBROUTINE L2CXFT ........................................ 500 C 501 RETURN 502 END L2CXFT Local Symbols NAME CLASS TYPE EXPLANATION MODE. . . . . . . . . . . argument IPRINT. . . . . . . . . . argument RES . . . . . . . . . . . argument PAROLD. . . . . . . . . . argument EC. . . . . . . . . . . . argument DC. . . . . . . . . . . . argument U . . . . . . . . . . . . argument Z . . . . . . . . . . . . argument D . . . . . . . . . . . . argument B . . . . . . . . . . . . argument E . . . . . . . . . . . . argument SXOLD . . . . . . . . . . argument ITER. . . . . . . . . . . argument PAR . . . . . . . . . . . argument NACT. . . . . . . . . . . argument IACT. . . . . . . . . . . argument NK. . . . . . . . . . . . argument COEF. . . . . . . . . . . argument IAKN. . . . . . . . . . . argument AKN . . . . . . . . . . . argument N . . . . . . . . . . . . argument I1. . . . . . . . . . . . argument SX. . . . . . . . . . . . argument F . . . . . . . . . . . . argument X . . . . . . . . . . . . argument AMINLA. . . . . . . . . . local REAL min Lagrange multiplier CKNOT . . . . . . . . . . local REAL temporary CVMAX . . . . . . . . . . local REAL max constraint violation CX1 . . . . . . . . . . . local REAL temporary CX2 . . . . . . . . . . . local REAL temporary CX3 . . . . . . . . . . . local REAL temporary I . . . . . . . . . . . . local INTEGER temporary I0. . . . . . . . . . . . local INTEGER I0=I1-1 ITREF . . . . . . . . . . local INTEGER see comments on working space ITEST . . . . . . . . . . local INTEGER see comments on working space INKNOT. . . . . . . . . . local INTEGER data index of inserted knot IPAR. . . . . . . . . . . local INTEGER temporary IAKNJ . . . . . . . . . . local INTEGER data index of KNOUT J . . . . . . . . . . . . local INTEGER temporary JJ. . . . . . . . . . . . local INTEGER temporary K . . . . . . . . . . . . local INTEGER temporary KL. . . . . . . . . . . . local INTEGER temporary KNTOUT. . . . . . . . . . local INTEGER index of deleted knot R . . . . . . . . . . . . local REAL residual of obj function RATIO . . . . . . . . . . local REAL temporary SUM2. . . . . . . . . . . local REAL temporary T . . . . . . . . . . . . local REAL abs residual of obj function THETA . . . . . . . . . . local REAL 0 ó THETA < 1 ZERO. . . . . . . . . . . local REAL ZERO=0. L2CXFT Global Symbols NAME CLASS NUMBER of lines of code including comments ADDKNT. . . . . . . . . . external CFLGRF. . . . . . . . . . external DELKNT. . . . . . . . . . external L2CXFT. . . . . . . . . . SUBROUTINE 502 L2ILS . . . . . . . . . . external MESSGE. . . . . . . . . . external Subroutine L2ILS ---------------- Line# Source Line 1 C...................................................................... 2 C GIVEN N NOISY FUNCTION MEASUREMENTS F(I): I=1,2,...N 3 C AT THE ABSCISSAE X(1) THE INNER PRODUCT OF TWO VECTORS AND, 63 C THE COEFFICIENTS OF SX(.) ARE DETERMINED BY SOLVING THE 64 C NORMAL EQUATIONS 65 C D(I)*COEF(I-1)+E(I)*COEF(I)+D(I+1)*COEF(I+1)=B(I) 66 C FOR I=I1-1,I1,...,J 67 C WHERE 68 C D(I)=, I=I1-1,I1,...,J 69 C E(I)=, I=I1-1,I1,...,J 70 C B(I)=, I=I1-1,I1,...,J. 71 C 72 C THEN PROCEED AS IN REFERENCE QUOTED. 73 C...................................................................... 74 C NO TEST ON THE CONDITIONS N THE INNER PRODUCT OF TWO VECTORS AND, 28 C THE SPLINE COEFFICIENTS COEF(.) ARE DETERMINED BY SOLVING THE 29 C NORMAL EQUATIONS 30 C D(J)*COEF(J-1)+E(J)*COEF(J)+D(J+1)*COEF(J+1)=B(J) 31 C FOR J=I1-1,I1,...,NK 32 C WHERE 33 C D(J)=, J=I1-1,I1,...,NK 34 C E(J)=, J=I1-1,I1,...,NK 35 C B(J)=, J=I1-1,I1,...,NK. 36 C...................................................................... 37 I0=I1-1 38 C 39 C.....DERIVE COEF(.) BY FACTORIZATION 40 C 41 UE(I0)=E(I0) 42 COEF(I0)=B(I0) 43 DO 10 I=I0+1,NK 44 UE(I)=E(I)-D(I)/UE(I-1)*D(I) 45 10 COEF(I)=B(I)-D(I)*COEF(I-1)/UE(I-1) 46 COEF(NK)=COEF(NK)/UE(NK) 47 DO 20 I=NK-1,I0,-1 48 20 COEF(I)=(COEF(I)-D(I+1)*COEF(I+1))/UE(I) 49 C 50 RETURN 51 END COEFF Local Symbols NAME CLASS TYPE EXPLANATION N . . . . . . . . . . . . argument I1. . . . . . . . . . . . argument NK. . . . . . . . . . . . argument UE. . . . . . . . . . . . argument D . . . . . . . . . . . . argument B . . . . . . . . . . . . argument E . . . . . . . . . . . . argument COEF. . . . . . . . . . . argument I . . . . . . . . . . . . local INTEGER temporary I0. . . . . . . . . . . . local INTEGER I0=I1-1 COEFF Global Symbols NAME CLASS NUMBER of lines of code including comments COEFF . . . . . . . . . . SUBROUTINE 51 Subroutine SXATXI ----------------- Line# Source Line 1 C....................................................................... 2 C GENERATE SPLINE SX(.) AT X(.) FROM COEF(.) OVER AKN(.). 3 C....................................................................... 4 SUBROUTINE SXATXI(SX,X,COEF,AKN,IAKN,I1,N,NK) 5 * IMPLICIT DOUBLE PRECISION (A-H,O-Z) 6 DIMENSION SX(I1:N),X(I1:N),COEF(I1-1:N),AKN(I1-2:N),IAKN(I1-1:N) 7 C 8 C.... I N P U T .... 9 C X,COEF,AKN,IAKN,I1,N,NK AS THEY ARE DEFINED IN SUBROUTINE L2CXFT. 10 C 11 C.... O U T P U T .... 12 C SX(I1:N) REAL ARRAY CONTAINING THE SPLINE COMPONENTS AT X(I) 13 C I=I1,I1+1,...,N. 14 C 15 C.... M E T H O D (GENERATE SPLINE COMPONENTS FROM COEFFICIENTS) .... 16 C SEE ALGORITHM 5.8 OF L.L. SCHUMAKER, SPLINE FUNCTIONS: BASIC 17 C THEORY, J. WILEY AND SONS, 1981, P.194 . 18 C....................................................................... 19 I0=I1-1 20 SX(I1)=COEF(I0) 21 SX(N)=COEF(NK) 22 I=I1+1 23 L=I0 24 50 IF (I.LT.N .AND. IAKN(L).LT.I .AND. I.LE.IAKN(L+1)) THEN 25 DENOM=AKN(L+1)-AKN(L) 26 A1=(X(I)-AKN(L))/DENOM 27 A2=1-A1 28 SX(I)=COEF(L)*A2+COEF(L+1)*A1 29 I=I+1 30 GO TO 50 31 END IF 32 L=L+1 33 IF (L.LT.NK) GO TO 50 34 C 35 RETURN 36 END SXATXI Local Symbols NAME CLASS TYPE EXPLANATION NK. . . . . . . . . . . . argument N . . . . . . . . . . . . argument I1. . . . . . . . . . . . argument IAKN. . . . . . . . . . . argument AKN . . . . . . . . . . . argument COEF. . . . . . . . . . . argument X . . . . . . . . . . . . argument SX. . . . . . . . . . . . argument I . . . . . . . . . . . . local INTEGER temporary, I1 ó I ó N L . . . . . . . . . . . . local INTEGER temporary, I0 ó L ó NK A1. . . . . . . . . . . . local REAL temporary A2. . . . . . . . . . . . local REAL temporary DENOM . . . . . . . . . . local REAL temporary I0. . . . . . . . . . . . local INTEGER I0=I1-1 SXATXI Global Symbols NAME CLASS NUMBER of lines of code including comments SXATXI. . . . . . . . . . SUBROUTINE 36 Subroutine LAGRNG ----------------- Line# Source Line 1 C....................................................................... 2 C THIS SUBROUTINE CALCULATES THE LAGRANGE MULTIPLIERS PAR(.) 3 C ASSOCIATED WITH THE SPLINE SX(.) IN SUBROUTINE L2CXFT. 4 C....................................................................... 5 SUBROUTINE LAGRNG(PAR,SX,F,IACT,NACT,EC,DC,UE,I1,N) 6 * IMPLICIT DOUBLE PRECISION (A-H,O-Z) 7 DIMENSION SX(I1:N),F(I1:N),DC(I1:N),EC(I1:N),UE(I1-1:N), 8 1 IACT(I1:N),PAR(I1:N) 9 C 10 C.... I N P U T .... 11 C SX,F,IACT,DC,EC,I1,N,NACT AS THEY ARE DEFINED IN SUBROUTINE 12 C L2CXFT. 13 C 14 C.... O U T P U T .... 15 C PAR(I1:N-2) REAL ARRAY CONTAINING THE LAGRANGE MULTIPLIERS SO 16 C THAT PAR(IACT(I)), I=I1,I1+1,...,NACT ARE THE MULTI- 17 C PLIERS OF THE ACTIVE CONSTRAINTS WHILE THE REMAINING 18 C MULTIPLIERS EQUAL TO ZERO. 19 C 20 C.... W O R K I N G S P A C E .... 21 C UE(I1-1:NK+1) REAL ARRAY THAT PROVIDES WORKING SPACE FOR THE 22 C FACTORIZATION THAT GIVES PAR(.). IT KEEPS THE 23 C DIAGONAL ELEMENTS OF THE FORWARD SWEEP OF THE 24 C FACTORIZATION. 25 C 26 C.... M E T H O D (CALCULATE THE LAGRANGE MULTIPLIERS) ..... 27 C THE FIRST ORDER OPTIMALITY CONDITIONS ARE 28 C SUM(PAR(II)*A(II,J))=2*(SX(II+1)-F(II+1)), 29 C FOR II=IACT(I), I=I1,...,NACT 30 C WHERE A(II,J) DENOTES THE J-TH COMPONENT OF THE II-TH SECOND 31 C DIVIDED DIFFERENCE. ONLY THE A(II,J) FOR J=II,II+1,II+2 ARE 32 C NON-ZERO. THESE CONDITIONS CONSTITUTE AN OVERDETERMINED SYSTEM 33 C OF EQUATIONS. CHOOSE THE EQUATIONS THAT CORRESPOND TO THE SECOND 34 C NON-ZERO COMPONENT OF EACH A(II,.) AND SOLVE BY GAUSS ELIMINATION 35 C WITHOUT PIVOTING THE DERIVED BLOCK DIAGONAL POSITIVE DEFINITE 36 C SYSTEM OF EQUATIONS. PROCEED AS IT IS EXPLAINED IN THE TEXTUAL 37 C PART OF THIS ALGORITHM. 38 C....................................................................... 39 IF (NACT.LE.I1-1) GO TO 170 40 C 41 C THE GRADIENT OF THE OBJECTIVE FUNCTION IS DEFINED BY, SAY, 42 C GRAD(I)=2*(SX(I+1)-F(I+1)), I=I1,I1+1,...,N-2. HOWEVER, INSTEAD 43 C OF AN ARRAY GRAD(.) THE SUBROUTINE USES DIRECTLY THE VALUES OF 44 C THE GRADIENT AT X(.) WHEN NECESSARY. 45 C 46 C CALCULATE SIZE K OF A BLOCK OF THE LAGRANGE COEFFICIENT 47 C MATRIX. 48 C 49 I=I1 50 110 K=1 51 120 IF (I.LT.NACT .AND. IACT(I).EQ.IACT(I+1)-1) THEN 52 I=I+1 53 K=K+1 54 GO TO 120 55 END IF 56 C 57 130 IF (K.EQ.1) THEN 58 II=IACT(I) 59 PAR(II)=(SX(II+1)-F(II+1))/EC(II) 60 ELSE 61 C 62 C FACTORIZE THE K*K 3/DIAG SYSTEM OF LAGRANGE COEF. 63 C 64 UE(I-K+1)=EC(IACT(I-K+1)) 65 DO 140 L=I-K+2,I 66 DI=DC(IACT(L)) 67 140 UE(L)=EC(IACT(L))-DI/UE(L-1)*DI 68 C 69 C USE THIS FACTORIZATION TO SOLVE FOR PAR(.). 70 C 71 GRAD=SX(IACT(I-K+1)+1)-F(IACT(I-K+1)+1) 72 PAR(IACT(I-K+1))=GRAD 73 DO 150 L=I-K+2,I 74 GRAD=SX(IACT(L)+1)-F(IACT(L)+1) 75 150 PAR(IACT(L))=GRAD-DC(IACT(L))*PAR(IACT(L-1))/UE(L-1) 76 PAR(IACT(I))=PAR(IACT(I))/UE(I) 77 DO 160 L=I-1,I-K+1,-1 78 160 PAR(IACT(L))=PAR(IACT(L))/UE(L)-DC(IACT(L+1))* 79 1 PAR(IACT(L+1))/UE(L) 80 END IF 81 I=I+1 82 IF (I-NACT) 110,110,170 83 C 84 170 RETURN 85 END LAGRNG Local Symbols NAME CLASS TYPE EXPLANATIONS N . . . . . . . . . . . . argument I1. . . . . . . . . . . . argument UE. . . . . . . . . . . . argument DC. . . . . . . . . . . . argument EC. . . . . . . . . . . . argument NACT. . . . . . . . . . . argument IACT. . . . . . . . . . . argument F . . . . . . . . . . . . argument SX. . . . . . . . . . . . argument PAR . . . . . . . . . . . argument I . . . . . . . . . . . . local INTEGER temporary, I1 ó I ó N K . . . . . . . . . . . . local INTEGER matrix size, 1 ó K ó NACT-I0 L . . . . . . . . . . . . local INTEGER temporary DI. . . . . . . . . . . . local REAL temporary II. . . . . . . . . . . . local INTEGER temporary GRAD. . . . . . . . . . . local REAL GRADient evaluation LAGRNG Global Symbols NAME CLASS NUMBER of lines of code including comments LAGRNG. . . . . . . . . . SUBROUTINE 85 Subroutine ADDKNT ----------------- Line# Source Line 1 C...................................................................... 2 C THIS SUBROUTINE INSERTS A KNOT, NAMELY THE DATA POINT X(INKNOT), 3 C INTO THE SPLINE KNOT SEQUENCE OF THE THE LINEAR SPLINE 4 C APPROXIMATION SX(.) DEFINED IN SUBROUTINE L2CXFT. 5 C 6 C CALLS SUBROUTINE PLUS. 7 C...................................................................... 8 SUBROUTINE ADDKNT(INKNOT,E,B,D,X,F,AKN,IAKN,IACT,NK,NACT,I1,N) 9 * IMPLICIT DOUBLE PRECISION (A-H,O-Z) 10 DIMENSION E(I1-1:N),B(I1-1:N),D(I1-1:N),X(I1:N),F(I1:N), 11 1 AKN(I1-2:N),IAKN(I1-1:N),IACT(I1:N) 12 C 13 C.... I N P U T .... 14 C INKNOT INTEGER VARIABLE, SUCH THAT X(INKNOT) IS INSERTED 15 C INTO THE KNOT SEQUENCE AKN(.). 16 C E,B,D,X,F,AKN,IAKN,IACT,I1,N,NK,NACT AS THEY ARE 17 C DEFINED IN SUBROUTINE L2CXFT. 18 C 19 C.... O U T P U T .... 20 C E,B,D,AKN,IAKN,IACT AFTER INSERTING X(INKNOT) INTO AKN(.). 21 C NK INCREASED BY ONE. 22 C NACT DECREASED BY ONE. 23 C 24 C.... M E T H O D .... 25 C FIND POSITION J IN THE KNOT SEQUENCE AKN(.) SUCH THAT 26 C AKN(J-1) < X(INKNOT) < AKN(J) BEFORE INSERTION AND 27 C AKN(J)=X(INKNOT) AFTER INSERTION. 28 C CALL SUBROUTINE PLUS TO MODIFY THE NORMAL EQUATIONS THAT ARE 29 C AFFECTED BY THE KNOT INSERTION. HENCE THE 'HAT' BASIS FUNCTIONS 30 C AFFECTED BY THE INSERTION HAVE TO BE DETERMINED. THE J-TH BASIS 31 C FUNCTION IS DEFINED FROM X(IP) TO X(IR), AND IS CENTERED AT 32 C X(IQ). ALSO THE (J-1)-TH BASIS FUNCTIONS IS CENTERED AT X(IP) AND 33 C IS DEFINED FROM X(IPL) TO X(IQ). THE (J+1)-TH BASIS FUNCTION 34 C IS CENTERED AT X(IR) AND IS DEFINED FROM X(IQ) TO X(IRR). 35 C REVISE IACT(.), THE INDICES OF THE ACTIVE CONSTRAINTS. 36 C...................................................................... 37 I0=I1-1 38 C 39 C GIVEN X(INKNOT), APPLY BISECTION IN [I0,NK) IN ORDER TO FIND IP 40 C SUCH THAT AKN(IP) < X(INKNOT) < AKN(IP+1). 41 C 42 IP=I0 43 IR=NK 44 20 IF (IR-IP.EQ.1) GO TO 30 45 IQ=(IP+IR)/2 46 IF (INKNOT .LT. IAKN(IQ)) THEN 47 IR=IQ 48 ELSE 49 IP=IQ 50 END IF 51 GO TO 20 52 30 CONTINUE 53 C 54 DO 100 I=NK,IP+1,-1 55 B(I+1)=B(I) 56 E(I+1)=E(I) 57 D(I+1)=D(I) 58 AKN(I+1)=AKN(I) 59 100 IAKN(I+1)=IAKN(I) 60 J=IP+1 61 E(J)=0. 62 B(J)=0. 63 D(J)=0. 64 NK=NK+1 65 AKN(J)=X(INKNOT) 66 IAKN(J)=INKNOT 67 C 68 C DEFINE SUPPORTS OF NEW BASIS FUNCTIONS THAT FOLLOW THE KNOT 69 C INSERTION. 70 C 71 IF (J.EQ.I1) THEN 72 IPL=I1 73 ELSE 74 IPL=IAKN(J-2) 75 END IF 76 C 77 IF (J.EQ.NK-1) THEN 78 IRR=N 79 ELSE 80 IRR=IAKN(J+2) 81 END IF 82 C 83 IP=IAKN(J-1) 84 IQ=IAKN(J) 85 IR=IAKN(J+1) 86 C 87 C MODIFY NORMAL EQUATIONS. 88 C 89 CALL PLUS(X,F,AKN,E,B,D,NK,I1,N,IPL,IP,IQ,IR,IRR,J) 90 C 91 C REVISE ACTIVE CONSTRAINTS. 92 C NACT>=I1 BECAUSE THE KNOT INSERTION SHOWS THAT THERE EXISTS 93 C AT LEAST ONE ACTIVE CONSTRAINT. 94 C 95 IF (NACT.EQ.I1) GO TO 140 96 C 97 C I1 I1 and / IAKN(J+2), if J < NK-1 IRR = \ N , if J = NK-1. In effect the (J-1)-th basis function is defined on [X(IPL), X(IQ)], the J-th basis function on [X(IP),X(IR)] and the (J+1)-th basis function on [X(IQ),X(IRR]. Therefore subroutine PLUS is safely applied to calculate E(ll) and B(ll) for ll=J-1,J,J+1 and D(ll) for ll=J,J+1. Finally the instructions that include label 110 revise IACT(.) and reduce NACT by one. The work required by subroutine ADDKNT amounts to at most O(N-I1) operations. SUBROUTINE PLUS --------------- This subroutine is called by subroutine ADDKNT to revise the equations (2.4) after inserting the data point X(INKNOT) into AKN(.). The arguments J, IPL, IP, IR and IRR have already been defined in subroutine ADDKNT. The calculations that are employed in this subroutine are based on the definitions of the arrays E, D, B that follow the equation (2.4). The work required by subroutine PLUS amounts to O(IR-IP) operations. SUBROUTINE DELKNT ----------------- This subroutine removes the knot AKN(KNTOUT), I1-1