***** 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