C monanova.f C This software implements the method called MONANOVA. C The author of this software is Joseph B Kruskal C Copyright (c) 1993 by AT&T. C Permission to use, copy, modify, and distribute this software for any C purpose without fee is hereby granted, provided that this entire C notice is included in all copies of any software which is or includes C a copy or modification of this software and in all copies of the C supporting documentation for such software. C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMP- C LIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. C This software comes from the FIRST MDS Package of AT&T Bell Laboratories C For explanation of the method this software implements, see C "Analysis of Factorial Experiments by Estimating Monotone C Transformations of the Data" (1965) by Joseph B. Kruskal, in Journal C of the Royal Statistical Society, Series B, 27(2), pp 251-263 C "MONANOVA, a Fortran-IV Program for Monotone Analysis of Variance C (Non-Metric Analysis of Factorial Experiments)" (1969) by Joseph B. C Kruskal and Frank J. Carmone in C Behavioral Science 14, pp 165-166 C "MONANOVA, A Fortran IV Program for Monotone Analysis of Variance" C (1969), Joseph B. Kruskal and Frank J. Carmone, in Journal of C Marketing Research 6, p 497 C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ C MONA 1 C M O N A N O V A FOUR(4) CHARACTER VERSION MONA 2 C NONMETRIC ANALYSIS OF FACTORIAL DESIGN MONA 3 C MONA A3 C DISTRIBUTED BY MONA B3 C MARKETING SCIENCE INSTITUTE MONA C3 C PHILADELPHIA,PENNA. 19104 MONA D3 C MONA 4 C THIS PROGRAM WAS WRITTEN BY MONA 5 C MONA 6 C DR. J. B. KRUSKAL MONA 7 C BELL TELEPHONE LABORATORIES MONA 8 C MURRAY HILL, N. J. MONA 9 C MONA 10 C MONA 11 C ADAPTED FOR IBM S/360 BY MONA 12 C MONA 13 C F. J. CARMONE, JR. MONA 14 C MARKETING SCIENCE INSTITUTE MONA 15 C PHILADELPHIA, PA. MONA 16 C NEW ADDRESS C DEPARTMENT OF ECONOMICS C UNIV OF WATERLOO C WATERLOO, ONTARIO CANADA C MONA 17 C JULY 1968 MONA 18 C MONA 19 C MONA 20 C REFERENCE MONA 21 C KRUSKAL, J. B. "ANALYSIS OF FACTORIAL EXPERIMENTS BY MONA 22 C ESTIMATING MONOTONE TRANSFORMATIONS OF THE DATA", MONA 23 C JOURNAL OF ROYAL STATISTICAL SOCIETY,SERIES B, VOL. 27, MONA 24 C NO. 2 (1965), PP. 251-63. MONA 25 C MONA 26 C MONA 27 C MONA 28 C GENERAL REMARKS MONA 29 C MONA 30 C THIS PROGRAM CONTAINS THE FOLLOWING MAJOR ROUTINES MONA 31 C MONA 32 C OPTEX MAIN ROUTINE MONA 33 C OPSOL CALCULATES STEP SIZE AND OTHER VARIABLES FOR TERMINATIOMONA 34 C CALC CALCULATES STRESS AND READS IN DATA MONA 35 C MFIT PERFORMS WEIGHTED LEAST SQUARES REGRESSION MONA 36 C CCACT READS AND INTERPERTS CONTROL CARDS MONA 37 C SORT SORTS ARRAYS BY USING POINTERS; SIMPLE, BUT VERY FAST MONA 38 C PLOTR PLOTS TWO MODELS VERSUS ORIGINAL DATA MONA 39 C MONA 40 C MONA 41 C NO USE IS MADE OF SPECIAL OR NON-STANDARD SOFTWARE. MONA 42 C MONA 43 C ALL NORMAL INPUT-OUTPUT HANDLED BY CALC AND CCACT. MONA 44 C MONA 45 C MFIT HAS EMERGENCY DIAGNOSTIC OUTPUT. MONA 46 C ALL INPUT AND OUTPUT IS ONTO FORTRAN LOGICAL UNITS CONTROLLED BY MONA 47 C THESE VARIABLES. LREAD, LPRINT, LPUNCH, LSCRAT. MONA 48 C UNIT NUMBERS 5, 6, 7, 2 HAVE BEEN USED RESPECTIVELY. MONA 49 C TO CHANGE THESE ASSIGNMENTS, IT SUFFICES TO CHANGE THE VALUES MONA 50 C FOR THE FOUR VARIABLES SET IN THE BLOCK DATA SUBPROGRAM. MONA 51 C MONA 52 C NOTE THAT THE SCRATCH UNIT IS USED IN A VERY MINOR WAY. MONA 53 C IT IS USED ONLY BY CCACT MONA 54 C MANY INSTALLATIONS WILL HAVE ALTERNATE METHODS OF DOING THE SAME THINGMONA 55 C MONA 56 CALL OPTEX MONA 57 CALL EXIT MONA 58 END MONA 59 C OPTS 1 SUBROUTINE OPTSOL OPTS 2 1 ( SRATST, STRMIN, NOIT, SRTAVW, COSAVW, ACSAVW, OPTS 3 2 LCONSW, STRESS ) OPTS 4 COMMON /OPTEX1/ LREAD,LPRINT,LPUNCH,LSCRAT OPTS 5 C INITIALIZE UNIVERSAL VARIABLES OPTS 6 100 ITNO=0 OPTS 7 GRMULT=0.0 OPTS 8 GRMAG=1.0 OPTS 9 SRATAV=1.2 OPTS 10 COSAV=0.0 OPTS 11 ACSAV=0.2 OPTS 12 LXMIN=1 OPTS 13 C INITIALIZE CONFIGURATION AND GRADIENTS OPTS 14 200 GO TO (210,220), LCONSW OPTS 15 C CREATE INITIAL CONFIGURATION OPTS 16 210 CALL CALC( 2, 0.0, 0.0, 0.0 ) OPTS 17 C NORMALIZE CONFIGURATION AND PRINT INPUT OPTS 18 220 CALL CALC( 6, 0.0, 0.0, 0.0 ) OPTS 19 C PERFORM ITERATION OPTS 20 300 STRLST=STRESS OPTS 21 CALL CALC( 4, STRESS, GRGR, GRGL ) OPTS 22 305 GLMAG=GRMAG OPTS 23 GRMAG = SQRT ( GRGR ) OPTS 24 IF(STRESS) 9999, 320, 310 OPTS 25 310 IF(GRGR ) 9999, 320, 330 OPTS 26 320 LXMIN=2 OPTS 27 GO TO 1000 OPTS 28 330 IF(ITNO) 9999, 340, 350 OPTS 29 340 SRAT=1.2 OPTS 30 CAGRGL=0.0 OPTS 31 GO TO 360 OPTS 32 350 SRAT=STRLST/STRESS OPTS 33 CAGRGL = GRGL / (GRMAG*GLMAG) OPTS 34 C GO TO UPDATING SECTION OPTS 35 360 GO TO 1000 OPTS 36 380 CALL CALC( 5, GRMULT, 0.0, 0.0 ) OPTS 37 ITNO=ITNO+1 OPTS 38 GO TO 300 OPTS 39 C TERMINATION HAS OCCURRED. PERFORM TERMINAL ACTIONS. OPTS 40 C NORMALIZE CONFIGURATION. OPTS 41 400 CALL CALC( 6, 0.0, 0.0,-50.0 ) OPTS 42 C AGAIN CALCULATE STRESS TO GUARANTEE AGAINST ROUNDING OPTS 43 C EFFECTS DURING NORMALIZATION OPTS 44 CALL CALC( 4, STRESS, GRGR, GRGL ) OPTS 45 RETURN OPTS 46 9999 CALL OPTERR( 6HOPTSOL ) OPTS 47 CALL EXIT OPTS 48 C UPDATE VARIOUS VARIABLES. DECIDE WHETHER TO TERMINATE. OPTS 49 1000 IF(ITNO) 9999, 1010, 1100 OPTS 50 1010 WRITE(LPRINT,1) OPTS 51 1 FORMAT( '0HISTORY OF COMPUTATION.'/1X) OPTS 52 WRITE(LPRINT,2) OPTS 53 2 FORMAT( ' ITERATION STRESS SRAT SRTAVG CAGRGL COSAV', OPTS 54 1 ' ACSAV GRMAG GRMULT STEP' ) OPTS 55 LFOOT=1 OPTS 56 1100 GO TO (1110, 1300 ),LXMIN OPTS 57 1110 COSAV = CAGRGL*COSAVW + COSAV*(1.0-COSAVW) OPTS 58 ACSAV = ABS (CAGRGL)*ACSAVW + ACSAV*(1.0-ACSAVW) OPTS 59 SRATAV = (SRAT**SRTAVW) * (SRATAV**(1.0-SRTAVW)) OPTS 60 IF(ITNO) 9999, 1200, 1250 OPTS 61 1200 GRMULT=25.0*STRESS OPTS 62 STEP=GRMULT*GRMAG OPTS 63 GO TO 1300 OPTS 64 1250 ANG=4.0**COSAV OPTS 65 TEMP1 = 1.0 + (AMIN1 (1.0,1.0/SRATAV) ) ** 5 OPTS 66 TEMP3=ABS (COSAV) OPTS 67 TEMP2= 1.0 + ( (ACSAV-TEMP3) / (1.0 - TEMP3 ) ) ** 3 OPTS 68 BIAS = 1.6 / (TEMP1 * TEMP2 ) OPTS 69 GOODLK=SQRT(AMIN1(1.0,SRAT)) OPTS 70 STEP = STEP * ANG * BIAS * GOODLK OPTS 71 GRMULT=STEP/GRMAG OPTS 72 1300 WRITE(LPRINT,3 ) ITNO,STRESS,SRAT,SRATAV,CAGRGL, OPTS 73 1COSAV,ACSAV,GRMAG,GRMULT,STEP OPTS 74 3 FORMAT(I10,F7.3,2F8.4,3F7.3,3F9.5) OPTS 75 2000 IF(STRESS) 9999, 2010, 2020 OPTS 76 2010 LSTOP=1 OPTS 77 GO TO 3000 OPTS 78 2020 IF(GRMAG) 9999, 2030, 2040 OPTS 79 2030 LSTOP=2 OPTS 80 GO TO 3000 OPTS 81 2040 TEMP1 = 0.5 * (1.0+SRATST) OPTS 82 TEMP2 = ABS (TEMP1 - 1.0) OPTS 83 IF( ABS (SRAT-TEMP1) - TEMP2 ) 2050, 2050, 2070 OPTS 84 2050 IF( ABS (SRATAV-TEMP1) - TEMP2 ) 2060, 2060, 2070 OPTS 85 2060 LSTOP=3 OPTS 86 GO TO 3000 OPTS 87 2070 IF(STRESS-STRMIN) 2080, 2080, 2090 OPTS 88 2080 LSTOP=4 OPTS 89 GO TO 3000 OPTS 90 2090 IF(NOIT-ITNO) 9999, 2100, 2110 OPTS 91 2100 LSTOP=5 OPTS 92 GO TO 3000 OPTS 93 2110 GO TO (380, 400), LFOOT OPTS 94 3000 GO TO ( 3010, 3020 ), LFOOT OPTS 95 3010 LFOOT=2 OPTS 96 WRITE(LPRINT, 2) OPTS 97 WRITE(LPRINT, 30) OPTS 98 30 FORMAT(1X) OPTS 99 3020 GO TO ( 3030, 3040, 3050, 3060, 3070), LSTOP OPTS 100 3030 WRITE(LPRINT, 31) OPTS 101 31 FORMAT(24H ZERO STRESS WAS REACHED ) OPTS 102 GO TO 2020 OPTS 103 3040 WRITE(LPRINT,32) OPTS 104 32 FORMAT(21H0MINIMUM WAS ACHIEVED ) OPTS 105 GO TO 2040 OPTS 106 3050 WRITE(LPRINT,33) OPTS 107 33 FORMAT( ' APPROXIMATE MINIMUM WAS REACHED.') OPTS 108 GO TO 2070 OPTS 109 3060 WRITE(LPRINT,34) OPTS 110 34 FORMAT(32H SATISFACTORY STRESS WAS REACHED ) OPTS 111 GO TO 2090 OPTS 112 3070 WRITE(LPRINT,35) OPTS 113 35 FORMAT(39H MAXIMUM NUMBER OF ITERATIONS WERE USED ) OPTS 114 GO TO 2110 OPTS 115 END OPTS 116 C CCAC 1 SUBROUTINE CCACT CCAC 2 CCCACT CCACT FOR MDSCAL CCAC 3 REAL*4 ATAB(4,22) CCAC 4 INTEGER*4 TAB(4,22),MTAB(4) CCAC 5 EQUIVALENCE (TAB,ATAB) CCAC 6 INTEGER LPAR( 25 ) CCAC 7 REAL PAR( 25 ) CCAC 8 EQUIVALENCE (LPAR,PAR) CCAC 9 INTEGER IPARAM(5) CCAC 10 REAL C(73),WORD(3) CCAC 11 REAL CHTAB(13) CCAC 12 REAL BLANK, DOT CCAC 13 INTEGER BLSW,NUMSW,DECSW,TYPE,XTYPE,T,TA, PARNO, TABNO CCAC 14 COMMON /OPTEX1/ LREAD, LPRINT, LPUNCH, LSCRAT CCAC 15 COMMON /OPTEX2/ LPAR CCAC 16 COMMON /OPTEX3/ MTAB, TAB CCAC 17 DATA BLANK, DOT, EQUALS, COMMA, C(73) /1H ,1H.,1H=,1H,,1H / CCAC 18 DATA CHTAB/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H.,1H+,1H-/ CCAC 19 DATA ASTER/1H*/ CCAC 20 1 FORMAT(72A1) CCAC 21 11 FORMAT(1H0,72A1) CCAC 22 2 FORMAT(1X,I1) CCAC 23 3 FORMAT(1X,18A4) CCAC 24 4 FORMAT(1X,I18) CCAC 25 5 FORMAT(1X,F18.3) CCAC 26 6 FORMAT(1X,72A1) CCAC 27 9 FORMAT(51H MDSCAL DIAGNOSTIC. CONTROL CARD ABOVE IS IMPROPER.) CCAC 28 99 FORMAT( 12H ITEM NUMBER, I5, 20H HAS EXPECTED TYPE , I5, CCAC 29 1 18H AND ACTUAL TYPE, I5, 16H . THIS ITEM IS', 3A4, 1H' ) CCAC 30 C READ AND PRINT CONTROL CARD CCAC 31 100 READ (LREAD,1) (C(I),I=1,72) CCAC 32 WRITE(LPRINT,11) (C(I),I=1,72) CCAC 33 C ANALYZE CONTROL CARD INTO TOKENS CCAC 34 C EACH 'TOKEN' IS DELIMITED BY BLANKS CCAC 35 C THERE ARE THREE TYPES OF TOKENS. CCAC 36 C ALPHABETIC, INTEGER, DECIMAL CCAC 37 C ALPHABETIC UNLESS FIRST CHARACTER IS DIGIT OR DEC POINT CCAC 38 C DECIMALS DISTINGUISHED BY DECIMAL POINT CCAC 39 REWIND LSCRAT CCAC 40 BLSW=1 CCAC 41 NUMSW=1 CCAC 42 DECSW=1 CCAC 43 K=0 CCAC 44 DO 300 I=1,73 CCAC 45 X=C(I) CCAC 46 GO TO (110,120), BLSW CCAC 47 110 IF(X.EQ.BLANK .OR. X.EQ.EQUALS .OR. X.EQ. COMMA ) GO TO 300CCAC 48 IF(X.EQ.ASTER) GO TO 310 CCAC 49 BLSW=2 CCAC 50 JA=I CCAC 51 DO 115 L=1,13 CCAC 52 115 IF(X.EQ.CHTAB(L)) NUMSW=2 CCAC 53 IF(X.EQ.DOT) DECSW=2 CCAC 54 GO TO 300 CCAC 55 120 IF(X.EQ.BLANK .OR. X.EQ.EQUALS .OR. X.EQ. COMMA ) GO TO 130CCAC 56 IF(X.EQ.ASTER) GO TO 130 CCAC 57 IF(X.EQ.DOT) DECSW=2 CCAC 58 GO TO 300 CCAC 59 130 K=K+1 CCAC 60 JB = MIN0 (I-1,JA+16) CCAC 61 JC=18-(JB-JA+1) CCAC 62 TYPE=1 CCAC 63 IF(NUMSW.EQ.2) TYPE=NUMSW+DECSW-1 CCAC 64 WRITE (LSCRAT,2) TYPE CCAC 65 GO TO (140,150),NUMSW CCAC 66 140 WRITE (LSCRAT,6) (C(J),J=JA,JB), (BLANK,J=1,JC) CCAC 67 GO TO 160 CCAC 68 150 WRITE (LSCRAT,6) (BLANK,J=1,JC), (C(J),J=JA,JB) CCAC 69 GO TO 160 CCAC 70 160 BLSW=1 CCAC 71 NUMSW=1 CCAC 72 DECSW=1 CCAC 73 IF(X.EQ.ASTER) GO TO 310 CCAC 74 GO TO 300 CCAC 75 300 CONTINUE CCAC 76 C ANALYZE TOKENS AND SET PARAMETER VALUES ACCORDINGLY CCAC 77 310 KB=K CCAC 78 IF(KB.EQ.0) RETURN CCAC 79 REWIND LSCRAT CCAC 80 XTYPE=1 CCAC 81 IPARAM(1) = 1 CCAC 82 DO 1000 J=1,KB CCAC 83 READ (LSCRAT,2) TYPE CCAC 84 IF(TYPE.NE.XTYPE) GO TO 999 CCAC 85 GO TO (400,410,420), XTYPE CCAC 86 400 READ (LSCRAT,3) (WORD(L),L=1,3) CCAC 87 GO TO 510 CCAC 88 410 READ (LSCRAT,4) INTPAR CCAC 89 LPAR(PARNO)=INTPAR CCAC 90 GO TO 430 CCAC 91 420 READ (LSCRAT,5) DECPAR CCAC 92 PAR(PARNO)=DECPAR CCAC 93 430 XTYPE=1 CCAC 94 GO TO 1000 CCAC 95 510 TABNO = IPARAM(1) CCAC 96 MA = MTAB(TABNO) CCAC 97 MB = MTAB(TABNO+1) - 1 CCAC 98 LMISSW = 0 CCAC 99 DO 700 M=MA,MB CCAC 100 IF(WORD(1).NE.ATAB(1,M)) GO TO 700 CCAC 101 LMISSW = 1 CCAC 102 600 XTYPE=1 CCAC 103 IPARAM(1) = 1 CCAC 104 PARNO=TAB(3,M) CCAC 105 LTEMP = TAB(2,M) CCAC 106 GO TO (610,620,630,640,650), LTEMP CCAC 107 C NAME OF INTEGER PARAMETER CCAC 108 610 XTYPE=2 CCAC 109 GO TO 700 CCAC 110 C NAME OF DECIMAL PARAMETER CCAC 111 620 XTYPE=3 CCAC 112 GO TO 700 CCAC 113 C IMPLICITLY SPECIFIED INTEGER PARAMETER CCAC 114 630 LPAR(PARNO) = TAB(4,M) CCAC 115 GO TO 700 CCAC 116 C IMPLICITLY SPECIFIED DECIMAL PARAMETER CCAC 117 640 PAR(PARNO) = ATAB(4,M) CCAC 118 GO TO 700 CCAC 119 C INTERNAL PARAMETER OF CCACT PROGRAM CCAC 120 650 IPARAM(PARNO) = TAB(4,M) CCAC 121 GO TO 700 CCAC 122 700 CONTINUE CCAC 123 IF(LMISSW.EQ.0) GO TO 999 CCAC 124 1000 CONTINUE CCAC 125 1001 RETURN CCAC 126 999 WRITE (LPRINT,9) CCAC 127 WRITE (LPRINT,99) K, XTYPE, TYPE, (WORD(L),L=1,3) CCAC 128 CALL EXIT CCAC 129 END CCAC 130 SUBROUTINE OPTERR( NAME ) OPTR 1 COMMON /OPTEX1/ LREAD,LPRINT,LPUNCH,LSCRAT OPTR 2 REAL*8 NAME OPTR 3 WRITE(LPRINT,1) NAME OPTR 4 1 FORMAT( '0MONANOVA ERROR STOP FROM ', A6 ) OPTR 5 CALL EXIT OPTR 6 END OPTR 7 C OPTX 1 SUBROUTINE OPTEX OPTX 2 DIMENSION CC(12), TITLE(20) OPTX 3 COMMON /OPTEX1/ LREAD,LPRINT,LPUNCH,LSCRAT OPTX 4 COMMON /OPTEX2/ SRATST,STRMIN,NOIT,SRTAVW,COSAVW,ACSAVW, OPTX 5 1 LCNTRL ,CUTOFF,MFITCO,LPUNSW,LPRINS,LCONSW OPTX 6 COMMON /TIL/ TITLE OPTX 7 DATA BLANK/1H / OPTX A7 C INITIALIZE OPTEX OPTX 8 10 CONTINUE OPTX 9 WRITE(LPRINT, 1) OPTX 10 1 FORMAT(1H1) OPTX 11 NOIT = 50 OPTX 12 STRMIN = 0.01 OPTX 13 SRATST=1.001 OPTX 14 SRTAVW=0.33334 OPTX 15 COSAVW=0.6666 OPTX 16 ACSAVW=0.6666 OPTX 17 LCONSW=1 OPTX 18 DO 11 K = 1,20 OPTX A18 11 TITLE(K) = BLANK OPTX B18 C PERFORM INITIAL ACTIONS OPTX 19 CALL CALC( 9, 0.0, 0.0, 0.0 ) OPTX 20 12 FORMAT(20A4) OPTX 21 13 FORMAT(1X,20A4) OPTX 22 100 LCNTRL=1 OPTX 23 CALL CCACT OPTX 24 GO TO (100, 200, 300, 400, 500, 600,700), LCNTRL OPTX 25 C READ DATA AND PERFORM INITIAL MANIPULATIONS OPTX 26 200 CALL CALC( 3, 0.0, 0.0, 0.0 ) OPTX 27 WRITE(LPRINT,2) OPTX 28 GO TO 100 OPTX 29 C READ CONFIGURATION OPTX 30 300 CALL CALC( 1, 0.0, 0.0, 0.0 ) OPTX 31 LCONSW=2 OPTX 32 WRITE(LPRINT,2) OPTX 33 GO TO 100 OPTX 34 C READ PARAMETERS OPTX 35 400 CALL CALC( 7, 0.0, 0.0, 0.0 ) OPTX 36 GO TO 100 OPTX 37 C SOLVE THE SYSTEM OPTX 38 500 CALL OPTSOL(SRATST, STRMIN, NOIT, SRTAVW, COSAVW, ACSAVW, OPTX 39 1 LCONSW, STRESS ) OPTX 40 C CREATE OUTPUT OPTX 41 WRITE(LPRINT,50) STRESS OPTX 42 50 FORMAT( '0FINAL CONFIGURATION HAS STRESS OF', 2PF7.1, OPTX 43 1 ' PERCENT.' ) OPTX 44 CALL CALC( 8, 0.0, 0.0, 0.0 ) OPTX 45 GO TO 10 OPTX 46 C PERFORM TERMINAL FUNCTIONS OPTX 47 600 CALL CALC( 10, 0.0, 0.0, 0.0 ) OPTX 48 C CONTROL MAY OR MAY NOT RETURN HERE. OPTX 49 RETURN OPTX 50 2 FORMAT(1H ) OPTX 51 700 READ(LREAD,12) TITLE OPTX 52 WRITE(LPRINT,13) TITLE OPTX 53 GO TO 100 OPTX 54 END OPTX 55 C MFIT 1 SUBROUTINE MFIT(DIST,DHAT,LBLOCK,MM,LFITSW,IPOINT) MFIT 2 C FIT PERFORMS WEIGHTED-LEAST-SQUARES MONOTONE REGRESSION MFIT 3 C THIS ROUTINE FINDS THE VALUES FOR DHAT WHICH ARE MONOTONIC MFIT 4 C AND WHICH BEST APPROXIMATE THE VALUES OF DIST, MFIT 5 C IN THE SENSE THAT THE SUM OF THE SQUARED DEVIATIONS, MFIT 6 C EACH ONE WEIGHTED BY THE CORRESPONDING VALUE IN WW , MFIT 7 C IS A MINIMUM. MFIT 8 C IT DIRECTLY USES SORT. MFIT 9 DIMENSION IPOINT(1), DIST(1),DHAT(1),LBLOCK(1) MFIT 10 COMMON /OPTEX1/ LREAD,LPRINT,LPUNCH,LSCRAT MFIT A10 C FORM FIRST APPROXIMATION TO CORRECT PARTITON MFIT 11 C IF LFITSW=0, DO PLAIN REGRESSION MFIT 12 C IF LFITSW GT 0, DO BLOCK REGRESSION MFIT 13 C IF LFITSW=1, USE PRIMARY APPROACH. THUS WE SORT EACH MFIT 14 C BLOCK INDICATED BY LBLOCK ACCORDING TO DIST VALUES. MFIT 15 C BLOCK OF EQUAL VALUES OF DATA ACCORDING TO DIST VALUES. MFIT 16 C THEN WE CREATE PARTITON INTO BLOCKS OF SIZE 1 TO START. MFIT 17 C IF LFITSW=2, USE SECONDARY APPROACH. THUS WE START WITHMFIT 18 C PARTITION INTO BLOCKS INDICATED BY LBLOCK. MFIT 19 C WITHIN EACH BLOCK OF WHATEVER SIZE, THE FIRST DHAT VALUEMFIT 20 C GIVES THE WEIGHTED TOTAL OF THE DIST VALUES IN THAT BLOCK, WHILE MFIT 21 C THE LAST DHAT VALUE GIVES THE SUM OF THE WEIGHTS FOR THATMFIT 22 C BLOCK. SIMILARLY, WITHIN EACH BLOCK, THE FIRST LBLOCK MFIT 23 C VALUE AND THE LAST LBLOCK VALUE BOTH GIVE THE SIZE OF THEMFIT 24 C BLOCK. MFIT 25 C BLOCKS OF SIZE 1 FORM A PARTIAL EXCEPTION. EVERYTHIN IS THE SAME MFIT 26 C EXCEPT THAT THE SUM OF THE W VALUES IS NOT FOUND IN THE MFIT 27 C LAST DHAT VALUE IN THE BLOCK. MFIT 28 IF(LFITSW.GT.0)GOTO100 MFIT 29 DO90M=1,MM MFIT 30 DHAT(M)=DIST(M) MFIT 31 90 LBLOCK(M)=1 MFIT 32 GOTO500 MFIT 33 100 MA=1 MFIT 34 110 MAP=IPOINT(MA) MFIT 35 K=LBLOCK(MAP) MFIT 36 MB=MA+K-1 MFIT 37 GO TO ( 200, 300 ), LFITSW MFIT 38 C PRIMARY APPROACH MFIT 39 200 CONTINUE MFIT 40 IF(K-1) 9999, 220, 210 MFIT 41 C IF K=1, SAVE SORTING TIME MFIT 42 210 CALL SORT(IPOINT(MA), K,DIST( 1),+1) MFIT 43 C 'SORT' WILL SORT THE K ELEMENTS OF 'DIST' IN ALGEBRAIC MFIT 44 C ORDER. MFIT 45 220 DO 230 M=MA,MB MFIT 46 MP=IPOINT(M) MFIT 47 DHAT(MP)=DIST(MP) MFIT 48 LBLOCK(MP)=1 MFIT 49 230 CONTINUE MFIT 50 GO TO 400 MFIT 51 C SECONDARY APPROACH MFIT 52 300 CONTINUE MFIT 53 TEMP1=0.0 MFIT 54 DO 310 M=MA,MB MFIT 55 MP=IPOINT(M) MFIT 56 TEMP1=TEMP1+DIST(MP) MFIT 57 310 CONTINUE MFIT 58 MAP=IPOINT(MA) MFIT 59 DHAT(MAP)=TEMP1 MFIT 60 C PROCEED TO NEXT BOCK. QUERY. IS IT THE LAST MFIT 61 400 MA=MA+K MFIT 62 IF(MA-MM-1) 110, 500, 9999 MFIT 63 C START MAIN COMPUTATION *************************************MFIT 64 500 MA=1 MFIT 65 510 LUD=2 MFIT 66 MAP=IPOINT(MA) MFIT 67 NSATIS=0 MFIT 68 520 K=LBLOCK(MAP) MFIT 69 MB=MA+K-1 MFIT 70 MBP=IPOINT(MB) MFIT 71 WT=FLOAT(K) MFIT 72 550 DAV = DHAT(MAP)/ WT MFIT 73 GO TO ( 600, 700 ), LUD MFIT 74 C IS BLOCK DOWN-SATISFIED. IF NOT, MERGE MFIT 75 600 CONTINUE MFIT 76 IF(MA-1) 9999, 630, 610 MFIT 77 610 MBD=MA-1 MFIT 78 MBDP=IPOINT(MBD) MFIT 79 KD=LBLOCK(MBDP) MFIT 80 MAD=MBD-KD+1 MFIT 81 MADP=IPOINT(MAD) MFIT 82 WTD=FLOAT(KD) MFIT 83 617 DAVD = DHAT(MADP)/ WTD MFIT 84 IF( DAVD-DAV ) 630, 620, 620 MFIT 85 620 KNEW=K+KD MFIT 86 LBLOCK(MADP)=KNEW MFIT 87 LBLOCK(MBP)=KNEW MFIT 88 DTONEW = DHAT(MADP ) + DHAT(MAP) MFIT 89 DHAT(MADP)= DTONEW MFIT 90 NSATIS=0 MFIT 91 MA=MAD MFIT 92 MAP=MADP MFIT 93 MAP=MADP MFIT 94 GO TO 800 MFIT 95 630 NSATIS=NSATIS+1 MFIT 96 GO TO 800 MFIT 97 C IS BLOCK UP-SATISFIED. IF NOT, MERGE MFIT 98 700 CONTINUE MFIT 99 IF(MB-MM) 710, 730, 9999 MFIT 100 710 MAU=MB+1 MFIT 101 MAUP=IPOINT(MAU) MFIT 102 KU=LBLOCK(MAUP) MFIT 103 MBU=MAU+KU-1 MFIT 104 MBUP=IPOINT(MBU) MFIT 105 WTU=FLOAT(KU) MFIT 106 717 DAVU = DHAT(MAUP)/ WTU MFIT 107 IF( DAV-DAVU ) 730, 720, 720 MFIT 108 720 KNEW=K+KU MFIT 109 LBLOCK(MAP)=KNEW MFIT 110 LBLOCK(MBUP)=KNEW MFIT 111 DTONEW = DHAT(MAP)+ DHAT(MAUP) MFIT 112 DHAT(MAP)= DTONEW MFIT 113 NSATIS=0 MFIT 114 GO TO 800 MFIT 115 730 NSATIS=NSATIS+1 MFIT 116 GO TO 800 MFIT 117 C PROCEED TO NEXT BLOCK IF READY. MFIT 118 800 LUD = 3-LUD MFIT 119 C QUERY. IS BLOCK BOTH UP AND DOWN SATISFIED. IF NOT, RETURMFIT 120 IF(NSATIS-1) 520, 520, 810 MFIT 121 C QUERY. IS THIS LAST BLOCK. IF NOT, GO ON TO NEXT BLOCK. MFIT 122 810 CONTINUE MFIT 123 IF(MB-MM) 820, 900, 9999 MFIT 124 820 MA=MB+1 MFIT 125 GO TO 510 MFIT 126 C MAIN COMPUTATION COMPLETE. PLACE ANSWERS IN DHAT. MFIT 127 900 MA=1 MFIT 128 910 MAP =IPOINT(MA) MFIT 129 K=LBLOCK(MAP) MFIT 130 MB=MA+K-1 MFIT 131 IF(K-1) 9999, 940, 920 MFIT 132 920 TEMP1=DHAT(MAP)/FLOAT(K) MFIT 133 DO 930 M=MA,MB MFIT 134 MP=IPOINT(M) MFIT 135 930 DHAT(MP)=TEMP1 MFIT 136 GO TO 945 MFIT 137 940 DHAT(MAP) = DIST(MAP) MFIT 138 945 MA = MB + 1 MFIT 139 LKK=7 MFIT 140 IF(MA-MM-1) 910, 950, 9999 MFIT 141 950 RETURN MFIT 142 C TROUBLE EXIT MFIT 143 9999 WRITE(LPRINT,99) MFIT 144 99 FORMAT(58H0MFIT DIAGNOSTIC. IMPOSSIBLE BRANCH TAKEN ON IF STATEMENMFIT 145 1T.) MFIT 146 RETURN MFIT 147 END MFIT 148 SUBROUTINE OPTPUN OPUN 1 DIMENSION TITLE(20) OPUN 2 COMMON /TIL/ TITLE OPUN 3 COMMON /OPTEX1/ LREAD,LPRINT,LPUNCH,LSCRAT OPUN 4 WRITE(LPUNCH,19) TITLE OPUN 5 19 FORMAT(20A4/'CONFIGURATION') OPUN 6 RETURN OPUN 7 END OPUN 8 SUBROUTINE OPTPRI OPRI 1 DIMENSION TITLE(20) OPRI 2 COMMON /OPTEX1/ LREAD,LPRINT,LPUNCH,LSCRAT OPRI 3 COMMON /TIL/ TITLE OPRI 4 WRITE(LPRINT,13) TITLE OPRI 5 13 FORMAT(1X,20A4) OPRI 6 RETURN OPRI 7 END OPRI 8 SUBROUTINE SORT(A,N,B,SWITCH) SORT 1 C SORT 2 C SORT ROUTINE MANIPULATES POINTERS IN "A" WITHOUT MOVING SORT 3 C VALUES IN "B" SORT 4 C THE N VALUES FROM A(1) TO A(N) ARE SORTED IN ALGEBRAIC ORDER. SORT 5 C IF SWITCH IS +, ASCENDING ORDER. SORT 6 C IF SWITCH IS -, DESCENDING ORDER. SORT 7 C **** NOTE THIS IS NOT THE SLOW SORT REFERRED TO IN THE WRITE UP BYSORT 8 C KRUSKAL AND CARMONE,BUT A MODIFIED VERSION OF THE SORT 9 C D. L. SHELL ALGORITHM. SORT 10 C SORT 11 DIMENSION A(1),B(1) SORT 12 INTEGER SWITCH,A,TEMP SORT 13 IF(N.LE.1) GO TO 999 SORT 14 M=1 SORT 15 106 M=M+M SORT 16 IF(M.LE.N) GO TO 106 SORT 17 M=M-1 SORT 18 994 M = M/2 SORT 19 IF(M.EQ.0) GO TO 999 SORT 20 KK = N-M SORT 21 J = 1 SORT 22 992 I = J SORT 23 II=A(I) SORT 24 996 IM = I + M SORT 25 IMM=A(IM) SORT 26 IF(SWITCH) 810,810,800 SORT 27 800 IF(B(II).GT.B(IMM)) GO TO 110 SORT 28 GO TO 995 SORT 29 810 IF(B(II).LT.B(IMM)) GO TO 110 SORT 30 995 J = J+1 SORT 31 IF(J.GT.KK) GO TO 994 SORT 32 GO TO 992 SORT 33 110 TEMP=A(I ) SORT 34 A(I )=A(IM ) SORT 35 A(IM )=TEMP SORT 36 140 I = I-M SORT 37 IF(I.LT.1) GO TO 995 SORT 38 II=A(I) SORT 39 GO TO 996 SORT 40 999 CONTINUE SORT 41 RETURN SORT 42 END SORT 43 C BDAT 1 BLOCK DATA BDAT 2 INTEGER LPAR( 25 ) BDAT 3 REAL PAR( 25 ) BDAT 4 EQUIVALENCE (LPAR,PAR) BDAT 5 REAL*4 ATAB(4,22) BDAT 6 INTEGER*4 TAB(4,22),MTAB(4) BDAT 7 EQUIVALENCE (TAB,ATAB) BDAT 8 COMMON /OPTEX1/ LREAD, LPRINT, LPUNCH, LSCRAT BDAT 9 COMMON /OPTEX2/ LPAR BDAT 10 COMMON /OPTEX3/ MTAB, TAB BDAT 11 DATA LREAD, LPRINT, LPUNCH, LSCRAT /5, 6,7,2/ BDAT 12 DATA MTAB /1,18,20,22/ BDAT 13 DATA TAB/4HSTRM,2,2,0,'SRAT',2,1,0,'ITER',1,3,0,'COSA',2,5,0, BDAT 14 1'ACSA',2,6,0, 'SRTA',2,4,0,'DATA',3,7,2,'CONF',3,7,3,'PARA',3,7,4,BDAT 15 2 'COMP',3,7,5,'STOP',3,7,6, 'TITL',3,7,7, 'CUTO',2,8,0,'TIES', BDAT 16 35,1,2,'CARD',3,10,1,'NOCA',3,10,2,'PRIN',5,1,3, BDAT 17 4 'PRIM',3,9,1,'SECO',3,9,2,'NOIN',3,11,1,'INPU',3,11,2,0,0,0,0/BDAT 18 END BDAT 19 C CALC 1 SUBROUTINE CALC( LLL, AAA, BBB, CCC ) CALC 2 C FOR CALC NONMETRIC ANALYSIS OF FACTORIAL DESIGCALC 3 C DIMENSION DATA(501), INDEX(500, 4), X(500), XHAT(500) CALC 4 DIMENSION DATA(2001), INDEX(2000,5), X(2000), XHAT(2000) C DIMENSION IPOINT(501),LTIES(501) CALC 5 DIMENSION IPOINT(2001), LTIES(2001) C DIMENSION CON(4,100), GR(4,100), GL(4,100), GR2(4,100) CALC 6 DIMENSION CON(20,10), GR(20,10), GL(20,10), GR2(20,10) C DIMENSION II(4), IT(4) CALC 7 DIMENSION II(20), IT(20) DIMENSION LPAR(25),FMTCON(20),FMTDAT(20) CALC 8 COMMON /OPTEX1/ LREAD,LPRINT,LPUNCH,LSCRAT CALC 9 COMMON /OPTEX2/ LPAR CALC 10 EQUIVALENCE(PARVAL,LPARVA) CALC 11 EQUIVALENCE (CUTOFF,LPAR(8)),(MFITCO,LPAR(9)),(LPUNSW,LPAR(10)) CALC 12 EQUIVALENCE (LPRINS,LPAR(11)),(LCONSW,LPAR(12)) CALC 13 GO TO ( 1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000, CALC 14 1 1, 11000 ), LLL CALC 150000 C ********************************************************************CALC 16 C READ CONFIGURATION CALC 17 1000 READ(LREAD,100) KK,(II(K),K=1,KK) CALC 18 WRITE(LPRINT,100) KK,(II(K),K=1,KK) CALC 19 100 FORMAT(24I3) CALC 20 READ(LREAD,101) FMTCON CALC 21 WRITE(LPRINT,105) FMTCON CALC 22 101 FORMAT(20A4) CALC 23 105 FORMAT(2X,20A4) CALC 24 DO 1010 K=1,KK CALC 25 IIK=II(K) CALC 26 1010 READ(LREAD,FMTCON)(CON(K,I),I=1,IIK) CALC 27 RETURN CALC 28 C ********************************************************************CALC 29 C CREATE INITIAL CONFIGURATION CALC 30 2000 DO 2005 K=1,KK CALC 31 IIK=II(K) CALC 32 DO 2005 I=1,IIK CALC 33 CON(K,I)=0 CALC 34 GR (K,I)=0 CALC 35 GL(K,I)=0 CALC 36 GR2(K,I)=0 CALC 37 2005 CONTINUE CALC 38 DO 2010 L=1,LDATA CALC 39 M=IPOINT(L) CALC 40 TEMP=DATA(M) CALC 41 DO 2010 K=1,KK CALC 42 I=INDEX(M,K) CALC 43 CON(K,I)=CON(K,I)+TEMP CALC 44 2010 GR(K,I)=GR(K,I)+1.0 CALC 45 DO 2020 K=1,KK CALC 46 IIK=II(K) CALC 47 DO 2020 I=1,IIK CALC 48 CON(K,I)=CON(K,I)/GR(K,I) CALC 49 GR (K,I)=0 CALC 50 2020 CONTINUE CALC 51 RETURN CALC 52 C ********************************************************************CALC 53 C READ DATA CALC 54 3000 READ(LREAD,100) KK,(II(K),K=1,KK),LREPLI CALC 55 WRITE(LPRINT,100) KK,(II(K),K=1,KK),LREPLI CALC 56 READ(LREAD,101) FMTDAT CALC 57 WRITE(LPRINT,105) FMTDAT CALC 58 IIS=0 CALC 59 IIP=1 CALC 60 DO 3010 K=1,KK CALC 61 IIS=IIS+II(K) CALC 62 IIP=IIP*II(K) CALC 63 3010 IT(K)=1 CALC 64 NOD=IIP*LREPLI CALC 65 3020 READ(LREAD,FMTDAT) (DATA(L), L=1,NOD ) CALC 66 L=0 CALC 67 LR=0 CALC 68 DATMAX=-1.0E+20 CALC 69 DATMIN= -DATMAX CALC 70 DO 3090 LA=1,NOD CALC 71 TDATA= DATA(LA) CALC 72 IF( TDATA - CUTOFF ) 3045, 3045, 3030 CALC 73 3030 L=L+1 CALC 74 DATMAX =AMAX1 ( DATMAX, TDATA ) CALC 75 DATMIN =AMIN1 (DATMIN, TDATA ) CALC 76 DATA(L) = TDATA CALC 77 DO 3040 K=1,KK CALC 78 3040 INDEX(L,K)=IT(K) CALC 79 3045 LR=LR+1 CALC 80 IF( LR-LREPLI ) 3090, 3050, 9993 CALC 81 3050 LR=0 CALC 82 IT(KK)=IT(KK)+1 CALC 83 DO 3080 KC=1,KK CALC 84 K=KK+1-KC CALC 85 IF( IT(K)-II(K)-1 ) 3090, 3060, 9993 CALC 86 3060 IT(K)=1 CALC 87 IF(K-1) 9993, 3100, 3070 CALC 88 3070 IT(K-1)=IT(K-1)+1 CALC 89 3080 CONTINUE CALC 90 3090 CONTINUE CALC 91 3100 LDATA=L CALC 92 LDATAP = LDATA + 1 CALC 93 DATA(LDATAP) = 1.0E+20 CALC 94 DO 3110 L=1,LDATAP CALC 95 3110 IPOINT(L) = L CALC 96 CALL SORT(IPOINT(1),LDATA,DATA(1),+1) CALC 97 RETURN CALC 98 9993 CALL OPTERR( 6HMOD1RD ) CALC 99 CALL EXIT CALC 100 C ********************************************************************CALC 101 C CALCULATE STRESS AND GRADIENT CALC 102 4000 XBAR=0.0 CALC 103 DO 4010 L=1,LDATA CALC 104 TEMP=0.0 CALC 105 M=IPOINT(L) CALC 106 DO 4005 K=1,KK CALC 107 I=INDEX(M,K) CALC 108 4005 TEMP=TEMP+CON(K,I) CALC 109 X(M)=TEMP CALC 110 4010 XBAR=XBAR+TEMP CALC 111 XBAR=XBAR/FLOAT (LDATA) CALC 112 K=1 CALC 113 L1=1 CALC 114 M=IPOINT(1) CALC 115 TEMP2=DATA(M) CALC 116 DO 50 L=1,LDATA CALC 117 TEMP1=TEMP2 CALC 118 M=IPOINT(L+1) CALC 119 TEMP2=DATA(M) CALC 120 IF(TEMP1-TEMP2) 30, 20, 30 CALC 121 20 K=K+1 CALC 122 GO TO 50 CALC 123 30 L2=L1+K-1 CALC 124 DO 40 LA=L1,L2 CALC 125 M=IPOINT(LA) CALC 126 LTIES(M)=K CALC 127 40 CONTINUE CALC 128 L1=L2+1 CALC 129 K=1 CALC 130 50 CONTINUE CALC 131 CALL MFIT(X(1),XHAT(1),LTIES(1),LDATA,MFITCO,IPOINT) CALC 132 U=0 CALC 133 T=0 CALC 134 DO 4030 K=1,KK CALC 135 IIK=II(K) CALC 136 DO 4030 I=1,IIK CALC 137 GR(K,I)=0 CALC 138 4030 GR2(K,I)=0 CALC 139 DO 4050 L=1,LDATA CALC 140 M=IPOINT(L) CALC 141 TEMP1=X(M)-XHAT(M) CALC 142 TEMP2=X(M)-XBAR CALC 143 U=U+TEMP1**2 CALC 144 T=T+TEMP2**2 CALC 145 DO 4040 K=1,KK CALC 146 I=INDEX(M,K) CALC 147 GR(K,I)=GR(K,I)+TEMP1 CALC 148 4040 GR2(K,I)=GR2(K,I)+TEMP2 CALC 149 4050 CONTINUE CALC 150 U=SQRT (U) CALC 151 T=SQRT (T) CALC 152 STRESS=U/T CALC 153 GRGR=0 CALC 154 GRGL=0 CALC 155 IF(U) 9994, 4080, 4060 CALC 156 4060 COEF1=1.0/(U*T) CALC 157 COEF2=-U/T**3 CALC 158 DO 4070 K=1,KK CALC 159 IIK=II(K) CALC 160 DO 4070 I=1,IIK CALC 161 GR(K,I)=-COEF1*GR(K,I) - COEF2*GR2(K,I) CALC 162 GRGR=GRGR+GR(K,I)**2 CALC 163 4070 GRGL=GRGL+GR(K,I)*GL(K,I) CALC 164 4080 CONTINUE CALC 165 AAA=STRESS CALC 166 BBB=GRGR CALC 167 CCC=GRGL CALC 168 RETURN CALC 169 9994 CALL OPTERR( 6HMOD1CA ) CALC 170 RETURN CALC 171 C ********************************************************************CALC 172 C MOVE TO NEW CONFIGURATION CALC 173 5000 GRMULT=AAA CALC 174 DO 5010 K=1,KK CALC 175 IIK=II(K) CALC 176 DO 5010 I=1,IIK CALC 177 CON(K,I)=CON(K,I)+GRMULT*GR(K,I) CALC 178 GL(K,I)=GR(K,I) CALC 179 GR (K,I)=0 CALC 180 5010 GR2(K,I)=0 CALC 181 RETURN CALC 182 C ********************************************************************CALC 183 C NORMALIZE CONFIGURATION CALC 184 6000 SMULT=0 CALC 185 IF(LPRINS.NE.2) GO TO 6040 CALC 186 IF(CCC.LT.0.0) GO TO 6040 CALC 187 WRITE(LPRINT,602) CALC 188 602 FORMAT(1H1,' SEQ. NO. DATA SUBSCRIPTS'/) CALC 189 DO 6050 L=1,LDATA CALC 190 WRITE(LPRINT,601) L,DATA(L),(INDEX(L,K),K=1,KK) CALC 191 601 FORMAT(I6,4X,F10.5,20I5) 6050 CONTINUE CALC 193 IF(LCONSW.EQ.1) GO TO 6040 CALC 194 WRITE(LPRINT,603) CALC 195 603 FORMAT(///5X,' CONFIGURATION'/) CALC 196 DO 8011 K=1,KK CALC 197 IIK=II(K) CALC 198 WRITE(LPRINT,8005) IIK,(CON(K,I),I=1,IIK) CALC 199 8011 CONTINUE CALC 200 6040 DO 6020 K=1,KK CALC 201 AV=0 CALC 202 IIK=II(K) CALC 203 DO 6005 I=1,IIK CALC 204 GL(K,I)=0 CALC 205 6005 AV=AV+CON(K,I) CALC 206 AV=AV/FLOAT (IIK) CALC 207 DO 6010 I=1,IIK CALC 208 CON(K,I)=CON(K,I)-AV CALC 209 6010 SMULT=SMULT+CON(K,I)**2 CALC 210 6020 CONTINUE CALC 211 SMULT=SQRT ( FLOAT (IIS)/SMULT) CALC 212 DO 6030 K=1,KK CALC 213 IIK=II(K) CALC 214 DO 6030 I=1,IIK CALC 215 6030 CON(K,I)=SMULT*CON(K,I) CALC 216 RETURN CALC 217 C ********************************************************************CALC 218 C READ PARAMETERS CALC 219 7000 READ(LREAD,103) CUTOFF,MFITCO CALC 220 103 FORMAT(2F10.4) CALC 221 RETURN CALC 222 C ********************************************************************CALC 223 C PRINT, PUNCH, AND DRAW OUTPUT CALC 224 8000 CALL OPTPRI CALC 225 DO 8010 K=1,KK CALC 226 IIK=II(K) CALC 227 WRITE(LPRINT, 8005) IIK,(CON(K,I),I=1,IIK) CALC 228 8005 FORMAT(1X,I3,10F7.3/ (4X,10F7.3) ) CALC 229 8010 CONTINUE CALC 230 WRITE(LPRINT,8001) CALC 231 8001 FORMAT(1H1) CALC 232 CALL PLOTR(DATA,XHAT,X,LDATA,6,+1) CALC 233 IF(LPUNSW.EQ.2) RETURN CALC 234 CALL OPTPUN CALC 235 WRITE(LPUNCH, 106 ) CALC 236 106 FORMAT(2X, ' (2X,10F7.3) ') CALC 237 WRITE(LPUNCH,8020) KK,(II(K),K=1,KK) CALC 238 8020 FORMAT(20I3) CALC 239 DO 8040 K=1,KK CALC 240 IIK=II(K) CALC 241 WRITE(LPUNCH,8030) (CON(K,I),I=1,IIK) CALC 242 8030 FORMAT(2X,10F7.3) CALC 243 8040 CONTINUE CALC 244 RETURN CALC 245 C ********************************************************************CALC 246 C INITIAL ACTIONS CALC 247 9000 CUTOFF= -1.23E+20 CALC 248 MFITCO=1 CALC 249 LPUNSW=2 CALC 250 RETURN CALC 251 C ********************************************************************CALC 252 1 CONTINUE CALC 2530000 RETURN CALC 254 C ********************************************************************CALC 255 C RECEIVE LOCAL PARAMETER FROM CONTROL CARD CALC 256 11000 CONTINUE CALC 257 C CALL OPTPAR(LPARAM,PARVAL) CALC 258 LPARAM = 3 CALCA259 GO TO (11010, 11020, 11030), LPARAM CALC 259 11010 CUTOFF=PARVAL CALC 260 RETURN CALC 261 11020 MFITCO=LPARVA CALC 262 RETURN CALC 263 11030 GO TO 9981 CALC 264 9981 CALL OPTERR( 6HMOD1SP ) CALC 265 CALL EXIT CALC 266 END CALC 267 C PLOT 1 SUBROUTINE PLOTR(X,Y,YP,NPOI,OUT,ID) PLOT 2 C PLOT 3 C ROUTINE TO GENERATE A ONE PAGE PLOT OF ARRAY -X- VS. -Y-. PLOT 4 C PLOT 5 C THE PARAMETERS -XMAX- -XMIN- YMAX- - YMIN- INDICATE THE UPPER PLOT 6 C AND LOWER BOUNDS FOR EACH AXIS. IF XMAX=XMIN THE ROUTINE GENERATES PLOT 7 C ITS OWN BOUNDS FOR THE X AXIS, AND SIMILARLY IF YMAX=YMIN. PLOT 8 C PLOT 9 C IT IS ASSUMED THAT THE ENTRIES HAVE -NPOI- ENTRIES PLOT 10 C PLOT 11 C THE PLOTTING IS DONE ON TAPE -OUT-. PLOT 12 C PLOT 13 C IF -ID- IS NEGATIVE, AXES WILL BE INCLUDED ON THE PLOT. PLOT 14 C IF -ID- IS POSITIVE, NO AXES WILL APPEAR. PLOT 15 C PLOT 16 C -LONG- IS THE LENGTH OF THE ARRAYS -X- AND -Y- IN THE CALLING PROGRAM.PLOT 17 C PLOT 18 C WRITTEN BY FORREST W. YOUNG, NOVEMBER, 1965 PLOT 19 C VERSION 3, APRIL 1967 PLOT 20 C PLOT 21 C--------ADAPTED BY F. J. CARMONE,JR. 5/1/67 PLOT 22 C PLOT 23 C PLOT 24 DIMENSION X(NPOI),Y(NPOI),ITEM(55,101),SMALL(21),HOLL(11),YP(NPOI)PLOT 25 REAL LABELS(55) PLOT 26 REAL ITEM PLOT 27 INTEGER OUT PLOT 28 DATA HOLL,DD,PLUS,BLANK/1H ,1HX,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9, PLOT 29 11HM,1H.,1H*,1H /,AIE,AMINUS/1HI,1H-/ PLOT 30 DATA LABELS/15*1H ,'T','H','E',' ','X',' ','A','R','E',' ', PLOT 31 1'T','H','E',' ','L','I','N','E','A','R',' ','M','O','D','E','L', PLOT 32 1 14*1H / PLOT 33 DATA AAB/1HX/,DOO/'0'/,PL/'+'/ PLOT 34 3001 FORMAT(14X,103H*.****.****.****.****.****.****.****.****.****.****PLOT 35 1.****.****.****.****.****.****.****.****.****.****.*) PLOT 36 3300 FORMAT(' ',A1,F9.2,1X,105A1,F11.2) PLOT 37 3301 FORMAT(15X,A1,10(F9.4,A1)) PLOT 38 3302 FORMAT(10X,11F10.4) PLOT 39 DO 115 I=1,55 PLOT 40 DO 115 J=1,101 PLOT 41 115 ITEM(I,J)=BLANK PLOT 42 116 CONTINUE PLOT 43 IF(ID.GT.0)GO TO 119 PLOT 44 DO 117 I=1,55 PLOT 45 117 ITEM(I,51)=AIE PLOT 46 DO 118 I=1,101 PLOT 47 118 ITEM(28,I)=AMINUS PLOT 48 119 CONTINUE PLOT 49 XMAX=X(1) PLOT 50 XMIN =X(2) PLOT 51 DO 121 I=1,NPOI PLOT 52 IF(X(I).GT.XMAX)XMAX=X(I) PLOT 53 IF(X(I).LT.XMIN)XMIN=X(I) PLOT 54 121 CONTINUE PLOT 55 YMAX=Y(1) PLOT 56 YMIN=Y(2) PLOT 57 DO 123 I=1,NPOI PLOT 58 IF(Y(I).GT.YMAX)YMAX=Y(I) PLOT 59 IF(YP(I).GT.YMAX)YMAX=YP(I) PLOT 60 IF(Y(I).LT.YMIN)YMIN=Y(I) PLOT 61 IF(YP(I).LT.YMIN)YMIN=YP(I) PLOT 62 123 CONTINUE PLOT 63 SPANX=XMAX-XMIN PLOT 64 SPANY=YMAX-YMIN PLOT 65 XMAX=XMAX+.05*SPANX PLOT 66 YMAX=YMAX+.05*SPANY PLOT 67 XMIN=XMIN-.05*SPANX PLOT A67 YMIN=YMIN-.05*SPANY PLOT B67 SPANX=XMAX-XMIN PLOT 68 SPANY=YMAX-YMIN PLOT 69 DELX=SPANX/100.0 PLOT 70 DELY=SPANY/54.0 PLOT 71 VALUE=YMAX+DELY PLOT 72 SMALL(1)=XMIN PLOT 73 DO 120 I=1,20 PLOT 74 120 SMALL(I+1)=SMALL(I)+5.0*DELX PLOT 75 DO 135 II=1,NPOI PLOT 76 I=(YMAX-Y(II))/DELY+1.5 PLOT 77 IP=(YMAX-YP(II))/DELY+1.5 PLOT 78 J=(X(II)-XMIN)/DELX+1.5 PLOT 79 IF(I.GT.56.OR.I.LT.1.OR.J.LT.1.OR.J.GT.101)GO TO 135 PLOT 80 IF(IP.GT.56.OR.IP.LT.1)GO TO 135 PLOT 81 201 IF(ITEM(IP,J).NE.BLANK) GO TO 202 PLOT 86 ITEM(IP,J)=AAB PLOT 87 GO TO 203 202 ITEM(IP,J)=PL PLOT 89 203 CONTINUE ITEM( I,J)=DOO PLOT 83 135 CONTINUE PLOT 90 WRITE(OUT,3301)DD,(SMALL(I),DD,I=2,20,2) PLOT 91 WRITE(OUT,3302)(SMALL(I),I=1,21,2) PLOT 92 WRITE(OUT,3001) PLOT 93 DO 330 I=1,55 PLOT 94 VALUE=VALUE-DELY PLOT 95 A=BLANK PLOT 96 B=PLUS PLOT 97 L=I+2 PLOT 98 IF(L/3*3-L)330,310,330 PLOT 99 310 B=DD PLOT 100 IF(L/2*2-L)330,320,330 PLOT 101 320 A=DD PLOT 102 330 WRITE(OUT,3300)LABELS(I),VALUE,A,B,(ITEM(I,J),J=1,101),B,A,VALUE PLOT 103 WRITE(OUT,3001) PLOT 104 WRITE(OUT,3301)DD,(SMALL(I),DD,I=2,20,2) PLOT 105 WRITE(OUT,3302)(SMALL(I),I=1,21,2) PLOT 106 CONTINUE PLOT 107 RETURN PLOT 108 END PLOT 109