home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
High Voltage Shareware
/
high1.zip
/
high1
/
DIR9
/
ACMALG01.ZIP
/
ACM600.FOR
< prev
next >
Wrap
Text File
|
1986-09-12
|
33KB
|
1,139 lines
C ALGORITHM 600, COLLECTED ALGORITHMS FROM ACM.
C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 2,
C JUN., 1983, P. 258-259.
SUBROUTINE QUINAT(N, X, Y, B, C, D, E, F)
C
INTEGER N
REAL X(N), Y(N), B(N), C(N), D(N), E(N), F(N)
C
C
C
C QUINAT COMPUTES THE COEFFICIENTS OF A QUINTIC NATURAL QUINTIC SPLI
C S(X) WITH KNOTS X(I) INTERPOLATING THERE TO GIVEN FUNCTION VALUES:
C S(X(I)) = Y(I) FOR I = 1,2, ..., N.
C IN EACH INTERVAL (X(I),X(I+1)) THE SPLINE FUNCTION S(XX) IS A
C POLYNOMIAL OF FIFTH DEGREE:
C S(XX) = ((((F(I)*P+E(I))*P+D(I))*P+C(I))*P+B(I))*P+Y(I) (*)
C = ((((-F(I)*Q+E(I+1))*Q-D(I+1))*Q+C(I+1))*Q-B(I+1))*Q+Y(I+1)
C WHERE P = XX - X(I) AND Q = X(I+1) - XX.
C (NOTE THE FIRST SUBSCRIPT IN THE SECOND EXPRESSION.)
C THE DIFFERENT POLYNOMIALS ARE PIECED TOGETHER SO THAT S(X) AND
C ITS DERIVATIVES UP TO S"" ARE CONTINUOUS.
C
C INPUT:
C
C N NUMBER OF DATA POINTS, (AT LEAST THREE, I.E. N > 2)
C X(1:N) THE STRICTLY INCREASING OR DECREASING SEQUENCE OF
C KNOTS. THE SPACING MUST BE SUCH THAT THE FIFTH POWER
C OF X(I+1) - X(I) CAN BE FORMED WITHOUT OVERFLOW OR
C UNDERFLOW OF EXPONENTS.
C Y(1:N) THE PRESCRIBED FUNCTION VALUES AT THE KNOTS.
C
C OUTPUT:
C
C B,C,D,E,F THE COMPUTED SPLINE COEFFICIENTS AS IN (*).
C (1:N) SPECIFICALLY
C B(I) = S'(X(I)), C(I) = S"(X(I))/2, D(I) = S"'(X(I))/6,
C E(I) = S""(X(I))/24, F(I) = S""'(X(I))/120.
C F(N) IS NEITHER USED NOR ALTERED. THE FIVE ARRAYS
C B,C,D,E,F MUST ALWAYS BE DISTINCT.
C
C OPTION:
C
C IT IS POSSIBLE TO SPECIFY VALUES FOR THE FIRST AND SECOND
C DERIVATIVES OF THE SPLINE FUNCTION AT ARBITRARILY MANY KNOTS.
C THIS IS DONE BY RELAXING THE REQUIREMENT THAT THE SEQUENCE OF
C KNOTS BE STRICTLY INCREASING OR DECREASING. SPECIFICALLY:
C
C IF X(J) = X(J+1) THEN S(X(J)) = Y(J) AND S'(X(J)) = Y(J+1),
C IF X(J) = X(J+1) = X(J+2) THEN IN ADDITION S"(X(J)) = Y(J+2).
C
C NOTE THAT S""(X) IS DISCONTINUOUS AT A DOUBLE KNOT AND, IN
C ADDITION, S"'(X) IS DISCONTINUOUS AT A TRIPLE KNOT. THE
C SUBROUTINE ASSIGNS Y(I) TO Y(I+1) IN THESE CASES AND ALSO TO
C Y(I+2) AT A TRIPLE KNOT. THE REPRESENTATION (*) REMAINS
C VALID IN EACH OPEN INTERVAL (X(I),X(I+1)). AT A DOUBLE KNOT,
C X(J) = X(J+1), THE OUTPUT COEFFICIENTS HAVE THE FOLLOWING VALUES:
C Y(J) = S(X(J)) = Y(J+1)
C B(J) = S'(X(J)) = B(J+1)
C C(J) = S"(X(J))/2 = C(J+1)
C D(J) = S"'(X(J))/6 = D(J+1)
C E(J) = S""(X(J)-0)/24 E(J+1) = S""(X(J)+0)/24
C F(J) = S""'(X(J)-0)/120 F(J+1) = S""'(X(J)+0)/120
C AT A TRIPLE KNOT, X(J) = X(J+1) = X(J+2), THE OUTPUT
C COEFFICIENTS HAVE THE FOLLOWING VALUES:
C Y(J) = S(X(J)) = Y(J+1) = Y(J+2)
C B(J) = S'(X(J)) = B(J+1) = B(J+2)
C C(J) = S"(X(J))/2 = C(J+1) = C(J+2)
C D(J) = S"'((X(J)-0)/6 D(J+1) = 0 D(J+2) = S"'(X(J)+0)/6
C E(J) = S""(X(J)-0)/24 E(J+1) = 0 E(J+2) = S""(X(J)+0)/24
C F(J) = S""'(X(J)-0)/120 F(J+1) = 0 F(J+2) = S""'(X(J)+0)/120
C
INTEGER I, M
REAL B1, P, PQ, PQQR, PR, P2, P3, Q, QR, Q2, Q3, R, R2, S, T, U, V
C
IF (N.LE.2) GO TO 190
C
C COEFFICIENTS OF A POSITIVE DEFINITE, PENTADIAGONAL MATRIX,
C STORED IN D,E,F FROM 2 TO N-2.
C
M = N - 2
Q = X(2) - X(1)
R = X(3) - X(2)
Q2 = Q*Q
R2 = R*R
QR = Q + R
D(1) = 0.
E(1) = 0.
D(2) = 0.
IF (Q.NE.0.) D(2) = 6.*Q*Q2/(QR*QR)
C
IF (M.LT.2) GO TO 40
DO 30 I=2,M
P = Q
Q = R
R = X(I+2) - X(I+1)
P2 = Q2
Q2 = R2
R2 = R*R
PQ = QR
QR = Q + R
IF (Q) 20, 10, 20
10 D(I+1) = 0.
E(I) = 0.
F(I-1) = 0.
GO TO 30
20 Q3 = Q2*Q
PR = P*R
PQQR = PQ*QR
D(I+1) = 6.*Q3/(QR*QR)
D(I) = D(I) + (Q+Q)*(15.*PR*PR+(P+R)*Q*(20.*PR+7.*Q2)+Q2*(8.*
* (P2+R2)+21.*PR+Q2+Q2))/(PQQR*PQQR)
D(I-1) = D(I-1) + 6.*Q3/(PQ*PQ)
E(I) = Q2*(P*QR+3.*PQ*(QR+R+R))/(PQQR*QR)
E(I-1) = E(I-1) + Q2*(R*PQ+3.*QR*(PQ+P+P))/(PQQR*PQ)
F(I-1) = Q3/PQQR
30 CONTINUE
C
40 IF (R.NE.0.) D(M) = D(M) + 6.*R*R2/(QR*QR)
C
C FIRST AND SECOND ORDER DIVIDED DIFFERENCES OF THE GIVEN FUNCTION
C VALUES, STORED IN B FROM 2 TO N AND IN C FROM 3 TO N
C RESPECTIVELY. CARE IS TAKEN OF DOUBLE AND TRIPLE KNOTS.
C
DO 60 I=2,N
IF (X(I).NE.X(I-1)) GO TO 50
B(I) = Y(I)
Y(I) = Y(I-1)
GO TO 60
50 B(I) = (Y(I)-Y(I-1))/(X(I)-X(I-1))
60 CONTINUE
DO 80 I=3,N
IF (X(I).NE.X(I-2)) GO TO 70
C(I) = B(I)*0.5
B(I) = B(I-1)
GO TO 80
70 C(I) = (B(I)-B(I-1))/(X(I)-X(I-2))
80 CONTINUE
C
C SOLVE THE LINEAR SYSTEM WITH C(I+2) - C(I+1) AS RIGHT-HAND SIDE.
C
IF (M.LT.2) GO TO 100
P = 0.
C(1) = 0.
E(M) = 0.
F(1) = 0.
F(M-1) = 0.
F(M) = 0.
C(2) = C(4) - C(3)
D(2) = 1./D(2)
C
IF (M.LT.3) GO TO 100
DO 90 I=3,M
Q = D(I-1)*E(I-1)
D(I) = 1./(D(I)-P*F(I-2)-Q*E(I-1))
E(I) = E(I) - Q*F(I-1)
C(I) = C(I+2) - C(I+1) - P*C(I-2) - Q*C(I-1)
P = D(I-1)*F(I-1)
90 CONTINUE
C
100 I = N - 1
C(N-1) = 0.
C(N) = 0.
IF (N.LT.4) GO TO 120
DO 110 M=4,N
C I = N-2, ..., 2
I = I - 1
C(I) = (C(I)-E(I)*C(I+1)-F(I)*C(I+2))*D(I)
110 CONTINUE
C
C INTEGRATE THE THIRD DERIVATIVE OF S(X).
C
120 M = N - 1
Q = X(2) - X(1)
R = X(3) - X(2)
B1 = B(2)
Q3 = Q*Q*Q
QR = Q + R
IF (QR) 140, 130, 140
130 V = 0.
T = 0.
GO TO 150
140 V = C(2)/QR
T = V
150 F(1) = 0.
IF (Q.NE.0.) F(1) = V/Q
DO 180 I=2,M
P = Q
Q = R
R = 0.
IF (I.NE.M) R = X(I+2) - X(I+1)
P3 = Q3
Q3 = Q*Q*Q
PQ = QR
QR = Q + R
S = T
T = 0.
IF (QR.NE.0.) T = (C(I+1)-C(I))/QR
U = V
V = T - S
IF (PQ) 170, 160, 170
160 C(I) = C(I-1)
D(I) = 0.
E(I) = 0.
F(I) = 0.
GO TO 180
170 F(I) = F(I-1)
IF (Q.NE.0.) F(I) = V/Q
E(I) = 5.*S
D(I) = 10.*(C(I)-Q*S)
C(I) = D(I)*(P-Q) + (B(I+1)-B(I)+(U-E(I))*P3-(V+E(I))*Q3)/PQ
B(I) = (P*(B(I+1)-V*Q3)+Q*(B(I)-U*P3))/PQ -
* P*Q*(D(I)+E(I)*(Q-P))
180 CONTINUE
C
C END POINTS X(1) AND X(N).
C
P = X(2) - X(1)
S = F(1)*P*P*P
E(1) = 0.
D(1) = 0.
C(1) = C(2) - 10.*S
B(1) = B1 - (C(1)+S)*P
C
Q = X(N) - X(N-1)
T = F(N-1)*Q*Q*Q
E(N) = 0.
D(N) = 0.
C(N) = C(N-1) + 10.*T
B(N) = B(N) + (C(N)-T)*Q
190 RETURN
END
C
C
SUBROUTINE QUINEQ(N, Y, B, C, D, E, F)
C
INTEGER N
REAL Y(N), B(N), C(N), D(N), E(N), F(N)
C
C
C
C QUINEQ COMPUTES THE COEFFICIENTS OF QUINTIC NATURAL QUINTIC SPLINE
C S(X) WITH EQUIDISTANT KNOTS X(I) INTERPOLATING THERE TO GIVEN
C FUNCTION VALUES:
C S(X(I)) = Y(I) FOR I = 1,2, ..., N.
C IN EACH INTERVAL (X(I),X(I+1)) THE SPLINE FUNCTION S(XX) IS
C A POLYNOMIAL OF FIFTH DEGREE:
C S(XX)=((((F(I)*P+E(I))*P+D(I))*P+C(I))*P+B(I))*P+Y(I) (*)
C =((((-F(I)*Q+E(I+1))*Q-D(I+1))*Q+C(I+1))*Q-B(I+1))*Q+Y(I+1)
C WHERE P = (XX - X(I))/X(I+1) - X(I))
C AND Q = (X(I+1) - XX)/(X(I+1) - X(I)).
C (NOTE THE FIRST SUBSCRIPT IN THE SECOND EXPRESSION.)
C THE DIFFERENT POLYNOMIALS ARE PIECED TOGETHER SO THAT S(X) AND
C ITS DERIVATIVES UP TO S"" ARE CONTINUOUS.
C
C INPUT:
C
C N NUMBER OF DATA POINTS, (AT LEAST THREE, I.E. N > 2)
C Y(1:N) THE PRESCRIBED FUNCTION VALUES AT THE KNOTS
C
C OUTPUT:
C
C B,C,D,E,F THE COMPUTED SPLINE COEFFICIENTS AS IN (*).
C (1:N) IF X(I+1) - X(I) = 1., THEN SPECIFICALLY:
C B(I) = S'X(I)), C(I) = S"(X(I))/2, D(I) = S"'(X(I))/6,
C E(I) = S""(X(I))/24, F(I) = S""'(X(I)+0)/120.
C F(N) IS NEITHER USED NOR ALTERED. THE ARRAYS
C Y,B,C,D MUST ALWAYS BE DISTINCT. IF E AND F ARE
C NOT WANTED, THE CALL QUINEQ(N,Y,B,C,D,D,D) MAY
C BE USED TO SAVE STORAGE LOCATIONS.
C
INTEGER I, M
REAL P, Q, R, S, T, U, V
C
IF (N.LE.2) GO TO 50
C
M = N - 3
P = 0.
Q = 0.
R = 0.
S = 0.
T = 0.
D(M+1) = 0.
D(M+2) = 0.
IF (M.LE.0) GO TO 30
DO 10 I=1,M
U = P*R
B(I) = 1./(66.-U*R-Q)
R = 26. - U
C(I) = R
D(I) = Y(I+3) - 3.*(Y(I+2)-Y(I+1)) - Y(I) - U*S - Q*T
Q = P
P = B(I)
T = S
S = D(I)
10 CONTINUE
C
I = N - 2
DO 20 M=4,N
C I = N-3, ..., 1
I = I - 1
D(I) = (D(I)-C(I)*D(I+1)-D(I+2))*B(I)
20 CONTINUE
C
30 M = N - 1
Q = 0.
R = D(1)
T = R
V = R
DO 40 I=2,M
P = Q
Q = R
R = D(I)
S = T
T = P - Q - Q + R
F(I) = T
U = 5.*(-P+Q)
E(I) = U
D(I) = 10.*(P+Q)
C(I) = 0.5*(Y(I+1)+Y(I-1)+S-T) - Y(I) - U
B(I) = 0.5*(Y(I+1)-Y(I-1)-S-T) - D(I)
40 CONTINUE
C
F(1) = V
E(1) = 0.
D(1) = 0.
C(1) = C(2) - 10.*V
B(1) = Y(2) - Y(1) - C(1) - V
E(N) = 0.
D(N) = 0.
C(N) = C(N-1) + 10.*T
B(N) = Y(N) - Y(N-1) + C(N) - T
50 RETURN
END
C
C
C
C
SUBROUTINE QUINDF(N, X, Y, B, C, D, E, F)
C
INTEGER N
REAL X(N), Y(N), B(N), C(N), D(N), E(N), F(N)
C
C
C
C QUINDF COMPUTES THE COEFFICIENTS OF A QUINTIC NATURAL QUINTIC
C SPLINE S(X) WITH KNOTS X(I) FOR WHICH THE FUNCTION VALUES Y(I) AND
C THE FIRST DERIVATIVES B(I) ARE SPECIFIED AT X(I), I = 1,2, ...,N.
C IN EACH INTERVAL (X(I),X(I+1)) THE SPLINE FUNCTION S(XX) IS
C A POLYNOMIAL OF FIFTH DEGREE:
C S(XX)=((((F(I)*P+E(I))*P+D(I))*P+C(I))*P+B(I))*P+Y(I) (*)
C WHERE P = XX - X(I).
C
C INPUT:
C
C N NUMBER OF DATA POINTS, (AT LEAST TWO, I.E. N > 1)
C X(1:N) THE STRICTLY INCREASING OR DECREASING SEQUENCE OF
C KNOTS. THE SPACING MUST BE SUCH THAT THE FIFTH
C POWER OF X(I+1) - X(I) CAN BE FORMED WITHOUT
C OVERFLOW OR UNDERFLOW OF EXPONENTS
C Y(1:N) THE PRESCRIBED FUNCTION VALUES AT THE KNOTS
C B(1:N) THE PRESCRIBED DERIVATIVE VALUES AT THE KNOTS
C
C OUTPUT:
C
C C,D,E,F THE COMPUTED SPLINE COEFFICIENTS AS IN (*).
C (1:N) E(N) AND F(N) ARG NEITHER USED NOR ALTERED.
C THE ARRAYS C,D,E,F MUST ALWAYS BE DISTINCT.
C
INTEGER I, M, N1
REAL CC, G, H, HH, H2, P, PP, Q, QQ, R, RR
C
IF (N.LE.1) GO TO 40
N1 = N - 1
CC = 0.
HH = 0.
PP = 0.
QQ = 0.
RR = 0.
G = 0.
DO 10 I=1,N1
H = 1./(X(I+1)-X(I))
H2 = H*H
D(I) = 3.*(HH+H) - G*HH
P = (Y(I+1)-Y(I))*H2*H
Q = (B(I+1)+B(I))*H2
R = (B(I+1)-B(I))*H2
CC = 10.*(P-PP) - 5.*(Q-QQ) + R + RR + G*CC
C(I) = CC
G = H/D(I)
HH = H
PP = P
QQ = Q
RR = R
10 CONTINUE
C
C(N) = (-10.*PP+5.*QQ+RR+G*CC)/(3.*HH-G*HH)
I = N
DO 20 M=1,N1
C I = N-1, ..., 1
I = I - 1
D(I+1) = 1./(X(I+1)-X(I))
C(I) = (C(I)+C(I+1)*D(I+1))/D(I)
20 CONTINUE
C
DO 30 I=1,N1
H = D(I+1)
P = (((Y(I+1)-Y(I))*H-B(I))*H-C(I))*H
Q = ((B(I+1)-B(I))*H-C(I)-C(I))*H
R = (C(I+1)-C(I))*H
G = Q - 3.*P
RR = R - 3.*(P+G)
QQ = -RR - RR + G
F(I) = RR*H*H
E(I) = QQ*H
D(I) = -RR - QQ + P
30 CONTINUE
C
D(N) = 0.
E(N) = 0.
F(N) = 0.
40 RETURN
END
C DRIVER PROGRAM FOR TEST OF QUINAT, QUINEQ AND QUINDF
C FOLLOWS.
C
INTEGER N,NM1,M,MM,MM1,I,K,J,JJ
REAL Z
REAL X(200),Y(200),B(200),BB(200),CC(200),DD(200),EE(200),
* FF(200),A(200,6),C(6),DIFF(5),COM(5)
C
C N NUMBER OF DATA POINTS.
C M 2*M-1 IS ORDER OF SPLINE.
C M = 3 ALWAYS FOR QUINTIC SPLINE.
C NN,NM1,MM,
C MM1,I,K,
C J,JJ TEMPORARY INTEGER VARIABLES.
C Z,P TEMPORARY REAL VARIABLES.
C X(1:N) THE SEQUENCE OF KNOTS.
C Y(1:N) THE PRESCRIBED FUNCTION VALUES AT THE KNOTS.
C B(1:N) THE PRESCRIBED DERIVATIVE VALUES AT THE KNOTS.
C BB,CC,DD,
C EE,FF(1:N) THE COMPUTED SPLINE COEFFICIENTS
C A(1:N,1:6) TWO DIMENSIONAL ARRAY WHOSE COLUMNS ARE
C Y, BB, CC, DD, EE, FF.
C DIFF(1:5) MAXIMUM VALUES OF DIFFERENCES OF VALUES AND
C DERIVATIVES TO RIGHT AND LEFT OF KNOTS.
C COM(1:5) MAXIMUM VALUES OF COEFFICIENTS.
C
C
C TEST OF QUINAT WITH NONEQUIDISTANT KNOTS AND
C EQUIDISTANT KNOTS FOLLOWS.
C
NOUT=6
WRITE(NOUT,10)
10 FORMAT(50H1 TEST OF QUINAT WITH NONEQUIDISTANT KNOTS)
N = 5
X(1) = -3.0
X(2) = -1.0
X(3) = 0.0
X(4) = 3.0
X(5) = 4.0
Y(1) = 7.0
Y(2) = 11.0
Y(3) = 26.0
Y(4) = 56.0
Y(5) = 29.0
M = 3
MM = 2*M
MM1 = MM - 1
WRITE(NOUT,15) N,M
15 FORMAT(5H-N = ,I3,7H M =,I2)
CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
DO 20 I = 1,N
A(I,1) = Y(I)
A(I,2) = BB(I)
A(I,3) = CC(I)
A(I,4) = DD(I)
A(I,5) = EE(I)
A(I,6) = FF(I)
20 CONTINUE
DO 30 I = 1,MM1
DIFF(I) = 0.0
COM(I) = 0.0
30 CONTINUE
DO 70 K = 1,N
DO 35 I = 1,MM
C(I) = A(K,I)
35 CONTINUE
WRITE(NOUT,40) K
40 FORMAT(40H ---------------------------------------,I3,
* 45H --------------------------------------------)
WRITE(NOUT,45) X(K)
45 FORMAT(F12.8)
IF (K .EQ. N) WRITE(NOUT,55) C(1)
IF (K .EQ. N) GO TO 75
WRITE(NOUT,55) (C(I),I=1,MM)
DO 50 I = 1,MM1
IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
50 CONTINUE
55 FORMAT(6F16.8)
Z = X(K+1) - X(K)
DO 60 I = 2,MM
DO 60 JJ = I,MM
J = MM + I - JJ
C(J-1) = C(J)*Z + C(J-1)
60 CONTINUE
WRITE(NOUT,55) (C(I),I=1,MM)
DO 65 I = 1,MM1
IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 65
Z = ABS(C(I) - A(K+1,I))
IF (Z .GT. DIFF(I)) DIFF(I) = Z
65 CONTINUE
70 CONTINUE
75 WRITE(NOUT,80)
80 FORMAT(41H MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
WRITE(NOUT,85) (DIFF(I),I=1,MM1)
85 FORMAT(5E18.9)
WRITE(NOUT,90)
90 FORMAT(42H MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
WRITE(NOUT,95) (COM(I),I=1,MM1)
95 FORMAT(5F16.8)
M = 3
DO 200 N = 10,100,10
MM = 2*M
MM1 = MM - 1
NM1 = N -1
DO 100 I = 1,NM1,2
X(I) = I
X(I+1) = I + 1
Y(I) = 1.
Y(I+1) = 0.
100 CONTINUE
IF (MOD(N,2) .EQ. 0) GOTO 105
X(N) = N
Y(N) = 1.
105 CONTINUE
WRITE(NOUT,110) N,M
110 FORMAT(5H-N = ,I3,7H M =,I2)
CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
DO 115 I = 1,N
A(I,1) = Y(I)
A(I,2) = BB(I)
A(I,3) = CC(I)
A(I,4) = DD(I)
A(I,5) = EE(I)
A(I,6) = FF(I)
115 CONTINUE
DO 120 I = 1, MM1
DIFF(I) = 0.0
COM(I) = 0.0
120 CONTINUE
DO 165 K = 1,N
DO 125 I = 1,MM
C(I) = A(K,I)
125 CONTINUE
IF (N .GT. 10) GOTO 140
WRITE(NOUT,130) K
130 FORMAT(40H ---------------------------------------,I3,
* 45H --------------------------------------------)
WRITE(NOUT,135) X(K)
135 FORMAT(F12.8)
IF (K .EQ. N) WRITE(NOUT,150) C(1)
140 CONTINUE
IF (K .EQ. N) GO TO 170
IF (N .LE. 10) WRITE(NOUT,150) (C(I), I=1,MM)
DO 145 I = 1,MM1
IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
145 CONTINUE
150 FORMAT(6F16.8)
Z = X(K+1) - X(K)
DO 155 I = 2,MM
DO 155 JJ = I,MM
J = MM + I - JJ
C(J-1) = C(J)*Z + C(J-1)
155 CONTINUE
IF (N .LE. 10) WRITE(NOUT,150) (C(I), I=1,MM)
DO 160 I = 1,MM1
IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 160
Z = ABS(C(I) - A(K+1,I))
IF (Z .GT. DIFF(I)) DIFF(I) = Z
160 CONTINUE
165 CONTINUE
170 WRITE(NOUT,175)
175 FORMAT(41H MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
WRITE(NOUT,180) (DIFF(I),I=1,MM1)
180 FORMAT(5E18.9)
WRITE(NOUT,185)
185 FORMAT(42H MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
WRITE(NOUT,190) (COM(I),I=1,MM1)
190 FORMAT(5E18.9)
200 CONTINUE
C
C
C TEST OF QUINEQ FOLLOWS.
C
WRITE(NOUT,210)
210 FORMAT(18H1 TEST OF QUINEQ)
M = 3
DO 400 N = 10,100,10
MM = 2*M
MM1 = MM - 1
NM1 = N -1
DO 300 I = 1,NM1,2
X(I) = I
X(I+1) = I + 1
Y(I) = 1.
Y(I+1) = 0.
300 CONTINUE
IF (MOD(N,2) .EQ. 0) GOTO 305
X(N) = FLOAT(N)
Y(N) = 1.
305 CONTINUE
WRITE(NOUT,310) N,M
310 FORMAT(5H-N = ,I3,7H M =,I2)
CALL QUINEQ(N,Y,BB,CC,DD,EE,FF)
DO 315 I = 1,N
A(I,1) = Y(I)
A(I,2) = BB(I)
A(I,3) = CC(I)
A(I,4) = DD(I)
A(I,5) = EE(I)
A(I,6) = FF(I)
315 CONTINUE
DO 320 I = 1, MM1
DIFF(I) = 0.0
COM(I) = 0.0
320 CONTINUE
DO 365 K = 1,N
DO 325 I = 1,MM
C(I) = A(K,I)
325 CONTINUE
IF (N .GT. 10) GOTO 340
WRITE(NOUT,330) K
330 FORMAT(40H ---------------------------------------,I3,
* 45H --------------------------------------------)
WRITE(NOUT,335) X(K)
335 FORMAT(F12.8)
IF (K .EQ. N) WRITE(NOUT,350) C(1)
340 CONTINUE
IF (K .EQ. N) GO TO 370
IF (N .LE. 10) WRITE(NOUT,350) (C(I), I=1,MM)
DO 345 I = 1,MM1
IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
345 CONTINUE
350 FORMAT(6F16.8)
Z = 1.
DO 355 I = 2,MM
DO 355 JJ = I,MM
J = MM + I - JJ
C(J-1) = C(J)*Z + C(J-1)
355 CONTINUE
IF (N .LE. 10) WRITE(NOUT,350) (C(I), I=1,MM)
DO 360 I = 1,MM1
IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 360
Z = ABS(C(I) - A(K+1,I))
IF (Z .GT. DIFF(I)) DIFF(I) = Z
360 CONTINUE
365 CONTINUE
370 WRITE(NOUT,375)
375 FORMAT(41H MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
WRITE(NOUT,380) (DIFF(I),I=1,MM1)
380 FORMAT(5E18.9)
WRITE(NOUT,385)
385 FORMAT(42H MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
WRITE(NOUT,390) (COM(I),I=1,MM1)
390 FORMAT(5E18.9)
400 CONTINUE
C
C
C TEST OF QUINDF WITH NONEQUIDISTANT KNOTS FOLLOWS.
C
WRITE(NOUT,410)
410 FORMAT(50H1 TEST OF QUINDF WITH NONEQUIDISTANT KNOTS)
N = 5
X(1) = -3.0
X(2) = -1.0
X(3) = 0.0
X(4) = 3.0
X(5) = 4.0
Y(1) = 7.0
Y(2) = 11.0
Y(3) = 26.0
Y(4) = 56.0
Y(5) = 29.0
B(1) = 2.0
B(2) = 15.0
B(3) = 10.0
B(4) =-27.0
B(5) =-30.0
M = 3
MM = 2*M
MM1 = MM - 1
WRITE(NOUT,415) N,M
415 FORMAT(5H-N = ,I3,7H M =,I2)
CALL QUINDF(N,X,Y,B,CC,DD,EE,FF)
DO 420 I = 1,N
A(I,1) = Y(I)
A(I,2) = B(I)
A(I,3) = CC(I)
A(I,4) = DD(I)
A(I,5) = EE(I)
A(I,6) = FF(I)
420 CONTINUE
DO 430 I = 1,MM1
DIFF(I) = 0.0
COM(I) = 0.0
430 CONTINUE
DO 470 K = 1,N
DO 435 I = 1,MM
C(I) = A(K,I)
435 CONTINUE
WRITE(NOUT,440) K
440 FORMAT(40H ---------------------------------------,I3,
* 45H --------------------------------------------)
WRITE(NOUT,445) X(K)
445 FORMAT(F12.8)
IF (K .EQ. N) WRITE(NOUT,455) C(1)
IF (K .EQ. N) GO TO 475
WRITE(NOUT,455) (C(I),I=1,MM)
DO 450 I = 1,MM1
IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
450 CONTINUE
455 FORMAT(6F16.8)
Z = X(K+1) - X(K)
DO 460 I = 2,MM
DO 460 JJ = I,MM
J = MM + I - JJ
C(J-1) = C(J)*Z + C(J-1)
460 CONTINUE
WRITE(NOUT,455) (C(I),I=1,MM)
DO 465 I = 1,MM1
IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 465
Z = ABS(C(I) - A(K+1,I))
IF (Z .GT. DIFF(I)) DIFF(I) = Z
465 CONTINUE
470 CONTINUE
475 WRITE(NOUT,480)
480 FORMAT(41H MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
WRITE(NOUT,485) (DIFF(I),I=1,MM1)
485 FORMAT(5E18.9)
WRITE(NOUT,490)
490 FORMAT(42H MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
WRITE(NOUT,495) (COM(I),I=1,MM1)
495 FORMAT(5F16.8)
M = 3
DO 600 N = 10,100,10
MM = 2*M
MM1 = MM - 1
P = 0.0
DO 500 I = 1,N
P = P + ABS(SIN((FLOAT(I)))) + 0.001*FLOAT(I)
X(I) = P
Y(I) = COS((FLOAT(I))) - 0.5
B(I) = COS((FLOAT(2*I))) - 0.5
500 CONTINUE
WRITE(NOUT,510) N,M
510 FORMAT(5H-N = ,I3,7H M =,I2)
CALL QUINDF(N,X,Y,B,CC,DD,EE,FF)
DO 515 I = 1,N
A(I,1) = Y(I)
A(I,2) = B(I)
A(I,3) = CC(I)
A(I,4) = DD(I)
A(I,5) = EE(I)
A(I,6) = FF(I)
515 CONTINUE
DO 520 I = 1, MM1
DIFF(I) = 0.0
COM(I) = 0.0
520 CONTINUE
DO 565 K = 1,N
DO 525 I = 1,MM
C(I) = A(K,I)
525 CONTINUE
IF (N .GT. 10) GOTO 540
WRITE(NOUT,530) K
530 FORMAT(40H ---------------------------------------,I3,
* 45H --------------------------------------------)
WRITE(NOUT,535) X(K)
535 FORMAT(F12.8)
IF (K .EQ. N) WRITE(NOUT,550) C(1)
540 CONTINUE
IF (K .EQ. N) GO TO 570
IF (N .LE. 10) WRITE(NOUT,550) (C(I), I=1,MM)
DO 545 I = 1,MM1
IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
545 CONTINUE
550 FORMAT(6F16.8)
Z = X(K+1) - X(K)
DO 555 I = 2,MM
DO 555 JJ = I,MM
J = MM + I - JJ
C(J-1) = C(J)*Z + C(J-1)
555 CONTINUE
IF (N .LE. 10) WRITE(NOUT,550) (C(I), I=1,MM)
DO 560 I = 1,MM1
IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 560
Z = ABS(C(I) - A(K+1,I))
IF (Z .GT. DIFF(I)) DIFF(I) = Z
560 CONTINUE
565 CONTINUE
570 WRITE(NOUT,575)
575 FORMAT(41H MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
WRITE(NOUT,580) (DIFF(I),I=1,MM1)
580 FORMAT(5E18.9)
WRITE(NOUT,585)
585 FORMAT(42H MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
WRITE(NOUT,590) (COM(I),I=1,MM1)
590 FORMAT(5E18.9)
600 CONTINUE
C
C
C TEST OF QUINAT WITH NONEQUIDISTANT DOUBLE KNOTS FOLLOWS.
C
WRITE(NOUT,610)
610 FORMAT(50H1 TEST OF QUINAT WITH NONEQUIDISTANT DOUBLE KNOTS)
N = 5
X(1) = -3.
X(2) = -3.
X(3) = -1.
X(4) = -1.
X(5) = 0.
X(6) = 0.
X(7) = 3.
X(8) = 3.
X(9) = 4.
X(10) = 4.
Y(1) = 7.
Y(2) = 2.
Y(3) = 11.
Y(4) = 15.
Y(5) = 26.
Y(6) = 10.
Y(7) = 56.
Y(8) = -27.
Y(9) = 29.
Y(10) = -30.
M = 3
NN = 2*N
MM = 2*M
MM1 = MM - 1
WRITE(NOUT,615) N,M
615 FORMAT(5H-N = ,I3,7H M =,I2)
CALL QUINAT(NN,X,Y,BB,CC,DD,EE,FF)
DO 620 I = 1,NN
A(I,1) = Y(I)
A(I,2) = BB(I)
A(I,3) = CC(I)
A(I,4) = DD(I)
A(I,5) = EE(I)
A(I,6) = FF(I)
620 CONTINUE
DO 630 I = 1,MM1
DIFF(I) = 0.0
COM(I) = 0.0
630 CONTINUE
DO 670 K = 1,NN
DO 635 I = 1,MM
C(I) = A(K,I)
635 CONTINUE
WRITE(NOUT,640) K
640 FORMAT(40H ---------------------------------------,I3,
* 45H --------------------------------------------)
WRITE(NOUT,645) X(K)
645 FORMAT(F12.8)
IF (K .EQ. NN) WRITE(NOUT,655) C(1)
IF (K .EQ. NN) GO TO 675
WRITE(NOUT,655) (C(I),I=1,MM)
DO 650 I = 1,MM1
IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
650 CONTINUE
655 FORMAT(6F16.8)
Z = X(K+1) - X(K)
DO 660 I = 2,MM
DO 660 JJ = I,MM
J = MM + I - JJ
C(J-1) = C(J)*Z + C(J-1)
660 CONTINUE
WRITE(NOUT,655) (C(I),I=1,MM)
DO 665 I = 1,MM1
IF ((K .GE. NN-1) .AND. (I .NE. 1)) GO TO 665
Z = ABS(C(I) - A(K+1,I))
IF (Z .GT. DIFF(I)) DIFF(I) = Z
665 CONTINUE
670 CONTINUE
675 WRITE(NOUT,680)
680 FORMAT(41H MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
WRITE(NOUT,685) (DIFF(I),I=1,MM1)
685 FORMAT(5E18.9)
WRITE(NOUT,690)
690 FORMAT(42H MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
WRITE(NOUT,695) (COM(I),I=1,MM1)
695 FORMAT(5F16.8)
M = 3
DO 800 N = 10,100,10
NN = 2*N
MM = 2*M
MM1 = MM - 1
P = 0.0
DO 700 I = 1,N
P = P + ABS(SIN(FLOAT(I)))
X(2*I-1) = P
X(2*I) = P
Y(2*I-1) = COS(FLOAT(I)) - 0.5
Y(2*I) = COS(FLOAT(2*I)) - 0.5
700 CONTINUE
WRITE(NOUT,710) N,M
710 FORMAT(5H-N = ,I3,7H M =,I2)
CALL QUINAT(NN,X,Y,BB,CC,DD,EE,FF)
DO 715 I = 1,NN
A(I,1) = Y(I)
A(I,2) = BB(I)
A(I,3) = CC(I)
A(I,4) = DD(I)
A(I,5) = EE(I)
A(I,6) = FF(I)
715 CONTINUE
DO 720 I = 1, MM1
DIFF(I) = 0.0
COM(I) = 0.0
720 CONTINUE
DO 765 K = 1,NN
DO 725 I = 1,MM
C(I) = A(K,I)
725 CONTINUE
IF (N .GT. 10) GOTO 740
WRITE(NOUT,730) K
730 FORMAT(40H ---------------------------------------,I3,
* 45H --------------------------------------------)
WRITE(NOUT,735) X(K)
735 FORMAT(F12.8)
IF (K .EQ. NN) WRITE(NOUT,750) C(1)
740 CONTINUE
IF (K .EQ. NN) GO TO 770
IF (N .LE. 10) WRITE(NOUT,750) (C(I), I=1,MM)
DO 745 I = 1,MM1
IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
745 CONTINUE
750 FORMAT(6F16.8)
Z = X(K+1) - X(K)
DO 755 I = 2,MM
DO 755 JJ = I,MM
J = MM + I - JJ
C(J-1) = C(J)*Z + C(J-1)
755 CONTINUE
IF (N .LE. 10) WRITE(NOUT,750) (C(I), I=1,MM)
DO 760 I = 1,MM1
IF ((K .GE. NN-1) .AND. (I .NE. 1)) GO TO 760
Z = ABS(C(I) - A(K+1,I))
IF (Z .GT. DIFF(I)) DIFF(I) = Z
760 CONTINUE
765 CONTINUE
770 WRITE(NOUT,775)
775 FORMAT(41H MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
WRITE(NOUT,780) (DIFF(I),I=1,MM1)
780 FORMAT(5E18.9)
WRITE(NOUT,785)
785 FORMAT(42H MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
IF (ABS(C(1)) .GT. COM(1)) COM(1) = ABS(C(1))
WRITE(NOUT,790) (COM(I),I=1,MM1)
790 FORMAT(5E18.9)
800 CONTINUE
C
C
C TEST OF QUINAT WITH NONEQUIDISTANT KNOTS, ONE DOUBLE KNOT,
C ONE TRIPLE KNOT, FOLLOWS.
C
WRITE(NOUT,805)
WRITE(NOUT,810)
805 FORMAT(51H1 TEST OF QUINAT WITH NONEQUIDISTANT KNOTS,)
810 FORMAT(40H ONE DOUBLE, ONE TRIPLE KNOT)
N = 8
X(1) = -3.
X(2) = -1.
X(3) = -1.
X(4) = 0.
X(5) = 3.
X(6) = 3.
X(7) = 3.
X(8) = 4.
Y(1) = 7.
Y(2) = 11.
Y(3) = 15.
Y(4) = 26.
Y(5) = 56.
Y(6) = -30.
Y(7) = -7.
Y(8) = 29.
M = 3
MM = 2*M
MM1 = MM - 1
WRITE(NOUT,815) N,M
815 FORMAT(5H-N = ,I3,7H M =,I2)
CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
DO 820 I = 1,N
A(I,1) = Y(I)
A(I,2) = BB(I)
A(I,3) = CC(I)
A(I,4) = DD(I)
A(I,5) = EE(I)
A(I,6) = FF(I)
820 CONTINUE
DO 830 I = 1,MM1
DIFF(I) = 0.0
COM(I) = 0.0
830 CONTINUE
DO 870 K = 1,N
DO 835 I = 1,MM
C(I) = A(K,I)
835 CONTINUE
WRITE(NOUT,840) K
840 FORMAT(40H ---------------------------------------,I3,
* 45H --------------------------------------------)
WRITE(NOUT,845) X(K)
845 FORMAT(F12.8)
IF (K .EQ. N) WRITE(NOUT,855) C(1)
IF (K .EQ. N) GO TO 875
WRITE(NOUT,855) (C(I),I=1,MM)
DO 850 I = 1,MM1
IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
850 CONTINUE
855 FORMAT(6F16.8)
Z = X(K+1) - X(K)
DO 860 I = 2,MM
DO 860 JJ = I,MM
J = MM + I - JJ
C(J-1) = C(J)*Z + C(J-1)
860 CONTINUE
WRITE(NOUT,855) (C(I),I=1,MM)
DO 865 I = 1,MM1
IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 865
Z = ABS(C(I) - A(K+1,I))
IF (Z .GT. DIFF(I)) DIFF(I) = Z
865 CONTINUE
870 CONTINUE
875 WRITE(NOUT,880)
880 FORMAT(41H MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
WRITE(NOUT,885) (DIFF(I),I=1,MM1)
885 FORMAT(5E18.9)
WRITE(NOUT,890)
890 FORMAT(42H MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
WRITE(NOUT,895) (COM(I),I=1,MM1)
895 FORMAT(5F16.8)
C
C
C TEST OF QUINAT WITH NONEQUIDISTANT KNOTS, TWO DOUBLE KNOTS,
C ONE TRIPLE KNOT,FOLLOWS.
C
WRITE(NOUT,905)
WRITE(NOUT,910)
905 FORMAT(51H1 TEST OF QUINAT WITH NONEQUIDISTANT KNOTS,)
910 FORMAT(40H TWO DOUBLE, ONE TRIPLE KNOT)
N = 10
X(1) = 0.
X(2) = 2.
X(3) = 2.
X(4) = 3.
X(5) = 3.
X(6) = 3.
X(7) = 5.
X(8) = 8.
X(9) = 9.
X(10)= 9.
Y(1) = 163.
Y(2) = 237.
Y(3) = -127.
Y(4) = 119.
Y(5) = -65.
Y(6) = 192.
Y(7) = 293.
Y(8) = 326.
Y(9) = 0.
Y(10)= -414.0
M = 3
MM = 2*M
MM1 = MM - 1
WRITE(NOUT,915) N,M
915 FORMAT(5H-N = ,I3,7H M =,I2)
CALL QUINAT(N,X,Y,BB,CC,DD,EE,FF)
DO 920 I = 1,N
A(I,1) = Y(I)
A(I,2) = BB(I)
A(I,3) = CC(I)
A(I,4) = DD(I)
A(I,5) = EE(I)
A(I,6) = FF(I)
920 CONTINUE
DO 930 I = 1,MM1
DIFF(I) = 0.0
COM(I) = 0.0
930 CONTINUE
DO 970 K = 1,N
DO 935 I = 1,MM
C(I) = A(K,I)
935 CONTINUE
WRITE(NOUT,940) K
940 FORMAT(40H ---------------------------------------,I3,
* 45H --------------------------------------------)
WRITE(NOUT,945) X(K)
945 FORMAT(F12.8)
IF (K .EQ. N) WRITE(NOUT,955) C(1)
IF (K .EQ. N) GO TO 975
WRITE(NOUT,955) (C(I),I=1,MM)
DO 950 I = 1,MM1
IF (ABS(A(K,I)) .GT. COM(I)) COM(I) = ABS(A(K,I))
950 CONTINUE
955 FORMAT(6F16.8)
Z = X(K+1) - X(K)
DO 960 I = 2,MM
DO 960 JJ = I,MM
J = MM + I - JJ
C(J-1) = C(J)*Z + C(J-1)
960 CONTINUE
WRITE(NOUT,955) (C(I),I=1,MM)
DO 965 I = 1,MM1
IF ((K .GE. N-1) .AND. (I .NE. 1)) GO TO 965
Z = ABS(C(I) - A(K+1,I))
IF (Z .GT. DIFF(I)) DIFF(I) = Z
965 CONTINUE
970 CONTINUE
975 WRITE(NOUT,980)
980 FORMAT(41H MAXIMUM ABSOLUTE VALUES OF DIFFERENCES )
WRITE(NOUT,985) (DIFF(I),I=1,MM1)
985 FORMAT(5E18.9)
WRITE(NOUT,990)
990 FORMAT(42H MAXIMUM ABSOLUTE VALUES OF COEFFICIENTS )
IF (ABS(C(1)) .GT. COM(1)) COM(1) =ABS(C(1))
WRITE(NOUT,995) (COM(I),I=1,MM1)
995 FORMAT(5F16.8)
STOP
END