Directory of image this file is from
This file as a plain text file
***** OS8 IS UP ********
TEST FPP2A FAILS BY HANGING AT LOC 1614 - 1615
ALTHOUGH FPP WORKS WITH SIGSYS-12, IT DOES NOT WORK
WITH FORTRAN-IV (FRTS).
RETURN TO OS8 BY TYPING ESC QUIT ESC ESC
ESC = ESCAPE KEY
ALL DIAGNOSTICS ARE ON MOUNTED DISK
DISK BOOT IS C(30) = 6743
C(31) = 5031 START LSW AT 30 IN 8 MODE
START FPP2A BY CALLING .R FPP2A
@ A8F PARAMETERS
C X - THE ARGUMENT VALUE SPECIFIED BY INPUT.
C ARG - THE INPUT VECTOR (DIMENSION NDIM) OF ARGUMENT
C VALUES OF THE TABLE (POSSIBLY DESTROYED).
C VAL - THE INPUT VECTOR (DIMENSION NDIM) OF FUNCTION
C VALUES OF THE TABLE (DESTROYED).
C Y - THE RESULTING INTERPOLATED FUNCTION VALUE.
C NDIM - AN INPUT VALUE WHICH SPECIFIES THE NUMBER OF
C POINTS IN TABLE (ARG,VAL).
C EPS - AN INPUT CONSTANT WHICH IS USED AS UPPER BOUND
C FOR THE ABSOLUTE ERROR.
C IER - A RESULTING ERROR PARAMETER.
C
C REMARKS
C (1) TABLE (ARG,VAL) SHOULD REPRESENT A SINGLE-VALUED
C FUNCTION AND SHOULD BE STORED IN SUCH A WAY, THAT THE
C DISTANCES ABS(ARG(I)-X) INCREASE WITH INCREASING
C SUBSCRIPT I. TO GENERATE THIS ORDER IN TABLE (ARG,VAL),
C SUBROUTINES ATSG, ATSM OR ATSE COULD BE USED IN A
C PREVIOUS STAGE.
C (2) NO ACTION BESIDES ERROR MESSAGE IN CASE NDIM LESS
C THAN 1.
C (3) INTERPOLATION IS TERMINATED EITHER IF THE DIFFERENCE
C BETWEEN TWO SUCCESSIVE INTERPOLATED VALUES IS
C ABSOLUTELY LESS THAN TOLERANCE EPS, OR IF THE ABSOLUTE
C VALUE OF THIS DIFFERENCE STOPS DIMINISHING, OR AFTER
C (NDIM-1) STEPS (THE NUMBER OF POSSIBLE STEPS IS
C DIMINISHED IF AT ANY STAGE INFINITY ELEMENT APPEARS IN
C THE DOWNWARD DIAGONAL OF INVERTED-DIFFERENCES-SCHEME
C AND IF IT IS IMPOSSIBLE TO ELIMINATE THIS INFINITY
C ELEMENT BY INTERCHANGING OF TABLE POINTS).
C FURTHER IT IS TERMINATED IF THE PROCEDURE DISCOVERS TWO
C ARGUMENT VALUES IN VECTOR ARG WHICH ARE IDENTICAL.
C DEPENDENT ON THESE FOUR CASES, ERROR PARAMETER IER IS
C CODED IN THE FOLLOWING FORM
C IER=0 - IT WAS POSSIBLE TO REACH THE REQUIRED
C ACCURACY (NO ERROR).
C IER=1 - IT WAS IMPOSSIBLE TO REACH THE REQUIRED
C ACCURACY BECAUSE OF ROUNDING ERRORS.
C IER=2 - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE
C NDIM IS LESS THAN 2, OR THE REQUIRED ACCURACY
C COULD NOT BE REACHED BY MEANS OF THE GIVEN
C TABLE. NDIM SHOULD BE INCREASED.
C IER=3 - THE PROCEDURE DISCOVERED TWO ARGUMENT VALUES
C IN VECTOR ARG WHICH ARE IDENTICAL.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C METHOD
C INTERPOLATION IS DONE BY CONTINUED FRACTIONS AND INVERTED-
C DIFFERENCES-SCHEME. ON RETURN Y CONTAINS AN INTERPOLATED
C FUNCTION VALUE AT POINT X, WHICH IS IN THE SENSE OF REMARK
C (3) OPTIMAL WITH RESPECT TO GIVEN TABLE. FOR REFERENCE, SEE
C F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,
C MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.395-406.
C
C ..................................................................
C
SUBROUTINE ACFI(X,ARG,VAL,Y,NDIM,EPS,IER)
C
C
DIMENSION ARG(1),VAL(1)
IER=2
IF(NDIM)20,20,1
1 Y=VAL(1)
DELT2=0.
IF(NDIM-1)20,20,2
C
C PREPARATIONS FOR INTERPOLATION LOOP
2 P2=1.
P3=Y
Q2=0.
Q3=1.
C
C
C START INTERPOLATION LOOP
DO 16 I=2,NDIM
II=0
P1=P2
P2=P3
Q1=Q2
Q2=Q3
Z=Y
DELT1=DELT2
JEND=I-1
C
C COMPUTATION OF INVERTED DIFFERENCES
3 AUX=VAL(I)
DO 10 J=1,JEND
H=VAL(I)-VAL(J)
IF(ABS(H)-1.E-6*ABS(VAL(I)))4,4,9
4 IF(ARG(I)-ARG(J))5,17,5
5 IF(J-JEND)8,6,6
C
C INTERCHANGE ROW I WITH ROW I+II
6 II=II+1
III=I+II
IF(III-NDIM)7,7,19
7 VAL(I)=VAL(III)
VAL(III)=AUX
AUX=ARG(I)
ARG(I)=ARG(III)
ARG(III)=AUX
GOTO 3
C
C COMPUTATION OF VAL(I) IN CASE VAL(I)=VAL(J) AND J LESS THAN I-1
8 VAL(I)=1.7E38 0
GOTO 10
C
C COMPUTATION OF VAL(I) IN CASE VAL(I) NOT EQUAL TO VAL(J)
9 VAL(I)=(ARG(I)-ARG(J))/H
10 CONTINUE
C INVERTED DIFFERENCES ARE COMPUTED
C
C COMPUTATION OF NEW Y
P3=VAL(I)*P2+(X-ARG(I-1))*P1
Q3=VAL(I)*Q2+(X-ARG(I-1))*Q1
IF(Q3)11,12,11
11 Y=P3/Q3
GOTO 13
12 Y=1.7E38 0
13 DELT2=ABS(Z-Y)
IF(DELT2-EPS)19,19,14
14 IF(I-8)16,15,15
15 IF(DELT2-DELT1)16,18,18
16 CONTINUE
C END OF INTERPOLATION LOOP
C
C
RETURN
C
C THERE ARE TWO IDENTICAL ARGUMENT VALUES IN VECTOR ARG
17 IER=3
RETURN
C
C TEST VALUE DELT2 STARTS OSCILLATING
18 Y=Z
IER=1
RETURN
C
C THERE IS SATISFACTORY ACCURACY WITHIN NDIM-1 STEPS
19 IER=0
20 RETURN
END
C
C ..................................................................
C
C SAMPLE MAIN PROGRAM FOR MATRIX ADDITION - ADSAM
C
C PURPOSE
C MATRIX ADDITION SAMPLE PROGRAM
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C MADD
C MATIN
C MXOUT
C LOC
C
C METHOD
C TWO INPUT MATRICES ARE READ FROM THE STANDARD INPUT DEVICE.
C THEY ARE ADDED AND THE RESULTANT MATRIX IS LISTED ON
C THE STANDARD OUTPUT DEVICE. THIS CAN BE REPEATED FOR ANY
C NUMBER OF PAIRS OF MATRICES UNTIL A BLANK CARD IS
C ENCOUNTERED
C
C ..................................................................
C
C MATRICES ARE DIMENSIONED FOR 1000 ELEMENTS. THEREFORE, PRODUCT
C OF NUMBER OF ROWS BY NUMBER OF COLUMNS CANNOT EXCEED 1000.
C
DIMENSION A(1000),B(1000),R(1000)
C
10 FORMAT(1H1,15HMATRIX ADDITION)
11 FORMAT(1H0,44HDIMENSIONED AREA TOO SMALL FOR INPUT MATRIX ,I4)
12 FORMAT(1H0,20HEXECUTION TERMINATED)
13 FORMAT(1H0,32HMATRIX DIMENSIONS NOT CONSISTENT)
14 FORMAT(1H0,42HINCORRECT NUMBER OF DATA CARDS FOR MATRIX ,I4)
15 FORMAT(1H0,18HGO ON TO NEXT CASE)
16 FORMAT(1H0,11HEND OF CASE)
C OPEN (UNIT=5, DEVICE='CDR', ACCESS='SEQIN')
C OPEN (UNIT=6, DEVICE='LPT', ACCESS='SEQOUT')
C
C ..................................................................
C
WRITE(6,10)
20 CALL MATIN(ICODA,A,1000,NA,MA,MSA,IER)
IF( NA ) 25,95,25
25 IF(IER-1) 40,30,35
30 WRITE(6,11) ICODA
GO TO 45
35 WRITE(6,14) ICODA
37 WRITE(6,12)
GO TO 95
40 CALL MXOUT(ICODA,A,NA,MA,MSA,60,120,2)
45 CALL MATIN(ICODB,B,1000,NB,MB,MSB,IER)
IF(IER-1) 60,50,55
50 WRITE(6,11) ICODB
WRITE(6,15)
GO TO 20
55 WRITE(6,14) ICODB
GO TO 37
60 IF(NA-NB) 75,70,75
70 IF(MA-MB) 75,80,75
75 WRITE(6,13)
WRITE(6,15)
GO TO 20
80 CALL MXOUT(ICODB,B,NB,MB,MSB,60,120,2)
ICODR=ICODA+ICODB
CALL MADD(A,B,R,NA,MA,MSA,MSB)
MSR=MSA
IF(MSA-MSB) 90,90,85
85 MSR=MSB
90 CALL MXOUT(ICODR,R,NA,MA,MSR,60,120,2)
WRITE(6,16)
GO TO 20
95 CONTINUE
END
C
C ..................................................................
C
C SUBROUTINE AHI
C
C PURPOSE
C TO INTERPOLATE FUNCTION VALUE Y FOR A GIVEN ARGUMENT VALUE
C X USING A GIVEN TABLE (ARG,VAL) OF ARGUMENT, FUNCTION, AND
C DERIVATIVE VALUES.
C
C USAGE
C CALL AHI (X,ARG,VAL,Y,NDIM,EPS,IER)
C
C DESCRIPTION OF PARAMETERS
C X - THE ARGUMENT VALUE SPECIFIED BY INPUT.
C ARG - THE INPUT VECTOR (DIMENSION NDIM) OF ARGUMENT
C VALUES OF THE TABLE (NOT DESTROYED).
C VAL - THE INPUT VECTOR (DIMENSION 2*NDIM) OF FUNCTION
C AND DERIVATIVE VALUES OF THE TABLE (DESTROYED).
C FUNCTION AND DERIVATIVE VALUES MUST BE STORED IN
C PAIRS, THAT MEANS BEGINNING WITH FUNCTION VALUE AT
C POINT ARG(1) EVERY FUNCTION VALUE MUST BE FOLLOWED
C BY THE VALUE OF DERIVATIVE AT THE SAME POINT.
C Y - THE RESULTING INTERPOLATED FUNCTION VALUE.
C NDIM - AN INPUT VALUE WHICH SPECIFIES THE NUMBER OF
C POINTS IN TABLE (ARG,VAL).
C EPS - AN INPUT CONSTANT WHICH IS USED AS UPPER BOUND
C FOR THE ABSOLUTE ERROR.
C IER - A RESULTING ERROR PARAMETER.
C
C REMARKS
C (1) TABLE (ARG,VAL) SHOULD REPRESENT A SINGLE-VALUED
C FUNCTION AND SHOULD BE STORED IN SUCH A WAY, THAT THE
C DISTANCES ABS(ARG(I)-X) INCREASE WITH INCREASING
C SUBSCRIPT I. TO GENERATE THIS ORDER IN TABLE (ARG,VAL),
C SUBROUTINES ATSG, ATSM OR ATSE COULD BE USED IN A
C PREVIOUS STAGE.
C (2) NO ACTION BESIDES ERROR MESSAGE IN CASE NDIM LESS
C THAN 1.
C (3) INTERPOLATION IS TERMINATED EITHER IF THE DIFFERENCE
C BETWEEN TWO SUCCESSIVE INTERPOLATED VALUES IS
C ABSOLUTELY LESS THAN TOLERANCE EPS, OR IF THE ABSOLUTE
C VALUE OF THIS DIFFERENCE STOPS DIMINISHING, OR AFTER
C (2*NDIM-2) STEPS. FURTHER IT IS TERMINATED IF THE
C PROCEDURE DISCOVERS TWO ARGUMENT VALUES IN VECTOR ARG
C WHICH ARE IDENTICAL. DEPENDENT ON THESE FOUR CASES,
C ERROR PARAMETER IER IS CODED IN THE FOLLOWING FORM
C IER=0 - IT WAS POSSIBLE TO REACH THE REQUIRED
C ACCURACY (NO ERROR).
C IER=1 - IT WAS IMPOSSIBLE TO REACH THE REQUIRED
C ACCURACY BECAUSE OF ROUNDING ERRORS.
C IER=2 - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE
C NDIM IS LESS THAN 2, OR THE REQUIRED ACCURACY
C COULD NOT BE REACHED BY MEANS OF THE GIVEN
C TABLE. NDIM SHOULD BE INCREASED.
C IER=3 - THE PROCEDURE DISCOVERED TWO ARGUMENT VALUES
C IN VECTOR ARG WHICH ARE IDENTICAL.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C METHOD
C INTERPOLATION IS DONE BY MEANS OF AITKENS SCHEME OF
C HERMITE INTERPOLATION. ON RETURN Y CONTAINS AN INTERPOLATED
C FUNCTION VALUE AT POINT X, WHICH IS IN THE SENSE OF REMARK
C (3) OPTIMAL WITH RESPECT TO GIVEN TABLE. FOR REFERENCE, SEE
C F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,
C MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.314-317, AND
C GERSHINSKY/LEVINE, AITKEN-HERMITE INTERPOLATION,
C JACM, VOL.11, ISS.3 (1964), PP.352-356.
C
C ..................................................................
C
SUBROUTINE AHI(X,ARG,VAL,Y,NDIM,EPS,IER)
C
C
DIMENSION ARG(1),VAL(1)
IER=2
H2=X-ARG(1)
IF(NDIM-1)2,1,3
1 Y=VAL(1)+VAL(2)*H2
2 RETURN
C
C VECTOR ARG HAS MORE THAN 1 ELEMENT.
C THE FIRST STEP PREPARES VECTOR VAL SUCH THAT AITKEN SCHEME CAN BE
C USED.
3 I=1
DO 5 J=2,NDIM
H1=H2
H2=X-ARG(J)
Y=VAL(I)
VAL(I)=Y+VAL(I+1)*H1
H=H1-H2
IF(H)4,13,4
4 VAL(I+1)=Y+(VAL(I+2)-Y)*H1/H
5 I=I+2
VAL(I)=VAL(I)+VAL(I+1)*H2
C END OF FIRST STEP
C
C PREPARE AITKEN SCHEME
DELT2=0.
IEND=I-1
C
C START AITKEN-LOOP
DO 9 I=1,IEND
DELT1=DELT2
Y=VAL(1)
M=(I+3)/2
H1=ARG(M)
DO 6 J=1,I
K=I+1-J
L=(K+1)/2
H=ARG(L)-H1
IF(H)6,14,6
6 VAL(K)=(VAL(K)*(X-H1)-VAL(K+1)*(X-ARG(L)))/H
DELT2=ABS(Y-VAL(1))
IF(DELT2-EPS)11,11,7
7 IF(I-5)9,8,8
8 IF(DELT2-DELT1)9,12,12
9 CONTINUE
C END OF AITKEN-LOOP
C
10 Y=VAL(1)
RETURN
C
C THERE IS SUFFICIENT ACCURACY WITHIN 2*NDIM-2 ITERATION STEPS
11 IER=0
GOTO 10
C
C TEST VALUE DELT2 STARTS OSCILLATING
12 IER=1
RETURN
C
C THERE ARE TWO IDENTICAL ARGUMENT VALUES IN VECTOR ARG
13 Y=VAL(1)
14 IER=3
RETURN
END
C
C ..................................................................
C
C SUBROUTINE ALI
C
C PURPOSE
C TO INTERPOLATE FUNCTION VALUE Y FOR A GIVEN ARGUMENT VALUE
C X USING A GIVEN TABLE (ARG,VAL) OF ARGUMENT AND FUNCTION
C VALUES.
C
C USAGE
C CALL ALI (X,ARG,VAL,Y,NDIM,EPS,IER)
C
C DESCRIPTION OF PARAMETERS
C X - THE ARGUMENT VALUE SPECIFIED BY INPUT.
C ARG - THE INPUT VECTOR (DIMENSION NDIM) OF ARGUMENT
C VALUES OF THE TABLE (NOT DESTROYED).
C VAL - THE INPUT VECTOR (DIMENSION NDIM) OF FUNCTION
C VALUES OF THE TABLE (DESTROYED).
C Y - THE RESULTING INTERPOLATED FUNCTION VALUE.
C NDIM - AN INPUT VALUE WHICH SPECIFIES THE NUMBER OF
C POINTS IN TABLE (ARG,VAL).
C EPS - AN INPUT CONSTANT WHICH IS USED AS UPPER BOUND
C FOR THE ABSOLUTE ERROR.
C IER - A RESULTING ERROR PARAMETER.
C
C REMARKS
C (1) TABLE (ARG,VAL) SHOULD REPRESENT A SINGLE-VALUED
C FUNCTION AND SHOULD BE STORED IN SUCH A WAY, THAT THE
C DISTANCES ABS(ARG(I)-X) INCREASE WITH INCREASING
C SUBSCRIPT I. TO GENERATE THIS ORDER IN TABLE (ARG,VAL),
C SUBROUTINES ATSG, ATSM OR ATSE COULD BE USED IN A
C PREVIOUS STAGE.
C (2) NO ACTION BESIDES ERROR MESSAGE IN CASE NDIM LESS
C THAN 1.
C (3) INTERPOLATION IS TERMINATED EITHER IF THE DIFFERENCE
C BETWEEN TWO SUCCESSIVE INTERPOLATED VALUES IS
C ABSOLUTELY LESS THAN TOLERANCE EPS, OR IF THE ABSOLUTE
C VALUE OF THIS DIFFERENCE STOPS DIMINISHING, OR AFTER
C (NDIM-1) STEPS. FURTHER IT IS TERMINATED IF THE
C PROCEDURE DISCOVERS TWO ARGUMENT VALUES IN VECTOR ARG
C WHICH ARE IDENTICAL. DEPENDENT ON THESE FOUR CASES,
C ERROR PARAMETER IER IS CODED IN THE FOLLOWING FORM
C IER=0 - IT WAS POSSIBLE TO REACH THE REQUIRED
C ACCURACY (NO ERROR).
C IER=1 - IT WAS IMPOSSIBLE TO REACH THE REQUIRED
C ACCURACY BECAUSE OF ROUNDING ERRORS.
C IER=2 - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE
C NDIM IS LESS THAN 3, OR THE REQUIRED ACCURACY
C COULD NOT BE REACHED BY MEANS OF THE GIVEN
C TABLE. NDIM SHOULD BE INCREASED.
C IER=3 - THE PROCEDURE DISCOVERED TWO ARGUMENT VALUES
C IN VECTOR ARG WHICH ARE IDENTICAL.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C METHOD
C INTERPOLATION IS DONE BY MEANS OF AITKENS SCHEME OF
C LAGRANGE INTERPOLATION. ON RETURN Y CONTAINS AN INTERPOLATED
C FUNCTION VALUE AT POINT X, WHICH IS IN THE SENSE OF REMARK
C (3) OPTIMAL WITH RESPECT TO GIVEN TABLE. FOR REFERENCE, SEE
C F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,
C MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.49-50.
C
C ..................................................................
C
SUBROUTINE ALI(X,ARG,VAL,Y,NDIM,EPS,IER)
C
C
DIMENSION ARG(1),VAL(1)
IER=2
DELT2=0.
IF(NDIM-1)9,7,1
C
C START OF AITKEN-LOOP
1 DO 6 J=2,NDIM
DELT1=DELT2
IEND=J-1
DO 2 I=1,IEND
H=ARG(I)-ARG(J)
IF(H)2,13,2
2 VAL(J)=(VAL(I)*(X-ARG(J))-VAL(J)*(X-ARG(I)))/H
DELT2=ABS(VAL(J)-VAL(IEND))
IF(J-2)6,6,3
3 IF(DELT2-EPS)10,10,4
4 IF(J-5)6,5,5
5 IF(DELT2-DELT1)6,11,11
6 CONTINUE
C END OF AITKEN-LOOP
C
7 J=NDIM
8 Y=VAL(J)
9 RETURN
C
C THERE IS SUFFICIENT ACCURACY WITHIN NDIM-1 ITERATION STEPS
10 IER=0
GOTO 8
C
C TEST VALUE DELT2 STARTS OSCILLATING
11 IER=1
12 J=IEND
GOTO 8
C
C THERE ARE TWO IDENTICAL ARGUMENT VALUES IN VECTOR ARG
13 IER=3
GOTO 12
END
C
C ..................................................................
C
C SAMPLE MAIN PROGRAM FOR ANALYSIS OF VARIANCE - ANOVA
C
C PURPOSE
C (1) READ THE PROBLEM PARAMETER CARD FOR ANALYSIS OF VARI-
C ANCE, (2) CALL THE SUBROUTINES FOR THE CALCULATION OF SUMS
C OF SQUARES, DEGREES OF FREEDOM AND MEAN SQUARE, AND
C (3) PRINT FACTOR LEVELS, GRAND MEAN AND ANALYSIS OF VARI-
C ANCE TABLE.
C
C REMARKS
C THE PROGRAM HANDLES ONLY COMPLETE FACTORIAL DESIGNS. THERE-
C FORE, OTHER EXPERIMENTAL DESIGN MUST BE REDUCED TO THIS FORM
C PRIOR TO THE USE OF THE PROGRAM.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C AVDAT
C AVCAL
C MEANQ
C
C METHOD
C THE METHOD IS BASED ON THE TECHNIQUE DISCUSSED BY H. O.
C HARTLEY IN 'MATHEMATICAL METHODS FOR DIGITAL COMPUTERS',
C EDITED BY A. RALSTON AND H. WILF, JOHN WILEY AND SONS,
C 1962, CHAPTER 20.
C
C ..................................................................
C
C THE FOLLOWING DIMENSION MUST BE GREATER THAN OR EQUAL TO THE
C CUMULATIVE PRODUCT OF EACH FACTOR LEVEL PLUS ONE (LEVEL(I)+1)
C FOR I=1 TO K, WHERE K IS THE NUMBER OF FACTORS..
C
DIMENSION X(3000)
C
C THE FOLLOWING DIMENSIONS MUST BE GREATER THAN OR EQUAL TO THE
C NUMBER OF FACTORS..
C
DIMENSION HEAD(6),LEVEL(6),ISTEP(6),KOUNT(6),LASTS(6)
C
C THE FOLLOWING DIMENSIONS MUST BE GREATER THAN OR EQUAL TO 2 TO
C THE K-TH POWER MINUS 1, ((2**K)-1)..
C
DIMENSION SUMSQ(63),NDF(63),SMEAN(63)
C
C THE FOLLOWING DIMENSION IS USED TO PRINT FACTOR LABELS IN ANALYSIS
C OF VARIANCE TABLE AND IS FIXED..
C
DIMENSION FMT(15)
C ..................................................................
C
C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C STATEMENT WHICH FOLLOWS.
C
C DOUBLE PRECISION X,GMEAN,SUMSQ,SMEAN,SUM
C
C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C ROUTINE.
C
C ...............................................................
C
1 FORMAT(A4,A2,I2,A4,3X,11(A1,I4)/(A1,I4,A1,I4,A1,I4,A1,I4,A1,I4))
2 FORMAT(26H1ANALYSIS OF VARIANCE.....A4,A2//)
3 FORMAT(18H0LEVELS OF FACTORS/(3X,A1,7X,I4))
4 FORMAT(1H0//11H GRAND MEANF20.5////)
5 FORMAT(10H0SOURCE OF18X,7HSUMS OF10X,10HDEGREES OF9X,4HMEAN/10H VA
1RIATION18X,7HSQUARES11X,7HFREEDOM10X,7HSQUARES/)
6 FORMAT(1H 15A1,F20.5,10X,I6,F20.5)
7 FORMAT(6H TOTAL10X,F20.5,10X,I6)
8 FORMAT(12F6.0)
C OPEN (UNIT=5, DEVICE='CDR', ACCESS='SEQIN')
C OPEN (UNIT=6, DEVICE='LPT', ACCESS='SEQOUT')
C
C ..................................................................
C
C READ PROBLEM PARAMETER CARD
C
LOGICAL EOF
CALL CHKEOF (EOF)
100 READ (5,1) PR,PR1,K,BLANK,(HEAD(I),LEVEL(I),I=1,K)
IF (EOF) GOTO 999
C PR.....PROBLEM NUMBER (MAY BE ALPHAMERIC)
C PR1....PROBLEM NUMBER (CONTINUED)
C K......NUMBER OF FACTORS
C BLANK..BLANK FIELD
C HEAD...FACTOR LABELS
C LEVEL..LEVELS OF FACTORS
C
C PRINT PROBLEM NUMBER AND LEVELS OF FACTORS
C
WRITE (6,2) PR,PR1
WRITE (6,3) (HEAD(I),LEVEL(I),I=1,K)
C
C CALCULATE TOTAL NUMBER OF DATA
C
N=LEVEL(1)
DO 102 I=2,K
102 N=N*LEVEL(I)
C
C READ ALL INPUT DATA
C
READ (5,8) (X(I),I=1,N)
C
CALL AVDAT (K,LEVEL,N,X,L,ISTEP,KOUNT)
CALL AVCAL (K,LEVEL,X,L,ISTEP,LASTS)
CALL MEANQ (K,LEVEL,X,GMEAN,SUMSQ,NDF,SMEAN,ISTEP,KOUNT,LASTS)
C
C PRINT GRAND MEAN
C
WRITE (6,4) GMEAN
C
C PRINT ANALYSIS OF VARIANCE TABLE
C
WRITE (6,5)
LL=(2**K)-1
ISTEP(1)=1
DO 105 I=2,K
105 ISTEP(I)=0
DO 110 I=1,15
110 FMT(I)=BLANK
NN=0
SUM=0.0
120 NN=NN+1
L=0
DO 140 I=1,K
FMT(I)=BLANK
IF(ISTEP(I)) 130, 140, 130
130 L=L+1
FMT(L)=HEAD(I)
140 CONTINUE
WRITE (6,6) (FMT(I),I=1,15),SUMSQ(NN),NDF(NN),SMEAN(NN)
SUM=SUM+SUMSQ(NN)
IF(NN-LL) 145, 170, 170
145 DO 160 I=1,K
IF(ISTEP(I)) 147, 150, 147
147 ISTEP(I)=0
GO TO 160
150 ISTEP(I)=1
GO TO 120
160 CONTINUE
170 N=N-1
WRITE (6,7) SUM,N
GO TO 100
999 STOP
END
C
C ..................................................................
C
C SUBROUTINE APCH
C
C PURPOSE
C SET UP NORMAL EQUATIONS OF LEAST SQUARES FIT IN TERMS OF
C CHEBYSHEV POLYNOMIALS FOR A GIVEN DISCRETE FUNCTION
C
C USAGE
C CALL APCH(DATI,N,IP,XD,X0,WORK,IER)
C
C DESCRIPTION OF PARAMETERS
C DATI - VECTOR OF DIMENSION 3*N (OR DIMENSION 2*N+1)
C CONTAINING THE GIVEN ARGUMENTS, FOLLOWED BY THE
C FUNCTION VALUES AND N (RESPECTIVELY 1) WEIGHT
C VALUES. THE CONTENT OF VECTOR DATI REMAINS
C UNCHANGED.
C N - NUMBER OF GIVEN POINTS
C IP - DIMENSION OF LEAST SQUARES FIT, I.E. NUMBER OF
C CHEBYSHEV POLYNOMIALS USED AS FUNDAMENTAL FUNCTIONS
C IP SHOULD NOT EXCEED N
C XD - RESULTANT MULTIPLICATIVE CONSTANT FOR LINEAR
C TRANSFORMATION OF ARGUMENT RANGE
C X0 - RESULTANT ADDITIVE CONSTANT FOR LINEAR
C TRANSFORMATION OF ARGUMENT RANGE
C WORK - WORKING STORAGE OF DIMENSION (IP+1)*(IP+2)/2
C ON RETURN WORK CONTAINS THE SYMMETRIC COEFFICIENT
C MATRIX OF THE NORMAL EQUATIONS IN COMPRESSED FORM
C FOLLOWED IMMEDIATELY BY RIGHT HAND SIDE
C AND SQUARE SUM OF FUNCTION VALUES
C IER - RESULTING ERROR PARAMETER
C IER =-1 MEANS FORMAL ERRORS IN DIMENSION
C IER = 0 MEANS NO ERRORS
C IER = 1 MEANS COINCIDING ARGUMENTS
C
C REMARKS
C NO WEIGHTS ARE USED IF THE VALUE OF DATI(2*N+1) IS
C NOT POSITIVE.
C EXECUTION OF SUBROUTINE APCH IS A PREPARATORY STEP FOR
C CALCULATION OF LEAST SQUARES FITS IN CHEBYSHEV POLYNOMIALS
C IT SHOULD BE FOLLOWED BY EXECUTION OF SUBROUTINE APFS
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C METHOD
C THE LEAST SQUARE FIT IS DETERMINED USING CHEBYSHEV
C POLYNOMIALS AS FUNDAMENTAL FUNCTION SYSTEM.
C THE METHOD IS DISCUSSED IN THE ARTICLE
C A.T.BERZTISS, LEAST SQUARES FITTING TO IRREGULARLY SPACED
C DATA, SIAM REVIEW, VOL.6, ISS.3, 1964, PP. 203-227.
C
C ..................................................................
C
SUBROUTINE APCH(DATI,N,IP,XD,X0,WORK,IER)
C
C
C DIMENSIONED DUMMY VARIABLES
DIMENSION DATI(1),WORK(1)
C
C CHECK FOR FORMAL ERRORS IN SPECIFIED DIMENSIONS
IF(N-1)19,20,1
1 IF(IP)19,19,2
C
C SEARCH SMALLEST AND LARGEST ARGUMENT
2 IF(IP-N)3,3,19
3 XA=DATI(1)
X0=XA
XE=0.
DO 7 I=1,N
XM=DATI(I)
IF(XA-XM)5,5,4
4 XA=XM
5 IF(X0-XM)6,7,7
6 X0=XM
7 CONTINUE
C
C INITIALIZE CALCULATION OF NORMAL EQUATIONS
XD=X0-XA
M=(IP*(IP+1))/2
IEND=M+IP+1
MT2=IP+IP
MT2M=MT2-1
C
C SET WORKING STORAGE AND RIGHT HAND SIDE TO ZERO
DO 8 I=1,IP
J=MT2-I
WORK(J)=0.
WORK(I)=0.
K=M+I
8 WORK(K)=0.
C
C CHECK FOR DEGENERATE ARGUMENT RANGE
IF(XD)20,20,9
C
C CALCULATE CONSTANTS FOR REDUCTION OF ARGUMENTS
9 X0=-(X0+XA)/XD
XD=2./XD
SUM=0.
C
C START GREAT LOOP OVER ALL GIVEN POINTS
DO 15 I=1,N
T=DATI(I)*XD+X0
J=I+N
DF=DATI(J)
C
C CALCULATE AND STORE VALUES OF CHEBYSHEV POLYNOMIALS
C FOR ARGUMENT T
XA=1.
XM=T
IF(DATI(2*N+1))11,11,10
10 J=J+N
XA=DATI(J)
XM=T*XA
11 T=T+T
SUM=SUM+DF*DF*XA
DF=DF+DF
J=1
12 K=M+J
WORK(K)=WORK(K)+DF*XA
13 WORK(J)=WORK(J)+XA
IF(J-MT2M)14,15,15
14 J=J+1
XE=T*XM-XA
XA=XM
XM=XE
IF(J-IP)12,12,13
15 CONTINUE
WORK(IEND)=SUM+SUM
C
C CALCULATE MATRIX OF NORMAL EQUATIONS
LL=M
KK=MT2M
JJ=1
K=KK
DO 18 J=1,M
WORK(LL)=WORK(K)+WORK(JJ)
LL=LL-1
IF(K-JJ)16,16,17
16 KK=KK-2
K=KK
JJ=1
GOTO 18
17 JJ=JJ+1
K=K-1
18 CONTINUE
IER=0
RETURN
C
C ERROR RETURN IN CASE OF FORMAL ERRORS
19 IER=-1
RETURN
C
C ERROR RETURN IN CASE OF COINCIDING ARGUMENTS
20 IER=1
RETURN
END
C
C ..................................................................
C
C SUBROUTINE APFS
C
C PURPOSE
C PERFORM SYMMETRIC FACTORIZATION OF THE MATRIX OF THE NORMAL
C EQUATIONS FOLLOWED BY CALCULATION OF THE LEAST SQUARES FIT
C OPTIONALLY
C
C USAGE
C CALL APFS(WORK,IP,IRES,IOP,EPS,ETA,IER)
C
C DESCRIPTION OF PARAMETERS
C WORK - GIVEN SYMMETRIC COEFFICIENT MATRIX, STORED
C COMPRESSED, I.E UPPER TRIANGULAR PART COLUMNWISE.
C THE GIVEN RIGHT HAND SIDE OCCUPIES THE NEXT IP
C LOCATIONS IN WORK. THE VERY LAST COMPONENT OF WORK
C CONTAINS THE SQUARE SUM OF FUNCTION VALUES E0
C THIS SCHEME OF STORAGE ALLOCATION IS PRODUCED E.G.
C BY SUBROUTINE APLL.
C THE GIVEN MATRIX IS FACTORED IN THE FORM
C TRANSPOSE(T)*T AND THE GIVEN RIGHT HAND SIDE IS
C DIVIDED BY TRANSPOSE(T).
C THE UPPER TRIANGULAR FACTOR T IS RETURNED IN WORK IF
C IOP EQUALS ZERO.
C IN CASE OF NONZERO IOP THE CALCULATED SOLUTIONS ARE
C STORED IN THE COLUMNS OF TRIANGULAR ARRAY WORK OF
C CORRESPONDING DIMENSION AND E0 IS REPLACED BY THE
C SQUARE SUM OF THE ERRORS FOR FIT OF DIMENSION IRES.
C THE TOTAL DIMENSION OF WORK IS (IP+1)*(IP+2)/2
C IP - NUMBER OF FUNDAMENTAL FUNCTIONS USED FOR LEAST
C SQUARES FIT
C IRES - DIMENSION OF CALCULATED LEAST SQUARES FIT.
C LET N1, N2, DENOTE THE FOLLOWING NUMBERS
C N1 = MAXIMAL DIMENSION FOR WHICH NO LOSS OF
C SIGNIFICANCE WAS INDICATED DURING FACTORIZATION
C N2 = SMALLEST DIMENSION FOR WHICH THE SQUARE SUM OF
C THE ERRORS DOES NOT EXCEED TEST=ABS(ETA*FSQ)
C THEN IRES=MINO(IP,N1) IF IOP IS NONNEGATIVE
C AND IRES=MINO(IP,N1,N2) IF IOP IS NEGATIVE
C IOP - INPUT PARAMETER FOR SELECTION OF OPERATION
C IOP = 0 MEANS TRIANGULAR FACTORIZATION, DIVISION OF
C THE RIGHT HAND SIDE BY TRANSPOSE(T) AND
C CALCULATION OF THE SQUARE SUM OF ERRORS IS
C PERFORMED ONLY
C IOP = +1 OR -1 MEANS THE SOLUTION OF DIMENSION IRES
C IS CALCULATED ADDITIONALLY
C IOP = +2 OR -2 MEANS ALL SOLUTIONS FOR DIMENSION ONE
C UP TO IRES ARE CALCULATED ADDITIONALLY
C EPS - RELATIVE TOLERANCE FOR TEST ON LOSS OF SIGNIFICANCE.
C A SENSIBLE VALUE IS BETWEEN 1.E-3 AND 1.E-6
C ETA - RELATIVE TOLERANCE FOR TOLERATED SQUARE SUM OF
C ERRORS. A REALISTIC VALUE IS BETWEEN 1.E0 AND 1.E-6
C IER - RESULTANT ERROR PARAMETER
C IER =-1 MEANS NONPOSITIVE IP
C IER = 0 MEANS NO LOSS OF SIGNIFICANCE DETECTED
C AND SPECIFIED TOLERANCE OF ERRORS REACHED
C IER = 1 MEANS LOSS OF SIGNIFICANCE DETECTED OR
C SPECIFIED TOLERANCE OF ERRORS NOT REACHED
C
C REMARKS
C THE ABSOLUTE TOLERANCE USED INTERNALLY FOR TEST ON LOSS OF
C SIGNIFICANCE IS TOL=ABS(EPS*WORK(1)).
C THE ABSOLUTE TOLERANCE USED INTERNALLY FOR THE SQUARE SUM OF
C ERRORS IS ABS(ETA*FSQ).
C IOP GREATER THAN 2 HAS THE SAME EFFECT AS IOP = 2.
C IOP LESS THAN -2 HAS THE SAME EFFECT AS IOP =-2.
C IRES = 0 MEANS THE ABSOLUTE VALUE OF EPS IS NOT LESS THAN
C ONE AND/OR WORK(1) IS NOT POSITIVE AND/OR IP IS NOT POSITIVE
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C METHOD
C CALCULATION OF THE LEAST SQUARES FITS IS DONE USING
C CHOLESKYS SQUARE ROOT METHOD FOR SYMMETRIC FACTORIZATION.
C THE INCORPORATED TEST ON LOSS OF SIGNIFICANCE MEANS EACH
C RADICAND MUST BE GREATER THAN THE INTERNAL ABSOLUTE
C TOLERANCE TOL=ABS(EPS*WORK(1)).
C IN CASE OF LOSS OF SIGNIFICANCE IN THE ABOVE SENSE ONLY A
C SUBSYSTEM OF THE NORMAL EQUATIONS IS SOLVED.
C IN CASE OF NEGATIVE IOP THE TRIANGULAR FACTORIZATION IS
C TERMINATED PREMATURELY EITHER IF THE SQUARE SUM OF THE
C ERRORS DOES NOT EXCEED ETA*FSQ OR IF THERE IS INDICATION
C FOR LOSS OF SIGNIFICANCE
C
C ..................................................................
C
SUBROUTINE APFS(WORK,IP,IRES,IOP,EPS,ETA,IER)
C
C
C DIMENSIONED DUMMY VARIABLES
DIMENSION WORK(1)
IRES=0
C
C TEST OF SPECIFIED DIMENSION
IF(IP)1,1,2
C
C ERROR RETURN IN CASE OF ILLEGAL DIMENSION
1 IER=-1
RETURN
C
C INITIALIZE FACTORIZATION PROCESS
2 IPIV=0
IPP1=IP+1
IER=1
ITE=IP*IPP1/2
IEND=ITE+IPP1
TOL=ABS(EPS*WORK(1))
TEST=ABS(ETA*WORK(IEND))
C
C START LOOP OVER ALL ROWS OF WORK
DO 11 I=1,IP
IPIV=IPIV+I
JA=IPIV-IRES
JE=IPIV-1
C
C FORM SCALAR PRODUCT NEEDED TO MODIFY CURRENT ROW ELEMENTS
JK=IPIV
DO 9 K=I,IPP1
SUM=0.
IF(IRES)5,5,3
3 JK=JK-IRES
DO 4 J=JA,JE
SUM=SUM+WORK(J)*WORK(JK)
4 JK=JK+1
5 IF(JK-IPIV)6,6,8
C
C TEST FOR LOSS OF SIGNIFICANCE
6 SUM=WORK(IPIV)-SUM
IF(SUM-TOL)12,12,7
7 SUM=SQRT(SUM)
WORK(IPIV)=SUM
PIV=1./SUM
GOTO 9
C
C UPDATE OFF-DIAGONAL TERMS
8 SUM=(WORK(JK)-SUM)*PIV
WORK(JK)=SUM
9 JK=JK+K
C
C UPDATE SQUARE SUM OF ERRORS
WORK(IEND)=WORK(IEND)-SUM*SUM
C
C RECORD ADDRESS OF LAST PIVOT ELEMENT
IRES=IRES+1
IADR=IPIV
C
C TEST FOR TOLERABLE ERROR IF SPECIFIED
IF(IOP)10,11,11
10 IF(WORK(IEND)-TEST)13,13,11
11 CONTINUE
IF(IOP)12,22,12
C
C PERFORM BACK SUBSTITUTION IF SPECIFIED
12 IF(IOP)14,23,14
13 IER=0
14 IPIV=IRES
15 IF(IPIV)23,23,16
16 SUM=0.
JA=ITE+IPIV
JJ=IADR
JK=IADR
K=IPIV
DO 19 I=1,IPIV
WORK(JK)=(WORK(JA)-SUM)/WORK(JJ)
IF(K-1)20,20,17
17 JE=JJ-1
SUM=0.
DO 18 J=K,IPIV
SUM=SUM+WORK(JK)*WORK(JE)
JK=JK+1
18 JE=JE+J
JK=JE-IPIV
JA=JA-1
JJ=JJ-K
19 K=K-1
20 IF(IOP/2)21,23,21
21 IADR=IADR-IPIV
IPIV=IPIV-1
GOTO 15
C
C NORMAL RETURN
22 IER=0
23 RETURN
END
C
C ..................................................................
C
C SUBROUTINE APLL
C
C PURPOSE
C SET UP NORMAL EQUATIONS FOR A LINEAR LEAST SQUARES FIT
C TO A GIVEN DISCRETE FUNCTION
C
C USAGE
C CALL APLL(FFCT,N,IP,P,WORK,DATI,IER)
C SUBROUTINE FFCT REQUIRES AN EXTERNAL STATEMENT
C
C DESCRIPTION OF PARAMETERS
C FFCT - USER CODED SUBROUTINE WHICH MUST BE DECLARED
C EXTERNAL IN THE MAIN PROGRAM. IT IS CALLED
C CALL FFCT(I,N,IP,P,DATI,WGT,IER) AND RETURNS
C THE VALUES OF THE FUNDAMENTAL FUNCTIONS FOR
C THE I-TH ARGUMENT IN P(1) UP TO P(IP)
C FOLLOWED BY THE I-TH FUNCTION VALUE IN P(IP+1)
C N IS THE NUMBER OF ALL POINTS
C DATI IS A DUMMY PARAMETER WHICH IS USED AS ARRAY
C NAME. THE GIVEN DATA SET MAY BE ALLOCATED IN DATI
C WGT IS THE WEIGHT FACTOR FOR THE I-TH POINT
C IER IS USED AS RESULTANT ERROR PARAMETER IN FFCT
C N - NUMBER OF GIVEN POINTS
C IP - NUMBER OF FUNDAMENTAL FUNCTIONS USED FOR LEAST
C SQUARES FIT
C IP SHOULD NOT EXCEED N
C P - WORKING STORAGE OF DIMENSION IP+1, WHICH
C IS USED AS INTERFACE BETWEEN APLL AND THE USER
C CODED SUBROUTINE FFCT
C WORK - WORKING STORAGE OF DIMENSION (IP+1)*(IP+2)/2.
C ON RETURN WORK CONTAINS THE SYMMETRIC COEFFICIENT
C MATRIX OF THE NORMAL EQUATIONS IN COMPRESSED FORM,
C I.E. UPPER TRINGULAR PART ONLY STORED COLUMNWISE.
C THE FOLLOWING IP POSITIONS CONTAIN THE RIGHT
C HAND SIDE AND WORK((IP+1)*(IP+2)/2) CONTAINS
C THE WEIGHTED SQUARE SUM OF THE FUNCTION VALUES
C DATI - DUMMY ENTRY TO COMMUNICATE AN ARRAY NAME BETWEEN
C MAIN LINE AND SUBROUTINE FFCT.
C IER - RESULTING ERROR PARAMETER
C IER =-1 MEANS FORMAL ERRORS IN SPECIFIED DIMENSIONS
C IER = 0 MEANS NO ERRORS
C IER = 1 MEANS ERROR IN EXTERNAL SUBROUTINE FFCT
C
C REMARKS
C TO ALLOW FOR EASY COMMUNICATION OF INTEGER VALUES
C BETWEEN MAINLINE AND EXTERNAL SUBROUTINE FFCT, THE ERROR
C PARAMETER IER IS TREATED AS A VECTOR OF DIMENSION 1 WITHIN
C SUBROUTINE APLL. ADDITIONAL COMPONENTS OF IER MAY BE
C INTRODUCED BY THE USER FOR COMMUNICATION BACK AND FORTH.
C IN THIS CASE, HOWEVER, THE USER MUST SPECIFY IER AS A
C VECTOR IN HIS MAINLINE.
C EXECUTION OF SUBROUTINE APLL IS A PREPARATORY STEP FOR
C CALCULATION OF THE LINEAR LEAST SQUARES FIT.
C NORMALLY IT IS FOLLOWED BY EXECUTION OF SUBROUTINE APFS
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C THE EXTERNAL SUBROUTINE FFCT MUST BE FURNISHED BY THE USER
C
C METHOD
C HANDLING OF THE GIVEN DATA SET (ARGUMENTS,FUNCTION VALUES
C AND WEIGHTS) IS COMPLETELY LEFT TO THE USER
C ESSENTIALLY HE HAS THREE CHOICES
C (1) THE I-TH VALUES OF ARGUMENT, FUNCTION VALUE AND WEIGHT
C ARE CALCULATED WITHIN SUBROUTINE FFCT FOR GIVEN I.
C (2) THE I-TH VALUES OF ARGUMENT, FUNCTION VALUE AND WEIGHT
C ARE DETERMINED BY TABLE LOOK UP. THE STORAGE LOCATIONS
C REQUIRED ARE ALLOCATED WITHIN THE DUMMY ARRAY DATI
C (POSSIBLY IN P TOO, IN EXCESS OF THE SPECIFIED IP + 1
C LOCATIONS).
C ANOTHER POSSIBILITY WOULD BE TO USE COMMON AS INTERFACE
C BETWEEN MAIN LINE AND SUBROUTINE FFCT AND TO ALLOCATE
C STORAGE FOR THE DATA SET IN COMMON.
C (3) THE I-TH VALUES OF ARGUMENT, FUNCTION VALUE AND WEIGHT
C ARE READ IN FROM AN EXTERNAL DEVICE. THIS MAY BE EASILY
C ACCOMPLISHED SINCE I IS USED STRICTLY INCREASING FROM
C ONE UP TO N WITHIN APLL
C
C ..................................................................
C
SUBROUTINE APLL(FFCT,N,IP,P,WORK,DATI,IER)
C
C
C DIMENSIONED DUMMY VARIABLES
DIMENSION P(1),WORK(1),DATI(1),IER(1)
C
C CHECK FOR FORMAL ERRORS IN SPECIFIED DIMENSIONS
IF(N)10,10,1
1 IF(IP)10,10,2
2 IF(N-IP)10,3,3
C
C SET WORKING STORAGE AND RIGHT HAND SIDE TO