C @(#)tdspline.for 17.1.1.1 (ES0-DMD) 01/25/02 17:47:17 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 13:19 - 15 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: J.D.PONZ C. Modif 30-Jul-90 P.Ballester C C 010123 last modif C C------------------------------------------------------------------- SUBROUTINE SMOOT(X,Y,W,Q,M,XB,XE,K,S,N,T,C,FP,IOPT,IER,NRDATA, + FP0,FPOLD,NPLUS) C C IMPLICIT NONE C C GIVEN THE SET OF DATA POINTS (X(I),Y(I)) AND THE SET OF POSITIVE C NUMBERS W(I),I=1,2,...M, SMOOT DETERMINES A SMOOTH SPLINE C APPROXIMATION OF DEGREE K FOR THE INTERVAL XB <= X <= XE. C THE NUMBER OF KNOTS OF THIS SPLINE S(X) AND THE POSITION C T(J),J=1,...N IS CHOSEN AUTOMATICALLY BY THE ROUTINE. C THE SMOOTHNESS OF S(X) IS ACHIEVED BY MINIMALIZING THE SUM(D(R)**2) C WHERE D(R) STANDS FOR THE DISCONTINUITY JUMP OF THE KTH DERIVATIVE C OF S(X) AT THE KNOT T(R),R=K+2,...,N-K-1. THE AMOUNT OF SMOOTH- C NESS IS DETERMINED BY THE CONDITION THAT F(P)=SUM(W(I)*(Y(I)- C S(X(I)))**2) BE <= S, WITH S A GIVEN NON-NEGATIVE CONSTANT. C THE FIT S(X) IS GIVEN IN ITS B-SPLINE REPRESENTATION ( B-SPLINE C COEFFICIENTS C(J),J=1,2,...N-K-1) AND CAN BE EVALUATED BY MEANS C OF FUNCTION DERIV. C C CALLING SEQUENCE: C CALL SMOOT(X,Y,W,Q,M,XB,XE,K,S,N,T,C,FP,IOPT,IER) C C INPUT PARAMETERS: C X : ARRAY,MINIMUM LENGTH M,CONTAINING THE DATA POSITIONS X(I). C Y : ARRAY,MINIMUM LENGTH M,CONTAINING THE DATA VALUES Y(I). C W : ARRAY,MINIMUM LENGTH M,CONTAINING THE WEIGHTS W(I). C M : INTEGER VALUE,CONTAINING THE NUMBER OF DATA POINTS. C XB,XE:REAL VALUES, CONTAINING THE BOUNDARIES OF THE APPR. INTERVAL C K : INTEGER VALUE,CONTAINING THE DEGREE OF S(X). C S : REAL VALUE,CONTAINING THE SMOOTHING FACTOR. C IOPT: INTEGER VALUE WHICH TAKES THE VALUE 0 OR 1. C IOPT=0: THE ROUTINE WILL RESTART ALL COMPUTATIONS. C IOPT=1: THE ROUTINE WILL START WITH THE KNOTS FOUND AT THE C LAST CALL OF THE ROUTINE. IF IOPT=1 THE OUTPUT C PARAMETERS T AND N ARE INPUT PARAMETERS AS WELL. C C OUTPUT PARAMETERS: C T : ARRAY,LENGTH NEST (SEE DATA INITIALIZATION STATEMENT), C WHICH CONTAINS THE POSITION OF THE KNOTS,I.E. THE POSITION C OF THE INTERIOR KNOTS T(K+2),...T(N-K-1), AS WELL AS THE C POSITION OF THE KNOTS T(1)=T(2)=...=T(K+1)=XB AND XE = C T(N-K)=...=T(N) WHICH ARE NEEDED FOR THE B-SPLINE REPRESENT. C C : ARRAY,LENGTH NEST,WHICH CONTAINS THE B-SPLINE COEFFICIENTS. C N : INTEGER VALUE,CONTAINING THE TOTAL NUMBER OF KNOTS. C FP : REAL VALUE,WHICH CONTAINS THE SUM(WI*(YI-S(XI))**2),I=1,...M C IER : ERROR CODE C IER=0: NORMAL RETURN. C IER=-1:NORMAL RETURN, S(X) IS AN INTERPOLATING SPLINE. C IER=-2:NORMAL RETURN, S(X) IS A POLYNOMIAL OF DEGREE K. C IER>0 :ABNORMAL TERMINATION. C IER=1:THE REQUIRED STORAGE SPACE EXCEEDS THE AVAILABLE C STORAGE SPACE,SPECIFIED BY THE PARAMETER NEST. C PROBABLY CAUSES:NEST OR S TOO SMALL. C IER=2:A THEORETICALLY IMPOSSIBLE RESULT WAS FOUND DURING C THE ITERATION PROCESS. C PROBABLY CAUSES:TOL TOO SMALL C IER=3:THE MAXIMAL NUMBER OF ITERATIONS MAXIT HAS BEEN C REACHED. C PROBABLY CAUSES:MAXIT OR TOL TOO SMALL. C IER=10:SOME OF THE INPUTDATA ARE INVALID(SEE RESTRICTIONS). C C C INPUT/OUTPUT PARAMETERS: C NRDATA : ARRAY, LENGTH NEST KEEP INFO BETWEEN ITERATIONS C FPO : REAL, KEEP INFO C FPOLD : REAL, KEEP INFO C NPLUS : REAL, KEEP INFO C C RESTRICTIONS: C 1) M > K > 0 C 2) XB <= X(R) < X(R+1) <= XE, R=1,2,...M-1. C 3) W(R) > 0, R=1,2,...M. C 4) S >= 0. C 5) NEST >= 2*K+2. C C OTHER SUBROUTINES REQUIRED: C BSPLIN,COSSIN,ROTATE,BACK,NKNOT,DISCO AND RATION. INTEGER MXP PARAMETER (MXP=1000) INTEGER M,N,IOPT,IER,K,I,J,L,NMAX,NEST,MAXIT,NEW,ICHECK,N8 INTEGER NPLUS,MK1,K2,K3,K1,NMIN,ITER,NRINT,NK1,IT,I1,I2,I3 INTEGER NPL1,L0,NRDATA(MXP) REAL X(M),Y(M),W(M),T(MXP),C(MXP),FPINT(MXP), + Z(MXP),DIAG(MXP),DPRIME(MXP),A(MXP,5),G(MXP,6),B(MXP,7), + Q(M,6),H(7),FP,XE,XB,S,FP0,FPOLD,TOL REAL ACC,XI,YI,WI,PIV,FPART,TERM,STORE,P1 REAL PINV,P,COS,SIN,FPMS,F1,F2,F3,P2,P3,RATION C C NRDATA: INTEGER ARRAY,LENGTH NEST,WHICH GIVES THE NUMBER OF C DATA POINTS INSIDE EACH KNOT INTERVAL. C FP0 : REAL VALUE,WHICH CONTAINS THE WEIGHTED SUM OF SQUARES C OF RESIDUAL RIGHT HAND SIDES, FOR THE LEAST SQUARES C POLYMOMIAL OF DEGREE K. C FPOLD : REAL VALUE,WHICH CONTAINS THE WEIGHTED SUM OF RESIDUAL C RIGHT HAND SIDES FOR THE LEAST SQUARES SPLINE WHICH C CORRESPONDS TO THE LAST FOUND SET OF KNOTS BUT ONE. C NPLUS : INTEGER VALUE,GIVING THE NUMBER OF KNOTS OF THE LAST C SET MINUS THE NUMBER OF THE LAST SET BUT ONE. C C DATA INITIALIZATION STATEMENT TO SPECIFY C TOL : THE REQUESTED RELATIVE ACCURACY FOR THE ROOT OF F(P) = S. C MAXIT: THE MAXIMAL NUMBER OF ITERATIONS ALLOWED. C NEST : AN OVER-ESTIMATE OF THE NUMBER OF KNOTS N. THIS PARAMETER C MUST BE SET BY THE USER TO INDICATE THE STORAGE SPACE C AVAILABLE TO THE SUBROUTINE. THE DIMENSION SPECIFICATIONS C OF THE ARRAYS T,C,NRDATA,FPINT,Z,DIAG,DPRIME(N),A(N,K), C Q(M,K+1),G(N,K+1),B(N,K+2) AND H(K+2) DEPEND C ON N,M AND K. SINCE N IS UNKNOWN AT THE TIME THE C USER SETS UP THE DIMENSION INFORMATION AN OVER-ESTIMATE C OF THESE ARRAYS WILL GENERALLY BE MADE. THE FOLLOWING C REMARKS ARE INTENDED TO HELP THE USER C 1) 2*K+2 <= N <= M+K+1 C 2) THE SMALLER THE VALUE OF S, THE GREATER N WILL BE. C 3) NORMALLY N = M/2 IS AN OVER-ESTIMATE. DATA TOL/0.001/,MAXIT/40/,NEST/MXP/ ! Old values TOL=0.001 and MAXIT=40 C BEFORE STARTING COMPUTATIONS A DATA CHECK IS MADE. IF THE INPUT C DATA ARE INVALID CONTROLE IS IMMEDIATELY REPASSED TO THE DRIVER C PROGRAM (IER=10). 10 IER = 0 K1 = K + 1 K2 = K1 + 1 NMIN = 2*K1 IF (K.LE.0) IER = 10 IF (M.LT.K1 .OR. NEST.LT.NMIN) IER = 10 IF (S.LT.0.) IER = 10 IF (XB.GT.X(1) .OR. XE.LT.X(M) .OR. W(1).LE.0.) IER = 10 IF (IER.NE.0) GO TO 400 C CALCULATION OF ACC, THE ABSOLUTE TOLERANCE FOR THE ROOT OF F(P)=S. ACC = TOL*S CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PART 1: DETERMINATION OF THE NUMBER OF KNOTS AND THEIR POSITION C C ************************************************************** C C GIVEN A SET OF KNOTS WE COMPUTE THE LEAST-SQUARES SPLINE SINF(X). C C IF THE SUM F(P=INF) <= S WE ACCEPT THE CHOICE OF KNOTS. C C OTHERWISE WE HAVE TO INCREASE THEIR NUMBER. C C THE INITIAL CHOICE OF KNOTS DEPENDS ON THE VALUE OF S AND IOPT. C C IF S=0 WE HAVE SPLINE INTERPOLATION; IN THAT CASE THE NUMBER OF C C KNOTS EQUALS NMAX = M+K+1. C C IF S > 0 AND C C IOPT=0 WE FIRST COMPUTE THE LEAST-SQUARES POLYNOMIAL OF C C DEGREE K; N = NMIN = 2*K+2 C C IOPT=1 WE START WITH THE SET OF KNOTS FOUND AT THE LAST C C CALL OF THE ROUTINE, EXCEPT FOR THE CASE THAT S > FP0; THEN C C WE COMPUTE DIRECTLY THE LEAST-SQUARES POLYNOMIAL OF DEGREE K. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C DETERMINE NMAX, THE NUMBER OF KNOTS FOR SPLINE INTERPOLATION. NMAX = M + K1 IF (S.GT.0.) GO TO 50 C IF S=0, S(X) IS AN INTERPOLATING SPLINE. N = NMAX C TEST WHETHER THE REQUIRED STORAGE SPACE EXCEEDS THE AVAILABLE ONE. IF (N.GT.NEST) GO TO 380 C FIND THE POSITION OF THE INTERIOR KNOTS IN CASE OF INTERPOLATION. MK1 = M - K1 IF (MK1.EQ.0) GO TO 70 K3 = K/2 I = K2 J = K3 + 2 IF (K3*2.EQ.K) GO TO 30 DO 20 L = 1,MK1 T(I) = X(J) I = I + 1 J = J + 1 20 CONTINUE GO TO 70 30 DO 40 L = 1,MK1 T(I) = (X(J)+X(J-1))*0.5 ! 30-JUl-90 I = I + 1 J = J + 1 40 CONTINUE GO TO 70 C IF S>0 OUR INITIAL CHOICE OF KNOTS DEPENDS ON THE VALUE OF IOPT. C IF IOPT=0 OR IOPT=1 AND S>=FP0, WE START COMPUTING THE LEAST-SQUARES C POLYNOMIAL OF DEGREE K WHICH IS A SPLINE WITHOUT INTERIOR KNOTS. C IF IOPT=1 AND FP0>S WE START COMPUTING THE LEAST SQUARES SPLINE C ACCORDING TO THE SET OF KNOTS FOUND AT THE LAST CALL OF THE ROUTINE. 50 IF (IOPT.LE.0) GO TO 60 IF (FP0.GT.S) GO TO 70 60 N = NMIN NRDATA(1) = M - 2 C MAIN LOOP FOR THE DIFFERENT SETS OF KNOTS. M IS A SAVE UPPER BOUND C FOR THE NUMBER OF TRIALS. 70 DO 220 ITER = 1,M IF (N.EQ.NMIN) IER = -2 C FIND NRINT, TNE NUMBER OF KNOT INTERVALS. NRINT = N - NMIN + 1 C FIND THE POSITION OF THE ADDITIONAL KNOTS WHICH ARE NEEDED FOR C THE B-SPLINE REPRESENTATION OF S(X). NK1 = N - K1 I = N DO 80 J = 1,K1 T(J) = XB T(I) = XE I = I - 1 80 CONTINUE C COMPUTE THE B-SPLINE COEFFICIENTS OF THE LEAST-SQUARES SPLINE C SINF(X). THE OBSERVATION MATRIX A IS BUILT UP ROW BY ROW AND C REDUCED TO UPPER TRIANGULAR FORM BY GIVENS TRANSFORMATIONS C WITHOUT SQUARE ROOTS. AT THE SAME TIME FP=F(P=INF) IS COMPUTED. FP = 0. C INITIALIZE THE OBSERVATION MATRIX A. DO 100 I = 1,NK1 DIAG(I) = 0. Z(I) = 0. DO 90 J = 1,K A(I,J) = 0. 90 CONTINUE 100 CONTINUE L = K1 DO 150 IT = 1,M C FETCH THE CURRENT DATA POINT X(IT),Y(IT). XI = X(IT) YI = Y(IT) WI = W(IT) C SEARCH FOR KNOT INTERVAL T(L) <= XI < T(L+1). IF (XI.GE.T(L+1) .AND. L.NE.NK1) L = L + 1 C EVALUATE THE (K+1) NON-ZERO B-SPLINES AT XI AND STORE THEM IN Q. CALL BSPLIN(T,N,K,XI,L,H) DO 110 I = 1,K1 IF (H(I).LT.0.1E-07) H(I) = 0. Q(IT,I) = H(I) 110 CONTINUE C ROTATE THE NEW ROW OF THE OBSERVATION MATRIX INTO TRIANGLE BY C GIVENS TRANSFORMATIONS WITHOUT SQUARE ROOTS. J = L - K1 DO 130 I = 1,K1 IF (WI.LE.1.E-30) GO TO 150 J = J + 1 PIV = H(I) IF (PIV.EQ.0.) GO TO 130 C CALCULATE THE PARAMETERS OF THE GIVENS TRANSFORMATION. CALL COSSIN(PIV,WI,DIAG(J),COS,SIN) C TRANSFORMATIONS TO RIGHT HAND SIDE. CALL ROTATE(PIV,COS,SIN,YI,Z(J)) IF (I.EQ.K1) GO TO 140 I2 = 0 I3 = I + 1 DO 120 I1 = I3,K1 I2 = I2 + 1 C TRANSFORMATIONS TO LEFT HAND SIDE. CALL ROTATE(PIV,COS,SIN,H(I1),A(J,I2)) 120 CONTINUE 130 CONTINUE C ADD CONTRIBUTION OF THIS ROW TO THE SUM OF SQUARES OF RESIDUAL C RIGHT HAND SIDES. 140 FP = FP + WI*YI**2 150 CONTINUE IF (IER.EQ.-2) FP0 = FP C BACKWARD SUBSTITUTION TO OBTAIN THE B-SPLINE COEFFICIENTS. CALL BACK(A,Z,NK1,K,C) C TEST WHETHER THE APPROXIMATION SINF(X) IS AN ACCEPTABLE SOLUTION. FPMS = FP - S IF (ABS(FPMS).LT.ACC) GO TO 400 C IF F(P=INF) < S ACCEPT THE CHOICE OF KNOTS. IF (FPMS.LT.0.) GO TO 230 C IF N = NMAX, SINF(X) IS AN INTERPOLATING SPLINE. IF (N.EQ.NMAX) GO TO 390 C INCREASE THE NUMBER OF KNOTS. C IF N=NEST WE CANNOT INCREASE THE NUMBER OF KNOTS BECAUSE OF C THE STORAGE CAPACITY LIMITATION. IF (N.EQ.NEST) GO TO 380 C DETERMINE THE NUMBER OF KNOTS NPLUS WE ARE GOING TO ADD. IF (IER.EQ.0) GO TO 160 NPLUS = 1 IER = 0 GO TO 170 160 NPL1 = NPLUS*2 IF (FPOLD-FP.GT.ACC) NPL1 = FLOAT(NPLUS)*FPMS/ (FPOLD-FP) NPLUS = MIN0(NPLUS*2,MAX0(NPL1,NPLUS/2,1)) 170 FPOLD = FP C COMPUTE THE SUM(W(I)*(Y(I)-S(X(I)))**2) FOR EACH KNOT INTERVAL C T(J+K) <= X(I) <= T(J+K+1) AND STORE IT IN FPINT(J),J=1,2,...NRINT. FPART = 0. I = 1 L = K2 NEW = 0 DO 200 IT = 1,M IF (X(IT).LT.T(L) .OR. L.GT.NK1) GO TO 180 NEW = 1 L = L + 1 180 TERM = 0. L0 = L - K2 DO 190 J = 1,K1 L0 = L0 + 1 TERM = TERM + C(L0)*Q(IT,J) 190 CONTINUE TERM = W(IT)* (TERM-Y(IT))**2 FPART = FPART + TERM IF (NEW.EQ.0) GO TO 200 STORE = TERM*0.5 FPINT(I) = FPART - STORE I = I + 1 FPART = STORE NEW = 0 200 CONTINUE FPINT(NRINT) = FPART DO 210 L = 1,NPLUS C ADD A NEW KNOT. CALL NKNOT(X,M,T,N,FPINT,NRDATA,NRINT) C TEST WHETHER WE CANNOT FURTHER INCREASE THE NUMBER OF KNOTS. IF (N.EQ.NMAX .OR. N.EQ.NEST) GO TO 220 210 CONTINUE C RESTART THE COMPUTATIONS WITH THE NEW SET OF KNOTS. 220 CONTINUE C TEST WHETHER THE LEAST-SQUARES KTH DEGREE POLYNOMIAL IS A SOLUTION C OF OUR APPROXIMATION PROBLEM. 230 IF (IER.EQ.-2) GO TO 400 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PART 2: DETERMINATION OF THE SMOOTHING SPLINE SP(X). C C *************************************************** C C WE HAVE DETERMINED THE NUMBER OF KNOTS AND THEIR POSITION. C C WE NOW COMPUTE THE B-SPLINE COEFFICIENTS OF THE SMOOTHING SPLINE C C SP(X). THE OBSERVATION MATRIX A IS EXTENDED BY THE ROWS OF MATRIX C C B EXPRESSING THAT THE KTH DERIVATIVE DISCONTINUITIES OF SP(X) AT C C THE INTERIOR KNOTS T(K+2),...T(N-K-1) MUST BE ZERO. THE CORRES- C C PONDING WEIGHTS OF THESE ADDITIONAL ROWS ARE SET TO 1/SQRT(P). C C ITERATIVELY WE THEN HAVE TO DETERMINE THE VALUE OF P SUCH THAT C C F(P)=SUM(W(I)*(Y(I)-SP(X(I)))**2) BE = S. WE ALREADY KNOW THAT C C THE LEAST-SQUARES KTH DEGREE POLYNOMIAL CORRESPONDS TO P=0, AND C C THAT THE LEAST-SQUARES SPLINE CORRESPONDS TO P=INFINITY. THE C C ITERATION PROCESS WHICH IS PROPOSED HERE, MAKES USE OF RATIONAL C C INTERPOLATION. SINCE F(P) IS A CONVEX AND STRICTLY DECREASING C C FUNCTION OF P, IT CAN BE APPROXIMATED BY A RATIONAL FUNCTION C C R(P) = (U*P+V)/(P+W). THREE VALUES OF P(P1,P2,P3) WITH CORRESPOND- C C ING VALUES OF F(P) (F1=F(P1)-S,F2=F(P2)-S,F3=F(P3)-S) ARE USED C C TO CALCULATE THE NEW VALUE OF P SUCH THAT R(P)=S. CONVERGENCE IS C C GUARANTEED BY TAKING F1>0 AND F3<0. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C EVALUATE THE DISCONTINUITY JUMP OF THE KTH DERIVATIVE OF THE C B-SPLINES AT THE KNOTS T(L),L=K+2,...N-K-1 AND STORE IN B. CALL DISCO(T,N,K2,B) C INITIAL VALUE FOR P. P1 = 0. P2 = 0. F1 = FP0 - S P3 = -1. F3 = FPMS P = -F1/F3 ICHECK = 0 N8 = N - NMIN C ITERATION PROCESS TO FIND THE ROOT OF F(P) = S. DO 350 ITER = 1,MAXIT C THE ROWS OF MATRIX B WITH WEIGHT 1/SQRT(P) ARE ROTATED INTO C THE TRIANGULARISED OBSERVATION MATRIX A WHICH IS STORED IN G. IF (ABS(P).LT.1.0E-30) THEN !P = 0 => PINV = 1/P is in trouble GOTO 360 !indicate that iteration failed ENDIF C PINV = 1.0/P DO 250 I = 1,NK1 DPRIME(I) = DIAG(I) C(I) = Z(I) G(I,K1) = 0. DO 240 J = 1,K G(I,J) = A(I,J) 240 CONTINUE 250 CONTINUE DO 290 IT = 1,N8 C*PL*ERROR* Comment line too long C THE ROW OF MATRIX B IS ROTATED INTO TRIANGLE BY GIVENS TRANSFORMATIONS DO 260 I = 1,K2 H(I) = B(IT,I) 260 CONTINUE YI = 0. WI = PINV DO 280 J = IT,NK1 IF (WI.LE.1E-30) GO TO 290 PIV = H(1) C CALCULATE THE PARAMETERS OF THE GIVENS TRANSFORMATION. CALL COSSIN(PIV,WI,DPRIME(J),COS,SIN) C TRANSFORMATIONS TO RIGHT HAND SIDE. CALL ROTATE(PIV,COS,SIN,YI,C(J)) IF (J.EQ.NK1) GO TO 290 I2 = K1 IF (J.GT.N8) I2 = NK1 - J DO 270 I = 1,I2 C TRANSFORMATIONS TO LEFT HAND SIDE. CALL ROTATE(PIV,COS,SIN,H(I+1),G(J,I)) H(I) = H(I+1) 270 CONTINUE H(I2+1) = 0. 280 CONTINUE 290 CONTINUE C BACKWARD SUBSTITUTION TO OBTAIN THE B-SPLINE COEFFICIENTS. CALL BACK(G,C,NK1,K1,C) C COMPUTATION OF F(P). FP = 0. L = K2 DO 320 IT = 1,M IF (X(IT).LT.T(L) .OR. L.GT.NK1) GO TO 300 L = L + 1 300 L0 = L - K2 TERM = 0. DO 310 J = 1,K1 L0 = L0 + 1 TERM = TERM + C(L0)*Q(IT,J) 310 CONTINUE FP = FP + W(IT)* (TERM-Y(IT))**2 320 CONTINUE C TEST WHETHER THE APPROXIMATION SP(X) IS AN ACCEPTABLE SOLUTION. FPMS = FP - S C TYPE*,ITER,FPMS,ACC IF (ABS(FPMS).LT.ACC) GO TO 400 C TEST WHETHER THE MAXIMAL NUMBER OF ITERATIONS IS REACHED. IF (ITER.EQ.MAXIT) GO TO 360 C CARRY OUT ONE MORE STEP OF THE ITERATION PROCESS. P2 = P F2 = FPMS IF (ICHECK.NE.0) GO TO 340 IF ((F2-F3).GT.ACC) GO TO 330 C OUR INITIAL CHOICE OF P IS TOO LARGE. P = P*0.1E-02 P3 = P2 F3 = F2 GO TO 350 330 IF ((F1-F2).GT.ACC) GO TO 340 C OUR INITIAL CHOICE OF P IS TOO SMALL P = P*0.1E+04 P1 = P2 F1 = F2 GO TO 350 C TEST WETHER THE ITERATION PROCESS PROCEEDS AS THEORETICALLY C EXPECTED. 340 IF (F2.GE.F1 .OR. F2.LE.F3) THEN TOL = TOL*4. IF (TOL.GT.1.5) GO TO 370 GO TO 10 END IF C ICHECK = 1 C FIND THE NEW VALUE FOR P. P = RATION(P1,F1,P2,F2,P3,F3) C 350 CONTINUE C ERROR CODES AND MESSAGES. 360 IER = 3 GO TO 400 370 IER = 2 GO TO 400 380 IER = 1 GO TO 400 390 IER = -1 400 RETURN C C END SUBROUTINE BSPLIN(T,N,K,X,L,H) C SUBROUTINE BSPLIN EVALUATES THE (K+1) NON-ZERO B-SPLINES OF C DEGREE K AT T(L) <= X < T(L+1) USING THE STABLE RECURRENCE C RELATION OF DE BOOR AND COX. C THE DIMENSION SPECIFICATIONS OF THE FOLLOWING ARRAYS MUST BE C AT LEAST H(K+1),HH(K). C IMPLICIT NONE INTEGER N,J,K,I,L,LJ,LI REAL T(N),H(6),HH(5),X,F H(1) = 1. DO 30 J = 1,K DO 10 I = 1,J HH(I) = H(I) 10 CONTINUE H(1) = 0. DO 20 I = 1,J LI = L + I LJ = LI - J F = HH(I)/ (T(LI)-T(LJ)) H(I) = H(I) + F* (T(LI)-X) H(I+1) = F* (X-T(LJ)) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE COSSIN(PIV,WI,WW,COS,SIN) C SUBROUTINE COSSIN CALCULATES THE PARAMETERS OF A GIVENS C TRANSFORMATION WITHOUT SQUARE ROOTS. C IMPLICIT NONE REAL COS,SIN,PIV,WI,WW,STORE,DD STORE = PIV*WI DD = WW + STORE*PIV COS = WW/DD SIN = STORE/DD WW = DD WI = COS*WI RETURN END SUBROUTINE ROTATE(PIV,COS,SIN,A,B) C SUBROUTINE ROTATE APPLIES A GIVENS ROTATION TO A AND B. C IMPLICIT NONE REAL PIV,COS,SIN,A,B,STORE STORE = B B = COS*STORE + SIN*A A = A - PIV*STORE RETURN END SUBROUTINE BACK(A,Z,N,K,C) C IMPLICIT NONE INTEGER MXP PARAMETER (MXP=1000) INTEGER I,J,I1,M,L,N,K REAL A,Z,C,STORE C SUBROUTINE BACK CALCULATES THE SOLUTION OF THE SYSTEM OF C EQUATIONS A*C = Z WITH A A N X N UNIT UPPER TRIANGULAR MATRIX C OF BANDWIDTH K+1. C ATTENTION: THE FIRST DIMENSION SPECIFICATION OF MATRIX A MUST C BE THE SAME AS IN THE CALLING PROGRAM. DIMENSION A(MXP,K),Z(N),C(N) C(N) = Z(N) I = N - 1 IF (I.EQ.0) GO TO 30 DO 20 J = 2,N STORE = Z(I) I1 = K IF (J.LE.K) I1 = J - 1 M = I DO 10 L = 1,I1 M = M + 1 STORE = STORE - C(M)*A(I,L) 10 CONTINUE C(I) = STORE I = I - 1 20 CONTINUE 30 RETURN END SUBROUTINE NKNOT(X,M,T,N,FPINT,NRDATA,NRINT) C IMPLICIT NONE INTEGER MXP PARAMETER (MXP=1000) INTEGER K,J,M,N,NRINT,JBEGIN,JPOINT,NUMBER,NRDATA(MXP) INTEGER MAXPT,MAXBEG,NEXT,JJ,JK,IHALF,NRX REAL X(M),T(MXP),FPINT(MXP),FPMAX C SUBROUTINE NKNOT LOCATES AN ADDITIONAL KNOT FOR A SPLINE OF DEGREE C K AND ADJUSTS THE CORRESPONDING PARAMETERS,I.E. C T : THE POSITION OF THE KNOTS. C N : THE NUMBER OF KNOTS. C NRINT : THE NUMBER OF KNOTINTERVALS. C FPINT : THE SUM OF SQUARES OF RESIDUAL RIGHT HAND SIDES C FOR EACH KNOT INTERVAL. C NRDATA: THE NUMBER OF DATA POINTS INSIDE EACH KNOT INTERVAL. C THE ARRAYS T,FPINT AND NRDATA MUST HAVE THE SAME DIMENSION C SPECIFICATIONS AS IN THE CALLING PROGRAM. K = (N-NRINT-1)/2 C SEARCH FOR KNOT INTERVAL T(NUMBER+K) <= X <= T(NUMBER+K+1) WHERE C FPINT(NUMBER) IS MAXIMAL ON THE CONDITION THAT NRDATA(NUMBER) C NOT EQUALS ZERO. FPMAX = 0. JBEGIN = 1 DO 20 J = 1,NRINT JPOINT = NRDATA(J) IF (FPMAX.GE.FPINT(J) .OR. JPOINT.EQ.0) GO TO 10 FPMAX = FPINT(J) NUMBER = J MAXPT = JPOINT MAXBEG = JBEGIN 10 JBEGIN = JBEGIN + JPOINT + 1 20 CONTINUE C LET COINCIDE THE NEW KNOT T(NUMBER+K+1) WITH A DATA POINT X(NRX) C INSIDE THE OLD KNOT INTERVAL T(NUMBER+K) <= X <= T(NUMBER+K+1). IHALF = MAXPT/2 + 1 NRX = MAXBEG + IHALF NEXT = NUMBER + 1 IF (NEXT.GT.NRINT) GO TO 40 C ADJUSTS THE DIFFERENT PARAMETERS. DO 30 J = NEXT,NRINT JJ = NEXT + NRINT - J FPINT(JJ+1) = FPINT(JJ) NRDATA(JJ+1) = NRDATA(JJ) JK = JJ + K T(JK+1) = T(JK) 30 CONTINUE 40 NRDATA(NUMBER) = IHALF - 1 NRDATA(NEXT) = MAXPT - IHALF FPINT(NUMBER) = FPMAX*FLOAT(NRDATA(NUMBER))/FLOAT(MAXPT) FPINT(NEXT) = FPMAX*FLOAT(NRDATA(NEXT))/FLOAT(MAXPT) JK = NEXT + K T(JK) = X(NRX) N = N + 1 NRINT = NRINT + 1 RETURN END SUBROUTINE DISCO(T,N,K2,B) C IMPLICIT NONE INTEGER MXP PARAMETER (MXP=1000) INTEGER N,K1,K2,K,L,NK1,J,LP,LMK,IK,LJ,LK,JK,I REAL T(N),B(MXP,K2),H(12),PROD C SUBROUTINE DISCO CALCULATES THE DISCONTINUITY JUMPS OF THE KTH C DERIVATIVE OF THE B-SPLINES OF DEGREE K AT THE KNOTS T(K+2)..T(N-K-1) C THE FIRST DIMENSION SPECIFICATION OF THE MATRIX B MUST BE THE SAME AS C IN THE CALLING PROGRAM; H MUST HAVE A DIMENSION SPECIFICATION AT C LEAST 2*K+2. K1 = K2 - 1 K = K1 - 1 NK1 = N - K1 DO 40 L = K2,NK1 LMK = L - K1 DO 10 J = 1,K1 IK = J + K1 LJ = L + J LK = LJ - K2 H(J) = T(L) - T(LK) H(IK) = T(L) - T(LJ) 10 CONTINUE LP = LMK DO 30 J = 1,K2 JK = J + K PROD = 1. DO 20 I = J,JK PROD = PROD*H(I) 20 CONTINUE LK = LP + K1 B(LMK,J) = (T(LK)-T(LP))/PROD LP = LP + 1 30 CONTINUE 40 CONTINUE RETURN END REAL FUNCTION RATION(P1,F1,P2,F2,P3,F3) REAL P,P1,F1,P2,F2,P3,F3,H1,H2,H3 C GIVEN THREE POINTS (P1,F1),(P2,F2) AND (P3,F3), FUNCTION RATION C GIVES THE VALUE OF P SUCH THAT THE RATIONAL INTERPOLATING FUNCTION C IMPLICIT NONE C OF THE FORM R(P) = (U*P+V)/(P+W) EQUALS ZERO AT P. IF (P3.GT.0.) GO TO 10 C VALUE OF P IN CASE P3 = INFINITY. P = (P1* (F1-F3)*F2-P2* (F2-F3)*F1)/ ((F1-F2)*F3) GO TO 20 C VALUE OF P IN CASE P3 ^= INFINITY. 10 H1 = F1* (F2-F3) H2 = F2* (F3-F1) H3 = F3* (F1-F2) P = - (P1*P2*H3+P2*P3*H1+P3*P1*H2)/ (P1*H1+P2*H2+P3*H3) C ADJUST THE VALUE OF P1,F1,P3 AND F3 SUCH THAT F1 > 0 AND F3 < 0. 20 IF (F2.LT.0.) GO TO 30 P1 = P2 F1 = F2 GO TO 40 30 P3 = P2 F3 = F2 40 RATION = P RETURN END REAL FUNCTION DERIV(T,N,C,NK1,NU,ARG,L) C FUNCTION DERIV EVALUATES A SPLINE S(X) OF DEGREE K WHICH IS C GIVEN IN ITS NORMALIZED B-SPLINE REPRESENTATION OR CALCULATES C IMPLICIT NONE C DERIVATIVES OF ANY SPECIFIED ORDER NU. C C CALLING SEQUENCE C VALUE = DERIV(T,N,C,NK1,NU,ARG,L) C C INPUT PARAMETERS: C T : ARRAY,MINIMUM LENGTH N, WHICH CONTAINS THE POSITION C OF THE KNOTS OF S(X),I.E. THE POSITION OF THE INTERIOR C KNOTS T(K+2),...T(N-K-1) AS WELL AS THE POSITION OF THE C KNOTS T(1),...T(K+1) AND T(N-K),...T(N) WHICH ARE NEEDED C FOR THE B-SPLINE REPRESENTATION. C N : INTEGER VALUE GIVING THE TOTAL NUMBER OF KNOTS OF S(X). C C : ARRAY,LENGTH NK1, CONTAINING THE B-SPLINE COEFFICIENTS. C NK1 : INTEGER VALUE,GIVING THE DIMENSION OF S(X),I.E.NK1 = N-K-1. C NU : INTEGER VALUE WHICH GIVES THE ORDER OF THE DERIVATIVE. C ARG : REAL VALUE,GIVING THE VALUE OF THE ARGUMENT. C L : INTEGER VALUE WHICH SPECIFIES THE POSITION OF THE ARGUMENT C I.E. T(L) <= ARG < T(L+1) OR C L = NK1 IF ARG = T(NK1+1). C C OUTPUT PARAMETER: C VALUE: REAL VALUE,GIVING THE VALUE OF THE NUTH DERIVATIVE OF C S(X) AT X = ARG. C C OTHER SUBROUTINES REQUIRED: NONE. C C RESTRICTIONS: C 1) NU >= 0 C 2) T(K+1) <= ARG <= T(NK1+1) C THE DIMENSION SPECIFICATION OF THE ARRAY H MUST BE AT LEAST K+1. INTEGER N,NK1,NU,L,I,J,K1,IK,JJ,NU1,NU2,LI,LJ REAL T(N),C(NK1),H(6),ARG DERIV = 0. K1 = N - NK1 IF (NU.LT.0 .OR. NU.GE.K1) RETURN DO 10 I = 1,K1 IK = L + I - K1 H(I) = C(IK) 10 CONTINUE IF (NU.EQ.0) GO TO 40 NU1 = NU + 1 DO 30 J = 2,NU1 DO 20 JJ = J,K1 I = J + K1 - JJ LI = L + I - K1 LJ = L + I - J + 1 H(I) = (H(I)-H(I-1))/ (T(LJ)-T(LI)) 20 CONTINUE 30 CONTINUE IF (NU.EQ.K1-1) GO TO 70 40 NU2 = NU + 2 DO 60 J = NU2,K1 DO 50 JJ = J,K1 I = J + K1 - JJ LI = L + I - K1 LJ = L + I - J + 1 H(I) = ((ARG-T(LI))*H(I)+ (T(LJ)-ARG)*H(I-1))/ + (T(LJ)-T(LI)) 50 CONTINUE 60 CONTINUE 70 DERIV = H(K1) IF (NU.EQ.0) RETURN DO 80 I = 1,NU DERIV = DERIV*FLOAT(K1-I) 80 CONTINUE RETURN END SUBROUTINE IERSPL (IER,LEVEL) C C. IDENTIFICATION C C Author : P. BALLESTER ESO-Garching C C Versions : 1.00 27.06.90 C C. PURPOSE C C This subroutine displays error messages using MIDAS interface C STTPUT, depending on the error code returned by SMOSPL (IER) C and a level (LEVEL) C C. INPUTS : C C IER : Error code returned by the subroutine SMOSPL C C LEVEL : = 0 No messages C = 1 Only severe error messages C = 2 Severe error messages and diagnosis C = 3 Maximum of information C C --- IMPLICIT NONE C INTEGER IER,LEVEL,STAT CHARACTER*80 STRING C IF (IER.GT.0) THEN IF (LEVEL.GT.0) THEN STRING = 'Error while interpolating:' CALL STTPUT (STRING,STAT) ENDIF IF (LEVEL.GT.1) THEN IF (IER.EQ.1) THEN STRING = '--> Insufficient storage space' CALL STTPUT (STRING,STAT) ELSE IF (IER.EQ.2) THEN STRING = '--> Theoretically impossible result' CALL STTPUT (STRING,STAT) STRING = + ' The requested relative accuracy is probably too small' CALL STTPUT (STRING,STAT) ELSE IF (IER.EQ.3) THEN STRING = '--> Maximal number of iterations has been reached' CALL STTPUT (STRING,STAT) ELSE IF (IER.EQ.10) THEN STRING = '--> Invalid input data' CALL STTPUT (STRING,STAT) ENDIF ENDIF ENDIF IF (IER.LE.0.AND.LEVEL.GT.2) THEN STRING = 'Warning : Normal return from spline interpolation' CALL STTPUT (STRING,STAT) IF (IER.EQ.-1) THEN STRING = 'S(x) is an interpolating spline' CALL STTPUT (STRING,STAT) ELSE IF (IER.EQ.-2) THEN STRING = 'S(x) is a polynomial of degree k' CALL STTPUT (STRING,STAT) ENDIF ENDIF RETURN END SUBROUTINE SORSPL (TAB,AUX1,AUX2,NP,SENSE) C C. IDENTIFICATION C C Author : P. BALLESTER ESO-Garching C C Versions : 1.00 27.06.90 C 1.10 17.07.90 C C. PURPOSE C C Sort values by increasing values in a 1D real array (TAB). C Elements of the auxiliary arrays AUX1 and AUX2 are also moved. C This subroutine is intended to be used before the SMOSPL routine C to sort the data, if necessary. The syntax is : C CALL SORSPL (X,Y,W,NP,SENS) C C. INPUT/OUTPUT : C C TAB Real array Address of the reference array C AUX1 Real array Address of the first auxiliary array C AUX2 Real array Address of the second auxiliary array C NP Integer Number of significant values in TAB,AUX1,AUX2 C SENSE Integer >0 = sort by increasing values C <0 = sort by decreasing values C C. VERSIONS C C 1.00 Sort reference and auxiliary arrays by increasing values C 1.10 Sort by increasing or decreasing values. C C ------------------------------------------------------------------------- C IMPLICIT NONE C INTEGER FLAG,NP,PTR,SENSE REAL TAB(NP),AUX1(NP),AUX2(NP),BUFFER,SS IF (SENSE.GT.0) SS = 1.0 ! Sort by increasing values IF (SENSE.LT.0) SS = -1.0 ! Sort by decreasing values 100 FLAG = 0 DO 200 PTR = 1,NP-1 IF ((SS*TAB(PTR)).GT.(SS*TAB(PTR+1))) THEN FLAG = 1 BUFFER = TAB (PTR) TAB (PTR) = TAB (PTR+1) TAB (PTR+1) = BUFFER BUFFER = AUX1 (PTR) AUX1 (PTR) = AUX1 (PTR+1) AUX1 (PTR+1) = BUFFER BUFFER = AUX2 (PTR) AUX2 (PTR) = AUX2 (PTR+1) AUX2 (PTR+1) = BUFFER ENDIF 200 CONTINUE IF (FLAG.NE.0) GOTO 100 RETURN END SUBROUTINE ARGSPL (ARG,KNOTS,NBKNOT,DEGREE,PREPOS,POSARG,NK1) IMPLICIT NONE INTEGER NBKNOT,DEGREE,POSARG REAL ARG,KNOTS(NBKNOT) INTEGER NK1,PREPOS,STAT CHARACTER STRING*80 NK1 = NBKNOT - DEGREE - 1 IF (ARG.EQ.KNOTS(NK1+1)) THEN POSARG = NK1 GOTO 550 ENDIF DO 500 POSARG = PREPOS , NK1 IF (KNOTS(POSARG).LE.ARG.AND.ARG.LT.KNOTS(POSARG+1)) THEN cd WRITE (STRING,5500) POSARG,KNOTS(POSARG),ARG,KNOTS(POSARG+1) cd CALL STTPUT (STRING,STAT) PREPOS = POSARG GOTO 550 ENDIF 500 CONTINUE C C --- 5.1.2. Cannot find argument position C STRING = '*** Argument position not defined' CALL STTPUT (STRING,STAT) WRITE (STRING,5500) POSARG,KNOTS(POSARG),ARG,KNOTS(POSARG+1) 5500 FORMAT ('Pos : ',I5,' Relation : ',F6.1,' <',F6.1,' <',F6.1) CALL STTPUT (STRING,STAT) 550 RETURN END