#!/bin/sh # To unpack, sh this message. For more explanation, read the next few lines. # --This message is a .shar file, i.e., a shell archive. # --To unpack the files it contains into some empty directory DIR, # --first cd (change directory) to DIR. # --Then put this message into a file in DIR. # --(Use a FILENAME which differs from those to be created!) # --Remove from file FILENAME the lines prior to the "#!/bin/sh" line above. # --Finally execute "sh FILENAME". PATH=/bin:/usr/bin echo processing file mapclus.f 1>&2 cat > mapclus.f <<'CUT HERE------------ mapclus.f' C mapclus.f C The authors of this software are Phipps Arabie and J Douglas Carroll 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 SECOND MDS Package of AT&T Bell Laboratories C For explanation of the method this software implements, see C "Three-Way Scaling and Clustering" by P Arabie, J D Carroll, and C W DeSarbo, a monograph published by Sage Publications, Beverly Hills, C Calif. in 1987 as series 07, item 065, in the Sage University Papers C Arabie,P., & Carroll,J.D. (1980). MAPCLUS: A mathematical programming C approach to fitting the ADCLUS model. Psychometrika, 45, 211-235. C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ C$ FORTY FORM,XREF MAPC 5 C MAPCLUS, JULY, 1980 (SOME OUTPUT FORMATS UPDATED 7/82) MAPC 10 C PROGRAM MAPCLS(PARAM,OUTPUT,PNCH,MATRIX,INIC,TAPE5=PARAM, MAPC 15 C 2 TAPE6=OUTPUT,TAPE7=PNCH,TAPE41=MATRIX,TAPE16=INIC) MAPC 20 C MAPC 25 C THE FOLLOWING NON-ANSI STANDARD FORTRAN CHARACTERS (: ;) ARE USED MAPC 30 C FOR COLON AND SEMICOLON, RESPECTIVELY. IF YOUR COMPILER OR MAPC 35 C OPERATING SYSTEM CANNOT DIGEST THESE CHARACTERS, SUBSTITUTING MAPC 40 C A BLANK IS THE BEST REMEDY. MAPC 45 C MAPC 50 DOUBLE PRECISION DCON MAPC 55 INTEGER FMT, TITL MAPC 60 DIMENSION FMT(72), SVAF(050), IDAT(2), IDIG(10), GL(030,2), MAPC 65 * GRLEN(2), GRMEN(2), OBSLF(2), ASTEP(2) MAPC 70 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), MAPC 75 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, MAPC 80 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC MAPC 85 COMMON /A2/ G(030), KT1(10), KT2(10), TK(10), LP(15), IPUN, IERR, MAPC 90 * LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), TRISE(025) MAPC 95 COMMON /A4/ DAPI(030), BCONST, TM, U, V MAPC 100 COMMON /A5/ AK, ALFRAY(025), BETRAY(025), IADJ(050), JADJ MAPC 105 COMMON /A6/ DCON, DATA(030,030), ALOS, IXV MAPC 110 C MAPC 115 C BOTH ARRAYS IN COMMON /A7/ MUST BE DIMENSIONED FOR MAX(N30,N2). MAPC 120 C IN THE ORIGINAL SOURCE VERSION OF THIS PROGRAM, THOSE MAPC 125 C RESPECTIVE VALUES ARE (30,25); BOTH ARE PRECEDED BY A LEAD ZERO MAPC 130 C (TO CREATE A UNIQUE STRING IN COL. 1-72) AND CORRESPOND TO MAPC 135 C (MAXIMUM POSSIBLE NUMBER OF STIMULI, MAXIMUM POSSIBLE NUMBER OF MAPC 140 C SUBSETS INCLUDING THE COMPLETE SUBSET FOR THE ADDITIVE CONSTANT). MAPC 145 C MAPC 150 C MAPC 155 COMMON /A7/ A(030), IZ(030) MAPC 160 COMMON /A8/ W(025), IPRE, ITPRE MAPC 165 COMMON /A9/ ILAB(030,6), ALD1 MAPC 170 DATA IDIG(1), IDIG(2), IDIG(3), IDIG(4), IDIG(5), IDIG(6), MAPC 175 * IDIG(7), IDIG(8), IDIG(9), IDIG(10), IBL, IL, IR /1H0,1H1,1H2, MAPC 180 * 1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H ,1H(,1H)/ MAPC 185 C IDAT(1) = IDATEZ(0) MAPC 190 C MAPC 195 C MAPC 200 C ASSUME CARD INPUT STREAM AS LOGICAL DEVICE FOR READING PROGRAM MAPC 205 C CONTROL PARAMETERS. MAPC 210 C MAPC 215 L5 = 5 MAPC 220 C MAPC 225 C MAPC 230 LOUT = 6 MAPC 235 C MAPC 240 C THE NEXT STATEMENT ASSUMES THAT THE DEFAULT DEVICE FOR PUNCHED MAPC 245 C OUTPUT IS 7. THERE IS NO PROVISION FOR REWINDING. MAPC 250 C MAPC 255 LPUN7 = 7 MAPC 260 C MAPC 265 C N2 IS THE MAXIMUM NUMBER OF CLUSTERS (INCLUDING COMPLETE MAPC 270 C SUBSET FOR ADDITIVE CONSTANT) ALLOWED BY CURRENT DIMENSION MAPC 275 C STATEMENTS. MAPC 280 C MAPC 285 N2 = (025) MAPC 290 N3 = N2 - 1 MAPC 295 N5 = (050) MAPC 300 C MAPC 305 C N30 IS THE MAXIMUM NUMBER OF STIMULI ALLOWED BY CURRENT MAPC 310 C DIMENSION STATEMENTS. MAPC 315 C MAPC 320 N30 = (030) MAPC 325 BCRIT = 1.0E-5 MAPC 330 JTPRE = 6 MAPC 335 C MAPC 340 C FOR INSTALLATIONS OTHER THAN BTL, LREAD SHOULD BE 5 MAPC 345 C MAPC 350 5 LREAD = 5 MAPC 355 IPRE = 0 MAPC 360 VAF = 0.0 MAPC 365 ICON = 0 MAPC 370 ISWIT = 0 MAPC 375 IPRN = 0 MAPC 380 IXV = 0 MAPC 385 DO 10 I=1,15 MAPC 390 LP(I) = 0 MAPC 395 10 CONTINUE MAPC 400 C MAPC 405 C MAPC 410 WRITE (LOUT,99999) MAPC 415 99999 FORMAT (1H1///11X, 28HRUN OF MAPCLUS VERSION 1 , 7X, 7H A MATH,MAPC 420 * 30HEMATICAL PROGRAMMING APPROACH , 26HBY P. ARABIE AND J. D. CAR,MAPC 425 * 4HROLL/46X, 43H TO FITTING THE SHEPARD-ARABIE ADCLUS MODEL//) MAPC 430 C WRITE (LOUT,8) IDAT(1) MAPC 435 C 8 FORMAT(30X,A6) MAPC 440 C WRITE (LOUT,8) IDAT(1) MAPC 445 C 8 FORMAT(30X,A10) MAPC 450 C MAPC 455 C MAPC 460 C MAPC 465 C MAPC 470 READ (L5,99998,END=551) LD1,ICON,KREAD,IREWIN,M15,IPRN,IPUN, MAPC 475 * LPUNCH ,(LP(I), I = 1, 15) MAPC 480 99998 FORMAT (8I3, 1X, 15I1) MAPC 485 WRITE (LOUT,99997) LD1, ICON, KREAD, IREWIN, M15, IPRN, IPUN, MAPC 490 * LPUNCH, (I,LP(I),I=1,15) MAPC 495 99997 FORMAT (//4X, 48H LD1 ICON KREAD IREWIN M15 IPRN IPUN LPU,MAPC 500 * 3HNCH, 54H ABBREVIATED -0- AND FULL-LENGTH -1- PRINTING DUR,MAPC 505 * 6HING 15/2X, 3I6, 2X, I6, 1X, I5, 3I6, 13X, 14H ITERATIONS OF, MAPC 510 * 27H COMBINATORIAL OPTIMIZATION//50X, 15(1X, I2, 1H:, I1)) MAPC 515 C MAPC 520 C LD1 IS THE NUMBER OF CLUSTERS TO BE FITTED, AS SPECIFIED BY MAPC 525 C THE USER, AND DOES NOT INCLUDE THE COMPLETE SET ASSOCIATED MAPC 530 C WITH THE ADDITIVE CONSTANT. MAPC 535 C MAPC 540 WRITE (LOUT,99996) LD1 MAPC 545 99996 FORMAT (//10X, 36HNUMBER OF CLUSTERS (LD1) SPECIFIED =, I3, MAPC 550 * 60H (NOT INCLUDING THE COMPLETE SET FOR THE ADDITIVE CONSTANT).) MAPC 555 IF (LD1.GT.0) GO TO 15 MAPC 560 WRITE (LOUT,99995) MAPC 565 99995 FORMAT (///46H THIS PROGRAM WILL NOT WORK CORRECTLY IF LD1, , MAPC 570 * 51HTHE NUMBER OF SUBSETS IS NOT DECLARED AS A POSITIVE, 6H INTEG,MAPC 575 * 3HER./22H EXECUTION TERMINATES.) MAPC 580 STOP MAPC 585 15 IF (LD1.LE.N3) GO TO 20 MAPC 590 WRITE (LOUT,99994) LD1, N3 MAPC 595 99994 FORMAT (///20H YOU HAVE SPECIFIED , I3, 19H CLUSTERS, BUT THIS, MAPC 600 * 36H VERSION OF THE PROGRAM ONLY ALLOWS , I3, 9H CLUSTERS/ MAPC 605 * 52H (PLUS ONE FOR THE ADDITIVE CONSTANT). INCREASE THE, MAPC 610 * 34H DIMENSION STATEMENTS ACCORDINGLY.) MAPC 615 STOP MAPC 620 C MAPC 625 C MAPC 630 C IF ICON IS NEGATIVE, A RANDOM INITIAL CONFIGURATION MAPC 635 C WILL BE GENERATED, AND THE FIRST (-ICON) NUMBERS MAPC 640 C FROM THE RANDOM NUMBER GENERATOR WILL BE DISCARDED. MAPC 645 C IF ICON IS ZERO, A RATIONAL INITIAL CONFIGURATION MAPC 650 C WILL BE PRODUCED BY THE PROGRAM. MAPC 655 C IF ICON IS POSITIVE, A USER-SUPPLIED INITIAL MAPC 660 C CONFIGURATION WILL BE READ FROM LOGICAL DEVICE MAPC 665 C NUMBER (ICON). MAPC 670 C MAPC 675 20 IF (LD1.LE.N30) GO TO 25 MAPC 680 WRITE (LOUT,99993) MAPC 685 99993 FORMAT (/////37H WARNING: WHENEVER NUMBER OF DECLARED, 8H CLUSTER,MAPC 690 * 28HS EXCEEDS NUMBER OF STIMULI,, 28H MAKE SURE THAT THE DIMENSIO,MAPC 695 * 26HN STATEMENTS IN COMMON/A7//31H CORRESPOND TO MAX(N30, N2) (TH,MAPC 700 * 23HAT STATEMENT --MAY-- BE, 19H CORRECT AS IT IS.)///) MAPC 705 25 IF (KREAD.GT.0) GO TO 30 MAPC 710 C MAPC 715 C KREAD IS THE LOGICAL DEVICE NUMBER FOR THE NUMERICAL MAPC 720 C PROXIMITIES DATA AND THE (OPTIONAL) LABELS. THE DEFAULT VALUE MAPC 725 C IS 5. MAPC 730 C MAPC 735 WRITE (LOUT,99992) KREAD, LREAD MAPC 740 99992 FORMAT (//10X, 40HSINCE THE DATA SET WAS DECLARED TO BE ON, MAPC 745 * 15H DEVICE NUMBER , I4, 30H (KREAD) THE DEFAULT VALUE OF , I3, MAPC 750 * 22H WILL INSTEAD BE USED.) MAPC 755 GO TO 35 MAPC 760 30 LREAD = KREAD MAPC 765 35 WRITE (LOUT,99991) LREAD MAPC 770 99991 FORMAT (//9X, 46H DATA WILL BE READ FROM LOGICAL DEVICE NUMBER , MAPC 775 * I3) MAPC 780 IF (IREWIN.EQ.1 .AND. LREAD.NE.L5) GO TO 40 MAPC 785 C MAPC 790 C IREWIN ALLOWS THE USER THE OPTION OF HAVING THE LOGICAL DEVICE MAPC 795 C ASSOCIATED WITH KREAD REWOUND, UNLESS LREAD .EQ. L5 MAPC 800 C MAPC 805 WRITE (LOUT,99990) MAPC 810 99990 FORMAT (9X, 27H WHICH WILL NOT BE REWOUND.) MAPC 815 GO TO 45 MAPC 820 40 REWIND LREAD MAPC 825 WRITE (LOUT,99989) MAPC 830 99989 FORMAT (9X, 24H WHICH HAS BEEN REWOUND.) MAPC 835 45 WRITE (LOUT,99988) LREAD, LOUT MAPC 840 99988 FORMAT (//9X, 43H I/O CHANNELS FOR READING DATA AND PRINTING, MAPC 845 * 31H OUTPUT: READ PRINT/58X, 12H LREAD =, I3, MAPC 850 * 9H LOUT =, I3) MAPC 855 C MAPC 860 C IF M15, THE NUMBER OF MAJOR ITERATIONS, IS ZERO, ONLY (GLOBAL) MAPC 865 C REGRESSION WILL BE PERFORMED ON AN INITIAL CONFIGURATION, MAPC 870 C PRESUMABLY SUPPLIED BY THE USER ON CHANNEL ICON. MAPC 875 C MAPC 880 C IF M15 IS NEGATIVE, GLOBAL REGRESSION IS PERFORMED (AS WHEN M15 MAPC 885 C .EQ. 0), FOLLOWED IMMEDIATELY BY COMBINATORIAL OPTIMIZATION, MAPC 890 C UP TO A MAXIMUM OF 15 ITERATIONS. MAPC 895 C MAPC 900 C IF M15 IS POSITIVE, THEN IT DEFINES THE MAXIMUM POSSIBLE MAPC 905 C NUMBER OF MAJOR ITERATIONS (TOTAL OF PRE-POLISHING AND MAPC 910 C POLISHING). MAPC 915 C MAPC 920 C MAPC 925 WRITE (LOUT,99987) M15 MAPC 930 99987 FORMAT (//10X, 47HMAXIMUM NUMBER OF MAJOR ITERATIONS (M15) SPECIF,MAPC 935 * 3HIED, 5H IS =, I3, 1H.) MAPC 940 IF (M15.LE.N5) GO TO 50 MAPC 945 WRITE (LOUT,99986) M15, N5 MAPC 950 99986 FORMAT (//20H YOU HAVE SPECIFIED , I4, 23H ITERATIONS, BUT THIS V,MAPC 955 * 6HERSION, 40H OF THE PROGRAM ALLOWS FOR AT MOST ONLY , I4, MAPC 960 * 14H ITERATIONS. /43H THE DIMENSION STATEMENT COVERING ARRAY KP ,MAPC 965 * 12HIN THE MAIN , 44H PROGRAM WILL HAVE TO BE INCREASED. EXECUTI,MAPC 970 * 14HON TERMINATES.) MAPC 975 STOP MAPC 980 C MAPC 985 C MAPC 990 C IF IPRN .EQ. 1 THE INPUT DATA WILL BE LISTED (A GOOD IDEA ON MAPC 995 C THE FIRST TIME A GIVEN DATA SET IS ENTERED). MAPC1000 C MAPC1005 C IF IPUN .EQ. 1, THE FINAL CONFIGURATION WILL BE WRITTEN ON MAPC1010 C DISK (OR PUNCHED), USING LOGIAL DEVICE (LPUNCH). MAPC1015 C MAPC1020 50 IF (IPUN.EQ.0) GO TO 55 MAPC1025 IF (LPUNCH.EQ.0) LPUNCH = LPUN7 MAPC1030 WRITE (LOUT,99985) LPUNCH MAPC1035 99985 FORMAT (/10X, 48HCONFIGURATION OF SUBSETS WILL BE WRITTEN TO DISK,MAPC1040 * 12H ON CHANNEL , I3, 35H AT END OF (MAJOR) POST-ITERATIONS.) MAPC1045 C MAPC1050 C MAPC1055 C GET THE INITIAL CONFIGURATION OF P(I). MAPC1060 C MAPC1065 55 IF (ICON) 60, 70, 75 MAPC1070 60 IX = -ICON MAPC1075 WRITE (LOUT,99984) IX MAPC1080 99984 FORMAT (//10X, 41HSINCE ICON IS NEGATIVE, A RANDOM INITIAL , MAPC1085 * 31HCONFIGURATION WILL BE GENERATED/10X, 13HAND THE FIRST, I5, MAPC1090 * 34H RANDOM NUMBERS WILL BE DISCARDED.//) MAPC1095 DO 65 I=1,IX MAPC1100 AMN = RUNIFV(0) MAPC1105 65 CONTINUE MAPC1110 GO TO 80 MAPC1115 70 WRITE (LOUT,99983) MAPC1120 99983 FORMAT (//10X, 47HSINCE ICON = 0, A RATIONAL INITIAL CONFIGURATIO,MAPC1125 * 1HN, 37H WILL BE GENERATED (DEFAULT OPTION).//) MAPC1130 GO TO 80 MAPC1135 75 WRITE (LOUT,99982) ICON MAPC1140 99982 FORMAT (//10X, 47HSINCE ICON IS POSITIVE, A USER-SUPPLIED INITIAL,MAPC1145 * 51H CONFIGURATION WILL BE READ FROM LOGICAL DEVICE NO., I4, 1H.//MAPC1150 * ) MAPC1155 80 IF (M15.GT.0) GO TO 90 MAPC1160 IF (ICON.GT.0) GO TO 85 MAPC1165 WRITE (LOUT,99981) M15, ICON MAPC1170 99981 FORMAT (//51H USER HAS SPECIFIED INCOMPATIBLE OPTIONS. ALTHOUGH, MAPC1175 * 47H COMPUTATION IS TO START WITH REGRESSION (M15 =, I3, 2H),/ MAPC1180 * 61H NO USER-SUPPLIED INITIAL CONFIGUATION HAS BEEN GIVEN (ICON =,MAPC1185 * I3, 25H). EXECUTION TERMINATES.) MAPC1190 STOP MAPC1195 85 IF (M15.EQ.0) GO TO 100 MAPC1200 WRITE (LOUT,99980) MAPC1205 99980 FORMAT (//10X, 40HSINCE M15 IS NEGATIVE, THERE WILL BE NO , MAPC1210 * 45HITERATIVE COMPUTATION OF WEIGHTS AND SUBSETS,/10X, 8HAND THE ,MAPC1215 * 41HPROGRAM WILL PROCEED TO DO COMBINATORIAL , 15HOPTIMIZATION IM,MAPC1220 * 27HMEDIATELY AFTER REGRESSION.) MAPC1225 GO TO 100 MAPC1230 C MAPC1235 C IF THERE ARE TO BE NO ITERATIONS OF THE MATHEMATICAL PROGRAMMING MAPC1240 C ALGORITHM, THEN DON"T READ IN THE REMAINING (IRRELEVANT) VARIABLESMAPC1245 C MAPC1250 C MAPC1255 C MAPC1260 C MAPC1265 C MAPC1270 90 READ (L5,99979) (KP(I),I=1,M15) MAPC1275 99979 FORMAT (72I1) MAPC1280 C MAPC1285 C KP IS A BINARY VECTOR THAT CONTROLS THE AMOUNT OF MAPC1290 C PRINTED OUTPUT FOR EACH MAJOR ITERATION. WHEN M15 .LE. 0, KP IS MAPC1295 C IRRELEVANT AND IS THEREFORE NOT READ IN. MAPC1300 C MAPC1305 WRITE (LOUT,99978) MAPC1310 99978 FORMAT (///6X, 47HABBREVIATED (0) AND FULL-LENGTH (1) PRINTOUT FO,MAPC1315 * 2HR , 17HMAJOR ITERATIONS:) MAPC1320 WRITE (LOUT,99977) (I, I = 1, M15) MAPC1325 99977 FORMAT (/6X, 35I3) MAPC1330 WRITE (LOUT,99977) (KP(I),I=1,M15) MAPC1335 C MAPC1340 C MAPC1345 C MAPC1350 C MAPC1355 C MAPC1360 C MAPC1365 READ (L5,99998) ITPRE, IUPPER, JLIM, JPOST MAPC1370 WRITE (LOUT,99976) ITPRE, IUPPER, JLIM, JPOST MAPC1375 99976 FORMAT (//6X, 26HITPRE IUPPER JLIM JPOST/5X, 2I6, 1X, I6, 1X, MAPC1380 * I6) MAPC1385 C MAPC1390 C ITPRE IS THE MAXIMUM NUMBER OF MAJOR ITERATIONS ALLOWED PRIOR TO MAPC1395 C POLISHING. (I. E., PRE-POLISHING ITERATIONS.) IF ITPRE IS MAPC1400 C NEGATIVE, ITERATIVE COMPUTATION BEGINS WITH DE NOVO MAJOR MAPC1405 C ITERATIONS. MAPC1410 C MAPC1415 C IF ITPRE IS ZERO, THE DEFAULT VALUE OF 6 IS SUPPLIED. MAPC1420 C THERE IS NO OPTION TO BEGIN COMPUTATION WITH THE FIRST MAPC1425 C POLISHING ITERATION, BUT ITPRE.LT.0 STARTS COMPUTATION WITH MAPC1430 C DE NOVO ITERATIONS. MAPC1435 C MAPC1440 C MAPC1445 C IUPPER IS THE MAXIMUM NUMBER OF INNER ITERATIONS IN THE OUTER MAPC1450 C ITERATIONS OF THE PRE-POLISHING MAJOR ITERATIONS. MAPC1455 C MAPC1460 C JLIM IS THE MAXIMUM NUMBER OF TIMES ALPHA AND BETA CAN BE MAPC1465 C ADJUSTED DURING AN OUTER ITERATION OF A PRE-POLISHING MAJOR MAPC1470 C ITERATION. MAPC1475 C MAPC1480 C JPOST IS THE MAXIMUM NUMBER OF INNER ITERATIONS OF AN OUTER MAPC1485 C ITERATION OF A POLISHING MAJOR ITERATION. MAPC1490 C MAPC1495 C MAPC1500 C MAPC1505 C THE NEXT STATEMENT ALLOWS THE OPTION OF STARTING A RUN WITH MAPC1510 C DE NOVO MAJOR ITERATIONS, AND PROBABLY SHOULD BE USED ONLY MAPC1515 C FOR RESUMING AN EARLIER ANALYSIS. MAPC1520 C MAPC1525 KTPRE = ITPRE MAPC1530 IF (ITPRE.GE.0) GO TO 95 MAPC1535 IPRE = 1 MAPC1540 WRITE (LOUT,99975) ITPRE MAPC1545 99975 FORMAT (//10X, 42HUSER HAS DECLARED THAT ITPRE IS NEGATIVE (, I3, MAPC1550 * 52H), SO THAT COMPUTATION WILL BEGIN WITH DE NOVO MAJOR, MAPC1555 * 12H ITERATIONS.) MAPC1560 IF (ICON.LT.0) WRITE (LOUT,99974) ICON MAPC1565 99974 FORMAT (//10X, 47HUSER HAS SPECIFIED RANDOM INITIAL CONFIGURATION,MAPC1570 * 8H (ICON =, I3, 45H) WHICH IS INCONSISTENT WITH A DE NOVO START;/MAPC1575 * 10X, 47HINITIAL CONFIGURATION WILL INSTEAD BE RATIONAL.) MAPC1580 95 IF (ITPRE.EQ.0) ITPRE = JTPRE MAPC1585 ITPRE1 = ITPRE + 1 MAPC1590 IF (IUPPER.LE.0) IUPPER = 50 MAPC1595 IF (ITPRE.GT.0) WRITE (LOUT,99973) ITPRE, IUPPER MAPC1600 99973 FORMAT (//10X, 34HMAXIMUM NUMBER OF MAJOR ITERATIONS, 9H PRIOR TO,MAPC1605 * 20H POLISHING (ITPRE) =, I4, 1H,/10X, 17HWITH A MAXIMUM OF, I5, MAPC1610 * 26H INNER ITERATIONS (IUPPER), 21H PER OUTER ITERATION.) MAPC1615 C MAPC1620 IF (JLIM.LE.0) JLIM = 4 MAPC1625 KLIM = JLIM MAPC1630 IF (JPOST.LE.0) JPOST = 75 MAPC1635 IF (ITPRE.GT.0) WRITE (LOUT,99972) KLIM MAPC1640 99972 FORMAT (//9X, 20H USER HAS SPECIFIED , 18HA LIMIT (JLIM) OF , I3, MAPC1645 * 28H CALL(S) TO ADJUST PER OUTER, 28H ITERATION PRIOR TO POLISHIN,MAPC1650 * 3HG. ) MAPC1655 WRITE (LOUT,99971) JPOST MAPC1660 99971 FORMAT (//10X, 14HTHERE WILL BE , I3, 24H INNER ITERATIONS (JPOST,MAPC1665 * 2H) , 37HPER OUTER ITERATION DURING POLISHING.) MAPC1670 C MAPC1675 C MAPC1680 C MAPC1685 C MAPC1690 READ (L5,99968) ALF0 MAPC1695 C MAPC1700 C ALF0 IS THE INITIAL ESTIMATE OF THE STEP SIZE (DEFAULT =1.0) MAPC1705 C MAPC1710 IF (ALF0.LE.0.0) ALF0 = 1.0 MAPC1715 WRITE (LOUT,99970) ALF0 MAPC1720 99970 FORMAT (//6X, 6HALF0 =, F12.8, 31H (INTERVAL FOR INITIAL ESTIMATE,MAPC1725 * 15H OF STEP SIZE).) MAPC1730 C MAPC1735 C MAPC1740 C MAPC1745 C MAPC1750 READ (L5,99968) ALPHAI, AK, ACRIT, BLCRIT MAPC1755 C MAPC1760 WRITE (LOUT,99969) ALPHAI, AK, ACRIT, BLCRIT MAPC1765 99969 FORMAT (//6X, 39HALPHAI AK ACRIT BLCRIT/3X, MAPC1770 * 4F11.8) MAPC1775 C MAPC1780 C ALPHAI IS THE WEIGHT FOR THE A-PART OF THE LOSS FUNCTION. MAPC1785 C MAPC1790 C AK IS THE CONSTANT IN THE FORMULA FOR ADJUSTING ALPHA AND BETA. MAPC1795 C MAPC1800 C ACRIT IS THE CRITERION OF CONVERGENCE FOR THE RELATIVE GRADIENT. MAPC1805 C MAPC1810 C BLCRIT IS THE CRITERION OF CONVERGENCE FOR THE B-PART OF THE MAPC1815 C LOSS FUNCTION. MAPC1820 C MAPC1825 99968 FORMAT (4F10.6) MAPC1830 C MAPC1835 C ESTABLISHING A DEFAULT VALUE FOR ALPHAI: MAPC1840 C MAPC1845 IF (ALPHAI.LE.0.0) ALPHAI = .40 MAPC1850 BETAI = 1.0 - ALPHAI MAPC1855 IF (AK.LE.0.0) AK = 2.0 MAPC1860 IF (ACRIT.LE.0.0) ACRIT = 0.005 MAPC1865 IF (BLCRIT.LE.0.0) BLCRIT = .05 MAPC1870 WRITE (LOUT,99967) ALPHAI, AK, ACRIT, BLCRIT MAPC1875 99967 FORMAT (//10X, 8HALPHAI =, F12.5, 28H (FOR WEIGHTING A-PART OF LO,MAPC1880 * 13HSS FUNCTION).//10X, 5HAK = , F12.5, 22H (FOR ADJUSTING ALPHA ,MAPC1885 * 10HAND BETA).//10X, 7HACRIT =, F12.5, 23H (CRITERION FOR CONVERG,MAPC1890 * 27HENCE OF RELATIVE GRADIENT).//10X, 8HBLCRIT =, F12.5, 6H (CRIT,MAPC1895 * 20HERION FOR B-PART OF , 15HLOSS FUNCTION).//) MAPC1900 C MAPC1905 C MAPC1910 C MAPC1915 C MAPC1920 C MAPC1925 100 READ (LREAD,99965) TITL MAPC1930 C MAPC1935 C MAPC1940 READ (LREAD,99966) N, ISWIT, LABL MAPC1945 99966 FORMAT (3I2) MAPC1950 C MAPC1955 C MAPC1960 READ (LREAD,99965) FMT MAPC1965 99965 FORMAT (72A1) MAPC1970 C MAPC1975 C N IS THE NUMBER OF STIMULI. MAPC1980 C MAPC1985 C ISWIT = 1 FOR DISSIMILARITIES DATA. MAPC1990 C ISWIT .LE. 0 FOR SIMILARITIES DATA. MAPC1995 C MAPC2000 C IF ILAB .EQ. 1, THE PROGRAM EXPECTS TO FIND 6-CHARACTER LABELS MAPC2005 C (12 TO A CARD, IN COLUMNS 1 - 72) IMMEDIATELY FOLLOWING THE MAPC2010 C NUMERICAL DATA; OTHERWISE, INTERNALLY GENERATED LABELS OF 1,...,N MAPC2015 C WILL BE USED FOR THE STIMULI. MAPC2020 C MAPC2025 C MAPC2030 IF (N.LE.N30) GO TO 105 MAPC2035 WRITE (LOUT,99964) N30, N MAPC2040 99964 FORMAT (47H IN THIS VERSION OF MAPCLUS, THE MAXIMUM NUMBER, MAPC2045 * 14H OF STIMULI IS, I3, 24H, BUT YOU HAVE SPECIFIED, I3, 1H.// MAPC2050 * 22H EXECUTION TERMINATES.) MAPC2055 STOP MAPC2060 105 WRITE (LOUT,99963) TITL, N, ISWIT, IPRN, FMT MAPC2065 99963 FORMAT (//1H1, 25HTITLE OF INPUT DATA SET: , 72A1///10H N ISWI,MAPC2070 * 6HT LABL/1X, I3, 2X, I3, 5X, I1//28H SPECIFIED FORMAT OF THE INP,MAPC2075 * 11HUT DATA IS://1X, 72A1/) MAPC2080 IF (LABL.EQ.0) WRITE (LOUT,99962) MAPC2085 99962 FORMAT (//1X, 46HTHE PROGRAM WILL GENERATE NUMERICAL LABELS FOR, MAPC2090 * 13H THE STIMULI.) MAPC2095 IF (LABL.EQ.1) WRITE (LOUT,99961) MAPC2100 99961 FORMAT (//1X, 48HUSER HAS DECLARED THAT LABELS ARE TO BE READ IN ,MAPC2105 * 33HIMMEDIATELY AFTER NUMERICAL DATA.) MAPC2110 LD2 = LD1 + 1 MAPC2115 ALD1 = LD1 MAPC2120 BN = N MAPC2125 BIN = 1.0/BN MAPC2130 N1 = N - 1 MAPC2135 AN1 = N1 MAPC2140 NDO = N*N1/2 MAPC2145 AB = NDO MAPC2150 ABNV = 1.0/AB MAPC2155 AMN = 1.0E20 MAPC2160 AMX = -1.0E20 MAPC2165 C MAPC2170 C READ IN LOWERHALFMATRIX BY ROWS MAPC2175 C MAPC2180 DO 115 I=2,N MAPC2185 I1 = I - 1 MAPC2190 READ (LREAD,FMT) (DATA(IV,I),IV=1,I1) MAPC2195 DO 110 IV=1,I1 MAPC2200 DT = DATA(IV,I) MAPC2205 IF (DT.GT.AMX) AMX = DT MAPC2210 IF (DT.LT.AMN) AMN = DT MAPC2215 110 CONTINUE MAPC2220 115 CONTINUE MAPC2225 C MAPC2230 C MAPC2235 C SET UP STIMULUS LABELS OF THE INTEGERS 1,...,N. THESE LABELS MAPC2240 C ARE CONSTRUCTED EVEN IF THE USER SUPPLIES LABELS, SO THAT IN MAPC2245 C THE EVENT OF AN END OF FILE WHILE ATTEMPTING TO READ THE LATTER MAPC2250 C LABELS, THE PROGRAM STILL HAS USABLE LABELS. MAPC2255 C MAPC2260 DO 120 I=1,N MAPC2265 ILAB(I,1) = IBL MAPC2270 ILAB(I,6) = IBL MAPC2275 ILAB(I,2) = IL MAPC2280 ILAB(I,5) = IR MAPC2285 K4 = MOD(I,10) MAPC2290 ILAB(I,4) = IDIG(K4+1) MAPC2295 K4 = 1 + (I-K4)/10 MAPC2300 ILAB(I,3) = IDIG(K4) MAPC2305 120 CONTINUE MAPC2310 C MAPC2315 C THE DATA WILL SOON BE IN THE LOWER TRIANGULAR HALF OF THE ARRAY MAPC2320 C DATA(I,J). THE UPPER HALF WILL BE USED FOR RESIDUALS. MAPC2325 C MAPC2330 C NOW TRANSFORM THEM TO BE SIMILARITIES (NOT DISSIMILARITIES) MAPC2335 C ON THE INTERVAL [0,1]. MAPC2340 C MAPC2345 IF (LABL.NE.1) GO TO 125 MAPC2350 READ (LREAD,99960,END=493) ((ILAB(I,J), J = 1, 6), I= 1, N) MAPC2355 99960 FORMAT (72A1) MAPC2360 GO TO 125 MAPC2365 493 WRITE (LOUT,99959) MAPC2370 99959 FORMAT (///42H ********PROGRAM READ AN END-OF-FILE WHILE, MAPC2375 * 40H ATTEMPTING TO READ USER-SUPPLIED LABELS/9X, 12HWILL INSTEAD, MAPC2380 * 44H USE INTEGER LABELING. EXECUTION CONTINUES.//) MAPC2385 125 IF (ISWIT.GT.0) GO TO 130 MAPC2390 BMN = AMN MAPC2395 C3 = -1.0 MAPC2400 GO TO 135 MAPC2405 130 BMN = AMX MAPC2410 C3 = 1.0 MAPC2415 135 BMX = AMX - AMN MAPC2420 IF (BMX.LT.1.0E-6) GO TO 530 MAPC2425 C1 = 1.0/BMX MAPC2430 C2 = BMN*C1 MAPC2435 DO 145 I=2,N MAPC2440 I1 = I - 1 MAPC2445 DO 140 J=1,I1 MAPC2450 DATA(I,J) = C3*C1*(BMN-DATA(J,I)) MAPC2455 140 CONTINUE MAPC2460 145 CONTINUE MAPC2465 IF (ISWIT.LE.0) GO TO 150 MAPC2470 WRITE (LOUT,99958) C2, C1 MAPC2475 99958 FORMAT (/42H THE INPUT DATA ARE DISSIMILARITIES WHICH , 7HHAVE BE,MAPC2480 * 35HEN LINEARLY TRANSFORMED AS FOLLOWS://40X, 10HY =, MAPC2485 * F10.5, 3H - , F10.5, 13H(Y )/41X, 7HNEW SIM, 27X, MAPC2490 * 10HRAW DISSIM) MAPC2495 GO TO 155 MAPC2500 150 WRITE (LOUT,99957) C1, C2 MAPC2505 99957 FORMAT (/39H THE INPUT DATA ARE SIMILARITIES WHICH , 9HHAVE BEEN, MAPC2510 * 33H LINEARLY TRANSFORMED AS FOLLOWS://40X, 10HY =, F10.5, MAPC2515 * 10H(Y ), 3H - , F10.5/41X, 7HNEW SIM, 14X, 7HRAW SIM) MAPC2520 155 IF (IPRN.NE.1) GO TO 170 MAPC2525 WRITE (LOUT,99956) MAPC2530 99956 FORMAT (44H1 NO TRANSFORMED RAW SUBSCRIPTS, 6H ,MAPC2535 * 9H LABELS/49H VALUES DATA OF STIMULI ,MAPC2540 * 12H OF STIMULI/44H I J//MAPC2545 * ) MAPC2550 K = 0 MAPC2555 DO 165 I=2,N MAPC2560 I1 = I - 1 MAPC2565 DO 160 J=1,I1 MAPC2570 K = K + 1 MAPC2575 IF (MOD(K,53).EQ.0) WRITE (LOUT,99956) MAPC2580 WRITE (LOUT,99955) K, DATA(I,J), DATA(J,I), I, J, MAPC2585 * (ILAB(I,L),L=1,6), (ILAB(J,L),L=1,6) MAPC2590 160 CONTINUE MAPC2595 165 CONTINUE MAPC2600 99955 FORMAT (1X, I5, F12.4, F14.4, 2I6, 4X, 6A1, 4X, 6A1) MAPC2605 C MAPC2610 170 DO 180 I=1,N MAPC2615 DO 175 IW=1,LD1 MAPC2620 P(I,IW) = 0.0 MAPC2625 175 CONTINUE MAPC2630 P(I,LD2) = 1.0 MAPC2635 180 CONTINUE MAPC2640 VARSIM = 0.0 MAPC2645 STOT = 0.0 MAPC2650 DO 190 I=2,N MAPC2655 I1 = I - 1 MAPC2660 DO 185 J=1,I1 MAPC2665 STOT = STOT + DATA(I,J) MAPC2670 185 CONTINUE MAPC2675 190 CONTINUE MAPC2680 STOTAV = STOT*ABNV MAPC2685 DO 200 I=2,N MAPC2690 I1 = I - 1 MAPC2695 DO 195 J=1,I1 MAPC2700 VARSIM = VARSIM + (STOTAV-DATA(I,J))**2 MAPC2705 195 CONTINUE MAPC2710 200 CONTINUE MAPC2715 DO 205 IW=1,LD2 MAPC2720 RISE(IW) = 0.0 MAPC2725 205 CONTINUE MAPC2730 IF (ICON.GT.0) CALL INICON(0, 0, ICON) MAPC2735 IF (ICON.LE.0 .AND. M15.GT.0) GO TO 215 MAPC2740 IF (ICON.GT.0 .OR. M15.NE.0) GO TO 210 MAPC2745 WRITE (LOUT,99954) ICON MAPC2750 99954 FORMAT (//45H YOU HAVE SPECIFIED THE OPTION FOR REGRESSION, MAPC2755 * 27H ONLY (M15 = 0 ITERATIONS),/30H BUT YOU HAVE DECLARED AN INFE,MAPC2760 * 23HASIBLE LOGICAL DEVICE #, 1H(, I3, 23H) AS THE VALUE OF ICON MAPC2765 * /46H ON THE PARAMETER CARD. EXECUTION TERMINATES.//) MAPC2770 STOP MAPC2775 C MAPC2780 C MAPC2785 C MAPC2790 C MAPC2795 C MAPC2800 C THIS IS THE OPTION FOR CALLING THE REGRESSION ROUTINE FOR A MAPC2805 C USER-SUPPLIED INIT. CONF., WHILE BYPASSING THE ITERATIVE MAPC2810 C FITTING PROCEDURE. MAPC2815 C MAPC2820 210 WRITE (LOUT,99944) TITL MAPC2825 NST = 1 MAPC2830 DENS(1) = 0.0 MAPC2835 KC = 1 MAPC2840 CALL REGR(0, 0, 0) MAPC2845 C MAPC2850 C THE FIRST ARGUMENT IN THE FOLLOWING CALL TO PRINTD IS 1 SO THAT MAPC2855 C DENSITY CAN BE COMPUTED AND PRINTED. MAPC2860 C MAPC2865 WRITE (LOUT,99942) MAPC2870 CALL PRINTD(1, 0) MAPC2875 C MAPC2880 C POSSIBLE BRANCH TO COMBINATORIAL OPTIMIZATION MAPC2885 C MAPC2890 IF (M15) 465, 520, 215 MAPC2895 215 DLAM = ALF0 MAPC2900 DO 220 IW=1,LD1 MAPC2905 ALFRAY(IW) = ALPHAI MAPC2910 BETRAY(IW) = BETAI MAPC2915 220 CONTINUE MAPC2920 C MAPC2925 C MAPC2930 C MAPC2935 C MAPC2940 C DO LOOP FOR "MAJOR" ITERATIONS MAPC2945 C MAPC2950 DO 450 MST=1,M15 MAPC2955 C MAPC2960 C MAPC2965 M24 = MST MAPC2970 BLMX(MST) = -1.0E20 MAPC2975 BLMN(MST) = 1.0E20 MAPC2980 IADJ(MST) = 0 MAPC2985 DENS(MST) = 0. MAPC2990 MP(MST) = 0 MAPC2995 NW(MST) = 0 MAPC3000 SVAF(MST) = 0.0 MAPC3005 MST1 = MST - 1 MAPC3010 KMOD = MOD(MST,2) MAPC3015 IF (IPRE.NE.1) GO TO 225 MAPC3020 C MAPC3025 C MECHANISM TO ALLOW A RESTART WITH DE NOVO ITERATIONS: MAPC3030 C MAPC3035 I = 0 MAPC3040 IF (ITPRE.LE.0) I = 1 MAPC3045 ITPRE1 = MST - I MAPC3050 ITPRE = ITPRE1 - 1 MAPC3055 C MAPC3060 C NO MORE PRE-ITERATIONS MAPC3065 C MAPC3070 IF (IPRE.EQ.1 .AND. KTPRE.GE.0 .AND. MST1.GT.0) WRITE MAPC3075 * (LOUT,99953) BLMX(MST1), BCRIT MAPC3080 99953 FORMAT (///40H SINCE THE MAXIMUM BLOS ACROSS SUBSETS (, F10.8, MAPC3085 * 53H) IS LESS THAN OR EQUAL TO THE STOPPING CRITERION OF , MAPC3090 * 8HBLCRIT =, F12.8/14H THERE WILL BE, 24H NO MORE PRE-POLISHING I,MAPC3095 * 10HTERATIONS.///) MAPC3100 C MAPC3105 C TAKING CARE OF NUMBER OF INNER ITERATIONS CONSTITUTING AN OUTER MAPC3110 C ITERATION, DEPENDING ON THE NATURE OF CURRENT MAJOR ITERATION. MAPC3115 C MAPC3120 225 IF (MST-ITPRE1) 230, 235, 240 MAPC3125 230 L15 = IUPPER MAPC3130 K15 = IUPPER MAPC3135 GO TO 250 MAPC3140 235 L15 = JPOST + 30 MAPC3145 WRITE (LOUT,99952) MST MAPC3150 99952 FORMAT (/49H 30 EXTRA INNER ITERATIONS FOR THE FIRST MAJOR, MAPC3155 * 12H ITERATION (, I3, 29H) WITH POLISHING (BETA = 1.0)//) MAPC3160 GO TO 245 MAPC3165 240 L15 = JPOST MAPC3170 245 K15 = JPOST MAPC3175 C MAPC3180 C NO EFFECTIVE LIMIT ON ADJUST CALLS FOR POST-POLISHING ITERATIONS MAPC3185 JLIM = 1000 MAPC3190 C MAPC3195 C MAPC3200 C MAPC3205 C THE FOLLOWING DO LOOP SETS UP THE OUTER ITERATIONS. MAPC3210 C MAPC3215 250 DO 440 IIW=1,LD1 MAPC3220 C MAPC3225 C KC DETERMINES PRINTING AT THE END OF A MAJOR ITERATION MAPC3230 C MAPC3235 KC = 0 MAPC3240 IF (IIW.EQ.LD1) KC = 1 MAPC3245 C MAPC3250 C PROCEDURE FOR REVERSING THE ORDER IN WHICH SUBSETS MAPC3255 C (AND THEIR WEIGHTS) ARE COMPUTED, FOR EVEN-NUMBERED MAPC3260 C MAJOR ITERATIONS. (ALSO USED BELOW IN COMBINATORIAL MAPC3265 C OPTIMIZATION ITERATIONS). MAPC3270 C MAPC3275 IF (KMOD.EQ.1) GO TO 255 MAPC3280 IW = LD2 - IIW MAPC3285 GO TO 260 MAPC3290 255 IW = IIW MAPC3295 C MAPC3300 C USER-SUPPLIED TRIAL STEP SIZE FOR FIRST INNER ITERATION MAPC3305 260 ALF0 = DLAM MAPC3310 IFLOP = 2 MAPC3315 OBSLF(2) = 1.0E20 MAPC3320 ASTEP(2) = 1.0E20 MAPC3325 RCORR = 0.0 MAPC3330 C MAPC3335 C GET THE RESIDUALS AND CENTER THEM AT THE BEGINNING OF OUTER IT. MAPC3340 C MAPC3345 CALL RESID(IW) MAPC3350 JADJ = 0 MAPC3355 C MAPC3360 C MAPC3365 C MAPC3370 C NOW SET UP THE INNER ITERATIONS. MAPC3375 C MAPC3380 C MAPC3385 C MAPC3390 DO 380 KL1=1,L15 MAPC3395 KL = KL1 MAPC3400 IFLIP = IFLOP MAPC3405 IFLOP = 3 - IFLOP MAPC3410 GRSUM = 0.0 MAPC3415 GRMEN(IFLOP) = 0.0 MAPC3420 C MAPC3425 C ON THE VERY FIRST POLISHING ITERATION, CALL ADJUST MAPC3430 C ON EVERY THIRD INNER ITERATION DURING THE EXTRA MAPC3435 C INNER ITERATIONS MAPC3440 C MAPC3445 IF (KL.GT.K15 .AND. MOD(KL,3).EQ.0 .AND. ALFRAY(IW).GT.1.0E-6) MAPC3450 * CALL ADJUST(MST, IW) MAPC3455 C MAPC3460 C NOW GET THE INITIAL CONFIGURATION FOR (ONLY) THE VERY FIRST MAPC3465 C INNER ITERATION OF SUBSET IW (DURING THE FIRST MAJOR ITERATION). MAPC3470 C MAPC3475 IF (KL+MST.EQ.2 .AND. ICON.LE.0 .AND. ITPRE.GT.0) CALL MAPC3480 * INICON(MST, IW, ICON) MAPC3485 IF (MST.LE.ITPRE .OR. KL.NE.1) GO TO 265 MAPC3490 C MAPC3495 C MAPC3500 C THE "DE NOVO" APPROACH TO COMPUTING P(I) AFTER INITIAL POLISHING MAPC3505 C MAPC3510 CALL INICON(MST, IW, 0) MAPC3515 ALFRAY(IW) = ALPHAI MAPC3520 BETRAY(IW) = BETAI MAPC3525 C MAPC3530 C MAPC3535 C MAPC3540 265 BCONST = BSUM(IW) MAPC3545 C MAPC3550 C MAPC3555 C NOW ESTIMATE THE VALUE OF RISE(IW) AND THE MAPC3560 C REGRESSION CONSTANT, AND THEN TRANSLATE THE RESIDUALS. MAPC3565 C MAPC3570 CALL COMPW(IW, 0) MAPC3575 C MAPC3580 C MAPC3585 C COMPUTE THE A PART OF THE LOSS FUNCTION (USED IN ESTIMATION MAPC3590 C OF OPTIMAL LAMBDA) AND THE B PART (EARLIER) VIA FUNCTION CALLS. MAPC3595 C (THE LATTER GIVES COEFFICIENT C FOR THE QUADRATIC INTERPOLATION MAPC3600 C FUNCTION.) MAPC3605 C MAPC3610 IF (KL.GT.1) GO TO 270 MAPC3615 CONST = ALFRAY(IW)*ASUM(IW) + BETRAY(IW)*BCONST MAPC3620 GO TO 275 MAPC3625 270 CONST = OBSLF(IFLIP) MAPC3630 C MAPC3635 C THEN EVALUATE THE PART OF THE EXPRESSIONS FOR EACH I = 1, N MAPC3640 C MAPC3645 275 CALL CDAPI(IW) MAPC3650 C MAPC3655 C MAPC3660 C NOTE THAT DB IS AVAILABLE AS A FUNCTION CALL MAPC3665 C MAPC3670 DO 280 KLM=1,N MAPC3675 GL(KLM,IFLOP) = BETRAY(IW)*DB(KLM,IW) + ALFRAY(IW)*DAPI(KLM) MAPC3680 GLA = GL(KLM,IFLOP) MAPC3685 IF (ABS(GLA).LE.1.0E-12) GO TO 280 MAPC3690 GRMEN(IFLOP) = GRMEN(IFLOP) + GLA MAPC3695 GRSUM = GRSUM + GLA**2 MAPC3700 280 CONTINUE MAPC3705 GRLEN(IFLOP) = GRSUM MAPC3710 GRSUM = SQRT(GRSUM) MAPC3715 IF (GRSUM.LE.1.0E-12) GO TO 385 MAPC3720 IF (KL.LE.1) GO TO 290 MAPC3725 RCORR = 0.0 MAPC3730 C = 0. MAPC3735 D = 0. MAPC3740 DO 285 I=1,N MAPC3745 C = C + GL(I,1) MAPC3750 D = D + GL(I,2) MAPC3755 RCORR = RCORR + GL(I,1)*GL(I,2) MAPC3760 285 CONTINUE MAPC3765 RCORR = (BN*RCORR-C*D)/(SQRT((BN*GRLEN(1)-GRMEN(1)**2)*(BN* MAPC3770 * GRLEN(2)-GRMEN(2)**2))) MAPC3775 290 CLAM = ALF0/GRSUM MAPC3780 ILOOP = 0 MAPC3785 DO 295 KLM=1,N MAPC3790 P(KLM,IW) = P(KLM,IW) - CLAM*GL(KLM,IFLOP) MAPC3795 295 CONTINUE MAPC3800 IF (KP(MST).EQ.1 .AND. KL.EQ.1) WRITE (LOUT,99951) MST, KL, IW, MAPC3805 * RISE(IW), ALFRAY(IW), BETRAY(IW), CLAM, ACRIT, GRSUM, RCORR MAPC3810 99951 FORMAT (6H MAJOR, I3, 7H, INNER, I4, 4H; W(, I2, 2H)=, F9.5, MAPC3815 * 7H ALPHA=, F6.4, 6H BETA=, F6.4, 8H LAMBDA=, F12.4, 7H ACRIT=, MAPC3820 * F9.7, 7H GRSUM=, F11.5, 5H COR=, F8.4) MAPC3825 C MAPC3830 C EVALUATE THE LOSS FUNCTION NOW THAT A STEP OF LENGTH ALF0 HAS MAPC3835 C BEEN TAKEN (USED FOR OPTIMAL LAMBDA). MAPC3840 C MAPC3845 FALPHA = ALFRAY(IW)*ASUM(IW) + BETRAY(IW)*BSUM(IW) MAPC3850 C MAPC3855 C THE FOLLOWING TWO STATEMENTS GIVE RISE TO BENIGN UNDERFLOW AND MAPC3860 C EXP DIAGNOSTICS MAPC3865 C MAPC3870 300 IF (ABS(ALF0).LE.1.0E-12) GO TO 305 MAPC3875 AA = (FALPHA-CONST+ALF0*GRSUM)/ALF0**2 MAPC3880 ALFMIN = GRSUM*0.5/AA MAPC3885 IF (AA.GT.0.001) GO TO 315 MAPC3890 305 IF (ABS(RCORR).GT.1.0E-12) GO TO 310 MAPC3895 ALFMIN = ALF0 MAPC3900 GO TO 315 MAPC3905 310 ALFMIN = ALF0*2.0**RCORR MAPC3910 315 ASTEP(IFLOP) = (ALFMIN-ALF0)/GRSUM MAPC3915 IF (AA.GT.0.001) GO TO 320 MAPC3920 IF (KP(MST).EQ.1) WRITE (LOUT,99950) AA, ALF0 MAPC3925 99950 FORMAT (27H QUADRATIC COEFFICIENT A =, F16.5, 15H; SUBSTITUTED V,MAPC3930 * 14HALUE OF ALF0 =, F11.5) MAPC3935 IF (KL.EQ.1) GO TO 330 MAPC3940 IF (KP(MST).EQ.1) WRITE (LOUT,99949) ALF0, ASTEP(IFLOP), ALFMIN MAPC3945 99949 FORMAT (34H USED COS MULTIPLIER OF OLD ALF0 =, F10.5, 9H TO GET N,MAPC3950 * 17HORMALIZED ASTEP =, F10.5, 4X, 13H NEW ALFMIN =, F10.5) MAPC3955 C MAPC3960 C IN CASE OVERFLOW CAUSES HAVOC ... MAPC3965 C MAPC3970 320 IF (ASTEP(IFLOP).GT.1.0E10) GO TO 385 MAPC3975 C MAPC3980 C MAPC3985 C NOW USE INTERMEDIATE STEP SIZE THAT SHOULD MINIMIZE THE LOSS MAPC3990 C FUNCTION. MAPC3995 C MAPC4000 DO 325 KLM=1,N MAPC4005 P(KLM,IW) = P(KLM,IW) - ASTEP(IFLOP)*GL(KLM,IFLOP) MAPC4010 325 CONTINUE MAPC4015 C MAPC4020 C COMPUTE THE OBSERVED VALUE OF THE LOSS FUNCTION MAPC4025 C MAPC4030 330 ALOS = ASUM(IW) MAPC4035 BCONST = BSUM(IW) MAPC4040 OBSLF(IFLOP) = ALFRAY(IW)*ALOS + BETRAY(IW)*BCONST MAPC4045 C MAPC4050 C GET THE ESTIMATED VALUE OF THE LOSS FUNCTION AT ALPHA MIN MAPC4055 C MAPC4060 PRELF = ALFMIN*(AA*ALFMIN-GRSUM) + CONST MAPC4065 C MAPC4070 C STRATEGY WHEN OBSLF IS GETTING WORSE: MAPC4075 C MAPC4080 IF (OBSLF(IFLOP).LE.OBSLF(IFLIP)) GO TO 335 MAPC4085 C CHECK TO SEE THAT PROCEDURE IS NOT CAUGHT IN A LOOP MAPC4090 ILOOP = ILOOP + 1 MAPC4095 IF (ILOOP.GT.3) GO TO 335 MAPC4100 ALF0 = ALFMIN MAPC4105 FALPHA = OBSLF(IFLOP) MAPC4110 GO TO 300 MAPC4115 C MAPC4120 C MAPC4125 C GET LENGTH OF "RELATIVE GRADIENT", WHICH IS THE PRODUCT OF MAPC4130 C LENGTH OF GRADIENT TIMES LENGTH OF P VECTOR. MAPC4135 C MAPC4140 335 RELGR = 0.0 MAPC4145 DO 340 KLM=1,N MAPC4150 RELGR = RELGR + P(KLM,IW)**2 MAPC4155 340 CONTINUE MAPC4160 RELGR = GRSUM*SQRT(RELGR) MAPC4165 IF (KP(MST).EQ.1 .AND. KL.EQ.1) WRITE (LOUT,99948) PRELF, MAPC4170 * OBSLF(IFLOP), ASTEP(IFLOP), RELGR, DCON, ALOS, BCONST MAPC4175 99948 FORMAT (14H F(ALPHA MIN)=, F11.5, 7H OBSLF=, F11.5, 7H ASTEP=, MAPC4180 * F12.5, 8H RELGR =, F11.4, 7H DCON =, F11.4, 7H ALOS =, F10.4, MAPC4185 * 6H BLOS=, F10.6) MAPC4190 C MAPC4195 C PUT IN NEW VALUE FOR LIMIT OF OPTIMAL STEP SIZE FOR MAPC4200 C NEXT INNER ITERATION. MAPC4205 C MAPC4210 ALF0 = ALFMIN MAPC4215 C MAPC4220 C CHECK FOR CONVERGENCE, UNLESS THIS IS THE LAST INNER ITERATION. MAPC4225 C MAPC4230 C IF THIS IS THE LAST (L15-TH) INNER ITERATION, NO TESTS ARE MAPC4235 C NECESSARY. MAPC4240 C MAPC4245 C MAPC4250 C MAPC4255 IF (KL.GT.K15) GO TO 360 MAPC4260 C MAPC4265 C MAPC4270 C IF THE "RELATIVE" GRADIENT IS NOT SHORT ENOUGH, DO NOT MAKE ANY MAPC4275 C CHANGES. MAPC4280 C MAPC4285 IF (RELGR.GT.ACRIT) GO TO 345 MAPC4290 C MAPC4295 C IF BLOS IS SMALL ENOUGH, AND THIS IS A PRE-ITERATION, MAPC4300 C TERMINATE THIS OUTER ITERATION. MAPC4305 C MAPC4310 IF (BCONST.GT.BLCRIT) GO TO 355 MAPC4315 C MAPC4320 IF (MST.GT.ITPRE) GO TO 350 MAPC4325 IF (KP(MST).EQ.1) WRITE (LOUT,99947) KL MAPC4330 99947 FORMAT (46H RELGR AND BCONST CRITERIA ARE BOTH SATISFIED., MAPC4335 * 46H THIS OUTER ITERATION OF A MAJOR PRE-ITERATION, 10H TERMINATE,MAPC4340 * 7HS AFTER, I3, 11H INNER ITS.) MAPC4345 GO TO 385 MAPC4350 C MAPC4355 345 IF (ABS(ASTEP(1)).LE.1.0E-5 .AND. ABS(ASTEP(2)).LE.1.0E-5) GO TO MAPC4360 * 355 MAPC4365 GO TO 365 MAPC4370 C MAPC4375 C MAPC4380 C SINCE THIS IS NOT A PRE-ITERATION, POSSIBLY CALL ADJUST AND MAPC4385 C CONTINUE THE INNER ITERATIONS. MAPC4390 C MAPC4395 350 IF (KP(MST).EQ.1 .AND. ALFRAY(IW).GT.1.0E-6) WRITE (LOUT,99946) MAPC4400 * RELGR, ACRIT, ALFRAY(IW), BETRAY(IW), BCONST, BLCRIT MAPC4405 99946 FORMAT (7H RELGR=, F10.8, 10H.LE.ACRIT=, F10.8, 2X, 10H CALLED AD,MAPC4410 * 12HJUST: ALPHA=, F9.6, 6H BETA=, F9.6, 8H BCONST=, F10.8, MAPC4415 * 8H BLCRIT=, F10.8) MAPC4420 355 IF (JADJ.EQ.JLIM) GO TO 385 MAPC4425 IF (ALFRAY(IW).GT.1.0E-6) CALL ADJUST(MST, IW) MAPC4430 C MAPC4435 C MAPC4440 C RETURN TO THE ORIGINAL (USER-SUPPLIED) ALPHA STEP SIZE AFTER MAPC4445 C CONVERGENCE HAS OCCURRED. MAPC4450 C MAPC4455 360 ALF0 = DLAM MAPC4460 C MAPC4465 C RE-CENTER THE RESIDUALS. MAPC4470 C MAPC4475 365 DO 375 I=2,N MAPC4480 I1 = I - 1 MAPC4485 DO 370 J=1,I1 MAPC4490 DATA(J,I) = DATA(J,I) + WCON MAPC4495 370 CONTINUE MAPC4500 375 CONTINUE MAPC4505 380 CONTINUE MAPC4510 C MAPC4515 385 IF (KP(MST).EQ.0) GO TO 390 MAPC4520 WRITE (LOUT,99951) MST, KL, IW, RISE(IW), ALFRAY(IW), BETRAY(IW), MAPC4525 * CLAM, ACRIT, GRSUM, RCORR MAPC4530 WRITE (LOUT,99948) PRELF, OBSLF(IFLOP), ASTEP(IFLOP), RELGR, MAPC4535 * DCON, ALOS, BCONST MAPC4540 C MAPC4545 C GET THE VAF BEFORE AS WELL AS AFTER POLISHING AT END OF OUTER MAPC4550 C ITERATION. MAPC4555 C MAPC4560 390 IF (KP(MST).EQ.1) CALL VARIAN(MST, 1, IW) MAPC4565 IF (BCONST.GT.BLMX(MST)) BLMX(MST) = BCONST MAPC4570 IF (BCONST.LT.BLMN(MST)) BLMN(MST) = BCONST MAPC4575 C MAPC4580 C MAPC4585 C SINCE BLOS IS USING PAIRWISE PRODUCTS OF THE P(I), THE SIGNS MAPC4590 C ARE POTENTIALLY REVERSED. SET THE SIGN OF THE LARGEST MAPC4595 C OF THE ABS(P(I)) TO BE POSITIVE, AND MULTIPLY THE REST BY MAPC4600 C THE SAME SIGN USED TO MAKE THE LARGEST POSITIVE. MAPC4605 C MAPC4610 I3 = 1 MAPC4615 AM = ABS(P(1,IW)) MAPC4620 DO 395 I=2,N MAPC4625 AN = ABS(P(I,IW)) MAPC4630 IF (AN.LE.AM) GO TO 395 MAPC4635 AM = AN MAPC4640 I3 = I MAPC4645 395 CONTINUE MAPC4650 IF (P(I3,IW).GE.0.0) GO TO 410 MAPC4655 DO 400 I=1,N MAPC4660 P(I,IW) = -P(I,IW) MAPC4665 400 CONTINUE MAPC4670 B3 = 0.0 MAPC4675 DO 405 KLM=1,N MAPC4680 B3 = B3 + P(KLM,IW) MAPC4685 405 CONTINUE MAPC4690 B3 = B3*BIN MAPC4695 B4 = TM/B3 MAPC4700 410 IF (MST.LE.ITPRE) GO TO 415 MAPC4705 IF (KP(MST).EQ.1) WRITE (LOUT,99945) IW, B3, TM, B4, MAPC4710 * (P(KLM,IW),KLM=1,N) MAPC4715 99945 FORMAT (17H PRE-POLISHED P (, I2, 15H): (MEAN P(J) =, F12.6, 1H), MAPC4720 * 17H MEAN(P(I)P(J)) =, F12.6, 8H RATIO =, F12.6, 11H; USED .707/ MAPC4725 * 64(1X, 16F8.4/)) MAPC4730 C MAPC4735 C FORCE THE T VECTOR TO BE BINARY. MAPC4740 C MAPC4745 CALL POLISH(MST, IW) MAPC4750 C MAPC4755 C IF SOME OF THE FOLLOWING WRITE STATEMENTS WERE DELETED, MAPC4760 C THE FOLLOWING CALL TO COMPW COULD BE ELIMINATED, SINCE MAPC4765 C THE SUBSEQUENT CALL TO REGR TAKES CARE OF ESTIMATING ALL MAPC4770 C THE WEIGHTS SIMULTANEOUSLY. MAPC4775 C MAPC4780 CALL COMPW(IW, 1) MAPC4785 C MAPC4790 C MAPC4795 C MAPC4800 415 IF (MST+IW.GT.2) CALL REGR(MST, IW, MST) MAPC4805 IF (MST.GT.ITPRE) GO TO 420 MAPC4810 IF (KP(MST).NE.1) GO TO 440 MAPC4815 IF (KC.EQ.1) WRITE (LOUT,99944) TITL MAPC4820 99944 FORMAT (/1X, 72A1) MAPC4825 WRITE (LOUT,99943) MST, IW, RISE(IW), IW, B3, TM, B4, MAPC4830 * (P(KLM,IW),KLM=1,N) MAPC4835 99943 FORMAT (/5H IT #, I3, 3H W(, I2, 2H)=, F12.6, 12H P MATRIX (, MAPC4840 * I2, 2H):, 6X, 13H( MEAN P(J) =, F12.6, 1H), 17H MEAN(P(I)P(J)) =,MAPC4845 * F12.6, 8H RATIO =, F12.6/64(1X, 16F8.4/)) MAPC4850 GO TO 440 MAPC4855 420 IF (KC.NE.1) GO TO 425 MAPC4860 99942 FORMAT (////) MAPC4865 GO TO 430 MAPC4870 425 IF (KP(MST).NE.1) GO TO 435 MAPC4875 430 CALL PRINTD(MST, IW) MAPC4880 435 IF (IPUN.EQ.1 .AND. KC.EQ.1) CALL PUN(MST) MAPC4885 440 CONTINUE MAPC4890 DO 445 IW=1,LD1 MAPC4895 IF (RISE(IW).LT.0.0) NW(MST) = NW(MST) + 1 MAPC4900 445 CONTINUE MAPC4905 WRITE (LOUT,99941) MST, IADJ(MST), NW(MST) MAPC4910 99941 FORMAT (/24H DURING MAJOR ITERATION , I4, 19H ADJUST WAS CALLED , MAPC4915 * I4, 22H TIMES AND THERE WERE , I4, 18H NEGATIVE WEIGHTS./////) MAPC4920 IF (IERR.NE.4) SVAF(MST) = VAF MAPC4925 C MAPC4930 C OPTION TO TERMINATE PRE-ITERATIONS, BASED ON VAF FROM MAPC4935 C ITERATIVE REGRESSION. MAPC4940 C MAPC4945 IF (BLMX(MST).LE.BCRIT .AND. MST.LE.ITPRE .AND. IPRE.NE.1) IPRE = MAPC4950 * 1 MAPC4955 C MAPC4960 IF (MST.LT.ITPRE1 .OR. MST1.LE.0) GO TO 450 MAPC4965 IF (SVAF(MST)-SVAF(MST1).GE.0.001 .OR. SVAF(MST)-SVAF(MST1) MAPC4970 * .LT.0.0 .OR. ABS(DENS(MST)-DENS(MST1)).GE.0.005) GO TO 450 MAPC4975 WRITE (LOUT,99940) MAPC4980 99940 FORMAT (///43H TERMINATED THIS SERIES OF MAJOR ITERATIONS, MAPC4985 * 34H BECAUSE OF NEGLIGIBLE IMPROVEMENT///) MAPC4990 GO TO 455 MAPC4995 450 CONTINUE MAPC5000 C MAPC5005 C MAPC5010 C MAPC5015 C MAPC5020 455 IF (M24.EQ.M15) WRITE (LOUT,99939) M24 MAPC5025 99939 FORMAT (//49H TERMINATED MAJOR ITERATIONS AFTER REACHING USER-, MAPC5030 * 19HSPECIFIED LIMIT OF , I3//) MAPC5035 WRITE (LOUT,99938) MAPC5040 99938 FORMAT (1H1) MAPC5045 WRITE (LOUT,99937) MAPC5050 99937 FORMAT (//50H SUMMARY OF ALTERNATING LEAST SQUARES COMPUTATION , MAPC5055 * 35HPRIOR TO COMBINATORIAL OPTIMIZATION//) MAPC5060 WRITE (LOUT,99944) TITL MAPC5065 M15 = MIN0(M24,M15) MAPC5070 WRITE (LOUT,99936) MAPC5075 99936 FORMAT (//20H VAF BY ITERATION, 19X, 13H DENSITY OF, 13X, MAPC5080 * 57HB-PART OF LOSS FUNCTION NUMBER OF PREMATURE/43X,MAPC5085 * 8H SUBSETS, 41H AT END OF OUTER ITS. , MAPC5090 * 30H ADJUSTMENTS POLISHING/41X, 13H(=0. PRIOR TO, 58X, MAPC5095 * 12H(=0 PRIOR TO/42X, 10HPOLISHING), 61X, 10HPOLISHING)) MAPC5100 WRITE (LOUT,99935) (MST,SVAF(MST),NW(MST),DENS(MST),BLMX(MST), MAPC5105 * BLMN(MST),IADJ(MST),MP(MST),MST=1,M15) MAPC5110 99935 FORMAT (//(2X, 4HVAF(, I3, 2H)=, F12.6, 2X, I3, 14H NEG WTS DEN,MAPC5115 * 3HS =, F10.5, 12H MAX BLOS =, F10.5, 12H MIN BLOS =, F10.7, MAPC5120 * 6H ADJ =, I4, 11H MID POL =, I5/)) MAPC5125 C MAPC5130 IF (IPRE.EQ.1 .OR. M24.GE.ITPRE1) GO TO 465 MAPC5135 WRITE (LOUT,99934) MAPC5140 99934 FORMAT (//47H SINCE THE SUBSETS HAVE NOT YET BEEN POLISHED, , MAPC5145 * 32HTHE P MATRIX IS NOT BINARY, AND /18H THERE WILL BE NO , MAPC5150 * 28H COMBINATORIAL OPTIMIZATION./////) MAPC5155 DO 460 IW=1,LD1 MAPC5160 WRITE (LOUT,99943) M24, IW, RISE(IW), IW, B3, TM, B4, MAPC5165 * (P(KLM,IW),KLM=1,N) MAPC5170 460 CONTINUE MAPC5175 GO TO 5 MAPC5180 C MAPC5185 C TIME FOR COMBINATORIAL OPTIMIZATION MAPC5190 C MAPC5195 465 IF (IERR.EQ.4) GO TO 470 MAPC5200 IF (M15.LE.0) GO TO 475 MAPC5205 IF (SVAF(M15).GT.0.0 .AND. SVAF(M15).LE.1.0) GO TO 475 MAPC5210 470 WRITE (LOUT,99933) MAPC5215 99933 FORMAT (///49H ******** SINCE LINEAR DEPENDENCY WAS PRESENT IN, MAPC5220 * 59H THE SOLUTION FROM THE ALTERNATING LEAST SQUARES PROCEDURE,/ MAPC5225 * 11X, 48HCOMBINATORIAL OPTIMIZATION WILL NOT BE EXECUTED.///) MAPC5230 GO TO 5 MAPC5235 475 KC = 1 MAPC5240 IXV = 1 MAPC5245 IC1 = 0 MAPC5250 DO 510 MST=1,15 MAPC5255 NST = MST MAPC5260 KMOD = MOD(NST,2) MAPC5265 IC0 = 0 MAPC5270 WRITE (LOUT,99932) NST MAPC5275 99932 FORMAT (///22H THIS IS ITERATION, I3, 18H OF COMBINATORIAL , MAPC5280 * 36HOPTIMIZATION, PHASE ONE (DOUBLETONS)///) MAPC5285 DO 490 IIW=1,LD1 MAPC5290 IF (KMOD.EQ.1) GO TO 480 MAPC5295 IW = LD2 - IIW MAPC5300 GO TO 485 MAPC5305 480 IW = IIW MAPC5310 485 CALL COMDUB(NST, IW) MAPC5315 IC0 = IC0 + ICOUN MAPC5320 490 CONTINUE MAPC5325 DENS(NST) = 0.0 MAPC5330 CALL PRINTD(NST, 0) MAPC5335 IF (IC0.EQ.LD1 .AND. IC1.EQ.LD1) GO TO 515 MAPC5340 IC1 = 0 MAPC5345 WRITE (LOUT,99931) NST MAPC5350 99931 FORMAT (///22H THIS IS ITERATION, I3, 18H OF COMBINATORIAL , MAPC5355 * 36HOPTIMIZATION, PHASE TWO (SINGLETONS)///) MAPC5360 DO 505 IIW=1,LD1 MAPC5365 IF (KMOD.EQ.1) GO TO 495 MAPC5370 IW = LD2 - IIW MAPC5375 GO TO 500 MAPC5380 495 IW = IIW MAPC5385 500 CALL COMSIN(NST, IW) MAPC5390 IC1 = IC1 + ICOUN MAPC5395 505 CONTINUE MAPC5400 DENS(NST) = 0.0 MAPC5405 CALL PRINTD(NST, 0) MAPC5410 IF (IC0.EQ.LD1 .AND. IC1.EQ.LD1) GO TO 515 MAPC5415 IXV = 0 MAPC5420 CALL REGR(0, 0, NST) MAPC5425 IXV = 1 MAPC5430 510 CONTINUE MAPC5435 515 IXV = 0 MAPC5440 CALL REGR(0, 0, NST) MAPC5445 IF (NST.EQ.15) WRITE (LOUT,99930) MAPC5450 99930 FORMAT (///42H PROGRAM REACHED LIMIT OF 15 ITERATIONS OF, MAPC5455 * 28H COMBINATORIAL OPTIMIZATION.///) MAPC5460 IF (NST.LT.15) WRITE (LOUT,99929) MAPC5465 99929 FORMAT (///42H SINCE THE PRECEDING ITERATION PRODUCED NO, MAPC5470 * 50H IMPROVEMENT IN VAF, COMBINATORIAL OPTIMIZATION IS, 7H TERMIN,MAPC5475 * 5HATED.///) MAPC5480 IF (IPUN.EQ.1) CALL PUN(NST) MAPC5485 C MAPC5490 C SORT THE SUBSETS BY WEIGHT AND PRINT THEM OUT. MAPC5495 C MAPC5500 520 DO 525 IW=1,LD1 MAPC5505 IZ(IW) = IW MAPC5510 A(IW) = RISE(IW) MAPC5515 525 CONTINUE MAPC5520 CALL SORT(LD1) MAPC5525 KC = 2 MAPC5530 WRITE (LOUT,99928) MAPC5535 99928 FORMAT (////45H1SUBSETS RANKED ACCORDING TO NUMERICAL WEIGHT///) MAPC5540 WRITE (LOUT,99944) TITL MAPC5545 CALL PRINTD(NST, IW) MAPC5550 CALL FINIS(BMN, C3, C1) MAPC5555 GO TO 5 MAPC5560 530 WRITE (LOUT,99927) MAPC5565 99927 FORMAT (///34H ******** TROUBLE EXIT **********//9X, 9H CHECK F,MAPC5570 * 21HORMAT OF INPUT DATA. ///) MAPC5575 551 CONTINUE MAPC5578 STOP MAPC5580 END MAPC5585 C FUNCTION IDATEZ(DUMMY) FINI 5 C CALL DATE(IDATEZ) FINI 10 C RETURN FINI 15 C END FINI 20 C$ FORTY FORM,XREF FINI 25 SUBROUTINE FINIS(BMN, C3, C1) FINI 30 DOUBLE PRECISION DCON FINI 35 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), FINI 40 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, FINI 45 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC FINI 50 COMMON /A6/ DCON, DATA(030,030), ALOS, IXV FINI 55 COMMON /A9/ ILAB(030,6), ALD1 FINI 60 WK = RISE(LD2) FINI 65 DO 15 I=1,N FINI 70 DO 10 J=1,I FINI 75 DATA(J,I) = 0. FINI 80 DO 5 K=1,LD1 FINI 85 DATA(J,I) = DATA(J,I) + RISE(K)*P(I,K)*P(J,K) FINI 90 5 CONTINUE FINI 95 DATA(J,I) = DATA(J,I) + WK FINI 100 10 CONTINUE FINI 105 15 CONTINUE FINI 110 WRITE (LOUT,99999) FINI 115 99999 FORMAT (1H1, 26H TRANSFORMED TRANSFORMED, 9X, 14HRAW UNTRANS,FINI 120 * 6HFORMED, 39H RESIDUAL SUBSCRIPTS LABELS/9H OBSE,FINI 125 * 17HRVED ESTIMATED, 36H OBSERVED ESTIMATED , FINI 130 * 35H OF STIMULI OF STIMULI/22H SIMILARITY SIMILA,FINI 135 * 32HRITY PROXIMITY PROXIMITY, 20X, 1HI, 5X, 1HJ, 6X, 1HI, FINI 140 * 8X, 1HJ//) FINI 145 C4 = 1.0/(C3*C1) FINI 150 K = 0 FINI 155 DO 30 I=1,N FINI 160 DO 25 J=1,I FINI 165 K = K + 1 FINI 170 IF (MOD(K,49).NE.0) GO TO 20 FINI 175 WRITE (LOUT,99997) FINI 180 WRITE (LOUT,99999) FINI 185 20 D1 = BMN - DATA(I,J)*C4 FINI 190 D2 = BMN - DATA(J,I)*C4 FINI 195 D3 = D1 - D2 FINI 200 IF (I.EQ.J) D1 = 0.0 FINI 205 WRITE (LOUT,99998) DATA(I,J), DATA(J,I), D1, D2, D3, I, J, FINI 210 * (ILAB(I,L),L=1,6), (ILAB(J,L),L=1,6) FINI 215 25 CONTINUE FINI 220 30 CONTINUE FINI 225 99998 FORMAT (2F12.4, 4X, 3F12.4, 5X, 2I6, 3X, 6A1, 3X, 6A1) FINI 230 WRITE (LOUT,99997) FINI 235 99997 FORMAT (//49H NOTE: SELF-SIMILARITY AND SELF-PROXIMITY VALUES, FINI 240 * 60H ARE ONLY PREDICTED BY, BUT NOT OBSERVED IN, THE INPUT DATA.) FINI 245 RETURN FINI 250 END FINI 255 C$ FORTY FORM HF00 5 FUNCTION HF(X) HF00 10 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), HF00 15 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, HF00 20 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC HF00 25 HF = X - ABNV*X*X HF00 30 RETURN HF00 35 END HF00 40 C$ FORTY FORM,XREF COMS 5 SUBROUTINE COMSIN(MST, IW) COMS 10 DOUBLE PRECISION DCON COMS 15 INTEGER TITL COMS 20 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), COMS 25 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, COMS 30 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC COMS 35 COMMON /A2/ G(030), KT1(10), KT2(10), TK(10), LP(15), IPUN, IERR, COMS 40 * LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), TRISE(025) COMS 45 COMMON /A6/ DCON, DATA(030,030), ALOS, IXV COMS 50 COMMON /A9/ ILAB(030,6), ALD1 COMS 55 ICOUN = 1 COMS 60 5 CALL RESID(IW) COMS 65 C COMS 70 C THAT GOT THE CENTERED RESIDUALS; NOW COMPUTE NUMERATOR (E) AND COMS 75 C DENOMINATOR (F) FOR V**2 COMS 80 C COMS 85 E = 0.0 COMS 90 F = 0.0 COMS 95 DO 15 I=2,N COMS 100 I1 = I - 1 COMS 105 PI = P(I,IW) COMS 110 DO 10 J=1,I1 COMS 115 PIJ = P(J,IW)*PI COMS 120 E = E + DATA(J,I)*PIJ COMS 125 F = F + PIJ COMS 130 10 CONTINUE COMS 135 15 CONTINUE COMS 140 VK2 = E*E/HF(F) COMS 145 99999 FORMAT (//47H VARIANCE IN THE RESIDUALS CURRENTLY ACCOUNTED , COMS 150 * 15H FOR BY SUBSET , I3, 2H =, F10.5) COMS 155 HTOT = 0.0 COMS 160 DO 20 I=1,N COMS 165 G(I) = 0. COMS 170 HTOT = HTOT + P(I,IW) COMS 175 20 CONTINUE COMS 180 IF (LP(MST).EQ.1) WRITE (LOUT,99999) IW, VK2 COMS 185 DO 30 I=1,N COMS 190 DO 25 J=1,N COMS 195 IF (J.EQ.I) GO TO 25 COMS 200 K1 = MIN0(I,J) COMS 205 K2 = MAX0(I,J) COMS 210 G(I) = G(I) + P(J,IW)*DATA(K1,K2) COMS 215 25 CONTINUE COMS 220 30 CONTINUE COMS 225 DO 35 K=1,LD2 COMS 230 TRISE(K) = RISE(K) COMS 235 35 CONTINUE COMS 240 IC10 = 0 COMS 245 RMAX = -1.0E20 COMS 250 DO 45 I=1,N COMS 255 PI2 = 2.0*P(I,IW) - 1.0 COMS 260 HT2 = HTOT - PI2 COMS 265 C COMS 270 C DO NOT CONSIDER DROPPING EITHER MEMBER OF A DYAD SUBSET COMS 275 C OR FORMING THE COMPLETE SUBSET. COMS 280 C COMS 285 IF (HT2.LT.1.98 .OR. HT2.GT.AN1) GO TO 45 COMS 290 R2 = 0.0 COMS 295 ANUM = E - PI2*G(I) COMS 300 IF (ABS(ANUM).LE.1.0E-12) GO TO 40 COMS 305 R2 = (ANUM**2)/HF(F-PI2*(HTOT-P(I,IW))) COMS 310 IF (R2.LE.RMAX) GO TO 40 COMS 315 C COMS 320 C UPDATE BEST V**2 COMS 325 C COMS 330 MR2 = I COMS 335 RMAX = R2 COMS 340 C COMS 345 C PRINT OUT (10 AT A TIME) THE ALLOWABLE SINGLETON CHANGES COMS 350 C AND RESULTING V**2 COMS 355 C COMS 360 40 CONTINUE COMS 365 IF (LP(MST).NE.1) GO TO 45 COMS 370 IC10 = IC10 + 1 COMS 375 KT1(IC10) = -FLOAT(I)*PI2 COMS 380 TK(IC10) = R2 COMS 385 IF ((IC10.LT.10 .AND. I.LT.N) .OR. (I.EQ.N .AND. IC10.EQ.0)) GO COMS 390 * TO 45 COMS 395 WRITE (LOUT,99998) (KT1(K),TK(K),K=1,IC10) COMS 400 99998 FORMAT (1X, 10(1X, 1H(, I3, F7.2, 1H))) COMS 405 IC10 = 0 COMS 410 45 CONTINUE COMS 415 IF (IC10.NE.0) WRITE (LOUT,99998) (KT1(K),TK(K),K=1,IC10) COMS 420 IF (RMAX.LE.VK2) RETURN COMS 425 P(MR2,IW) = 1.0 - P(MR2,IW) COMS 430 UVAF = TVAF COMS 435 CALL REGR(0, 0, MST) COMS 440 IF (IERR.EQ.4) GO TO 60 COMS 445 C COMS 450 C FIRST CHECK TO SEE IF THE WEIGHT OF ANY SUBSET WOULD COMS 455 C BECOME NEGATIVE. COMS 460 C COMS 465 DO 50 IK=1,LD1 COMS 470 LK = IK COMS 475 IF (RISE(IK).LE.0.0 .AND. TRISE(IK).GT.0.0) GO TO 55 COMS 480 50 CONTINUE COMS 485 GO TO 70 COMS 490 55 WRITE (LOUT,99997) (ILAB(MR2,I),I=1,6), IW, TVAF, UVAF, LK, COMS 495 * RISE(LK) COMS 500 99997 FORMAT (//19H REVERSING ELEMENT , 1X, 6A1, 11H OF SUBSET , I3, COMS 505 * 31H WOULD INCREASE OVERALL VAF TO , F10.6, 12H, REPLACING , COMS 510 * F10.6/55H BUT WAS NOT DONE SINCE THE CURRENTLY POSITIVE WEIGHT F,COMS 515 * 2HOR, 7H SUBSET, I3, 30H WOULD HAVE BECOME NEGATIVE (, F12.7, COMS 520 * 2H).) COMS 525 C COMS 530 C UNSUCCESSFUL CHANGE COMS 535 C COMS 540 60 P(MR2,IW) = 1.0 - P(MR2,IW) COMS 545 TVAF = UVAF COMS 550 DO 65 K=1,LD2 COMS 555 RISE(K) = TRISE(K) COMS 560 65 CONTINUE COMS 565 RETURN COMS 570 C COMS 575 C FOR A SUCCESSFUL CHANGE COMS 580 C COMS 585 70 ICOUN = 2 COMS 590 WRITE (LOUT,99996) (ILAB(MR2,I),I=1,6), IW, TVAF, UVAF COMS 595 99996 FORMAT (//19H REVERSING ELEMENT , 1X, 6A1, 11H OF SUBSET , I3, COMS 600 * 33H GAVE AN IMPROVED OVERALL VAF OF , F10.6, 11H REPLACING , COMS 605 * F10.6//) COMS 610 UVAF = TVAF COMS 615 C COMS 620 C NOW GO BACK AND LOOK FOR ANY OTHER SUCCESSFUL CHANGES WITHIN COMS 625 C SUBSET IW. COMS 630 C COMS 635 GO TO 5 COMS 640 END COMS 645 C$ FORTY FORM,XREF COMD 5 SUBROUTINE COMDUB(MST, IW) COMD 10 DOUBLE PRECISION DCON COMD 15 INTEGER TITL COMD 20 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), COMD 25 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, COMD 30 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC COMD 35 COMMON /A2/ G(030), KT1(10), KT2(10), TK(10), LP(15), IPUN, IERR, COMD 40 * LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), TRISE(025) COMD 45 COMMON /A6/ DCON, DATA(030,030), ALOS, IXV COMD 50 COMMON /A9/ ILAB(030,6), ALD1 COMD 55 ICOUN = 1 COMD 60 5 CALL RESID(IW) COMD 65 C COMD 70 C THAT GOT THE CENTERED RESIDUALS; NOW COMPUTE NUMERATOR (E) AND COMD 75 C DENOMINATOR (F) FOR V**2 COMD 80 C COMD 85 E = 0.0 COMD 90 F = 0.0 COMD 95 DO 15 I=2,N COMD 100 I1 = I - 1 COMD 105 PI = P(I,IW) COMD 110 DO 10 J=1,I1 COMD 115 PIJ = P(J,IW)*PI COMD 120 E = E + DATA(J,I)*PIJ COMD 125 F = F + PIJ COMD 130 10 CONTINUE COMD 135 15 CONTINUE COMD 140 VK2 = E*E/HF(F) COMD 145 99999 FORMAT (//51H VARIANCE IN THE RESIDUALS CURRENTLY ACCOUNTED FOR , COMD 150 * 10HBY SUBSET , I3, 2H =, F10.5) COMD 155 HTOT = 0.0 COMD 160 DO 20 I=1,N COMD 165 G(I) = 0. COMD 170 HTOT = HTOT + P(I,IW) COMD 175 20 CONTINUE COMD 180 IF (LP(MST).EQ.1) WRITE (LOUT,99999) IW, VK2 COMD 185 DO 30 I=1,N COMD 190 DO 25 J=1,N COMD 195 IF (J.EQ.I) GO TO 25 COMD 200 K1 = MIN0(I,J) COMD 205 K2 = MAX0(I,J) COMD 210 G(I) = G(I) + P(J,IW)*DATA(K1,K2) COMD 215 25 CONTINUE COMD 220 30 CONTINUE COMD 225 DO 35 K=1,LD2 COMD 230 TRISE(K) = RISE(K) COMD 235 35 CONTINUE COMD 240 IC10 = 0 COMD 245 RMAX = -1.0E20 COMD 250 DO 50 I=2,N COMD 255 I1 = I - 1 COMD 260 PI = P(I,IW) COMD 265 PI2 = 2.0*PI - 1.0 COMD 270 DO 45 J=1,I1 COMD 275 PJ = P(J,IW) COMD 280 PJ2 = 2.0*PJ - 1.0 COMD 285 C COMD 290 C DO NOT CONSIDER DROPPING EITHER MEMBER OF A DYAD SUBSET COMD 295 C OR FORMING THE COMPLETE SUBSET. COMD 300 C COMD 305 HT2 = HTOT - PI2 - PJ2 COMD 310 IF (HT2.LT.1.98 .OR. HT2.GT.AN1) GO TO 45 COMD 315 C COMD 320 C COMPUTE VARIANCE ACCOUNTED FOR BY SUBSET IW, WHEN A COMD 325 C DOUBLETON CHANGE IS MADE, REVERSING ELEMENTS I AND J OF SUBSET COMD 330 C IW. COMD 335 C COMD 340 R3 = 0. COMD 345 ANUM = E - PI2*G(I) - PJ2*G(J) + DATA(J,I)*PI2*PJ2 COMD 350 IF (ABS(ANUM).LE.1.0E-12) GO TO 40 COMD 355 R3 = (ANUM**2)/HF(F-PI2*(HTOT-PI)-PJ2*(HTOT-PJ)+PI2*PJ2) COMD 360 IF (R3.LE.RMAX) GO TO 40 COMD 365 C COMD 370 C UPDATE BEST V**2 COMD 375 C COMD 380 MR2 = I COMD 385 MR3 = J COMD 390 RMAX = R3 COMD 395 C COMD 400 C PRINT OUT (8 AT A TIME) THE PERMISSIBLE DOUBLETON CHANGES AND COMD 405 C RESULTING V**2. COMD 410 C COMD 415 40 CONTINUE COMD 420 IF (LP(MST).NE.1) GO TO 45 COMD 425 IC10 = IC10 + 1 COMD 430 KT2(IC10) = -FLOAT(J)*PJ2 COMD 435 KT1(IC10) = -FLOAT(I)*PI2 COMD 440 TK(IC10) = R3 COMD 445 IF ((IC10.LT.8 .AND. (I.LT.N .OR. J.LT.N1)) .OR. (I.EQ.N .AND. COMD 450 * J.EQ.N1 .AND. IC10.EQ.0)) GO TO 45 COMD 455 WRITE (LOUT,99998) (KT1(K),KT2(K),TK(K),K=1,IC10) COMD 460 99998 FORMAT (1X, 8(1X, 1H(, 2I3, F7.2, 1H))) COMD 465 IC10 = 0 COMD 470 45 CONTINUE COMD 475 50 CONTINUE COMD 480 IF (IC10.NE.0) WRITE (LOUT,99998) (KT1(K),KT2(K),TK(K),K=1,IC10) COMD 485 IF (RMAX.LE.VK2) RETURN COMD 490 P(MR3,IW) = 1.0 - P(MR3,IW) COMD 495 P(MR2,IW) = 1.0 - P(MR2,IW) COMD 500 UVAF = TVAF COMD 505 CALL REGR(0, 0, MST) COMD 510 IF (IERR.EQ.4) GO TO 65 COMD 515 C COMD 520 C FIRST CHECK TO SEE IF THE WEIGHT OF ANY SUBSET WOULD COMD 525 C BECOME NEGATIVE. COMD 530 C COMD 535 DO 55 IK=1,LD1 COMD 540 LK = IK COMD 545 IF (RISE(IK).LE.0.0 .AND. TRISE(IK).GT.0.0) GO TO 60 COMD 550 55 CONTINUE COMD 555 GO TO 75 COMD 560 60 WRITE (LOUT,99997) (ILAB(MR2,I),I=1,6), (ILAB(MR3,I),I=1,6), IW, COMD 565 * TVAF, UVAF, LK, RISE(LK) COMD 570 99997 FORMAT (/20H REVERSING ELEMENTS , 1X, 6A1, 1X, 6A1, 10H OF SUBSET,COMD 575 * 1H , I3, 31H WOULD INCREASE OVERALL VAF TO , F10.6, 9H, REPLACI, COMD 580 * 3HNG , F10.6/48H BUT WAS NOT DONE SINCE THE CURRENTLY POSITIVE W,COMD 585 * 9HEIGHT FOR, 7H SUBSET, I3, 29H WOULD HAVE BECOME NEGATIVE (, COMD 590 * F12.7, 2H).) COMD 595 C COMD 600 C UNSUCCESSFUL CHANGE COMD 605 C COMD 610 65 P(MR2,IW) = 1.0 - P(MR2,IW) COMD 615 P(MR3,IW) = 1.0 - P(MR3,IW) COMD 620 TVAF = UVAF COMD 625 DO 70 K=1,LD2 COMD 630 RISE(K) = TRISE(K) COMD 635 70 CONTINUE COMD 640 RETURN COMD 645 C COMD 650 C FOR A SUCCESSFUL CHANGE COMD 655 C COMD 660 75 ICOUN = 2 COMD 665 WRITE (LOUT,99996) (ILAB(MR2,I),I=1,6), (ILAB(MR3,I),I=1,6), IW, COMD 670 * TVAF, UVAF COMD 675 99996 FORMAT (/20H REVERSING ELEMENTS , 6A1, 1X, 6A1, 11H OF SUBSET , COMD 680 * I3, 33H GAVE AN IMPROVED OVERALL VAF OF , F10.6, 11H REPLACING , COMD 685 * F10.6//) COMD 690 UVAF = TVAF COMD 695 C COMD 700 C NOW GO BACK AND LOOK FOR ANY OTHER SUCCESSFUL CHANGES WITHIN COMD 705 C SUBSET IW. COMD 710 C COMD 715 GO TO 5 COMD 720 END COMD 725 C$ FORTY FORM PUN0 5 SUBROUTINE PUN(MST) PUN0 10 INTEGER FMT4 PUN0 15 DIMENSION LIST(030), FMT4(12) PUN0 20 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), PUN0 25 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, PUN0 30 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC PUN0 35 COMMON /A2/ G(030), KT1(10), KT2(10), TK(10), LP(15), IPUN, IERR, PUN0 40 * LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), TRISE(025) PUN0 45 DATA FMT4(1), FMT4(2), FMT4(3), FMT4(4), FMT4(5), FMT4(6), PUN0 50 * FMT4(7), FMT4(8), FMT4(9), FMT4(10), FMT4(11), FMT4(12) /1H(,1H1,PUN0 55 * 1H4,1HX,1H,,1H5,1H8,1HF,1H1,1H.,1H0,1H)/ PUN0 60 C PUN0 65 C SUBROUTINE TO PUNCH (TO DISK) POST-POLISHED SUBSETS PUN0 70 C PUN0 75 C PUN0 80 C NOTE THAT IF THE NUMBER OF STIMULI EVER EXCEEDS 58, OR IF PUN0 85 C CHANGES ARE MADE IN OUTPUT FORMAT STATEMENT 2, THEN FORMAT PUN0 90 C STATEMENT 1 AND ARRAY FMT4 MUST BE MODIFIED. PUN0 95 C PUN0 100 WRITE (LPUNCH,99999) MST, (TITL(I),I=1,61) PUN0 105 99999 FORMAT (4H IT., I3, 12H CONFIG FOR , 61A1) PUN0 110 WRITE (LPUNCH,99998) FMT4 PUN0 115 99998 FORMAT (1X, 12A1) PUN0 120 DO 10 I=1,LD1 PUN0 125 DO 5 J=1,N PUN0 130 LIST(J) = 0 PUN0 135 IF (P(J,I).GT.0.98 .AND. P(J,I).LT.1.02) LIST(J) = 1 PUN0 140 5 CONTINUE PUN0 145 WRITE (LPUNCH,99997) I, RISE(I), (LIST(J),J=1,N) PUN0 150 10 CONTINUE PUN0 155 99997 FORMAT (1X, I3, 1X, F8.5, 1X, 58I1) PUN0 160 RETURN PUN0 165 END PUN0 170 C$ FORTY FORM,XREF REGR 5 SUBROUTINE REGR(MST, IW, MRT) REGR 10 DOUBLE PRECISION DCON, Y REGR 15 DIMENSION IPVT(025), B(025,025), R(025,025) REGR 20 DIMENSION Y(025) REGR 25 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), REGR 30 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, REGR 35 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC REGR 40 COMMON /A2/ G(030), KT1(10), KT2(10), TK(10), LP(15), IPUN, IERR, REGR 45 * LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), TRISE(025) REGR 50 COMMON /A6/ DCON, DATA(030,030), ALOS, IXV REGR 55 COMMON /A8/ W(025), IPRE, ITPRE REGR 60 C REGR 65 C REGR 70 C REGR 75 C NOTE THAT SUBROUTINE REGR USES THE TRANSFORMED DATA REGR 80 C (LOWERHALF OF ARRAY DATA), UNLIKE MOST OF THE OTHER REGR 85 C SUBROUTINES, WHICH USE THE RESIDUALS (UPPERHALF). REGR 90 C REGR 95 IF (MST.NE.1 .OR. IW.EQ.LD1 .OR. ICON.GT.0) GO TO 10 REGR 100 C REGR 105 C FOR GLOBAL REGRESSION DURING THE FIRST PRE-POLISHING REGR 110 C MAJOR ITERATION WHEN LESS THAN ALL LD1 SUBSETS HAVE BEEN REGR 115 C FORMED, THE ADDITIVE CONSTANT HAS TO BE TACKED ON AS SUBSET REGR 120 C (IW + 1), UNLESS A USER-SUPPLIED INITIAL CONFIGURATION WAS GIVEN. REGR 125 C REGR 130 C REGR 135 C REGR 140 LD0 = IW + 1 REGR 145 DO 5 J=1,N REGR 150 P(J,LD0) = 1.0 REGR 155 5 CONTINUE REGR 160 GO TO 15 REGR 165 10 LD0 = LD2 REGR 170 15 DO 30 I=1,LD0 REGR 175 Y(I) = 0.0D0 REGR 180 DO 25 J=2,N REGR 185 PJI = P(J,I) REGR 190 J1 = J - 1 REGR 195 DO 20 L=1,J1 REGR 200 Y(I) = Y(I) + PJI*P(L,I)*DATA(J,L) REGR 205 20 CONTINUE REGR 210 25 CONTINUE REGR 215 30 CONTINUE REGR 220 DO 40 I=1,N2 REGR 225 DO 35 J=1,N2 REGR 230 B(I,J) = 0. REGR 235 R(I,J) = 0. REGR 240 35 CONTINUE REGR 245 40 CONTINUE REGR 250 DO 60 I=1,LD0 REGR 255 DO 55 J=1,I REGR 260 DO 50 L=2,N REGR 265 PLI = P(L,I) REGR 270 PLJ = P(L,J) REGR 275 L1 = L - 1 REGR 280 DO 45 M1=1,L1 REGR 285 R(I,J) = R(I,J) + PLI*PLJ*P(M1,I)*P(M1,J) REGR 290 45 CONTINUE REGR 295 50 CONTINUE REGR 300 55 CONTINUE REGR 305 60 CONTINUE REGR 310 CALL INVS2(LD0, R, N2, B, N2, IPVT, IERR) REGR 315 IF (IERR.EQ.0) GO TO 65 REGR 320 WRITE (LOUT,99999) REGR 325 99999 FORMAT (///51H ******** TROUBLE ENCOUNTERED WHILE INVERTING MATRI,REGR 330 * 1HX, 24H FOR REGRESSION ********/10X, 23HLINEAR DEPENDENCIES LIK,REGR 335 * 18HELY TO BE PRESENT.//) REGR 340 IF (LD0.NE.LD2) GO TO 105 REGR 345 RETURN REGR 350 65 DO 75 I=1,LD0 REGR 355 W(I) = 0.0 REGR 360 DO 70 J=1,LD0 REGR 365 W(I) = W(I) + Y(J)*B(J,I) REGR 370 70 CONTINUE REGR 375 75 CONTINUE REGR 380 JB = LD1 REGR 385 IF (LD0.NE.LD2) JB = IW REGR 390 IF (MRT.EQ.0) GO TO 85 REGR 395 IF (IXV.EQ.1) GO TO 80 REGR 400 IF ((KP(MRT).EQ.1 .OR. KC.EQ.1) .AND. MST.LE.ITPRE) WRITE REGR 405 * (LOUT,99998) W(LD0), (I,W(I),I=1,JB) REGR 410 99998 FORMAT (/50H WEIGHTS FOR SUBSETS, FITTED BY GLOBAL REGRESSION:, REGR 415 * 10X, 31H( PLUS AN ADDITIVE CONSTANT OF , F12.6, 1H)/64(2H (, I2, REGR 420 * 1H), F9.5, 2H (, I2, 1H), F9.5, 2H (, I2, 1H), F9.5, 2H (, I2, REGR 425 * 1H), F9.5, 2H (, I2, 1H), F9.5, 2H (, I2, 1H), F9.5, 2H (, I2, REGR 430 * 1H), F9.5, 2H (, I2, 1H), F9.5, 2H (, I2, 1H), F9.5/)) REGR 435 GO TO 85 REGR 440 80 IF (LP(MRT).EQ.1) WRITE (LOUT,99997) W(LD0), (I,W(I),I=1,JB) REGR 445 99997 FORMAT (47H -- TENTATIVE -- WEIGHTS FOR SUBSETS, FITTED BY, REGR 450 * 19H GLOBAL REGRESSION:, 10X, 31H( PLUS AN ADDITIVE CONSTANT OF , REGR 455 * F12.6, 1H)/64(2H (, I2, 1H), F9.6, 2H (, I2, 1H), F9.6, 2H (, REGR 460 * I2, 1H), F9.5, 2H (, I2, 1H), F9.5, 2H (, I2, 1H), F9.5, 2H (, REGR 465 * I2, 1H), F9.5, 2H (, I2, 1H), F9.5, 2H (, I2, 1H), F9.5, 2H (, REGR 470 * I2, 1H), F9.5/)) REGR 475 85 DO 90 I=1,JB REGR 480 RISE(I) = W(I) REGR 485 90 CONTINUE REGR 490 RISE(LD2) = W(LD0) REGR 495 IF (LD0.EQ.LD2) GO TO 100 REGR 500 DO 95 J=LD0,LD1 REGR 505 RISE(J) = 0.0 REGR 510 95 CONTINUE REGR 515 100 MRT1 = MST REGR 520 IF (MRT1.EQ.0) MRT1 = MRT REGR 525 IF (LD0.NE.LD2) GO TO 105 REGR 530 CALL VARIAN(MRT1, 2, IW) REGR 535 RETURN REGR 540 105 DO 110 J=1,N REGR 545 P(J,LD0) = 0. REGR 550 110 CONTINUE REGR 555 RETURN REGR 560 END REGR 565 C$ FORTY FORM,XREF INVS 5 SUBROUTINE INVS2(N, A, IA, B, IB, IPVT, IERR) INVS 10 INTEGER N, IA, IB, IPVT(N), IERR INVS 15 LOGICAL FAIL INVS 20 REAL A(IA,IA), B(IA,IA) INVS 25 C INVS 30 C THIS SUBROUTINE INVERTS THE MATRIX A WHERE A IS A SYMMETRIC INVS 35 C MATRIX STORED IN A 2 DIMENSIONAL ARRAY INVS 40 C INVS 45 C INVS 50 C INPUT PARAMETERS INVS 55 C INVS 60 C N THE ORDER OF THE PROBLEM INVS 65 C A A 2 DIMENSIONAL ARRAY CONTAINING INVS 70 C THE COEFFICIENT MATRIX. ONLY THE INVS 75 C LOWER TRIANGULAR PORTION NEED BE INVS 80 C SPECIFIED INVS 85 C IA THE ROW DIMENSION OF THE A MATRIX INVS 90 C IB THE ROW DIMENSION OF THE B MATRIX,INTO INVS 95 C WHICH THE INVERSE OF A WILL BE PLACED INVS 100 C INVS 105 C OUTPUT PARAMETERS INVS 110 C B AN ARRAY CONTAINING THE INVERSE OF A INVS 115 C IERR INTEGER VARIABLE GIVING ERROR CONDITION INVS 120 C 0 NO ERROR INVS 125 C 1 N.LT.1 INVS 130 C 2 IA.LT.N INVS 135 C 3 IB .LT.N INVS 140 C 4 SINGULAR MATRIX INVS 145 C INVS 150 C SCRATCH SPACE INVS 155 C IPVT INTEGER,LENGTH N INVS 160 IERR = 0 INVS 165 IF (N.LT.1) IERR = 1 INVS 170 IF (IA.LT.N) IERR = 2 INVS 175 IF (IB.LT.N) IERR = 3 INVS 180 IF (IERR.NE.0) RETURN INVS 185 DO 10 I=1,N INVS 190 DO 5 J=1,N INVS 195 B(I,J) = 0.0 INVS 200 5 CONTINUE INVS 205 B(I,I) = 1.0 INVS 210 10 CONTINUE INVS 215 C CALL LEQS2(N,A,IA,B,IB,N) INVS 220 CALL DSP2(N, A, IA, IPVT) INVS 225 DO 15 I=1,N INVS 230 CALL SSP2(N, A, IA, B(1,I), IPVT, FAIL) INVS 235 IF (.NOT.FAIL) GO TO 15 INVS 240 IERR = 4 INVS 245 RETURN INVS 250 15 CONTINUE INVS 255 RETURN INVS 260 END INVS 265 C$ FORTY FORM,XREF DSP2 5 SUBROUTINE DSP2(N, A, IA, INTER) DSP2 10 C DSP2 15 C GIVEN A REAL SYMMETRIC MATRIX OF ORDER N,THIS SUBROUTINE DSP2 20 C DETERMINES ITS DECOMPOSITION INTO PMDM(TRANSPOSE)P(TRANSPOSE) DSP2 25 C WHERE P IS A PERMUTATION MATRIX,M IS A UNIT LOWER TRIANGULAR DSP2 30 C MATRIX, AND D IS A BLOCK DIAGONAL MATRIX WITH BLOCKS OF ORDER DSP2 35 C 1 OR 2 WHERE D(I+1,I) IS NONZERO WHENEVER M(I+1,I) IS ZERO. DSP2 40 C ONLY THE LOWER TRIANGULAR PORTION OF A IS USED. THE DECOMPOSITION DSP2 45 C IS PLACED IN THE LOWER TRIANGULAR PORTION. THUS IF ALL THE ELEMENTS DSP2 50 C OF A ARE SPECIFIED,THE STRICT UPPER TRIANGLE IS NOT DESTROYED BUT DSP2 55 C THE DIAGONAL IS DESTROYED. ON OUTPUT THE VECTOR INTER OF DSP2 60 C LENGTH N WILL CONTAIN A RECORD OF THE PERMUTATIONS GENERATED. THE DSP2 65 C INTEGER VARIABLE IA GIVES THE ROW DIMENSION OF THE A MATRIX. DSP2 70 C DSP2 75 REAL A(IA,IA), TEMP, SAVE, DENOM, AII, AIP1I, AIP1 DSP2 80 REAL ALPHA, AAII, SIGMA, ALFLAM, LAMBDA DSP2 85 INTEGER INTER(N) DSP2 90 INTER(N) = N DSP2 95 ALPHA = (1.0+SQRT(17.0))/8. DSP2 100 I = 1 DSP2 105 5 IF (I.GE.N) RETURN DSP2 110 AAII = ABS(A(I,I)) DSP2 115 INTER(I) = I DSP2 120 C DSP2 125 C FIND THE LARGEST OFF DIAGONAL ELEMENT IN THE I-TH COLUMN DSP2 130 C DSP2 135 J = I + 1 DSP2 140 IP1 = I + 1 DSP2 145 LAMBDA = ABS(A(IP1,I)) DSP2 150 IP2 = I + 2 DSP2 155 IF (IP2.GT.N) GO TO 15 DSP2 160 DO 10 K=IP2,N DSP2 165 IF (ABS(A(K,I)).LE.LAMBDA) GO TO 10 DSP2 170 LAMBDA = ABS(A(K,I)) DSP2 175 J = K DSP2 180 10 CONTINUE DSP2 185 15 ALFLAM = ALPHA*LAMBDA DSP2 190 IF (AAII.GE.ALFLAM) GO TO 65 DSP2 195 C DSP2 200 C FIND THE LARGEST OFF DIAGONAL ELEMENT IN THE J-TH COLUMN DSP2 205 C DSP2 210 SIGMA = LAMBDA DSP2 215 JM1 = J - 1 DSP2 220 IF (IP1.GT.JM1) GO TO 25 DSP2 225 DO 20 K=IP1,JM1 DSP2 230 IF (ABS(A(J,K)).GT.SIGMA) SIGMA = ABS(A(J,K)) DSP2 235 20 CONTINUE DSP2 240 25 JP1 = J + 1 DSP2 245 IF (JP1.GT.N) GO TO 35 DSP2 250 DO 30 K=JP1,N DSP2 255 IF (ABS(A(K,J)).GT.SIGMA) SIGMA = ABS(A(K,J)) DSP2 260 30 CONTINUE DSP2 265 35 IF (AAII*SIGMA.GE.ALFLAM*LAMBDA) GO TO 65 DSP2 270 IF (ABS(A(J,J)).GE.ALPHA*SIGMA) GO TO 60 DSP2 275 C DSP2 280 C PERFORM A 2 BY 2 PIVOT STEP DSP2 285 C DSP2 290 INTER(I) = J DSP2 295 IF (IP2.GT.N) GO TO 55 DSP2 300 IF (J.EQ.IP1) GO TO 40 DSP2 305 CALL ISP2(A, IA, N, J, IP1) DSP2 310 TEMP = A(J,I) DSP2 315 A(J,I) = A(IP1,I) DSP2 320 A(IP1,I) = TEMP DSP2 325 40 AIP1I = A(IP1,I) DSP2 330 AII = A(I,I)/AIP1I DSP2 335 AIP1 = A(IP1,IP1) DSP2 340 DENOM = AII*AIP1 - AIP1I DSP2 345 DO 50 J=IP2,N DSP2 350 TEMP = (A(J,I)-AII*A(J,IP1))/DENOM DSP2 355 SAVE = -(AIP1*TEMP+A(J,IP1))/AIP1I DSP2 360 DO 45 K=J,N DSP2 365 A(K,J) = A(K,J) + A(K,I)*SAVE + A(K,IP1)*TEMP DSP2 370 45 CONTINUE DSP2 375 A(J,I) = SAVE DSP2 380 A(J,IP1) = TEMP DSP2 385 50 CONTINUE DSP2 390 55 INTER(IP1) = -1 DSP2 395 I = IP2 DSP2 400 GO TO 5 DSP2 405 C DSP2 410 C INTERCHANGE THE I-TH AND J-TH ROWS AND COLUMNS DSP2 415 C DSP2 420 60 INTER(I) = J DSP2 425 CALL ISP2(A, IA, N, J, I) DSP2 430 C DSP2 435 C PERFORM A 1 X 1 PIVOT DSP2 440 C DSP2 445 65 IF (A(I,I).EQ.0.0) GO TO 80 DSP2 450 AII = A(I,I) DSP2 455 DO 75 J=IP1,N DSP2 460 SAVE = -A(J,I)/AII DSP2 465 IF (SAVE.EQ.0.0) GO TO 75 DSP2 470 DO 70 K=J,N DSP2 475 A(K,J) = A(K,J) + A(K,I)*SAVE DSP2 480 70 CONTINUE DSP2 485 A(J,I) = SAVE DSP2 490 75 CONTINUE DSP2 495 80 I = IP1 DSP2 500 GO TO 5 DSP2 505 END DSP2 510 C$ FORTY FORM,XREF ISP2 5 SUBROUTINE ISP2(A, IA, N, J, I) ISP2 10 REAL A(IA,IA), TEMP ISP2 15 C INTERCHANGE THE ELEMENTS BELOW BOTH DIAGONALS ISP2 20 JP1 = J + 1 ISP2 25 IF (JP1.GT.N) GO TO 10 ISP2 30 DO 5 K=JP1,N ISP2 35 TEMP = A(K,J) ISP2 40 A(K,J) = A(K,I) ISP2 45 A(K,I) = TEMP ISP2 50 5 CONTINUE ISP2 55 10 IF (I+1.GT.J-1) GO TO 20 ISP2 60 IP1 = I + 1 ISP2 65 JM1 = J - 1 ISP2 70 DO 15 K=IP1,JM1 ISP2 75 TEMP = A(K,I) ISP2 80 A(K,I) = A(J,K) ISP2 85 A(J,K) = TEMP ISP2 90 15 CONTINUE ISP2 95 C INTERCHANGE THE DIAGONAL ELEMENTS ISP2 100 20 TEMP = A(I,I) ISP2 105 A(I,I) = A(J,J) ISP2 110 A(J,J) = TEMP ISP2 115 RETURN ISP2 120 END ISP2 125 C$ FORTY FORM,XREF SSP2 5 SUBROUTINE SSP2(N, A, IA, B, INTER, FAIL) SSP2 10 C SSP2 15 C THIS SUBROUTINE SOLVES THE LINEAR SYSTEM AX=B WHERE A IS A SSP2 20 C REAL SYMMETRIC MATRIX OF ORDER N AND ROW DIMENSION IA SSP2 25 C WHOSE DECOMPOSITION HAS BEEN COMPUTED BY THE SUBROUTINE SSP2 30 C DSP2 AND LEFT IN THE LOWER TRIANGULAR PORTION OF A AND B SSP2 35 C IS AN N VECTOR. THE SOLUTION X IS RETURNED IN THE VECTOR B. SSP2 40 C THE VECTOR INTER OF LENGTH N IS GENERATED BY DSP2 AND CONTAINS SSP2 45 C INFORMATION ABOUT THE PERMUTATIONS PERFORMED ON THE A MATRIX SSP2 50 C IN DSP2. SSP2 55 C SSP2 60 C IF THE MATRIX IS SINGULAR,THE LOGICAL VARIABLE FAIL WILL SSP2 65 C BE TRUE, AND THE SOLUTION WILL NOT BE FOUND. SSP2 70 REAL A(IA,IA), B(IA), SAVE, TEMP, DENOM SSP2 75 INTEGER INTER(N) SSP2 80 LOGICAL FAIL SSP2 85 FAIL = .FALSE. SSP2 90 I = 1 SSP2 95 C SSP2 100 C SOLVE M D Y=B FOR Y AND STORE Y IN B SSP2 105 C SSP2 110 5 IF (I.GE.N) GO TO 50 SSP2 115 IP1 = I + 1 SSP2 120 ICH = INTER(I) SSP2 125 SAVE = B(ICH) SSP2 130 IF (INTER(IP1).LT.0) GO TO 35 SSP2 135 B(ICH) = B(I) SSP2 140 C SSP2 145 C SSP2 150 C IN THE STATEMENTS THAT FOLLOW, THE NORM OF A() MULTIPLIED SSP2 155 C BY 1.0E-20 IS USED FOR A TEST OF A ZERO PIVOTAL ELEMENT. SSP2 160 C THIS CRITERION IS FAIRLY CONSERVATIVE AND SOMEWHAT MACHINE- SSP2 165 C DEPENDENT. SSP2 170 C SSP2 175 EPS = 0.0 SSP2 180 DO 15 L=1,N SSP2 185 DO 10 M=1,N SSP2 190 EPS = EPS + A(L,M)**2 SSP2 195 10 CONTINUE SSP2 200 15 CONTINUE SSP2 205 EPS = SQRT(EPS)*1.0E-20 SSP2 210 IF (ABS(A(I,I)).GE.EPS) GO TO 20 SSP2 215 FAIL = .TRUE. SSP2 220 RETURN SSP2 225 20 B(I) = SAVE/A(I,I) SSP2 230 IF (SAVE.EQ.0.0) GO TO 30 SSP2 235 DO 25 J=IP1,N SSP2 240 B(J) = B(J) + A(J,I)*SAVE SSP2 245 25 CONTINUE SSP2 250 30 I = IP1 SSP2 255 GO TO 5 SSP2 260 35 TEMP = B(I) SSP2 265 B(ICH) = B(IP1) SSP2 270 DENOM = A(I,I)*A(IP1,IP1)/A(IP1,I) - A(IP1,I) SSP2 275 B(IP1) = (SAVE*A(I,I)/A(IP1,I)-TEMP)/DENOM SSP2 280 B(I) = (SAVE-A(IP1,IP1)*B(IP1))/A(IP1,I) SSP2 285 IP2 = I + 2 SSP2 290 IF (IP2.GT.N) GO TO 45 SSP2 295 DO 40 J=IP2,N SSP2 300 B(J) = B(J) + A(J,I)*TEMP + A(J,IP1)*SAVE SSP2 305 40 CONTINUE SSP2 310 45 I = IP2 SSP2 315 GO TO 5 SSP2 320 C SSP2 325 C SSP2 330 50 IF (I.NE.N) GO TO 60 SSP2 335 IF (ABS(A(I,I)).GE.EPS) GO TO 55 SSP2 340 FAIL = .TRUE. SSP2 345 RETURN SSP2 350 55 B(I) = B(I)/A(I,I) SSP2 355 C NOW SOLVE M(TRANSPOSE) X = Y FOR X, WHERE Y IS STORED SSP2 360 C IN THE VECTOR B, AND STORE X IN B SSP2 365 C SSP2 370 60 I = N - 1 SSP2 375 IF (INTER(N).LT.0) I = N - 2 SSP2 380 65 IF (I.LE.0) RETURN SSP2 385 I1 = I SSP2 390 IF (INTER(I).LT.0) I1 = I - 1 SSP2 395 IP1 = I + 1 SSP2 400 DO 75 K=I1,I SSP2 405 SAVE = B(K) SSP2 410 DO 70 J=IP1,N SSP2 415 SAVE = SAVE + A(J,K)*B(J) SSP2 420 70 CONTINUE SSP2 425 B(K) = SAVE SSP2 430 75 CONTINUE SSP2 435 ICH = INTER(I1) SSP2 440 B(I) = B(ICH) SSP2 445 B(ICH) = SAVE SSP2 450 I = I1 - 1 SSP2 455 GO TO 65 SSP2 460 END SSP2 465 C$ FORTY FORM ASUM 5 FUNCTION ASUM(IW) ASUM 10 DOUBLE PRECISION DCON ASUM 15 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), ASUM 20 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, ASUM 25 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC ASUM 30 COMMON /A6/ DCON, DATA(030,030), ALOS, IXV ASUM 35 C ASUM 40 C COMPUTE THE A PART OF THE LOSS FUNCTION, WHICH IS A CONSTANT ASUM 45 C FOR ALL I. ASUM 50 C ASUM 55 ASUM = 0.0 ASUM 60 WIW = RISE(IW) ASUM 65 DO 10 I=2,N ASUM 70 PI = P(I,IW) ASUM 75 I1 = I - 1 ASUM 80 DO 5 J=1,I1 ASUM 85 TEST1 = DATA(J,I) - WIW*PI*P(J,IW) ASUM 90 IF (ABS(TEST1).LE.1.0E-12) GO TO 5 ASUM 95 ASUM = ASUM + TEST1**2 ASUM 100 5 CONTINUE ASUM 105 10 CONTINUE ASUM 110 ASUM = ASUM*DCON ASUM 115 RETURN ASUM 120 END ASUM 125 C$ FORTY FORM,XREF BSUM 5 FUNCTION BSUM(IW) BSUM 10 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), BSUM 15 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, BSUM 20 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC BSUM 25 COMMON /A4/ DAPI(030), BCONST, TM, U, V BSUM 30 C BSUM 35 C COMPUTE THE B PART OF THE LOSS FUNCTION, WHICH IS A CONSTANT BSUM 40 C FOR ALL I. BSUM 45 C BSUM 50 TM = 0.0 BSUM 55 DO 10 I=2,N BSUM 60 PI = P(I,IW) BSUM 65 I1 = I - 1 BSUM 70 DO 5 J=1,I1 BSUM 75 TM = TM + PI*P(J,IW) BSUM 80 5 CONTINUE BSUM 85 10 CONTINUE BSUM 90 TM = TM*ABNV BSUM 95 C BSUM 100 C TM IS THE MEAN (P(I,IW)*P(J,IW)) I .GT. J BSUM 105 C BSUM 110 U = 0.0 BSUM 115 V = 0.0 BSUM 120 DO 25 I=2,N BSUM 125 PI = P(I,IW) BSUM 130 I1 = I - 1 BSUM 135 DO 20 J=1,I1 BSUM 140 PIPJ = PI*P(J,IW) BSUM 145 TEST1 = (PIPJ-1.0)*PIPJ BSUM 150 IF (ABS(TEST1).LE.1.0E-12) GO TO 15 BSUM 155 U = U + TEST1**2 BSUM 160 15 CONTINUE BSUM 165 TEST2 = PIPJ - TM BSUM 170 IF (ABS(TEST2).LE.1.0E-12) GO TO 20 BSUM 175 V = V + TEST2**2 BSUM 180 20 CONTINUE BSUM 185 25 CONTINUE BSUM 190 C NOW ADD THE DIAGONAL ELEMENTS TO U BSUM 195 DO 30 I=1,N BSUM 200 PI2 = P(I,IW)**2 BSUM 205 TEST3 = (PI2-1.0)*PI2 BSUM 210 IF (ABS(TEST3).LE.1.0E-12) GO TO 30 BSUM 215 U = U + 0.5*TEST3**2 BSUM 220 30 CONTINUE BSUM 225 BSUM = U/V BSUM 230 RETURN BSUM 235 END BSUM 240 C$ FORTY FORM,XREF VARI 5 SUBROUTINE VARIAN(MST, IND, IW) VARI 10 DOUBLE PRECISION DCON VARI 15 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), VARI 20 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, VARI 25 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC VARI 30 COMMON /A2/ G(030), KT1(10), KT2(10), TK(10), LP(15), IPUN, IERR, VARI 35 * LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), TRISE(025) VARI 40 COMMON /A6/ DCON, DATA(030,030), ALOS, IXV VARI 45 COMMON /A8/ W(025), IPRE, ITPRE VARI 50 C VARI 55 C VARI 60 C VARI 65 C VARI 70 IF (IND.NE.1) GO TO 5 VARI 75 VAF = 1.0 - 4.0*ALOS*ABNV VARI 80 WRITE (LOUT,99999) MST, VAF, IW VARI 85 99999 FORMAT (/17H DURING MAJOR IT., I4, 25H VARIANCE ACCOUNTED FOR =, VARI 90 * F12.6, 29H IN THE RESIDUALS BY SUBSET (, I2, 16H), WITHOUT POLIS,VARI 95 * 4HHING) VARI 100 RETURN VARI 105 C VARI 110 C IF REGRESSION HAS JUST FAILED, THE FOLLOWING DEFINITION VARI 115 C IS INVALID. VARI 120 C VARI 125 5 CONS = RISE(LD2) VARI 130 SSQ = 0.0 VARI 135 DO 20 I=2,N VARI 140 I1 = I - 1 VARI 145 DO 15 J=1,I1 VARI 150 DATAJI = 0.0 VARI 155 DO 10 K=1,LD1 VARI 160 DATAJI = DATAJI + P(I,K)*P(J,K)*RISE(K) VARI 165 10 CONTINUE VARI 170 SSQ = SSQ + (DATA(I,J)-CONS-DATAJI)**2 VARI 175 15 CONTINUE VARI 180 20 CONTINUE VARI 185 VAF = (VARSIM-SSQ)/VARSIM VARI 190 TVAF = VAF VARI 195 IF (IXV.EQ.1) GO TO 30 VARI 200 IF (MST.LE.0) GO TO 25 VARI 205 IF (KP(MST).NE.1 .AND. KC.NE.1) RETURN VARI 210 25 WRITE (LOUT,99998) MST, VAF VARI 215 99998 FORMAT (/24H DURING MAJOR ITERATION , I4, 19H, OVERALL VARIANCE , VARI 220 * 15HACCOUNTED FOR =, F12.6, 31H AFTER USING REGRESSION WEIGHTS) VARI 225 RETURN VARI 230 30 WRITE (LOUT,99997) MST, VAF VARI 235 99997 FORMAT (17H DURING ITERATION, I4, 19H, OVERALL VARIANCE , VARI 240 * 31H--TENTATIVELY-- ACCOUNTED FOR =, F12.6, 18H AFTER USING REGRE,VARI 245 * 13HSSION WEIGHTS) VARI 250 RETURN VARI 255 END VARI 260 C$ FORTY FORM POLI 5 SUBROUTINE POLISH(MST, IW) POLI 10 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), POLI 15 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, POLI 20 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC POLI 25 COMMON /A4/ DAPI(030), BCONST, TM, U, V POLI 30 COMMON /A7/ A(030), IZ(030) POLI 35 C POLI 40 C SUBROUTINE TO FORCE P(,IW) TO BE ZEROES AND ONES. POLI 45 C POLI 50 C POLI 55 M1 = 0 POLI 60 DO 10 KLM=1,N POLI 65 A(KLM) = P(KLM,IW) POLI 70 IF (A(KLM).GE.0.25 .AND. A(KLM).LE.0.75) MP(MST) = MP(MST) + 1 POLI 75 IF (A(KLM).GE.0.707) GO TO 5 POLI 80 P(KLM,IW) = 0.0 POLI 85 GO TO 10 POLI 90 5 M1 = M1 + 1 POLI 95 P(KLM,IW) = 1.0 POLI 100 10 CONTINUE POLI 105 IF (M1.LE.1 .OR. M1.EQ.N) GO TO 15 POLI 110 RETURN POLI 115 C POLI 120 C IF A VECTOR OF ALL ONES OR ZEROES HAS BEEN GENERATED, PRINT POLI 125 C A DIAGNOSTIC. POLI 130 C POLI 135 15 IF (KP(MST).EQ.1) WRITE (LOUT,99999) IW, M1 POLI 140 99999 FORMAT (/49H ******** WARNING FROM SUBROUTINE POLISH ********/ POLI 145 * 7H SUBSET, I3, 24H IS TRIVIAL SINCE IT HAS, I3, 7H ONE(S), POLI 150 * 39H ON THE PRESENT POLISHING; FIX-UP TAKEN//) POLI 155 DO 20 KLM=1,N POLI 160 IZ(KLM) = KLM POLI 165 20 CONTINUE POLI 170 CALL SORT(N) POLI 175 IF (M1.EQ.N) GO TO 25 POLI 180 C POLI 185 C FOR A NULL SUBSET, CHANGE THE TWO LARGEST ENTRIES TO UNITIES. POLI 190 C POLI 195 IZ1 = IZ(1) POLI 200 P(IZ1,IW) = 1.0 POLI 205 IZ2 = IZ(2) POLI 210 P(IZ2,IW) = 1.0 POLI 215 RETURN POLI 220 C POLI 225 C OR IF THE SUBSET CONTAINS THE WHOLE SET OF STIMULI, CHANGE THE POLI 230 C SMALLEST ENTRY IN THE VECTOR TO ZERO. POLI 235 C POLI 240 25 IZN = IZ(N) POLI 245 P(IZN,IW) = 0.0 POLI 250 RETURN POLI 255 END POLI 260 C$ FORTY FORM SORT 5 SUBROUTINE SORT(N) SORT 10 COMMON /A7/ A(030), IZ(030) SORT 15 C SORT 20 C PUTS ELEMENTS OF VECTOR IN DESCENDING ORDER SORT 25 C SORT 30 C SORT 35 M = 1 SORT 40 5 M = M + M SORT 45 IF (M.LE.N) GO TO 5 SORT 50 M = M - 1 SORT 55 10 M = M/2 SORT 60 IF (M.EQ.0) GO TO 35 SORT 65 KK = N - M SORT 70 J = 1 SORT 75 15 I = J SORT 80 20 IM = I + M SORT 85 IF (A(I).LT.A(IM)) GO TO 30 SORT 90 25 J = J + 1 SORT 95 IF (J.GT.KK) GO TO 10 SORT 100 GO TO 15 SORT 105 30 TEMP = A(I) SORT 110 A(I) = A(IM) SORT 115 A(IM) = TEMP SORT 120 ITEMP = IZ(I) SORT 125 IZ(I) = IZ(IM) SORT 130 IZ(IM) = ITEMP SORT 135 I = I - M SORT 140 IF (I.LT.1) GO TO 25 SORT 145 GO TO 20 SORT 150 35 RETURN SORT 155 END SORT 160 C$ FORTY FORM,XREF INIC 5 SUBROUTINE INICON(MST, IW, ICON) INIC 10 DOUBLE PRECISION DCON INIC 15 DIMENSION U(030), FMT2(72) INIC 20 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), INIC 25 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, INIC 30 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC INIC 35 COMMON /A6/ DCON, DATA(030,030), ALOS, IXV INIC 40 COMMON /A8/ W(025), IPRE, ITPRE INIC 45 C INIC 50 C INIC 55 IF (ICON) 80, 5, 55 INIC 60 C INIC 65 C FOR A RATIONAL INITIAL CONFIGURATION INIC 70 C INIC 75 5 NP = 0 INIC 80 NN = 0 INIC 85 SUM = 0.0 INIC 90 SUMP = 0.0 INIC 95 SUMN = 0.0 INIC 100 DO 15 I=1,N INIC 105 U(I) = 0.0 INIC 110 DO 10 J=1,N INIC 115 IF (I.EQ.J) GO TO 10 INIC 120 K1 = MIN0(I,J) INIC 125 K2 = MAX0(I,J) INIC 130 U(I) = U(I) + DATA(K1,K2) INIC 135 10 CONTINUE INIC 140 15 CONTINUE INIC 145 DO 20 J=1,N INIC 150 SUM = SUM + U(J) INIC 155 20 CONTINUE INIC 160 SUM = SUM*BIN INIC 165 DO 35 J=1,N INIC 170 U(J) = U(J) - SUM INIC 175 IF (U(J)) 25, 35, 30 INIC 180 25 NN = NN + 1 INIC 185 SUMN = SUMN + U(J) INIC 190 GO TO 35 INIC 195 30 NP = NP + 1 INIC 200 SUMP = SUMP + U(J) INIC 205 35 CONTINUE INIC 210 IF (NN*NP.EQ.0) GO TO 45 INIC 215 SUMP = SUMP/FLOAT(NP) INIC 220 SUMN = SUMN/FLOAT(NN) INIC 225 BB = 1.0/(SUMP-SUMN) INIC 230 A = -SUMN/(SUMP-SUMN) INIC 235 C WRITE (LOUT,200) NN,NP,SUMP,SUMN,A,BB,(U(K),K=1, N) INIC 240 C 200 FORMAT(" NN =",I3," NP =",I3," SUMP =",F12.6," SUMN =",F12.6, INIC 245 C 2" A =",F12.6," BB =",F12.6/(1X,10F12.6/)) INIC 250 DO 40 J=1,N INIC 255 P(J,IW) = A + BB*U(J) INIC 260 40 CONTINUE INIC 265 C INIC 270 IF (MST.EQ.0) RETURN INIC 275 IF (MST.EQ.1 .OR. KP(MST).EQ.1) WRITE (LOUT,99999) IW, INIC 280 * (P(KLM,IW),KLM=1,N) INIC 285 99999 FORMAT (/53H RATIONAL INITIAL CONFIGURATION P VECTOR FOR SUBSET ,INIC 290 * 1H(, I2, 2H):/64(1X, 16F8.4/)) INIC 295 RETURN INIC 300 C INIC 305 C INIC 310 C ISSUE DIAGNOSTIC IF NN (NUMBER OF NEG. U(J)) OR NP (POS.) INIC 315 C IS ZERO. THAT IS, DATA(J,I) WAS A DOUBLE-CENTERED MATRIX. INIC 320 C INIC 325 45 WRITE (LOUT,99998) NN, NP INIC 330 99998 FORMAT (///50H ********DIAGNOSTIC FROM SUBROUTINE INICON********//INIC 335 * 19H NO. OF NEG U(J) =, I3, 15H AND POS U(J) =, I3//20X, INIC 340 * 30HLISTING OF MATRIX OF RESIDUALS//) INIC 345 DO 50 I=2,N INIC 350 I1 = I - 1 INIC 355 WRITE (LOUT,99997) (I,J,DATA(J,I),J=1,I1) INIC 360 50 CONTINUE INIC 365 99997 FORMAT (2H (, 2I3, 1H), F9.4, 2H (, 2I3, 1H), F9.4, 2H (, 2I3, INIC 370 * 1H), F9.4, 2H (, 2I3, 1H), F9.4, 2H (, 2I3, 1H), F9.4, 2H (, INIC 375 * 2I3, 1H), F9.4/) INIC 380 STOP INIC 385 C INIC 390 C INIC 395 C INIC 400 C IF (ICON .GT.0) READ IN A USER-SUPPLIED INITIAL CONFIGURATION FOR INIC 405 C P(,IW) VIA CHANNEL ICON. INIC 410 C INIC 415 55 WRITE (LOUT,99996) LD1, ICON INIC 420 99996 FORMAT (/////49H USER HAS SPECIFIED THAT INITIAL CONFIGURATION OF,INIC 425 * I4, 44H CLUSTERS IS TO BE READ FROM LOGICAL DEVICE , I3/6H WITH ,INIC 430 * 32HTHE TITLE AND FORMAT AS FOLLOWS:///) INIC 435 C INIC 440 C READ THE TITLE CARD AND PRINT IT; THEN DO THE SAME FOR THE FORMAT.INIC 445 C INIC 450 READ (ICON,99995) FMT2 INIC 455 WRITE (LOUT,99994) FMT2 INIC 460 99995 FORMAT (72A1) INIC 465 99994 FORMAT (//1X, 46HTITLE FOR USER-SUPPLIED INITIAL CONFIGUATION: , INIC 470 * 72A1//) INIC 475 99993 FORMAT (//1X, 33HUSER-SPECIFIED FORMAT OF INITIAL , 10HCONFIGURAT,INIC 480 * 8HION IS: , 72A1//) INIC 485 READ (ICON,99995) FMT2 INIC 490 WRITE (LOUT,99993) FMT2 INIC 495 WRITE (LOUT,99992) (I, I = 1, N) INIC 500 99992 FORMAT (///64(1X, I6, 15I8/)) INIC 505 DO 60 IX=1,LD1 INIC 510 READ (ICON,FMT2,END=67) (P(KLM,IX), KLM = 1, N) INIC 515 WRITE (LOUT,99991) IX, (P(KLM,IX),KLM=1,N) INIC 520 60 CONTINUE INIC 525 99991 FORMAT (//34H INITIAL CONFIGURATION P VECTOR (, I2, 2H):/64(1X, INIC 530 * 16F8.4/)) INIC 535 C INIC 540 C INIC 545 C THIS CHECK SPOTS SINGLETONS AND LUMPER, AND ALSO LOOKS INIC 550 C TO SEE IF USER-SUPPLIED P MATRIX WERE 0,1. INIC 555 C INIC 560 IDANG = 0 INIC 565 DO 75 J=1,LD1 INIC 570 IK = 0 INIC 575 DO 70 I=1,N INIC 580 PIJ = P(I,J) INIC 585 IF (PIJ.LT.1.01 .AND. PIJ.GT.0.99) GO TO 65 INIC 590 IF (-1.0E-04.LE.PIJ .AND. PIJ.LT.1.0E-04) GO TO 70 INIC 595 WRITE (LOUT,99990) I, J, P(I,J) INIC 600 99990 FORMAT (/35H WARNING: IN USER-SUPPLIED INITIAL , 13HCONFIGURATION,INIC 605 * 4H, P(, I3, 1H,, I3, 4H) = , F10.5, 25H INSTEAD OF BEING A BINAR,INIC 610 * 14HY COEFFICIENT.) INIC 615 65 IK = IK + 1 INIC 620 70 CONTINUE INIC 625 IF (IK.GT.1 .AND. IK.LT.N) GO TO 75 INIC 630 IDANG = 1 INIC 635 WRITE (LOUT,99989) J, IK INIC 640 99989 FORMAT (///17H ******** SUBSET , I3, 24H HAS TOO FEW (0 OR 1) OR, INIC 645 * 30H TOO MANY (N) ELEMENTS, NAMELY, I3, 15H AND WILL CAUSE/10X, INIC 650 * 47HTHE MULTIPLE LINEAR REGRESSION TO FAIL ********//) INIC 655 75 CONTINUE INIC 660 IF (IDANG.NE.0) STOP INIC 665 RETURN INIC 670 67 WRITE (LOUT,99988) ICON INIC 675 99988 FORMAT (///49H ******** READ AN END OF FILE WHILE ATTEMPTING TO, INIC 680 * 26H READ INIT. CONF. FOR T(J)/10X, 26HIN SUBROUTINE INICON, USIN,INIC 685 * 1HG, 13H CHANNEL NO. , I3, 9H ********) INIC 690 RETURN INIC 695 C INIC 700 C FOR A RANDOM INITIAL CONFIGURATION INIC 705 C INIC 710 80 DO 85 KLM=1,N INIC 715 P(KLM,IW) = RUNIFV(0) INIC 720 85 CONTINUE INIC 725 WRITE (LOUT,99987) IW, (P(KLM,IW),KLM=1,N) INIC 730 99987 FORMAT (//33H RANDOM INITIAL CONFIGURATION P (, I2, 2H):/64(1X, INIC 735 * 16F8.4/)) INIC 740 RETURN INIC 745 END INIC 750 C$ FORTY FORM RUNI 5 FUNCTION RUNIFV(L) RUNI 10 DATA MODULO, FLMOD, K /8192,8192.0,1/ RUNI 15 C RUNI 20 C RANDOM NUMBERS MODULO 2**13 RUNI 25 C BASED ON A MODIFICATION OF A ROUTINE SUGGESTED BY J. B. KRUSKAL RUNI 30 C (CACM, 1969, 12, 93-94), THIS SUBROUTINE IS SIMILAR TO THE ONE RUNI 35 C USED IN THE MULTIDIMENSIONAL SCALING PROGRAM KYST. RUNI 40 C RUNI 45 DO 5 I=1,3 RUNI 50 K = MOD(5*K,MODULO) RUNI 55 5 CONTINUE RUNI 60 RUNIFV = FLOAT(K)/FLMOD RUNI 65 RETURN RUNI 70 END RUNI 75 C$ FORTY FORM,XREF RESI 5 SUBROUTINE RESID(IW) RESI 10 DOUBLE PRECISION DCON, SUMA RESI 15 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), RESI 20 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, RESI 25 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC RESI 30 COMMON /A6/ DCON, DATA(030,030), ALOS, IXV RESI 35 C RESI 40 C RESI 45 C IN THE ALTERNATING LEAST SQUARES PROCEDURE, COMPUTE RESI 50 C THE RESIDUALS FOR SUBSET IW, AND CENTER THEM RESI 55 C RESI 60 C RESI 65 DO 10 I=2,N RESI 70 I1 = I - 1 RESI 75 DO 5 J=1,I1 RESI 80 DATA(J,I) = 0.0 RESI 85 5 CONTINUE RESI 90 10 CONTINUE RESI 95 DO 25 K=1,LD1 RESI 100 IF (K.EQ.IW) GO TO 25 RESI 105 WIW = RISE(K) RESI 110 DO 20 I=2,N RESI 115 I1 = I - 1 RESI 120 PIK = P(I,K) RESI 125 DO 15 J=1,I1 RESI 130 DATA(J,I) = DATA(J,I) + PIK*P(J,K)*WIW RESI 135 15 CONTINUE RESI 140 20 CONTINUE RESI 145 25 CONTINUE RESI 150 SUM = 0.0 RESI 155 DO 35 I=2,N RESI 160 I1 = I - 1 RESI 165 DO 30 J=1,I1 RESI 170 DATA(J,I) = DATA(I,J) - DATA(J,I) RESI 175 SUM = SUM + DATA(J,I) RESI 180 30 CONTINUE RESI 185 35 CONTINUE RESI 190 SUM = SUM*ABNV RESI 195 SUMA = 0.0 RESI 200 DO 45 I=2,N RESI 205 I1 = I - 1 RESI 210 DO 40 J=1,I1 RESI 215 DATA(J,I) = DATA(J,I) - SUM RESI 220 C RESI 225 C SINCE THE UPPER TRIANGULAR HALF OF DATA CONSISTS OF CENTERED RESI 230 C RESIDUALS, RESI 235 C THE SUM OF SQUARES EQUALS THE NUMERATOR OF THE VARIANCE. RESI 240 C RESI 245 SUMA = SUMA + DATA(J,I)**2 RESI 250 40 CONTINUE RESI 255 45 CONTINUE RESI 260 C RESI 265 C COMPUTE THE CONSTANT THAT NORMALIZES THE SUM OF SQUARED ERROR. RESI 270 C RESI 275 DCON = 0.25*AB/SUMA RESI 280 RETURN RESI 285 END RESI 290 C$ FORTY FORM ADJU 5 SUBROUTINE ADJUST(MST, IW) ADJU 10 COMMON /A5/ AK, ALFRAY(025), BETRAY(025), IADJ(050), JADJ ADJU 15 IADJ(MST) = IADJ(MST) + 1 ADJU 20 JADJ = JADJ + 1 ADJU 25 ALFRAY(IW) = AMAX1(0.000001,ALFRAY(IW)/(ALFRAY(IW)+AK*BETRAY(IW)))ADJU 30 BETRAY(IW) = 1.0 - ALFRAY(IW) ADJU 35 RETURN ADJU 40 END ADJU 45 C$ FORTY FORM,XREF CDAP 5 SUBROUTINE CDAPI(IW) CDAP 10 DOUBLE PRECISION DCON CDAP 15 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), CDAP 20 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, CDAP 25 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC CDAP 30 COMMON /A4/ DAPI(030), BCONST, TM, U, V CDAP 35 COMMON /A6/ DCON, DATA(030,030), ALOS, IXV CDAP 40 C CDAP 45 C COMPUTE THE PARTIAL OF A WITH RESPECT TO P(I,IW) CDAP 50 C CDAP 55 DO 10 I=1,N CDAP 60 WIW = RISE(IW) CDAP 65 PI = P(I,IW) CDAP 70 AS = 0. CDAP 75 DO 5 J=1,N CDAP 80 IF (J.EQ.I) GO TO 5 CDAP 85 PJ = P(J,IW) CDAP 90 K1 = MIN0(I,J) CDAP 95 K2 = MAX0(I,J) CDAP 100 AS = AS + PJ*(DATA(K1,K2)-WIW*PJ*PI) CDAP 105 5 CONTINUE CDAP 110 DAPI(I) = -2.0*WIW*AS*DCON CDAP 115 10 CONTINUE CDAP 120 RETURN CDAP 125 END CDAP 130 C$ FORTY FORM,XREF DB00 5 FUNCTION DB(KLM, IW) DB00 10 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), DB00 15 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, DB00 20 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC DB00 25 COMMON /A4/ DAPI(030), BCONST, TM, U, V DB00 30 C DB00 35 C FUNCTION TO ASSEMBLE ALL THE EXPRESSIONS FOR DB00 40 C EVALUATING DB/DP(I) DB00 45 C DB00 50 B1 = 0.0 DB00 55 B2 = 0.0 DB00 60 DO 5 I=1,N DB00 65 PI = P(I,IW) DB00 70 B1 = B1 + PI DB00 75 B2 = B2 + PI*PI DB00 80 5 CONTINUE DB00 85 EIK = ABNV*(B1-P(KLM,IW)) DB00 90 C DB00 95 C TM IS THE MEAN (P(I,IW)*P(J,IW)) I .GT. J DB00 100 C DB00 105 TM = ABNV*0.5*(B1*B1-B2) DB00 110 DUDPI = 0.0 DB00 115 DVDPI = 0.0 DB00 120 PI = P(KLM,IW) DB00 125 DO 15 J=1,N DB00 130 PJ = P(J,IW) DB00 135 PIPJ = PI*PJ DB00 140 IF (J.EQ.KLM) GO TO 10 DB00 145 DVDPI = DVDPI + (PJ-EIK)*(PIPJ-TM) DB00 150 10 DUDPI = DUDPI + PJ*PJ*(PIPJ*(2.0*PIPJ-3.0)+1.0) DB00 155 15 CONTINUE DB00 160 DUDPI = 2.0*PI*DUDPI DB00 165 DVDPI = 2.0*DVDPI DB00 170 DB = (V*DUDPI-U*DVDPI)/(V*V) DB00 175 RETURN DB00 180 END DB00 185 C$ FORTY FORM,XREF COMP 5 SUBROUTINE COMPW(IW, IAV) COMP 10 DOUBLE PRECISION DCON, DELBAR COMP 15 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), COMP 20 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, COMP 25 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC COMP 30 COMMON /A4/ DAPI(030), BCONST, TM, U, V COMP 35 COMMON /A6/ DCON, DATA(030,030), ALOS, IXV COMP 40 C COMP 45 C THIS SUBROUTINE COMPUTES THE ESTIMATE OF THE WEIGHT FOR COMP 50 C SUBSET IW AS WELL AS THE INNER ITERATION CONSTANT. THEN COMP 55 C TRANSLATE THE RESIDUALS BY THE CONSTANT FOR THE FIRST COMP 60 C CALL WHEN IAV = 0 COMP 65 C COMP 70 C COMP 75 C TM IS THE MEAN (P(I,IW)*P(J,IW)) I .GT. J COMP 80 C COMP 85 C COMP 90 C COMPUTE THE MEAN AND SQUARIANCE OF PAIRWISE PRODUCTS COMP 95 C IN STRAIGHTFORWARD MANNER TO AVOID TRUNCATION. COMP 100 C COMP 105 TM = 0.0 COMP 110 AW = 0.0 COMP 115 DO 10 I=2,N COMP 120 TA = P(I,IW) COMP 125 I1 = I - 1 COMP 130 DO 5 J=1,I1 COMP 135 TN = TA*P(J,IW) COMP 140 TM = TM + TN COMP 145 AW = AW + DATA(J,I)*TN COMP 150 5 CONTINUE COMP 155 10 CONTINUE COMP 160 C COMP 165 C DENOMINATOR OF RISE(IW) IS THE SUM OF SQUARES ABOUT TM. COMP 170 C COMP 175 TM = TM*ABNV COMP 180 SQA = 0.0 COMP 185 DO 20 I=2,N COMP 190 PI = P(I,IW) COMP 195 I1 = I - 1 COMP 200 DO 15 J=1,I1 COMP 205 SQA = SQA + (TM-PI*P(J,IW))**2 COMP 210 15 CONTINUE COMP 215 20 CONTINUE COMP 220 RISE(IW) = AW/SQA COMP 225 C COMP 230 C THE REST OF THE COMPUTATION IS NOT NECESSARY FOR THE COMP 235 C CLEANUP CALL TO COMPW. COMP 240 C COMP 245 IF (IAV.EQ.1) RETURN COMP 250 C COMP 255 C COMP 260 C NOW FOR COMPUTING THE REGRESSION CONSTANT FOR SUBSET IW. COMP 265 C COMP 270 DELBAR = 0.0D0 COMP 275 DO 30 I=2,N COMP 280 I1 = I - 1 COMP 285 DO 25 J=1,I1 COMP 290 DELBAR = DELBAR + DATA(J,I) COMP 295 25 CONTINUE COMP 300 30 CONTINUE COMP 305 DELBAR = ABNV*DELBAR COMP 310 WCON = -RISE(IW)*TM + DELBAR COMP 315 C COMP 320 C TRANSLATE THE RESIDUALS. COMP 325 C COMP 330 DO 40 I=2,N COMP 335 I1 = I - 1 COMP 340 DO 35 J=1,I1 COMP 345 DATA(J,I) = DATA(J,I) - WCON COMP 350 35 CONTINUE COMP 355 40 CONTINUE COMP 360 RETURN COMP 365 END COMP 370 C$ FORTY FORM,XREF PRIN 5 SUBROUTINE PRINTD(MST, IW) PRIN 10 COMMON /A1/ P(030,025), DENS(050), NW(050), MP(050), BLMX(050), PRIN 15 * BLMN(050), KP(050), AB, WCON, BIN, N2, N, LD2, N1, AN1, LOUT, PRIN 20 * VARSIM, VAF, ABNV, NDO, RISE(025), LD1, IFLOP, KMOD, KC PRIN 25 COMMON /A7/ A(030), IZ(030) PRIN 30 COMMON /A9/ ILAB(030,6), ALD1 PRIN 35 DIMENSION LSTIM(030) PRIN 40 C PRIN 45 C THIS SUBROUTINE PRINTS OUT THE STIMULI CONTAINED IN SUBSETS. PRIN 50 C PRIN 55 C IF KC .EQ. 0, SUBSET IW IS PRINTED. IF KC .EQ. 1, SUBSETS 1,...NPRIN 60 C ARE PRINTED. IF KC .EQ. 2, THE SUBSETS ARE PRINTED AFTER PRIN 65 C BEING RANKED ACCORDING TO WEIGHT. PRIN 70 C PRIN 75 WRITE (LOUT,99999) PRIN 80 99999 FORMAT (/1X, 6HSUBSET, 4X, 6HWEIGHT, 3X, 21HSTIMULI CONTAINED IN ,PRIN 85 * 6HSUBSET) PRIN 90 IF (KC.NE.0) GO TO 5 PRIN 95 IA = IW PRIN 100 IB = IW PRIN 105 IZ(IW) = IW PRIN 110 GO TO 15 PRIN 115 5 IA = 1 PRIN 120 IB = LD1 PRIN 125 IF (KC.EQ.2) GO TO 15 PRIN 130 DO 10 KX=1,LD1 PRIN 135 IZ(KX) = KX PRIN 140 10 CONTINUE PRIN 145 15 DO 25 KX=IA,IB PRIN 150 LX = IZ(KX) PRIN 155 LCNT = 0 PRIN 160 DO 20 I=1,N PRIN 165 IF (P(I,LX).LE.0.00001) GO TO 20 PRIN 170 LCNT = LCNT + 1 PRIN 175 LSTIM(LCNT) = I PRIN 180 20 CONTINUE PRIN 185 C PRIN 190 C IF YOUR COMPILER CANNOT HANDLE THE ILAB PART OF THE FOLLOWING PRIN 195 C WRITE STATEMENT, THEN WRITING PRIN 200 C PRIN 205 C (LSTIM(KK), KK = 1, LCNT) PRIN 210 C PRIN 215 C IN I FORMAT WILL LIST THE ORDINAL NUMBERS (1, ..., N) OF THE PRIN 220 C STIMULI (INSTEAD OF THEIR USER-SUPPLIED LABELS) FOR EACH PRIN 225 C SUBSET. ALTERNATIVELY, A SCRATCH ARRAY COULD BE CREATED FOR PRIN 230 C THE LINE OF LABELS BEING PRINTED IN THE PRESENT ARRANGEMENT. PRIN 235 C PRIN 240 WRITE (LOUT,99998) LX,RISE(LX), PRIN 245 * ((ILAB(LSTIM(KK),L), L = 1, 6), KK = 1, LCNT) PRIN 250 99998 FORMAT (/1X, 1H(, I4, 1H), F10.5, 2X, 13(6A1, 1X)/(19X, 6A1, 1X, PRIN 255 * 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1, 1X, PRIN 260 * 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1/)) PRIN 265 IF (KC.EQ.1 .AND. MST.GT.0) DENS(MST) = DENS(MST) + FLOAT(LCNT) PRIN 270 25 CONTINUE PRIN 275 IF (KC.EQ.0) RETURN PRIN 280 IF (MST.GT.0) GO TO 30 PRIN 285 WRITE (LOUT,99997) RISE(LD2) PRIN 290 99997 FORMAT (/31H PLUS AN ADDITIVE CONSTANT OF , F10.5) PRIN 295 RETURN PRIN 300 30 IF (KC.EQ.1) DENS(MST) = DENS(MST)*BIN/ALD1 PRIN 305 WRITE (LOUT,99996) RISE(LD2), DENS(MST) PRIN 310 99996 FORMAT (/30H PLUS AN ADDITIVE CONSTANT OF , F10.5, 15X, 7HDENSITY,PRIN 315 * 2H =, F10.5) PRIN 320 RETURN PRIN 325 END PRIN 330 CUT HERE------------ mapclus.f echo processing file mapclus.data 1>&2 cat > mapclus.data <<'CUT HERE------------ mapclus.data' # mapclus.data Sample data. C The authors of this software are Phipps Arabie and J Douglas Carroll 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 SECOND MDS Package of AT&T Bell Laboratories C Remove up to and including the following line before use. C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ 012 16 41 0 0 1 0 0 BLANK CONFUSION FREQUENCIES, ADAPTED FROM GIBSON, ET AL. (1963) 26 1 (25F3.0) 0. 0. 1. 0. 0. 1. 0. 4. 0. 0. 0. 1. 2. 1.15. 0. 6.16. 1. 1. 2. 2. 1. 1. 4. 8. 2. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 5. 1. 2. 1. 4. 0. 1. 0. 0. 0. 1. 0. 2. 0. 3. 0. 0. 0. 0. 1. 0. 2. 2. 0. 0.10. 1. 2. 0. 1. 0. 1. 1. 0. 0. 5. 0. 0. 1. 0. 1. 0. 1. 2. 1. 3. 2. 3. 0. 2. 6. 0.21. 1. 1. 6. 8. 0. 0. 2. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 2. 2. 2. 1. 1. 2. 0. 0. 0. 0. 0. 0. 0. 0. 3. 5. 1. 1. 4. 1. 0. 4. 1. 0. 0. 0.14. 5. 1. 6. 1. 0. 1. 3. 5. 3. 0. 0. 2. 1. 1. 0. 1.11. 0. 2. 1. 3. 0. 1. 1. 4. 1. 0. 3. 1. 1. 0. 1. 0. 1. 0. 1. 0. 0. 0. 1. 1. 0. 1. 1. 5. 0. 0. 9. 0. 1. 0. 1. 0. 0. 0. 0. 0. 7. 0. 0. 0. 0. 1. 0. 5. 0. 1. 1. 0. 2. 0. 0. 0. 1. 2. 4. 1. 1. 0. 0. 0. 0. 1. 0. 1. 3. 1. 0. 1. 0. 0. 2. 0. 1. 1. 5. 0. 0. 0. 0. 2. 0. 0. 1. 0. 0. 2. 0.21.12. 0. 1. 1. 0. 0. 0. 0. 1. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0.15. 1. 2. 2. 0. 1. 0. 0. 1. 0. 0. 6. 3. 1. 0. 0. 0. 0. 1. 0. 1. 0. 0. 5. 0. 0. 0. 2. 1. 0. 2. 0. 3. 2.14. 0.10. 2. 1. 1. 1. 0. 0. 1. 1. 0. 2. 1. 4. 2. 0. 1. 0. 0. 0. 2. 0. 1. 0. 0. 4. 3. A B C D E F G H I J K L M N O P Q R S T U V W X Y Z E. J. GIBSON'S PROPOSED FEATURES FOR UPPER CASE LETTERS (26F1.0) X 1 1111 1 1 1 HORIZONTAL 1 111 11 1111 1 1 1 1 VERTICAL 1 1 1 11111 DIAGONAL (FORWARD) 1 1 11 11 1111 DIAGONAL (BACK) 1 1 1111 CLOSED CURVE 1 1 OPEN-VERTICAL 1 1 1 1 OPEN-HORIZONTAL 11 11 1 1 111 1 1 INTERSECTION 1 1 1 1 1 CYCLIC CHANGE 11111 11 1 1 1 111111 SYMMETRY 1 1 11 1 11 1 1 1 1 VERTICAL 11 1 1 1 HORIZONTAL 08 41 35 1 1 17 000000000000000000000000000000000000000000000000000000000000000000000000000000 .60 MILLER-NICELY DATA FOR CONFUSIONS BETWEEN 16 CONSONANTS 160001 (15F5.3) X .229 01 .432 .241 02 .101 .057 .077 03 .124 .079 .084 .423 04 .052 .050 .063 .066 .157 05 .038 .050 .047 .030 .048 .115 06 .022 .013 .018 .046 .045 .024 .012 07 .025 .022 .020 .025 .041 .031 .033 .058 08 .013 .016 .030 .015 .039 .033 .021 .069 .342 09 .016 .022 .020 .035 .040 .023 .020 .210 .059 .054 10 .028 .016 .018 .032 .031 .026 .018 .145 .094 .120 .338 11 .025 .023 .025 .018 .033 .035 .017 .055 .106 .139 .080 .161 12 .019 .017 .019 .007 .017 .022 .012 .027 .089 .125 .029 .033 .136 13 .025 .022 .021 .016 .019 .017 .012 .038 .024 .032 .030 .034 .021 .016 14 .017 .018 .020 .012 .018 .013 .011 .024 .032 .030 .022 .028 .016 .030 .151 15 PA TA KA FA THETA SA SHA BA DA GA VA THAT ZA ZHA MA NA CUT HERE------------ mapclus.data #define END