Trailing-Edge
-
PDP-10 Archives
-
decuslib20-04
-
decus/20-0101/nln20c.pgm
There are 2 other files named nln20c.pgm in the archive. Click here to see a list.
C --- NONLINWOOD NONLINEAR LEAST-SQUARES CURVE-FITTING PROGRAM 01CT0010
C 01CT0020
C --- 1/79 01CT0030
C 01CT0040
C *** PROGRAM CURRENTLY DIMENSIONED FOR 20 MAX VARIABLES, 20 MAX 01CT0050
C --- COEFFICIENTS, 170 OBSERVATIONS AND 2 DELETE OBSERVATIONS CARDS. 01CT0060
C 01CT0070
C GAUSHAUS PROGRAM WAS ORIGINALLY WRITTEN BY D. A. MEETER, 01CT0080
C UNIV. OF WISCONSIN, USING D. W. MARQUART'S MAXIMUM NEIGHBORHOOD 01CT0090
C METHOD. PROGRAM WAS REVISED AT STANDARD OIL (INDIANA) BY 01CT0100
C F. S. WOOD TO PROVIDE NOMENCLATURE, PRINTOUT AND PLOTS CONSISTENT 01CT0110
C WITH THE LINEAR LEAST-SQUARES CURVE FITTING PROGRAM. IT WAS 01CT0120
C MODIFIED TO ALLOW FOR OVERLAYING BY R. E. HENSCHEL. THE PROGRAM 01CT0130
C IS IN DOUBLE PRECISION FOR IBM 360 AND 370 COMPUTERS (SHARE 01CT0140
C LIBRARY PROGRAM NUMBER 360-13.6.007). 01CT0150
C 01CT0160
C PROGRAM CONVERTED AT UNIV. CAL(BERKELEY) BY N. H. TIMM TO 01CT0170
C SINGLE PRECISION FOR CDC 6000 SERIES COMPUTERS (VIM LIBRARY 01CT0180
C PROGRAM NUMBER G2-CAL-NLWOOD). UPDATED BY TIMM, WOOD, AND 01CT0190
C B. E. FOSTER AND ELI COHEN, NORTHWESTERN UNIVERSITY. 01CT0191
C 01CT0200
C PROGRAM CONVERTED AT ALCOA R&D LABORATORIES BY R. F. KOHM 01CT0210
C FOR DECSYSTEM-10 AND 20 COMPUTERS (DECUS LIBRARY PROGRAM NUMBER 01CT0220
C NONLIN-10-258). PROGRAM ADAPTED FOR INTERACTIVE USE BY 01CT0221
C E. R. ZIEGEL AT STANDARD OIL (2ND). 01CT0222
C 01CT0230
C PROGRAM CONVERTED AT UNIV. WISCONSIN BY J. D. MURAT AND 01CT0240
C T. R. ZEISLER FOR BURROUGHS 4700 AND 6700 SYSTEM COMPUTERS (CUBE 01CT0250
C LIBRARY PROGRAM NUMBER WIS/NONLINWOOD). THE PROJECT WAS 01CT0260
C INITIATED BY D. S. RUMSEY OF BURROUGHS. 01CT0261
C 01CT0270
C PROGRAM CONVERTED AT U.S. NAVAL AVIONICS FACILITY, 01CT0271
C INDIANAPOLIS, IN. BY D. F. ZARNOW FOR HONEYWELL 6000/600/66 01CT0272
C COMPUTERS (HLSUA LIBRARY PROGRAM NUMBER GES-1207). 01CT0273
C 01CT0274
C PROGRAM CONVERTED AT JOHNSON SPACE CENTER BY J. E. KEITH 01CT0280
C FOR UNIVAC 1108 AND 1110 COMPUTERS. 01CT0290
C 01CT0300
C 01CT0310
C 01CT0320
C FOR GLOSSARY OF TERMS, DEFINITIONS OF STATISTICS AND 01CT0330
C INTERPRETATION OF RESULTS SEE USER'S MANUAL IN BOOK "FITTING 01CT0340
C EQUATIONS TO DATA" BY C. DANIEL AND F. S. WOOD WITH ASSISTANCE 01CT0350
C OF J. W. GORMAN, WILEY, 2ND EDITION, 1978. 01CT0360
C 01CT0370
C 01CT0380
C 01CT0390
C --- AUXILIARY EQUIPMENT NEEDED TO RUN PROGRAM 01CT0400
C --- CARD READER (UNIT 5), AND PRINTER (UNIT 6). 01CT0410
C 01CT0420
C 01CT0430
C 01CT0440
C 01CT0450
C 01CT0460
C --- USER'S INSTRUCTIONS 01CT0470
C --- 01CT0480
C --- ORDER OF CARDS FOR EACH PROBLEM 01CT0490
C --- 1 CONTROL CARD. 01CT0500
C --- 2 FORMAT CARD (72COL). E.G. (A6,F6.0,(NOIND*F6.0)) TO READ DATA. 01CT0510
C --- 3 DELETE-OBSERVATIONS CARD(S), IF ANY 01CT0511
C --- (1ST 6 OF 10 COLS. / OBSERVATION DELETED, 7 / CARD). 01CT0512
C --- 4 STARTING VALUES (GUESSES) OF COEFFICIENTS (10COL/COEF, 7/CARD). 01CT0520
C --- 5 INFORMATION CARDS FOR PRINTOUT,IF ANY(72COL/CARD, 12 CARDS MAX).01CT0530
C --- 6 NAMES OF COEFFICIENTS CARD(S) FOR PRINTOUT, IF ANY 01CT0540
C --- (1ST 6 OF 10 COLS. / COEFFICIENT, 7 / CARD). 01CT0550
C --- 7 DATA CARDS (IF NOT READ FROM A FILE). 01CT0560
C --- 8 END CARD (END IN FIRST 3 POSITIONS OF IDENTIFICATION). 01CT0570
C --- NOTE: FORMAT, DELETE, DATA AND END CARDS ARE NOT NEEDED IN 01CT0580
C --- SUBSEQUENT PROBLEMS IF SAME INPUT DATA ARE REUSED. 01CT0590
C --- 01CT0600
C --- CONTROL CARD 01CT0610
C --- COL. 1-20 NPROB :IDENTIFICATION OF PROBLEM, 5A4. 01CT0620
C --- 21 JSAME :0 OR BLANK, OBSERVATIONS READ FROM CARDS. 01CT0630
C --- 1, REUSE DATA FROM PREVIOUS PROBLEM. 01CT0640
C --- 23-24 NC :NO. OF COEFFICIENTS TO BE ESTIMATED, I2. 01CT0650
C --- 25-26 NFILE :FILE NUMBER IF DATA ARE TO BE READ FROM 01CT0660
C --- SEPARATE FILE, NO END CARD(S) IF ONLY ONE 01CT0670
C --- SET OF DATA IS ON EACH FILE. 01CT0680
C --- 27-28 NOIND :NO. OF INDEPENDENT VARIABLES TO BE READ IN, I2 01CT0690
C --- 31-32 MODEL :NO. OF THE EQUATION TO BE USED, I2. 01CT0700
C --- 01CT0710
C --- 33-36 FLAM :STARTING VALUE FOR LAMBDA, F4.2, (E.G. 0.1), 01CT0720
C --- USED AS A MULTIPLIER TO SCALE THE SPACE OR 01CT0730
C --- SIZE STEPS TAKEN. 01CT0740
C --- 37-40 FNU :VALUE OF NU, F4.0, (E.G. 10.), 01CT0750
C --- DIVISOR AND MULTIPLIER TO CHANGE SIZE OF 01CT0760
C --- LAMBDA DEPENDING ON WHETHER SUM OF SQUARES OF 01CT0770
C --- ITERATION IS NEAR OR FAR FROM MINIMUM. 01CT0780
C --- 43-44 MIT :MAX NUMBER OF ITERATIONS, I2, (E.G. 20). 01CT0790
C --- 45-48 DIFF :MULTIPLIER USED TO INCREMENT VALUE OF 01CT0800
C --- COEFFICIENTS, F4.3 (E.G. 0.01). 01CT0810
C --- 01CT0811
C --- NOTE: IF FLAM, FNU, MIT, AND DIFF ARE NOT DEFINED 01CT0820
C --- ON THE CONTROL CARD, THEIR LEVEL WILL BE SET 01CT0830
C --- AUTOMATICALLY TO THE ABOVE E.G. VALUES. 01CT0840
C --- 01CT0850
C --- STOP :CRITERIA FOR ENDING CONVERGING ITERATIONS. 01CT0860
C --- 49-56 STOPSS:SUM OF SQUARES CRITERION, F8.7 (E.G. 0.0001, 01CT0870
C --- A CHANGE OF LESS THAN 0.0001 IN THE RESIDUAL 01CT0880
C --- SUM OF SQUARES). 01CT0890
C --- 57-64 STOPRC:RATIO OF COEFFICIENTS CRITERION, F8.7, 01CT0900
C --- (E.G. 0.001 A CHANGE OF LESS THAN 0.001 IN 01CT0910
C --- THE RATIOS OF ALL COMPARABLE COEFFICIENTS). 01CT0920
C --- 01CT0921
C --- NOTE: STOPSS AND STOPRC CAN BE SET AT 0.0 IF CONTROL 01CT0930
C --- OF EITHER OR BOTH IS NOT DESIRED. 01CT0940
C --- 01CT0950
C --- 65-66 NEQU :NO OF INFORMATION CARDS TO BE READ FOR DISPLAY 01CT0960
C --- ON PRINTOUT IF DESIRED, 72 COL.EACH, 12 MAX. 01CT0970
C --- 68 NAMEBI:1 READ NAMES OF COEFFICIENTS FROM CARDS FOR 01CT0980
C --- DISPLAY ON PRINTOUT, 1ST 6 OF 10 COLS. / 01CT0990
C --- COEFFICIENT, 7 / CARD. 01CT1000
C --- 69 NPLOT :3 PLOT RESIDUALS VS. EACH INDEPENDENT VARIABLE.01CT1010
C --- 70 NOMITC:NO OF DELETE OBSERVATION CARDS, OBSERVATION 01CT1011
C --- IDENTIFICATION IN 1ST 6 OF 10 COLS., 7 / CARD. 01CT1012
C -- 01CT1020
C --- FORMAT CARD 01CT1030
C --- COL. 1-72 E.G. (A6, F6.0, (NOIND*F6.0) 01CT1040
C --- IDENT IDENTIFICATION OF OBSERVATION, A6 01CT1050
C --- Y(J) DEPENDENT VARIABLE, J-TH OBSERVATION, F6.0 01CT1060
C --- X(I,J) INDEPENDENT VARIABLE, I-TH VARIABLE, J-TH 01CT1070
C --- OBSERVATION, NOIND*(F6.0) 01CT1080
C -- 01CT1090
C *******************************************************************01CT1100
C -- 01CT1110
C --- EQUATIONS TO BE USED WILL BE WRITTEN IN FORTRAN AND PLACED 01CT1120
C --- IN SUBROUTINE MODEL1, 2, 3, 4, OR 5 AS FOLLOWS: 01CT1130
C --- 01CT1140
C --- SUBROUTINE MODEL1 (NPROB, B, FY, NOB, NC, X, NVARX, NOBMAX, NCMAX)01CT1150
C IMPLICIT REAL*8(A-H,O-Z) 01CT1160
C --- DIMENSION B(NCMAX), FY(NOBMAX), X(NVARX,NOBMAX) 01CT1170
C 01CT1180
C --- DO 10 J=1,NOB 01CT1190
C --- EXAMPLE OF EQUATION 01CT1200
C --- FY(J)=B(1)*X(1,J)**B(2) + B(3) 01CT1210
C --- 01CT1220
C -10 CONTINUE 01CT1230
C --- RETURN 01CT1240
C --- END 01CT1250
C --- NOTE: ALL EQUATIONS ARE TO BE WRITTEN SO THAT THE 01CT1260
C --- COEFFICIENTS ESTIMATED AND CALCULATED WILL BE 01CT1270
C --- POSITIVE. IN SINGLE PRECISION PROGRAMS DELETE 01CT1280
C --- IMPLICIT REAL*(A-H,O-Z) STATEMENTS. 01CT1281
C 01CT1290
C *******************************************************************01CT1300
C * *01CT1310
C * TO PERFORM TRANSFORMATIONS, ADD SUBROUTINE BLOCK DATA AND *01CT1320
C * APPROPRIATE CARDS IN MODEL SUBROUTINE. E.G. *01CT1330
C * *01CT1340
C BLOCK DATA *01CT1350
C IMPLICIT REAL*8(A-H,O-Z) *01CT1360
C COMMON /DATA05/ IOPEN *01CT1370
C DATA IOPEN /0/ *01CT1380
C END *01CT1390
C * *01CT1400
C SUBROUTINE MODEL2 (NPROB, B, FY, NOB, NC, X, NVARX, NOBMAX, NCMAX)01CT1410
C *01CT1420
C IN THIS EXAMPLE THE VALUES OF Y ARE LOGGED AND THE RECIPROCAL *01CT1430
C OF X IS TAKEN INITIALLY AND NOT AGAIN. *01CT1440
C *01CT1450
C IMPLICIT REAL*8(A-H,O-Z) *01CT1460
C COMMON /DATA01/ IDENT(170), Y(170) *01CT1470
C COMMON /DATA05/ IOPEN *01CT1480
C DIMENSION B(NCMAX), FY(NOBMAX), X(NVARX,NOBMAX) *01CT1490
C REAL*8 IDENT *01CT1500
C *01CT1510
C IF (IOPEN .NE. 0) GO TO 7 *01CT1520
C DO 5 J = 1, NOB *01CT1530
C Y(J) = DLOG(Y(J)) *01CT1540
C 5 X(1,J) = (1./X(1,J)) *01CT1550
C *01CT1560
C 7 DO 10 J = 1, NOB *01CT1570
C FY(J) = B(1)*X(1,J)**B(2) + B(3) *01CT1580
C 10 CONTINUE *01CT1590
C *01CT1600
C IOPEN = 1 *01CT1610
C RETURN *01CT1620
C END *01CT1630
C * *01CT1640
C *******************************************************************01CT1650
C 01CT1660
C --- NVARX AND NOBMAX ARE USED TO DIMENSION SUBSEQUENT SUBROUTINES. 01CT1670
C --- THE MAX NUMBER OF VARIABLES(NVARX) AND THE MAX NUMBER OF 01CT1680
C --- OBSERVATIONS(NOBMAX) ARE DEPENDENT ON THE COMPUTER CORE SIZE 01CT1690
C --- AVAILABLE. HOWEVER, ONE CAN BE TRADED FOR THE OTHER BY THE 01CT1700
C --- FOLLOWING EQUATIONS: (UNITS IN WORDS) 01CT1710
C --- CORE AVAILABLE = (CORE LESS BUFFER - PROGRAM EX VARS. AND OBSVS.) 01CT1720
C --- E.G. 48,100CA = ( 61,300CLB - 13,200PROGRAM EX V AND O) 01CT1730
C --- NOBMAX = (CORE AVAILABLE - 8NVARX - 3(NVARX)SQRD) / (4 + 2NVARX) 01CT1740
C --- E.G.(50VAR,386NOB), (60VAR,296NOB), (70VAR,228NOB), (80VAR,172NOB)01CT1741
C --- 01CT1742
C --- DELZ MUST TOTAL 2756+NOBMAX TO PROVIDE SPACE FOR PLOTS IN PITCHA 01CT1743
C 01CT1744
C * * * * * * * * * * * * * * * * * * * * * * * 01CT1750
C * 01CT1751
C * SUGGESTED OVERLAY CARDS TO BE INCLUDED IN THE LINK STEP 01CT1752
C INSERT DATA01,DATA02,DATA04,MINMAX 01CT1753
C OVERLAY AAA 01CT1754
C INSERT READIN 01CT1755
C OVERLAY AAA 01CT1756
C INSERT GAUSHS,TTEST,EIGENJ,MATINV 01CT1757
C OVERLAY BBB 01CT1758
C INSERT MODEL1 01CT1759
C OVERLAY BBB 01CT1760
C INSERT MODEL2 01CT1761
C OVERLAY BBB 01CT1762
C INSERT MODEL3 01CT1763
C OVERLAY BBB 01CT1764
C INSERT MODEL4 01CT1765
C OVERLAY BBB 01CT1766
C INSERT MODEL5 01CT1767
C OVERLAY AAA 01CT1768
C INSERT SORT 01CT1769
C OVERLAY AAA 01CT1770
C INSERT PITCHA,GRID,PACK,FORCE 01CT1771
C ENTRY MAIN 01CT1772
C * 01CT1773
C * * * * * * * * * * * * * * * * * * * * * * * 01CT1774
C 01CT1780
IMPLICIT REAL*8 (A-H,O-Z) 01CT1790
COMMON /DATA01/ IDENT(170), Y(170), X(20,170), B(20), DIFZ(20), 01CT1800
1 SIGNS(20), F(170), R(170), DELZ(170,20) 01CT1810
C 01CT1820
COMMON /DATA02/ STOPSS, STOPCR, FLAM, FNU, MIT, NC, NEQU, 01CT1830
1 NOMIT, NFOUND 01CT1840
C 01CT1850
COMMON /DATA04/ EQU(216), BI( 40), NPROB(5) 01CT1860
REAL*8 IDENT, JOMIT(14) 01CT1870
EXTERNAL MODEL1, MODEL2, MODEL3, MODEL4, MODEL5 01CT1880
DIMENSION LSORT(1) 01CT1890
EQUIVALENCE (LSORT(1), DELZ(1,1)) 01CT1900
C 01CT1910
C *******************************************************************01CT1920
C * *01CT1930
C * TO CHANGE THE NUMBER OF VARIABLES, COEFFICIENTS AND/OR 01CT1940
C * OBSERVATIONS: 01CT1941
C * 1) REPLACE LABEL COMMON DATA01 AND DATA04 IN MAIN *01CT1950
C * 2) INITIALIZE NVARX, NCMAX AND NOBMAX WHERE: *01CT1960
C * NVARX = MAXIMUM NUMBER OF VARIABLES *01CT1970
C * NCMAX = MAXIMUM NUMBER OF COEFFICIENTS *01CT1971
C * NOBMAX= MAXIMUM NUMBER OF OBSERVATIONS *01CT1980
C * 3) REPLACE DIMENSION STATEMENT IN SUBROUTINE GAUSHS *01CT1990
C * IF NCMAX IS CHANGED. *01CT2000
C * *01CT2010
C * EXAMPLE: NUMBER OF VARIABLES = 10, COEFFICIENTS = 15 AND *01CT2020
C * OBSERVATIONS = 200 *01CT2021
C * *01CT2030
C * REPLACE DATA01 COMMON TO *01CT2040
C COMMON /DATA01/ IDENT(200), Y(200), X(10,200), B(15), DIFZ(15), *01CT2050
C 1 SIGNS(15), F(200), R(200), DELZ(200,15) *01CT2060
C * REPLACE DATA04 COMMON WITH DIMENSION OF BI = 2*NCMAX *01CT2061
C COMMON /DATA04/ EQU(216), BI( 30), NPROB(5) *01CT2062
C * *01CT2070
C * AND REPLACE *01CT2080
NVARX = 20 01CT2090
NCMAX = 20 01CT2091
NOBMAX = 170 01CT2100
C * WITH *01CT2110
C NVARX = 10 *01CT2120
C NCMAX = 15 *01CT2121
C NOBMAX = 200 *01CT2130
C * *01CT2140
C * REPLACE DIMENSION STATEMENT IN SUBROUTINE GAUSHS TO: *01CT2150
C DIMENSION E(15), P(15), PHI(15), Q(15), TB(15), A(15,15), *01CT2160
C 1 D(15,15), U(15,15) *01CT2170
C * *01CT2180
C *******************************************************************01CT2190
C * *01CT2191
C * TO CHANGE THE MAXIMUM NUMBER OF DELETE OBSERVATION CARDS *01CT2192
C * READ AND HENCE THE MAXIMUM NUMBER OF DELETE OBSERVATIONS *01CT2193
C * ( 7 PER CARD): CHANGE MNOMTC, MAXIMUM NUMBER OF DELETE CARDS *01CT2194
C * AND DIMENSION OF JOMIT(MNOMTC*7) IN MAIN. *01CT2195
C * *01CT2196
C *******************************************************************01CT2197
NNPROB = 5 01CT2200
C --- MXNEQU = MAXIMUM NUMBER OF INFORMATION (EQUATION) CARDS. 01CT2201
4 MXNEQU = 12 01CT2202
C --- NEQUNO = DIMENSION OF NEQU. 01CT2203
NEQUNO = MXNEQU*18 01CT2210
C --- MNOMTC = MAXIMUM NUMBER OF DELETE OBSERVATION CARDS. 01CT2211
MNOMTC = 2 01CT2212
C --- MNOMIT = MNOMTC*7 OBSERVATIONS/CARD, 7(A6,4X), IDENT OF 01CT2213
C --- OBSERVATIONS TO BE DELETED, THE DIMENSIONS OF JOMIT. 01CT2214
MNOMIT = MNOMTC*7 01CT2215
C --- NBINO = DIMENSIONS OF BI, THE NAMES OF THE COEFFICIENTS. 01CT2216
NBINO = 2*NCMAX 01CT2220
C --- KTIN = FILE TO READ CONTROL CARDS 01CT2221
KTIN = 5 01CT2222
C --- NFILE = FILE TO READ DATA 01CT2223
C --- NFILE = KTIN UNLESS SPECIFIED ON CONTROL CARD, COLUMN 25-26 01CT2224
C --- KTOU = FILE FOR PRINTOUT OF PLOTS 01CT2225
KTOU = 6 01CT2226
KTINED = 0 01CT2230
100 IERGAU = 0 01CT2240
C 01CT2250
C --- CALL SUBROUTINE READIN TO READ CONTROL CARDS AND DATA. 01CT2260
C 01CT2270
CALL OPEN 01CT2271
200 CALL READIN (KTOU,IDENT, Y, X, EQU, BI, NPROB, B, DIFZ, SIGNS, 01CT2280
1 MODEL, KTINED, KTIN, JOMIT, MNOMIT, NVARX, NOBMAX, NEQUNO, 01CT2290
2 NBINO, NNPROB, NOB, NOIND, NPLOT, NCMAX) 01CT2291
C 01CT2300
IF (KTINED .NE. 0) GO TO 9999 01CT2310
C 01CT2320
C --- DETERMINE WHICH MODEL AND CALL SUBROUTINE GAUSHS. 01CT2330
C 01CT2340
IF (MODEL .GT. 0 .AND. MODEL .LT. 6) GO TO 250 01CT2341
WRITE(KTOU, 240) MODEL 01CT2342
240 FORMAT(1H1,38H***CONTROL CARD ERROR*** MODEL NUMBER ,I2, 01CT2343
1 10H REQUESTED ) 01CT2344
GO TO 200 01CT2345
250 IF (MODEL .EQ. 1) CALL GAUSHS (IDENT, Y, X, EQU, BI, NPROB, B, 01CT2350
1 KTOU, DIFZ, SIGNS, F, R, DELZ, MODEL1, IERGAU, MNOMIT, JOMIT, 01CT2360
2 NVARX, NOBMAX, NEQUNO, NBINO, NNPROB, NOB, NOIND, NPLOT, NCMAX) 01CT2361
IF (MODEL .EQ. 2) CALL GAUSHS (IDENT, Y, X, EQU, BI, NPROB, B, 01CT2370
1 KTOU, DIFZ, SIGNS, F, R, DELZ, MODEL2, IERGAU, MNOMIT, JOMIT, 01CT2380
2 NVARX, NOBMAX, NEQUNO, NBINO, NNPROB, NOB, NOIND, NPLOT, NCMAX) 01CT2381
IF (MODEL .EQ. 3) CALL GAUSHS (IDENT, Y, X, EQU, BI, NPROB, B, 01CT2390
1 KTOU, DIFZ, SIGNS, F, R, DELZ, MODEL3, IERGAU, MNOMIT, JOMIT, 01CT2400
2 NVARX, NOBMAX, NEQUNO, NBINO, NNPROB, NOB, NOIND, NPLOT, NCMAX) 01CT2401
IF (MODEL .EQ. 4) CALL GAUSHS (IDENT, Y, X, EQU, BI, NPROB, B, 01CT2410
1 KTOU, DIFZ, SIGNS, F, R, DELZ, MODEL4, IERGAU, MNOMIT, JOMIT, 01CT2420
2 NVARX, NOBMAX, NEQUNO, NBINO, NNPROB, NOB, NOIND, NPLOT, NCMAX) 01CT2421
IF (MODEL .EQ. 5) CALL GAUSHS (IDENT, Y, X, EQU, BI, NPROB, B, 01CT2430
1 KTOU, DIFZ, SIGNS, F, R, DELZ, MODEL5, IERGAU, MNOMIT, JOMIT, 01CT2440
2 NVARX, NOBMAX, NEQUNO, NBINO, NNPROB, NOB, NOIND, NPLOT, NCMAX) 01CT2441
C 01CT2450
IF (IERGAU .NE. 0) GO TO 100 01CT2460
C 01CT2470
C --- CALL SUBROUTINE SORT TO PRINT AND SORT RESIDUALS IN 01CT2480
C --- CHRONOLOGICAL AND ASCENDING ORDER. 01CT2490
C 01CT2500
CALL SORT (IDENT, Y, NPROB, F, R, DELZ, YCMIN, YCMAX, NOBMAX, 01CT2510
1 NNPROB, NOB, KTOU) 01CT2511
C 01CT2520
C --- CALL SUBROUTINE PITCHA TO PLOT RESIDUALS. 01CT2530
C 01CT2540
CALL PITCHA (X, EQU, BI, F, R, DELZ, NPROB, YCMIN, YCMAX, 01CT2550
1 LSORT(NOBMAX+1), LSORT(NOBMAX+1379), KTOU, 01CT2560
2 NVARX, NOBMAX, NEQUNO, NBINO, NNPROB, NOB, NOIND, NPLOT) 01CT2561
C 01CT2570
C CALL OPTIONAL OUTPUT SUBROUTINE 01CT2571
C 01CT2572
CALL OUTPUT (IDENT,Y,X,F,B,NOB,NOIND,NC,NOBMAX,NVARX,NCMAX,KTOU) 01CT2573
C 01CT2580
TYPE 9000 01CT2581
9000 FORMAT(//1H0,' *** PRINT OUTPUT FILE FOR COMPLETE INFORMATION *** 01CT2582
1 '////) 01CT2583
WRITE(KTOU, 300) NPROB 01CT2584
300 FORMAT(1H1,18HEND OF PROBLEM ,5A4) 01CT2590
GO TO 200 01CT2600
C 01CT2610
9999 STOP 01CT2620
END 01CT2630
SUBROUTINE GAUSHS (IDENT, Y, X, EQU, BI, NPROB, TH, KTOU, DIFZ, 03GH0010
1 SIGNS, F, R, DELZ, MODEL, IERGAU, MNOMIT, JOMIT, 03GH0020
2 NVARX, NOBMAX, NEQUNO, NBINO, NNPROB, NOB, NOIND, NPLOT, NCMAX) 03GH0021
IMPLICIT REAL*8(A-H,O-Z) 03GH0030
COMMON /DATA02/ STOPSS, STOPCR, FLAM, FNU, MIT, NC, NEQU, 03GH0040
1 NOMIT, NFOUND 03GH0050
C 03GH0060
REAL*8 IDENT, JOMIT(MNOMIT) 03GH0070
DIMENSION IDENT(NOBMAX), Y(NOBMAX), X(NVARX,NOBMAX), EQU(NEQUNO), 03GH0080
1 BI(NBINO), NPROB(NNPROB), TH(NCMAX), DIFZ(NCMAX), 03GH0090
2 SIGNS(NCMAX), F(NOBMAX), R(NOBMAX), DELZ(NOBMAX,NCMAX) 03GH0100
C *******************************************************************03GH0110
C * *03GH0120
C * REPLACE THE FOLLOWING DIMENSION STATEMENT WHEN THE MAX NUMBER *03GH0130
C * OF COEFFICIENTS (NCMAX) IS CHANGED. *03GH0140
C * *03GH0150
C *******************************************************************03GH0160
DIMENSION E(20), P(20), PHI(20), Q(20), TB(20), A(20,20), 03GH0170
1 D(20,20), U(20,20) 03GH0180
C 03GH0190
C --- NCMAX DIMENSIONS SUBROUTINES EIGENJ AND MATINV IN THE ARGUMENTS 03GH0200
C --- OF THEIR CALL STATEMENTS. 03GH0210
NP = NC 03GH0220
TYPE 3999 03GH0221
3999 FORMAT(//30H0NONLINWOOD - DEC VERSION 1979/17H0NONLINEAR LEAST-, 03GH0222
1 29HSQUARES CURVE-FITTING PROGRAM //) 03GH0223
TYPE 3998,NPROB,NOB, NP, FLAM, FNU, MIT 03GH0224
3998 FORMAT(31H NONLINEAR ESTIMATION, , 5A4/23H NUMBER OF OBS03GH0225
1ERVATIONS, I5/23H NUMBER OF COEFFICIENTS, I5/16H STARTING LAMBDA, 03GH0226
2 F16.3/15H STARTING NU , F17.3/22H MAX NO. OF ITERATIONS, I6// )03GH0227
WRITE(KTOU, 1000)NPROB,NOB,NP,FLAM,FNU,MIT 03GH0230
1000 FORMAT(31H1NONLINEAR ESTIMATION, , 5A4/23H NUMBER OF OBS03GH0240
1ERVATIONS, I5/23H NUMBER OF COEFFICIENTS, I5/16H STARTING LAMBDA, 03GH0250
2 F16.3/15H STARTING NU , F17.3/22H MAX NO. OF ITERATIONS, I6// )03GH0260
TYPE 1001 03GH0261
WRITE(KTOU, 1001) 03GH0270
1001 FORMAT(35H0INITIAL VALUES OF THE COEFFICIENTS ) 03GH0280
TYPE 4000, (I,I=1,NP) 03GH0281
4000 FORMAT(8(4X,I2,6X)) 03GH0282
WRITE(KTOU, 2000)(I,I=1,NP) 03GH0290
2000 FORMAT(10(4X,I2,6X)) 03GH0300
TYPE 4001, (TH(I),I=1,NP) 03GH0301
4001 FORMAT(1H ,1P6E12.5) 03GH0302
WRITE(KTOU, 2001)(TH(I),I=1,NP) 03GH0310
2001 FORMAT(1H0,1P10E12.5) 03GH0320
WRITE(KTOU, 1002 ) 03GH0330
1002 FORMAT(53H0PROPORTIONS USED IN CALCULATING DIFFERENCE QUOTIENTS ) 03GH0340
WRITE(KTOU, 2000)(I,I=1,NP) 03GH0350
WRITE(KTOU, 2001)(DIFZ(I),I=1,NP) 03GH0360
IF(NOIND - NVARX) 14,14,991 03GH0361
14 IF(NP-NCMAX)15,15,992 03GH0370
15 IF(NP-2)993,16,16 03GH0380
16 GA=FLAM 03GH0390
NIT=1 03GH0400
ASSIGN 225 TO IRAN 03GH0410
ASSIGN 265 TO JORDAN 03GH0420
ASSIGN 180 TO KUWAIT 03GH0430
DELT=1.0E-08*FLOAT(NP*NP/2) 03GH0440
IF(NOB-NOBMAX)18,18,994 03GH0450
18 IF(NOB-NP)995,20,20 03GH0460
20 IF(MIT-100)10,10,996 03GH0470
10 IF(STOPCR)997,40,30 03GH0480
40 IF(STOPSS)997,60,50 03GH0490
60 ASSIGN 270 TO IRAN 03GH0500
GO TO 70 03GH0510
50 ASSIGN 265 TO IRAN 03GH0520
GO TO 70 03GH0530
30 IF(STOPSS)997,80,70 03GH0540
80 ASSIGN 270 TO JORDAN 03GH0550
70 SSQ=0 03GH0560
CALL MODEL(NPROB, TH, F, NOB, NP, X, NVARX, NOBMAX, NCMAX, KTOU) 03GH0570
DO 90 I=1,NOB 03GH0580
R(I)=Y(I)-F(I) 03GH0590
90 SSQ=SSQ+R(I)*R(I) 03GH0600
TYPE 1003, SSQ 03GH0601
WRITE(KTOU, 1003)SSQ 03GH0610
1003 FORMAT(24H0INITIAL SUM OF SQUARES:, 1PE12.4) 03GH0620
C --- BEGIN ITERATION 03GH0630
100 WRITE(KTOU, 1004)NIT 03GH0640
TYPE 1004, NIT 03GH0641
1004 FORMAT(///35X,13HITERATION NO., I4) 03GH0650
GA=GA/FNU 03GH0660
DO 130 J=1,NP 03GH0670
TEMP=TH(J) 03GH0680
TH(J)=TH(J)+DIFZ(J)*TH(J) 03GH0690
Q(J)=0 03GH0700
CALL MODEL(NPROB,TH,DELZ(1,J),NOB,NP,X,NVARX,NOBMAX, NCMAX, KTOU) 03GH0710
DO 120 I=1,NOB 03GH0720
DELZ(I,J)=DELZ(I,J)-F(I) 03GH0730
120 Q(J)=Q(J)+DELZ(I,J)*R(I) 03GH0740
Q(J)=Q(J)/(DIFZ(J)*TH(J)) 03GH0750
C --- Q=XT*R (STEEPEST DESCENT) 03GH0760
130 TH(J)=TEMP 03GH0770
DO 150 I=1,NP 03GH0780
DO 151 J=1,I 03GH0790
SUM=0 03GH0800
DO 160 K=1,NOB 03GH0810
160 SUM=SUM+DELZ(K,I)*DELZ(K,J) 03GH0820
TEMP =SUM/(DIFZ(I)*TH(I)*DIFZ(J)*TH(J)) 03GH0830
C --- D=XT*X (MOMENT MATRIX) 03GH0840
D(J,I)=TEMP 03GH0850
151 D(I,J)=TEMP 03GH0860
150 E(I)=DSQRT(D(I,I)) 03GH0870
GO TO KUWAIT,(225,265,180,270,666) 03GH0880
C --- -ITERATION 1 ONLY- 03GH0890
180 DO 200 I=1,NP 03GH0900
DO 200 J=1,I 03GH0910
SUM =D(I,J) 03GH0920
A(J,I)=SUM 03GH0930
200 A(I,J)=SUM 03GH0940
C --- KOK AND LOL ARE ALTERNATE DIMENSIONS AMD LIMIT VALUES OF 03GH0950
C --- THE MATRIX INVERSION SUBROUTINE MATINV. 03GH0960
KOK=1 03GH0970
LOL=0 03GH0980
CALL EIGENJ(A,U,DELT,NP,NCMAX) 03GH0990
TYPE 1006 03GH0991
WRITE(KTOU, 1006 ) 03GH1000
1006 FORMAT(52H0EIGENVALUES OF MOMENT MATRIX - PRELIMINARY ANALYSIS ) 03GH1010
TYPE 4001, (A(I,I),I=1,NP) 03GH1011
WRITE(KTOU, 2001)(A(I,I),I=1,NP) 03GH1020
ASSIGN 666 TO KUWAIT 03GH1030
C --- -END ITERATION 1- 03GH1040
666 DO 153 I=1,NP 03GH1050
DO 153 J=1,I 03GH1060
C --- A= SCALED MOMENT MATRIX 03GH1070
A(I,J)=D(I,J)/(E(I)*E(J)) 03GH1080
153 A(J,I)=A(I,J) 03GH1090
DO 155 I=1,NP 03GH1100
P(I)=Q(I)/E(I) 03GH1110
PHI(I)=P(I) 03GH1120
155 A(I,I)=A(I,I)+GA 03GH1130
CALL MATINV(A,NP,P,KOK,DET,NCMAX) 03GH1140
C --- P/E = CORRECTION VECTOR 03GH1150
WRITE(KTOU, 1005)DET 03GH1160
1005 FORMAT(13H0DETERMINANT:, 1PE12.4) 03GH1170
STEP=1.0 03GH1180
170 DO 220 I=1,NP 03GH1190
220 TB(I)=P(I)*STEP/E(I)+TH(I) 03GH1200
DO 401 I=1,NP 03GH1210
IF(SIGNS(I))998, 401, 402 03GH1220
402 IF(TH(I)*TB(I))403,403,401 03GH1230
401 CONTINUE 03GH1240
SUMB=0 03GH1250
CALL MODEL(NPROB, TB, F, NOB, NP, X, NVARX, NOBMAX, NCMAX, KTOU) 03GH1260
DO 230 I=1,NOB 03GH1270
R(I)=Y(I)-F(I) 03GH1280
230 SUMB=SUMB+R(I)*R(I) 03GH1290
403 SUM1=0.0 03GH1300
SUM2=0.0 03GH1310
SUM3=0.0 03GH1320
DO 231 I=1,NP 03GH1330
SUM1=SUM1+P(I)*PHI(I) 03GH1340
SUM2=SUM2+P(I)*P(I) 03GH1350
231 SUM3=SUM3+PHI(I)*PHI(I) 03GH1360
ANGLE=57.2957795*DARCOS(SUM1/DSQRT(SUM2*SUM3)) 03GH1370
C --- PRINT ANGLE IN SCALED COORDINATES 03GH1380
WRITE(KTOU, 1041)ANGLE 03GH1390
1041 FORMAT(25H0ANGLE IN SCALED. COORD,=, F5.2, 8H DEGREES ) 03GH1400
DO 2401 I=1,NP 03GH1410
IF(SIGNS(I))998,2401,2402 03GH1420
2402 IF(TH(I)*TB(I))663,663,2401 03GH1430
2401 CONTINUE 03GH1440
IF(SUMB/SSQ-1.0)662,662,663 03GH1450
663 IF(ANGLE-30.0)665,665,664 03GH1460
665 STEP=STEP/2.0 03GH1470
GO TO 170 03GH1480
664 GA=GA*FNU 03GH1490
GO TO 666 03GH1500
C --- PRINT COEFFICIENTS VIA FIT 03GH1510
662 WRITE(KTOU, 1007) 03GH1520
TYPE 1007 03GH1521
1007 FORMAT(23H0VALUES OF COEFFICIENTS ) 03GH1530
DO 669 I=1,NP 03GH1540
669 TH(I)=TB(I) 03GH1550
TYPE 4000,(I,I=1,NP) 03GH1551
WRITE(KTOU, 2000)(I,I=1,NP) 03GH1560
TYPE 4001, (TH(I),I=1,NP) 03GH1561
WRITE(KTOU, 2006)(TH(I),I=1,NP) 03GH1570
2006 FORMAT(1H , 1P10E12.5) 03GH1580
TYPE 4040, SUMB 03GH1581
4040 FORMAT(16H0SUM OF SQUARES:, 1PE12.4) 03GH1582
WRITE(KTOU, 1040)GA,SUMB 03GH1590
1040 FORMAT(8H0LAMBDA:,1PE10.3, 20X,27HSUM OF SQUARES AFTER LEAST-, 03GH1600
1 12HSQUARES FIT:, 1PE12.4) 03GH1610
GO TO IRAN,(225,265,180,270,666) 03GH1620
225 DO 240 I=1,NP 03GH1630
IF(DABS(P(I)*STEP/E(I))/(1.0E-20+DABS(TH(I)))-STOPCR)240,240,250 03GH1640
C --- ITERATION STOPS---RELATIVE CHANGE IN EACH COEFFICIENT. 03GH1650
240 CONTINUE 03GH1660
TYPE 1009, STOPCR 03GH1661
WRITE(KTOU, 1009)STOPCR 03GH1670
1009 FORMAT(54H0ITERATION STOPS - RELATIVE CHANGE IN EACH COEFFICIENT, 03GH1680
1 11H LESS THAN:, 1PE12.4) 03GH1690
GO TO 280 03GH1700
250 GO TO JORDAN,(225,265,180,270,666) 03GH1710
265 IF(SSQ-SUMB-STOPSS)260,260,270 03GH1720
C --- ITERATION STOPS---RELATIVE CHANGE IN SUM OF SQUARES. 03GH1730
260 WRITE(KTOU, 1010)STOPSS 03GH1740
TYPE 1010, STOPSS 03GH1741
1010 FORMAT(63H0ITERATION STOPS ---------- CHANGE IN SUM OF SQUARES LES03GH1750
1S THAN:, 1PE12.4) 03GH1760
GO TO 280 03GH1770
270 SSQ=SUMB 03GH1780
NIT=NIT+1 03GH1790
IF(NIT-MIT)100,100,279 03GH1800
C --- END ITERATION 03GH1810
279 TYPE 4011, MIT 03GH1811
4011 FORMAT(40H0NUMBER OF ITERATIONS EXCEEDS MAXIMUM OF,I4) 03GH1812
280 SSQ=SUMB 03GH1820
IDF=NOB-NP 03GH1830
C --- PRINT CORRELATION MATRIX 03GH1840
TYPE 1015 03GH1841
WRITE(KTOU, 1015 ) 03GH1850
1015 FORMAT(19H0CORRELATION MATRIX ) 03GH1860
CALL MATINV(D,NP,P,LOL,DET,NCMAX) 03GH1870
TYPE 1005, DET 03GH1871
WRITE(KTOU, 1005)DET 03GH1880
DO 7692 I=1,NP 03GH1890
7692 E(I)=DSQRT(D(I,I)) 03GH1900
TYPE 4005, (I,I=1,NP) 03GH1901
4005 FORMAT(3X,7(4X,I2,5X)) 03GH1902
WRITE(KTOU, 2005)(I,I=1,NP) 03GH1910
2005 FORMAT(3X,10(4X,I2,5X)) 03GH1920
DO 390 I=1,NP 03GH1930
DO 340 J=I,NP 03GH1940
A(J,I)=D(J,I)/(E(I)*E(J)) 03GH1950
D(J,I)=D(J,I)/(DIFZ(I)*TH(I)*DIFZ(J)*TH(J)) 03GH1960
D(I,J)=D(J,I) 03GH1970
340 A(I,J)=A(J,I) 03GH1980
TYPE 4009, I,(A(J,I),J=1,I) 03GH1981
4009 FORMAT(I3,7(2X,F7.4,2X)) 03GH1982
390 WRITE(KTOU, 2009)I,(A(J,I),J=1,I) 03GH1990
2009 FORMAT(I3,10(2X,F7.4,2X)) 03GH2000
IF(IDF)995,410,7058 03GH2010
7058 RMS = SSQ/FLOAT(IDF) 03GH2020
SDEV = DSQRT(RMS) 03GH2030
YYMIN = 1E32 03GH2040
YYMAX = -YYMIN 03GH2050
DO 500 I = 1,NOB 03GH2060
CALL MINMAX(YYMIN, YYMAX, Y(I) ) 03GH2070
500 CONTINUE 03GH2080
YW = YYMAX - YYMIN 03GH2090
WRITE(KTOU, 1050) NPROB, YYMIN, YYMAX, YW 03GH2100
1050 FORMAT(1H1,28X,45HNONLINEAR LEAST-SQUARES CURVE-FITTING PROGRAM// 03GH2110
1 1X, 5A4, 2X, 9HDEP.VAR.:, 1X, 7HMIN Y =, 1PE10.3, 03GH2120
2 2X, 8HMAX Y = , 1PE10.3, 2X,10HRANGE Y = , 1PE10.3 / ) 03GH2130
IF(NEQU.EQ.0) GO TO 502 03GH2140
WRITE(KTOU, 1051) (EQU(I), I=1,NEQU ) 03GH2150
1051 FORMAT(1H , 18X, 18A4 ) 03GH2160
502 IF(NOMIT.EQ.0) GO TO 944 03GH2161
IF(NFOUND.EQ.0) GO TO 909 03GH2162
WRITE(KTOU, 907) (JOMIT(JZ),JZ=1,NFOUND) 03GH2163
907 FORMAT(19X,22HOBSERVATIONS DELETED: ,7(A6,2X),7(/41X,7(A6,2X))) 03GH2164
IF((NOMIT-NFOUND).EQ.0) GO TO 944 03GH2165
J = NFOUND + 1 03GH2166
909 WRITE(KTOU, 910) (JOMIT(JZ),JZ=J,NOMIT) 03GH2167
910 FORMAT(19X,24HOBSERVATIONS NOT FOUND: ,7(A6,2X),7(/43X,7(A6,2X))) 03GH2168
944 WRITE(KTOU, 1052 ) 03GH2170
1052 FORMAT(1H0, 62X,21H95% CONFIDENCE LIMITS / 03GH2180
1 1X,10HIND.VAR(I), 3X, 4HNAME, 5X, 9HCOEF.B(I), 6X,10HS.E. COEF , 03GH2190
2 4X, 7HT-VALUE, 4X,22HLOWER UPPER ) 03GH2200
TVAR=TTEST(IDF) 03GH2210
DO 550 I = 1,NP 03GH2220
SECOEF= E(I)*SDEV 03GH2230
TVALUE= TH(I)/SECOEF 03GH2240
TSEC = TVAR*SECOEF 03GH2250
TMCOE = TH(I) - TSEC 03GH2260
TPCOE = TH(I) + TSEC 03GH2270
K = 2*I 03GH2280
J = K - 1 03GH2290
WRITE(KTOU, 2010) I, BI(J), BI(K), TH(I), SECOEF, TVALUE, TMCOE, 03GH2300
1 TPCOE 03GH2310
2010 FORMAT(1H , 4X, I2, 6X, A4, A2, 2X, 1PE12.5, 5X, 1PE9.2, 4X, 03GH2320
1 0PF6.1, 5X, 1PE9.2, 5X, 1PE9.2) 03GH2330
550 CONTINUE 03GH2340
WRITE(KTOU, 1054) NOB 03GH2350
1054 FORMAT(20H0NO. OF OBSERVATIONS, 11X, I5) 03GH2360
WRITE(KTOU, 1056) NP 03GH2370
1056 FORMAT(20H NO. OF COEFFICIENTS, 11X, I5) 03GH2380
NDFR = IDF + 0.01 03GH2390
WRITE(KTOU, 1058) NDFR 03GH2400
1058 FORMAT(28H RESIDUAL DEGREES OF FREEDOM, 4X, I4) 03GH2410
WRITE(KTOU, 1062) SDEV 03GH2420
1062 FORMAT(26H RESIDUAL ROOT MEAN SQUARE, 3X, F16.8) 03GH2430
WRITE(KTOU, 1064) RMS 03GH2440
1064 FORMAT(21H RESIDUAL MEAN SQUARE, 8X, F16.8) 03GH2450
WRITE(KTOU, 1066) SSQ 03GH2460
1066 FORMAT(24H RESIDUAL SUM OF SQUARES, 5X, F16.8) 03GH2470
GO TO 99999 03GH2480
410 WRITE(KTOU, 1033)NPROB 03GH2490
1033 FORMAT(1H1,18HEND OF PROBLEM ,5A4) 03GH2500
GO TO 88888 03GH2510
991 WRITE(KTOU, 3034) NVARX 03GH2511
3034 FORMAT(1H0,11H***ERROR***/32H0EXCESSIVE NUMBER OF VARIABLES. , 03GH2512
1 36H PROGRAM DIMENSIONED FOR MAXIMUM OF , I3 ) 03GH2513
GO TO 410 03GH2514
992 WRITE(KTOU, 3035) NCMAX 03GH2520
3035 FORMAT(1H0,11H***ERROR***/35H0EXCESSIVE NUMBER OF COEFFICIENTS. , 03GH2530
1 36H PROGRAM DIMENSIONED FOR MAXIMUM OF , I3 ) 03GH2540
GO TO 410 03GH2550
993 WRITE(KTOU, 3036 ) 03GH2560
3036 FORMAT(1H0,11H***ERROR***/26H0LESS THAN 2 COEFFICIENTS. ) 03GH2570
GO TO 410 03GH2580
994 WRITE(KTOU, 3037) NOBMAX 03GH2590
3037 FORMAT(1H0,11H***ERROR***/35H0EXCESSIVE NUMBER OF OBSERVATIONS. , 03GH2600
1 36H PROGRAM DIMENSIONED FOR MAXIMUM OF , I4 ) 03GH2610
GO TO 410 03GH2620
995 WRITE(KTOU, 3038 ) 03GH2630
3038 FORMAT(1H0,11H***ERROR***/55H0NUMBER OF COEFFICIENTS EXCEEDS NUMBE03GH2640
1R OF OBSREVATIONS. ) 03GH2650
GO TO 410 03GH2660
996 WRITE(KTOU, 3039 ) 03GH2670
3039 FORMAT(1H0,11H***ERROR***/33H0EXCESSIVE NUMBER OF ITERATIONS. , 03GH2680
1 44H PROGRAM LIMIT SWITCH SET FOR MAXIMUM OF 100 ) 03GH2690
GO TO 410 03GH2700
997 WRITE(KTOU, 3040) STOPSS, STOPCR 03GH2710
3040 FORMAT(1H0,11H***ERROR***/36H0STOP CRITERIA IS NEGATIVE. STOPSS =,03GH2720
1 F12.7,10H, STOPRC =, F12.7 ) 03GH2730
GO TO 410 03GH2740
998 WRITE(KTOU, 3041) I 03GH2750
3041 FORMAT(1H0,11H***ERROR***/12H0COEFFICIENT ,I4,18H HAS NEGATIVE SIG03GH2760
1N ) 03GH2770
GO TO 410 03GH2780
88888 IERGAU = 1 03GH2790
99999 RETURN 03GH2800
END 03GH2810