#!/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 indclus.f 1>&2 cat > indclus.f <<'CUT HERE------------ indclus.f' C indclus.f C The authors of this software are J.Douglas Carroll and Phipps Arabie. 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 W DeSarbo, C a monograph published by Sage Publications, Beverly Hills, Calif. C in 1987 as series 07, item 065, in the Sage University Papers, and C "Additive clustering: Representation of similarities as combinations of C discrete overlapping properties" by R N Shepard and P Arabie in C Psychological Review, 1979, vol 86, pages 87-123. C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ CC$ FORTY FORM,XREF INDC 5 C INDCLUS, DECEMBER 1992 (MAINFRAME VERSION) INDC 10 C PROGRAM INDCLUS(PARAM,OUTPUT,PNCH,MATRIX,INIC,TAPE5=PARAM, INDC 15 C 2 TAPE6=OUTPUT,TAPE7=PNCH,TAPE41=MATRIX,TAPE16=INIC) INDC 20 C INDC 25 C THE FOLLOWING NON-ANSI STANDARD FORTRAN CHARACTERS (: ;) ARE USED INDC 30 C FOR COLON AND SEMICOLON, RESPECTIVELY. IF YOUR COMPILER OR INDC 35 C OPERATING SYSTEM CANNOT DIGEST THESE CHARACTERS, SUBSTITUTING INDC 40 C A BLANK IS THE BEST REMEDY. INDC 45 C INDC 50 DOUBLE PRECISION DCON INDC 55 INTEGER FMT, TITL INDC 60 DIMENSION FMT(72), SVAF(050), IDAT(2), IDIG(10), ING(2,4), INDC 65 * GL(030,2), GRLEN(2), GRMEN(2), OBSLF(2), ASTEP(2) INDC 70 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, INDC 75 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, INDC 80 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, INDC 85 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON INDC 90 COMMON /A2/ G(018,030), KT1(10), KT2(10), TK(10), LP(15), IPUN, INDC 95 * IERR, LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), INDC 100 * TRISE(018,017), E(018) INDC 105 COMMON /A4/ DAPI(030), BCONST, TM, U, V INDC 110 COMMON /A5/ AK, ALFRAY(017), BETRAY(017), IADJ(050), JADJ INDC 115 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), INDC 120 * ALOS, IXV INDC 125 COMMON /A7/ AZ(030), IZA(030), A(017), IZ(017) INDC 130 COMMON /A8/ W(018,017), IPRE, NNEG, NDOT, NLD1 INDC 135 COMMON /A9/ ILAB(030,6), MLAB(018,6), ALD1, SIGN INDC 140 DATA IDIG(1), IDIG(2), IDIG(3), IDIG(4), IDIG(5), IDIG(6), INDC 145 * IDIG(7), IDIG(8), IDIG(9), IDIG(10), IBL, IL, IR, IHY /1H0,1H1, INDC 150 * 1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1H ,1H(,1H),1H-/ INDC 155 DATA ING(1,1), ING(1,2), ING(1,3), ING(1,4), ING(2,1), ING(2,2), INDC 160 * ING(2,3), ING(2,4) /1H ,1HN,1HE,1HG,1HZ,1HE,1HR,1HO/ INDC 165 CX IDAT(1) = IDATEZ(0) INDC 170 CX XX = ITIMEZ(0)/100000. INDC 175 C INDC 180 C INDC 185 C ASSUME CARD INPUT STREAM AS LOGICAL DEVICE FOR READING PROGRAM INDC 190 C CONTROL PARAMETERS. INDC 195 C INDC 200 L5 = 5 INDC 205 C INDC 210 C INDC 215 LOUT = 6 INDC 220 C INDC 225 C THE NEXT STATEMENT ASSUMES THAT THE DEFAULT DEVICE FOR PUNCHED INDC 230 C OUTPUT IS 7. THERE IS NO PROVISION FOR REWINDING. INDC 235 C INDC 240 LPUN7 = 7 INDC 245 C INDC 250 C INDC 255 C INDC 260 C INDC 265 C INDC 270 C INDC 275 C NDOT IS THE UPPER BOUND FOR THE NUMBER OF STIMULI TAKEN TWO INDC 280 C AT A TIME (LATER COMPUTED AS NDO, BASED ON THE USER'S INDC 285 C DECLARED UPPER BOUND FOR N), AND NLD1 IS THE UPPER BOUND INDC 290 C FOR THE NUMBER OF CLUSTERS NOT INCLUDING THE COMPLETE SET INDC 295 C FOR THE ADDITIVE CONSTANT. THESE TWO VALUES ARE NECESSARY INDC 300 C TO CHECK FOR DIMENSIONS OF ARRAYS USED IN SUBROUTINE REGR IF INDC 305 C THE NON-NEGATIVE LEAST SQUARES ROUTINE IS CALLED BY REGR, INDC 310 C WHICH CAN HAPPEN ONLY IF NNEG .GT. 0 INDC 315 C INDC 320 C INDC 325 NDOT = (0435) INDC 330 NLD1 = (016) INDC 335 C INDC 340 C NLD1 SHOULD BE .GE. (N2) - 1 IF NON-NEGATIVE REGRESSION IS TO BE INDC 345 C USED. INDC 350 C INDC 355 C N2 IS THE MAXIMUM NUMBER OF CLUSTERS (INCLUDING COMPLETE INDC 360 C SUBSET FOR ADDITIVE CONSTANT) ALLOWED BY CURRENT DIMENSION INDC 365 C STATEMENTS. INDC 370 C INDC 375 N2 = (017) INDC 380 C INDC 385 C INDC 390 C M3TOP IS THE MAXIMUM NUMBER OF SOURCES OF DATA (E.G., SUBJECTS) INDC 395 C ALLOWED BY CURRENT DIMENSION STATEMENTS INDC 400 C INDC 405 M3TOP = (018) INDC 410 N5 = (050) INDC 415 C INDC 420 C N30 IS THE MAXIMUM NUMBER OF STIMULI (N) ALLOWED BY CURRENT INDC 425 C DIMENSION STATEMENTS. INDC 430 C INDC 435 N30 = (030) INDC 440 C INDC 445 C INDC 450 C INDC 455 C INDC 460 BCRIT = 1.0E-5 INDC 465 JTPRE = 6 INDC 470 C INDC 475 C FOR INSTALLATIONS OTHER THAN BTL, LREAD SHOULD BE 5 INDC 480 C INDC 485 5 LREAD = 5 INDC 490 IPRE = 0 INDC 495 VAF = 0.0 INDC 500 ICON = 0 INDC 505 ISWIT = 0 INDC 510 IPRN = 0 INDC 515 IXV = 0 INDC 520 DO 10 I=1,15 INDC 525 LP(I) = 0 INDC 530 10 CONTINUE INDC 535 C INDC 540 C INDC 545 WRITE (LOUT,99999) INDC 550 99999 FORMAT (1H1///11X, 27HRUN OF INDCLUS VERSION 1.2, 8X, 8H A MATHE,INDC 555 * 28HMATICAL PROGRAMMING APPROACH, 28H BY J. D. CARROLL AND P. ARA,INDC 560 * 19HBIE FOR FITTING THE/46X, 33H INDIVIDUAL DIFFERENCES (3 - WAY),INDC 565 * 37H GENERALIZATION OF THE SHEPARD-ARABIE, 13H ADCLUS MODEL//) INDC 570 CX WRITE (LOUT,8) IDAT(1),XX INDC 575 CX 8 FORMAT(30X,A6,5X,F8.6) INDC 580 C INDC 585 C INDC 590 C INDC 595 C INDC 600 READ (L5,99998,END=655) LD1,ICON,KREAD,IREWIN,M15,NNEG,IPRN,IPUN, INDC 605 * LPUNCH,(LP(I), I = 1, 15) INDC 610 99998 FORMAT (9I3, 1X, 15I1) INDC 615 WRITE (LOUT,99997) LD1, ICON, KREAD, IREWIN, M15, NNEG, IPRN, INDC 620 * IPUN, LPUNCH, (I,LP(I),I=1,15) INDC 625 99997 FORMAT (//4X, 48H LD1 ICON KREAD IREWIN M15 NNEG IPRN IPU,INDC 630 * 9HN LPUNCH, 48H ABBREVIATED -0- AND FULL-LENGTH -1- PRINTI,INDC 635 * 12HNG DURING 15/2X, 3I6, 2X, I6, 1X, I5, 4I6, 13X, 10H ITERATION,INDC 640 * 4HS OF, 27H COMBINATORIAL OPTIMIZATION//50X, 15(1X, I2, 1H:, I1))INDC 645 C INDC 650 C LD1 IS THE NUMBER OF CLUSTERS TO BE FITTED, AS SPECIFIED BY INDC 655 C THE USER, AND DOES NOT INCLUDE THE COMPLETE SET ASSOCIATED INDC 660 C WITH THE ADDITIVE CONSTANT. INDC 665 C INDC 670 WRITE (LOUT,99996) LD1 INDC 675 99996 FORMAT (//10X, 36HNUMBER OF CLUSTERS (LD1) SPECIFIED =, I3, INDC 680 * 60H (NOT INCLUDING THE COMPLETE SET FOR THE ADDITIVE CONSTANT).) INDC 685 IF (LD1.GT.0) GO TO 15 INDC 690 WRITE (LOUT,99995) INDC 695 99995 FORMAT (///46H THIS PROGRAM WILL NOT WORK CORRECTLY IF LD1, , INDC 700 * 51HTHE NUMBER OF SUBSETS IS NOT DECLARED AS A POSITIVE, 6H INTEG,INDC 705 * 3HER./22H EXECUTION TERMINATES.) INDC 710 99994 FORMAT (//50H CONSULT DOCUMENTATION FOR INDCLUS FOR DETAILS ON , INDC 715 * 30HCHANGING DIMENSION STATEMENTS.//) INDC 720 STOP INDC 725 C INDC 730 C INDC 735 C IF ICON IS NEGATIVE, A RANDOM INITIAL CONFIGURATION INDC 740 C WILL BE GENERATED, AND THE FIRST (-ICON) NUMBERS INDC 745 C FROM THE RANDOM NUMBER GENERATOR WILL BE DISCARDED. INDC 750 C IF ICON IS ZERO, A RATIONAL INITIAL CONFIGURATION INDC 755 C WILL BE PRODUCED BY THE PROGRAM. INDC 760 C IF ICON IS POSITIVE, A USER-SUPPLIED INITIAL INDC 765 C CONFIGURATION WILL BE READ FROM LOGICAL DEVICE INDC 770 C NUMBER (ICON). INDC 775 C INDC 780 15 IF (KREAD.GT.0) GO TO 20 INDC 785 C INDC 790 C KREAD IS THE LOGICAL DEVICE NUMBER FOR THE NUMERICAL INDC 795 C PROXIMITIES DATA AND THE (OPTIONAL) LABELS. THE DEFAULT VALUE INDC 800 C IS 5. INDC 805 C INDC 810 WRITE (LOUT,99993) KREAD, LREAD INDC 815 99993 FORMAT (//10X, 40HSINCE THE DATA SET WAS DECLARED TO BE ON, INDC 820 * 15H DEVICE NUMBER , I4, 30H (KREAD) THE DEFAULT VALUE OF , I3, INDC 825 * 22H WILL INSTEAD BE USED.) INDC 830 GO TO 25 INDC 835 20 LREAD = KREAD INDC 840 25 WRITE (LOUT,99992) LREAD INDC 845 99992 FORMAT (//9X, 46H DATA WILL BE READ FROM LOGICAL DEVICE NUMBER , INDC 850 * I3) INDC 855 IF (IREWIN.EQ.1 .AND. LREAD.NE.L5) GO TO 30 INDC 860 C INDC 865 C IREWIN ALLOWS THE USER THE OPTION OF HAVING THE LOGICAL DEVICE INDC 870 C ASSOCIATED WITH KREAD REWOUND, UNLESS LREAD .EQ. L5 INDC 875 C INDC 880 WRITE (LOUT,99991) INDC 885 99991 FORMAT (9X, 27H WHICH WILL NOT BE REWOUND.) INDC 890 GO TO 35 INDC 895 30 REWIND LREAD INDC 900 WRITE (LOUT,99990) INDC 905 99990 FORMAT (9X, 24H WHICH HAS BEEN REWOUND.) INDC 910 35 WRITE (LOUT,99989) LREAD, LOUT INDC 915 99989 FORMAT (//9X, 43H I/O CHANNELS FOR READING DATA AND PRINTING, INDC 920 * 31H OUTPUT: READ PRINT/58X, 12H LREAD =, I3, INDC 925 * 9H LOUT =, I3) INDC 930 C INDC 935 C IF M15, THE NUMBER OF MAJOR ITERATIONS, IS ZERO, ONLY (GLOBAL) INDC 940 C REGRESSION WILL BE PERFORMED ON AN INITIAL CONFIGURATION, INDC 945 C PRESUMABLY SUPPLIED BY THE USER ON CHANNEL ICON. INDC 950 C INDC 955 C IF M15 IS NEGATIVE, GLOBAL REGRESSION IS PERFORMED (AS WHEN M15 INDC 960 C .EQ. 0), FOLLOWED IMMEDIATELY BY COMBINATORIAL OPTIMIZATION, INDC 965 C UP TO A MAXIMUM OF 15 ITERATIONS. INDC 970 C INDC 975 C IF M15 IS POSITIVE, THEN IT DEFINES THE MAXIMUM POSSIBLE INDC 980 C NUMBER OF MAJOR ITERATIONS (TOTAL OF PRE-POLISHING AND INDC 985 C POLISHING). INDC 990 C INDC 995 C INDC1000 WRITE (LOUT,99988) M15 INDC1005 99988 FORMAT (//10X, 47HMAXIMUM NUMBER OF MAJOR ITERATIONS (M15) SPECIF,INDC1010 * 3HIED, 5H IS =, I3, 1H.) INDC1015 IF (M15.LE.N5) GO TO 40 INDC1020 WRITE (LOUT,99987) M15, N5 INDC1025 99987 FORMAT (//20H YOU HAVE SPECIFIED , I4, 23H ITERATIONS, BUT THIS V,INDC1030 * 6HERSION, 40H OF THE PROGRAM ALLOWS FOR AT MOST ONLY , I4, INDC1035 * 14H ITERATIONS. /43H THE DIMENSION STATEMENT COVERING ARRAY KP ,INDC1040 * 12HIN THE MAIN , 44H PROGRAM WILL HAVE TO BE INCREASED. EXECUTI,INDC1045 * 14HON TERMINATES.) INDC1050 WRITE (LOUT,99994) INDC1055 STOP INDC1060 C INDC1065 C NNEG DETERMINES THE POLICY ON NEGATIVE WEIGHTS. IF NNEG=O, INDC1070 C SUCH WEIGHTS ARE CONSIDERED NO MORE OBJECTIONABLE THAN POSITIVE INDC1075 C ONES. IF NNEG=1, NEGATIVE WEIGHTS ARE DEALT WITH AT THE END OF INDC1080 C EACH MAJOR ITERATION. IF NNEG=2, NEGATIVE WEIGHTS GO TO INDC1085 C THE ROUTINE FROM LAWSON AND HANSON WHENEVER GLOBAL REGRESSION IS INDC1090 C CALLED. INDC1095 C INDC1100 C INDC1105 C INDC1110 40 IF (NNEG-1) 45, 50, 55 INDC1115 45 NNEG = 0 INDC1120 WRITE (LOUT,99986) INDC1125 99986 FORMAT (//10X, 33HSINCE NNEG = 0, THERE WILL BE NO , 9HNON-NEGAT, INDC1130 * 33HIVITY CONSTRAINTS ON THE WEIGHTS.) INDC1135 GO TO 60 INDC1140 50 WRITE (LOUT,99985) INDC1145 99985 FORMAT (//10X, 36HSINCE NNEG = 1, THE WEIGHTS WILL BE , 7HCONSTRA,INDC1150 * 42HINED TO BE NON-NEGATIVE AT THE END OF EACH/10X, 10HMAJOR ITER,INDC1155 * 6HATION.) INDC1160 GO TO 60 INDC1165 55 NNEG = 2 INDC1170 WRITE (LOUT,99984) INDC1175 99984 FORMAT (//10X, 36HSINCE NNEG = 2, THE WEIGHTS WILL BE , 7HCONSTRA,INDC1180 * 42HINED TO BE NON-NEGATIVE AT THE END OF EACH/10X, 10HOUTER ITER,INDC1185 * 6HATION.) INDC1190 C INDC1195 C INDC1200 C INDC1205 C IF IPRN .EQ. 1 THE INPUT DATA WILL BE LISTED (A GOOD IDEA ON INDC1210 C THE FIRST TIME A GIVEN DATA SET IS ENTERED). INDC1215 C INDC1220 C IF IPUN .EQ. 1, THE FINAL CONFIGURATION WILL BE WRITTEN ON INDC1225 C DISK (OR PUNCHED), USING LOGIAL DEVICE (LPUNCH). INDC1230 C INDC1235 60 IF (IPUN.EQ.0) GO TO 65 INDC1240 IF (LPUNCH.EQ.0) LPUNCH = LPUN7 INDC1245 WRITE (LOUT,99983) LPUNCH INDC1250 99983 FORMAT (//10X, 47HCONFIGURATION OF SUBSETS WILL BE WRITTEN TO DIS,INDC1255 * 1HK, 12H ON CHANNEL , I3, 35H AT END OF (MAJOR) POST-ITERATIONS.)INDC1260 C INDC1265 C INDC1270 C GET THE INITIAL CONFIGURATION OF P(I). INDC1275 C INDC1280 65 IF (ICON) 70, 80, 85 INDC1285 70 IX = -ICON INDC1290 WRITE (LOUT,99982) IX INDC1295 99982 FORMAT (//10X, 41HSINCE ICON IS NEGATIVE, A RANDOM INITIAL , INDC1300 * 31HCONFIGURATION WILL BE GENERATED/10X, 13HAND THE FIRST, I5, INDC1305 * 34H RANDOM NUMBERS WILL BE DISCARDED.//) INDC1310 DO 75 I=1,IX INDC1315 AMN = RUNIFV(0) INDC1320 75 CONTINUE INDC1325 GO TO 90 INDC1330 80 WRITE (LOUT,99981) INDC1335 99981 FORMAT (//10X, 47HSINCE ICON = 0, A RATIONAL INITIAL CONFIGURATIO,INDC1340 * 1HN, 37H WILL BE GENERATED (DEFAULT OPTION).//) INDC1345 GO TO 90 INDC1350 85 WRITE (LOUT,99980) ICON INDC1355 99980 FORMAT (//10X, 47HSINCE ICON IS POSITIVE, A USER-SUPPLIED INITIAL,INDC1360 * 51H CONFIGURATION WILL BE READ FROM LOGICAL DEVICE NO., I4, 1H.//INDC1365 * ) INDC1370 90 IF (M15.GT.0) GO TO 100 INDC1375 C INDC1380 IF (ICON.GT.0) GO TO 95 INDC1385 WRITE (LOUT,99979) M15, ICON INDC1390 99979 FORMAT (//51H USER HAS SPECIFIED INCOMPATIBLE OPTIONS. ALTHOUGH, INDC1395 * 47H COMPUTATION IS TO START WITH REGRESSION (M15 =, I3, 2H),/ INDC1400 * 61H NO USER-SUPPLIED INITIAL CONFIGUATION HAS BEEN GIVEN (ICON =,INDC1405 * I3, 25H). EXECUTION TERMINATES.) INDC1410 STOP INDC1415 95 IF (M15.EQ.0) GO TO 110 INDC1420 WRITE (LOUT,99978) INDC1425 99978 FORMAT (//10X, 40HSINCE M15 IS NEGATIVE, THERE WILL BE NO , INDC1430 * 45HITERATIVE COMPUTATION OF WEIGHTS AND SUBSETS,/10X, 8HAND THE ,INDC1435 * 41HPROGRAM WILL PROCEED TO DO COMBINATORIAL , 15HOPTIMIZATION IM,INDC1440 * 27HMEDIATELY AFTER REGRESSION.) INDC1445 GO TO 110 INDC1450 C INDC1455 C IF THERE ARE TO BE NO ITERATIONS OF THE MATHEMATICAL PROGRAMMING INDC1460 C ALGORITHM, THEN DON'T READ IN THE REMAINING (IRRELEVANT) VARIABLESINDC1465 C INDC1470 C INDC1475 C INDC1480 C INDC1485 C INDC1490 100 READ (L5,99977) (KP(I),I=1,M15) INDC1495 99977 FORMAT (72I1) INDC1500 C INDC1505 C KP IS A BINARY VECTOR THAT CONTROLS THE AMOUNT OF INDC1510 C PRINTED OUTPUT FOR EACH MAJOR ITERATION. WHEN M15 .LE. 0, KP IS INDC1515 C IRRELEVANT AND IS THEREFORE NOT READ IN. INDC1520 C INDC1525 WRITE (LOUT,99976) INDC1530 99976 FORMAT (///6X, 47HABBREVIATED (0) AND FULL-LENGTH (1) PRINTOUT FO,INDC1535 * 2HR , 17HMAJOR ITERATIONS:) INDC1540 WRITE (LOUT,99975) (I, I = 1, M15) INDC1545 99975 FORMAT (/6X, 35I3) INDC1550 WRITE (LOUT,99975) (KP(I),I=1,M15) INDC1555 C INDC1560 C INDC1565 C INDC1570 C INDC1575 C INDC1580 C INDC1585 READ (L5,99998) ITPRE, IUPPER, JLIM, JPOST INDC1590 WRITE (LOUT,99974) ITPRE, IUPPER, JLIM, JPOST INDC1595 99974 FORMAT (//6X, 27HITPRE IUPPER JLIM JPOST/5X, 2I6, 3(1X, I6)) INDC1600 C INDC1605 C ITPRE IS THE MAXIMUM NUMBER OF MAJOR ITERATIONS ALLOWED PRIOR TO INDC1610 C POLISHING. (I. E., PRE-POLISHING ITERATIONS.) IF ITPRE IS INDC1615 C NEGATIVE, ITERATIVE COMPUTATION BEGINS WITH DE NOVO MAJOR INDC1620 C ITERATIONS. INDC1625 C INDC1630 C IF ITPRE IS ZERO, THE DEFAULT VALUE OF 6 IS SUPPLIED. INDC1635 C THERE IS NO OPTION TO BEGIN COMPUTATION WITH THE FIRST INDC1640 C POLISHING ITERATION, BUT ITPRE.LT.0 STARTS COMPUTATION WITH INDC1645 C DE NOVO ITERATIONS. INDC1650 C INDC1655 C INDC1660 C IUPPER IS THE MAXIMUM NUMBER OF INNER ITERATIONS IN THE OUTER INDC1665 C ITERATIONS OF THE PRE-POLISHING MAJOR ITERATIONS. INDC1670 C INDC1675 C JLIM IS THE MAXIMUM NUMBER OF TIMES ALPHA AND BETA CAN BE INDC1680 C ADJUSTED DURING AN OUTER ITERATION OF A PRE-POLISHING MAJOR INDC1685 C ITERATION. INDC1690 C INDC1695 C JPOST IS THE MAXIMUM NUMBER OF INNER ITERATIONS OF AN OUTER INDC1700 C ITERATION OF A POLISHING MAJOR ITERATION. INDC1705 C INDC1710 C INDC1715 C INDC1720 C INDC1725 C THE NEXT STATEMENT ALLOWS THE OPTION OF STARTING A RUN WITH INDC1730 C DE NOVO MAJOR ITERATIONS, AND PROBABLY SHOULD BE USED ONLY INDC1735 C FOR RESUMING AN EARLIER ANALYSIS. INDC1740 C INDC1745 KTPRE = ITPRE INDC1750 IF (ITPRE.GE.0) GO TO 105 INDC1755 IPRE = 1 INDC1760 WRITE (LOUT,99973) ITPRE INDC1765 99973 FORMAT (//10X, 42HUSER HAS DECLARED THAT ITPRE IS NEGATIVE (, I3, INDC1770 * 52H), SO THAT COMPUTATION WILL BEGIN WITH DE NOVO MAJOR, INDC1775 * 12H ITERATIONS.) INDC1780 IF (ICON.LT.0) WRITE (LOUT,99972) ICON INDC1785 99972 FORMAT (//10X, 47HUSER HAS SPECIFIED RANDOM INITIAL CONFIGURATION,INDC1790 * 8H (ICON =, I3, 45H) WHICH IS INCONSISTENT WITH A DE NOVO START;/INDC1795 * 10X, 47HINITIAL CONFIGURATION WILL INSTEAD BE RATIONAL.) INDC1800 105 IF (ITPRE.EQ.0) ITPRE = JTPRE INDC1805 ITPRE1 = ITPRE + 1 INDC1810 IF (IUPPER.LE.0) IUPPER = 50 INDC1815 IF (ITPRE.GT.0) WRITE (LOUT,99971) ITPRE, IUPPER INDC1820 99971 FORMAT (//10X, 34HMAXIMUM NUMBER OF MAJOR ITERATIONS, 9H PRIOR TO,INDC1825 * 20H POLISHING (ITPRE) =, I4, 1H,/10X, 17HWITH A MAXIMUM OF, I5, INDC1830 * 26H INNER ITERATIONS (IUPPER), 21H PER OUTER ITERATION.) INDC1835 C INDC1840 IF (JLIM.LE.0) JLIM = 4 INDC1845 KLIM = JLIM INDC1850 IF (JPOST.LE.0) JPOST = 75 INDC1855 IF (ITPRE.GT.0) WRITE (LOUT,99970) KLIM INDC1860 99970 FORMAT (//9X, 20H USER HAS SPECIFIED , 18HA LIMIT (JLIM) OF , I3, INDC1865 * 28H CALL(S) TO ADJUST PER OUTER, 28H ITERATION PRIOR TO POLISHIN,INDC1870 * 3HG. ) INDC1875 WRITE (LOUT,99969) JPOST INDC1880 99969 FORMAT (//10X, 14HTHERE WILL BE , I3, 24H INNER ITERATIONS (JPOST,INDC1885 * 2H) , 37HPER OUTER ITERATION DURING POLISHING.) INDC1890 C INDC1895 C INDC1900 READ (L5,99966) ALF0 INDC1905 C INDC1910 C ALF0 IS THE INITIAL ESTIMATE OF THE STEP SIZE (DEFAULT =1.0) INDC1915 C INDC1920 IF (ALF0.LE.0.0) ALF0 = 1.0 INDC1925 WRITE (LOUT,99968) ALF0 INDC1930 99968 FORMAT (//6X, 6HALF0 =, F12.8, 31H (INTERVAL FOR INITIAL ESTIMATE,INDC1935 * 15H OF STEP SIZE).) INDC1940 C INDC1945 C INDC1950 C INDC1955 C INDC1960 READ (L5,99966) ALPHAI, AK, ACRIT, BLCRIT INDC1965 C INDC1970 WRITE (LOUT,99967) ALPHAI, AK, ACRIT, BLCRIT INDC1975 99967 FORMAT (//6X, 39HALPHAI AK ACRIT BLCRIT/3X, INDC1980 * 4F11.8) INDC1985 C INDC1990 C ALPHAI IS THE WEIGHT FOR THE A-PART OF THE LOSS FUNCTION. INDC1995 C INDC2000 C AK IS THE CONSTANT IN THE FORMULA FOR ADJUSTING ALPHA AND BETA. INDC2005 C INDC2010 C ACRIT IS THE CRITERION OF CONVERGENCE FOR THE RELATIVE GRADIENT. INDC2015 C INDC2020 C BLCRIT IS THE CRITERION OF CONVERGENCE FOR THE B-PART OF THE INDC2025 C LOSS FUNCTION. INDC2030 C INDC2035 99966 FORMAT (4F10.6) INDC2040 C INDC2045 C ESTABLISHING A DEFAULT VALUE FOR ALPHAI: INDC2050 C INDC2055 IF (ALPHAI.LE.0.0) ALPHAI = .40 INDC2060 BETAI = 1.0 - ALPHAI INDC2065 IF (AK.LE.0.0) AK = 2.0 INDC2070 IF (ACRIT.LE.0.0) ACRIT = 0.005 INDC2075 IF (BLCRIT.LE.0.0) BLCRIT = .05 INDC2080 WRITE (LOUT,99965) ALPHAI, AK, ACRIT, BLCRIT INDC2085 99965 FORMAT (//10X, 8HALPHAI =, F12.5, 28H (FOR WEIGHTING A-PART OF LO,INDC2090 * 13HSS FUNCTION).//10X, 5HAK = , F12.5, 22H (FOR ADJUSTING ALPHA ,INDC2095 * 10HAND BETA).//10X, 7HACRIT =, F12.5, 23H (CRITERION FOR CONVERG,INDC2100 * 27HENCE OF RELATIVE GRADIENT).//10X, 8HBLCRIT =, F12.5, 6H (CRIT,INDC2105 * 20HERION FOR B-PART OF , 15HLOSS FUNCTION).//) INDC2110 C INDC2115 C INDC2120 C INDC2125 C INDC2130 C INDC2135 C INDC2140 110 READ (LREAD,99963) TITL INDC2145 C INDC2150 C INDC2155 READ (LREAD,99964) M3, MCON, N, ISWIT, LABL INDC2160 99964 FORMAT (5I2) INDC2165 C INDC2170 C INDC2175 READ (LREAD,99963) FMT INDC2180 99963 FORMAT (72A1) INDC2185 C INDC2190 C INDC2195 C M3 IS THE NUMBER OF SUBJECTS, OR OTHER SOURCES OF DATA, INDC2200 C EACH GIVING RISE TO A 2-WAY MATRIX OF PROXIMITIES INDC2205 C FOR THE N STIMULI. INDC2210 C INDC2215 C IF MCON.EQ.0 THE DATA ARE VIEWED AS BEING FROM A INDC2220 C MATRIX 'UNCONDITIONAL' DATA SOURCE. INDC2225 C INDC2230 C IF MCON.EQ.1 THE DATA ARE VIEWED AS BEING FROM A MATRIX INDC2235 C 'CONDITIONAL' DATA SOURCE. INDC2240 C INDC2245 C INDC2250 C N IS THE NUMBER OF STIMULI. INDC2255 C INDC2260 C ISWIT = 1 FOR DISSIMILARITIES DATA. INDC2265 C ISWIT .LE. 0 FOR SIMILARITIES DATA. INDC2270 C INDC2275 C IF ILAB .EQ. 1, THE PROGRAM EXPECTS TO FIND 6-CHARACTER LABELS INDC2280 C (12 TO A CARD, IN COLUMNS 1 - 72) IMMEDIATELY FOLLOWING THE INDC2285 C NUMERICAL DATA; OTHERWISE, INTERNALLY GENERATED LABELS OF 1,...,N INDC2290 C WILL BE USED FOR THE STIMULI. INDC2295 C INDC2300 C INDC2305 IF (N.LE.N30) GO TO 115 INDC2310 WRITE (LOUT,99962) N30, N INDC2315 99962 FORMAT (47H IN THIS VERSION OF INDCLUS, THE MAXIMUM NUMBER, INDC2320 * 14H OF STIMULI IS, I3, 24H, BUT YOU HAVE SPECIFIED, I3, 1H.// INDC2325 * 22H EXECUTION TERMINATES.) INDC2330 WRITE (LOUT,99994) INDC2335 STOP INDC2340 115 IF (M3.GT.0) GO TO 120 INDC2345 WRITE (LOUT,99961) M3 INDC2350 99961 FORMAT (54H THE DATA PARAMETER M3 IS REQUIRED TO BE POSITIVE, BUT,INDC2355 * 30H YOU HAVE SPECIFIED A VALUE OF, I4, 22H. EXECUTION TERMINATES/INDC2360 * /) INDC2365 STOP INDC2370 120 IF (M3.LE.M3TOP) GO TO 125 INDC2375 WRITE (LOUT,99960) M3, M3TOP INDC2380 99960 FORMAT (35H USER HAS SPECIFIED A VALUE OF M3 (, I4, 10H) EXCEEDIN,INDC2385 * 5HG THE, 34H MAXIMUM VALUE CURRENTLY ALLOWED (, I4, 2H)./ INDC2390 * 43H INCREASE DIMENSION STATEMENTS ACCORDINGLY;, 13H EXECUTION TE,INDC2395 * 9HRMINATES.//) INDC2400 WRITE (LOUT,99994) INDC2405 STOP INDC2410 C INDC2415 C INDC2420 C INDC2425 125 WRITE (LOUT,99959) TITL, M3, MCON, N, ISWIT, IPRN, FMT INDC2430 99959 FORMAT (//1H0, 25HTITLE OF INPUT DATA SET: , 72A1///10H M3 MCO,INDC2435 * 16HN N ISWIT LABL/3(2X, I3), I3, 5X, I1//18H SPECIFIED FORMAT ,INDC2440 * 21HOF THE INPUT DATA IS://1X, 72A1) INDC2445 IF (MCON.NE.1) GO TO 130 INDC2450 WRITE (LOUT,99958) INDC2455 99958 FORMAT (//39H SINCE MCON = 1, THE INPUT DATA WILL BE, 9H TREATED ,INDC2460 * 41HAS MATRIX CONDITIONAL. THUS, THE WEIGHTS /14H ARE POSSIBLY ,INDC2465 * 38HON A DIFERENT SCALE FROM THE ONE USED , 18HFOR AN UNCONDITION,INDC2470 * 12HAL ANALYSIS.) INDC2475 GO TO 135 INDC2480 130 MCON = 0 INDC2485 WRITE (LOUT,99957) INDC2490 99957 FORMAT (//41H THE INPUT DATA WILL BE TREATED AS MATRIX, 7H UNCOND,INDC2495 * 8HITIONAL.) INDC2500 135 IF (LABL.EQ.0) WRITE (LOUT,99956) INDC2505 99956 FORMAT (//1X, 46HTHE PROGRAM WILL GENERATE NUMERICAL LABELS FOR, INDC2510 * 13H THE STIMULI.) INDC2515 IF (LABL.EQ.1) WRITE (LOUT,99955) INDC2520 99955 FORMAT (//1X, 48HUSER HAS DECLARED THAT LABELS ARE TO BE READ IN ,INDC2525 * 33HIMMEDIATELY AFTER NUMERICAL DATA.) INDC2530 C INDC2535 LD2 = LD1 + 1 INDC2540 BN = N INDC2545 BIN = 1.0/BN INDC2550 N1 = N - 1 INDC2555 AN1 = N1 INDC2560 NDO = N*N1/2 INDC2565 RM = 1.0/FLOAT(M3) INDC2570 ABM = M3*NDO INDC2575 ABIN = RM/BN INDC2580 ALD1 = BIN/FLOAT(LD1) INDC2585 AMV = 1.0/ABM INDC2590 AB = NDO INDC2595 ABNV = 1.0/AB INDC2600 AMN = 1.0E20 INDC2605 AMX = -1.0E20 INDC2610 C INDC2615 C READ IN LOWERHALFMATRIX BY ROWS INDC2620 DO 150 M=1,M3 INDC2625 DO 145 I=2,N INDC2630 I1 = I - 1 INDC2635 READ (LREAD,FMT) (DATA(M,IV,I),IV=1,I1) INDC2640 DO 140 IV=1,I1 INDC2645 DT = DATA(M,IV,I) INDC2650 IF (DT.GT.AMX) AMX = DT INDC2655 IF (DT.LT.AMN) AMN = DT INDC2660 140 CONTINUE INDC2665 145 CONTINUE INDC2670 150 CONTINUE INDC2675 C INDC2680 C SET UP STIMULUS LABELS OF THE INTEGERS 1,...,N. THESE LABELS INDC2685 C ARE CONSTRUCTED EVEN IF THE USER SUPPLIES LABELS, SO THAT IN INDC2690 C THE EVENT OF AN END OF FILE WHILE ATTEMPTING TO READ THE LATTER INDC2695 C LABELS, THE PROGRAM STILL HAS USABLE LABELS. IT IS ASSUMED INDC2700 C THAT NEITHER N NOR M3 WILL EXCEED 99. INDC2705 C INDC2710 DO 155 I=1,N INDC2715 ILAB(I,1) = IBL INDC2720 ILAB(I,6) = IBL INDC2725 ILAB(I,2) = IL INDC2730 ILAB(I,5) = IR INDC2735 K4 = MOD(I,10) INDC2740 ILAB(I,4) = IDIG(K4+1) INDC2745 K4 = 1 + (I-K4)/10 INDC2750 ILAB(I,3) = IDIG(K4) INDC2755 155 CONTINUE INDC2760 C INDC2765 C NOW GET LABELS FOR THE 'SUBJECTS' INDC2770 C INDC2775 DO 160 M=1,M3 INDC2780 MLAB(M,1) = IHY INDC2785 MLAB(M,2) = IHY INDC2790 MLAB(M,5) = IHY INDC2795 MLAB(M,6) = IHY INDC2800 K4 = MOD(M,10) INDC2805 MLAB(M,4) = IDIG(K4+1) INDC2810 K4 = 1 + (M-K4)/10 INDC2815 MLAB(M,3) = IDIG(K4) INDC2820 160 CONTINUE INDC2825 IF (LABL.NE.1) GO TO 165 INDC2830 READ (LREAD,99954,END=162) ((ILAB(I,J), J = 1, 6), I= 1, N) INDC2835 READ (LREAD,99954,END=162) ((MLAB(M,J), J = 1, 6), M= 1, M3) INDC2840 99954 FORMAT (72A1) INDC2845 GO TO 165 INDC2850 162 WRITE (LOUT,99953) INDC2855 99953 FORMAT (///42H ********PROGRAM READ AN END-OF-FILE WHILE, INDC2860 * 40H ATTEMPTING TO READ USER-SUPPLIED LABELS/9X, 12HWILL INSTEAD, INDC2865 * 44H USE INTEGER LABELING. EXECUTION CONTINUES.//) INDC2870 C INDC2875 C INDC2880 C THE DATA WILL SOON BE IN THE LOWER TRIANGULAR HALF OF THE ARRAY INDC2885 C DATA(M,I,J). THE UPPER HALF WILL BE USED FOR RESIDUALS. INDC2890 C INDC2895 C NOW TRANSFORM THEM TO BE SIMILARITIES (NOT DISSIMILARITIES) INDC2900 C ON THE INTERVAL [0,1]. INDC2905 C INDC2910 C INDC2915 C INDC2920 165 IF (MCON.EQ.0) GO TO 170 INDC2925 C INDC2930 C NORMALIZE MATRIX CONDITIONAL DATA. INDC2935 C INDC2940 CALL NRMLZE INDC2945 GO TO 205 INDC2950 C INDC2955 C INDC2960 C THE FOLLOWING NORMALIZATION IS APPROPRIATE ONLY FOR MATRIX INDC2965 C UNCONDITIONAL DATA. INDC2970 C INDC2975 C INDC2980 170 IF (ISWIT.GT.0) GO TO 175 INDC2985 BMN = AMN INDC2990 C3 = -1.0 INDC2995 GO TO 180 INDC3000 175 BMN = AMX INDC3005 C3 = 1.0 INDC3010 180 BMX = AMX - AMN INDC3015 IF (BMX.LT.1.0E-6) GO TO 650 INDC3020 C1 = 1.0/BMX INDC3025 C2 = BMN*C1 INDC3030 DO 195 M=1,M3 INDC3035 DO 190 I=2,N INDC3040 I1 = I - 1 INDC3045 DO 185 J=1,I1 INDC3050 DATA(M,I,J) = C3*C1*(BMN-DATA(M,J,I)) INDC3055 185 CONTINUE INDC3060 190 CONTINUE INDC3065 195 CONTINUE INDC3070 IF (ISWIT.LE.0) GO TO 200 INDC3075 WRITE (LOUT,99952) C2, C1 INDC3080 99952 FORMAT (/42H THE INPUT DATA ARE DISSIMILARITIES WHICH , 7HHAVE BE,INDC3085 * 35HEN LINEARLY TRANSFORMED AS FOLLOWS://40X, 10HY =, INDC3090 * F10.5, 3H - , F10.5, 13H(Y )/41X, 7HNEW SIM, 27X, INDC3095 * 10HRAW DISSIM) INDC3100 GO TO 205 INDC3105 200 WRITE (LOUT,99951) C1, C2 INDC3110 99951 FORMAT (/39H THE INPUT DATA ARE SIMILARITIES WHICH , 9HHAVE BEEN, INDC3115 * 33H LINEARLY TRANSFORMED AS FOLLOWS://40X, 10HY =, F10.5, INDC3120 * 10H(Y ), 3H - , F10.5/41X, 7HNEW SIM, 14X, 7HRAW SIM) INDC3125 205 IF (IPRN.NE.1) GO TO 225 INDC3130 WRITE (LOUT,99950) INDC3135 99950 FORMAT (44H1 NO TRANSFORMED RAW SUBSCRIPTS, 6H ,INDC3140 * 9H LABELS/49H VALUES DATA OF STIMULI ,INDC3145 * 12H OF STIMULI/44H I J//INDC3150 * ) INDC3155 K = 0 INDC3160 K1 = 0 INDC3165 DO 220 M=1,M3 INDC3170 DO 215 I=2,N INDC3175 I1 = I - 1 INDC3180 DO 210 J=1,I1 INDC3185 K = K + 1 INDC3190 K1 = K1 + 1 INDC3195 IF (MOD(K1,50).EQ.0) WRITE (LOUT,99950) INDC3200 WRITE (LOUT,99949) K, DATA(M,I,J), DATA(M,J,I), I, J, INDC3205 * (ILAB(I,L),L=1,6), (ILAB(J,L),L=1,6), (MLAB(M,L),L=1,6) INDC3210 210 CONTINUE INDC3215 215 CONTINUE INDC3220 99949 FORMAT (1X, I5, F12.4, F13.4, 1X, 2I6, 4X, 6A1, 4X, 6A1, 10X, INDC3225 * 14H DATA SOURCE =, 2X, 6A1) INDC3230 K1 = K1 + 4 INDC3235 WRITE (LOUT,99935) INDC3240 220 CONTINUE INDC3245 C INDC3250 225 DO 235 I=1,N INDC3255 DO 230 IW=1,LD1 INDC3260 P(I,IW) = 0.0 INDC3265 230 CONTINUE INDC3270 P(I,LD2) = 1.0 INDC3275 235 CONTINUE INDC3280 VARSIM = 0.0 INDC3285 IF (MCON.EQ.1) GO TO 270 INDC3288 DO 250 M=1,M3 INDC3290 DM = 0.0 INDC3295 DO 245 I=2,N INDC3300 I1 = I - 1 INDC3305 DO 240 J=1,I1 INDC3310 DIJ = DATA(M,I,J) INDC3315 DM = DM + DIJ INDC3320 240 CONTINUE INDC3325 245 CONTINUE INDC3330 DMEAN(M) = DM*ABNV INDC3335 250 CONTINUE INDC3340 C INDC3350 C COMPUTE MEAN AND S.D.'S FOR MATRIX UNCONDITIONAL DATA. INDC3355 C INDC3360 DO 265 M=1,M3 INDC3365 BMX = 0.0 INDC3370 AMX = 0.0 INDC3375 SDDAT(M) = 0.0 INDC3380 DO 260 I=2,N INDC3385 I1 = I - 1 INDC3390 DO 255 J=1,I1 INDC3395 AMX = DATA(M,I,J) INDC3400 BMX = BMX + AMX INDC3405 SDDAT(M) = SDDAT(M) + AMX*AMX INDC3410 255 CONTINUE INDC3415 260 CONTINUE INDC3420 SDDAT(M) = SQRT((SDDAT(M)-BMX*BMX*ABNV)*ABNV) INDC3425 265 CONTINUE INDC3430 270 DO 285 M=1,M3 INDC3435 DM = 0.0 INDC3436 IF (MCON.EQ.0) DM = DMEAN(M) INDC3438 DO 280 I=2,N INDC3440 I1 = I - 1 INDC3445 DO 275 J=1,I1 INDC3450 DATA(M,I,J) = DATA(M,I,J) - DM INDC3455 VARSIM = VARSIM + DATA(M,I,J)**2 INDC3460 275 CONTINUE INDC3465 280 CONTINUE INDC3470 285 CONTINUE INDC3475 DO 295 M=1,M3 INDC3480 DO 290 IW=1,LD2 INDC3485 RISE(M,IW) = 0.0 INDC3490 290 CONTINUE INDC3495 295 CONTINUE INDC3500 IF (ICON.GT.0) CALL INICON(0, 0, ICON) INDC3505 IF (ICON.LE.0 .AND. M15.GT.0) GO TO 310 INDC3510 IF (ICON.GT.0 .OR. M15.NE.0) GO TO 300 INDC3515 WRITE (LOUT,99948) ICON INDC3520 99948 FORMAT (//45H YOU HAVE SPECIFIED THE OPTION FOR REGRESSION, INDC3525 * 27H ONLY (M15 = 0 ITERATIONS),/30H BUT YOU HAVE DECLARED AN INFE,INDC3530 * 23HASIBLE LOGICAL DEVICE #, 1H(, I3, 23H) AS THE VALUE OF ICON INDC3535 * /46H ON THE PARAMETER CARD. EXECUTION TERMINATES.//) INDC3540 STOP INDC3545 C INDC3550 C THIS IS THE OPTION FOR CALLING THE REGRESSION ROUTINE FOR A INDC3555 C USER-SUPPLIED INIT. CONF., WHILE BYPASSING THE ITERATIVE INDC3560 C FITTING PROCEDURE. INDC3565 C INDC3570 300 WRITE (LOUT,99937) TITL INDC3575 NST = 1 INDC3580 DENS(1) = 0.0 INDC3585 KC = 1 INDC3590 CALL REGR(0, 0, 0) INDC3595 C INDC3600 C THE FIRST ARGUMENT IN THE FOLLOWING CALL TO PRINTD IS 1 SO THAT INDC3605 C DENSITY CAN BE COMPUTED AND PRINTED. INDC3610 C INDC3615 WRITE (LOUT,99935) INDC3620 DO 305 M=1,M3 INDC3625 CALL PRINTD(1, 0, M) INDC3630 305 CONTINUE INDC3635 C INDC3640 C POSSIBLE BRANCH TO COMBINATORIAL OPTIMIZATION INDC3645 C INDC3650 IF (M15) 575, 635, 310 INDC3655 310 DLAM = ALF0 INDC3660 DO 315 IW=1,LD1 INDC3665 ALFRAY(IW) = ALPHAI INDC3670 BETRAY(IW) = BETAI INDC3675 315 CONTINUE INDC3680 C INDC3685 C INDC3690 C INDC3695 C INDC3700 C DO LOOP FOR 'MAJOR' ITERATIONS INDC3705 C INDC3710 DO 565 MST=1,M15 INDC3715 C INDC3720 C INDC3725 M24 = MST INDC3730 BLMX(MST) = -1.0E20 INDC3735 BLMN(MST) = 1.0E20 INDC3740 IADJ(MST) = 0 INDC3745 DENS(MST) = 0. INDC3750 MP(MST) = 0 INDC3755 NW(MST) = 0 INDC3760 SVAF(MST) = 0.0 INDC3765 MST1 = MST - 1 INDC3770 KMOD = MOD(MST,2) INDC3775 IF (IPRE.NE.1) GO TO 320 INDC3780 C INDC3785 C MECHANISM TO ALLOW A RESTART WITH DE NOVO ITERATIONS: INDC3790 C INDC3795 I = 0 INDC3800 IF (ITPRE.LE.0) I = 1 INDC3805 ITPRE1 = MST - I INDC3810 ITPRE = ITPRE1 - 1 INDC3815 C INDC3820 C NO MORE PRE-ITERATIONS INDC3825 C INDC3830 IF (IPRE.EQ.1 .AND. KTPRE.GE.0 .AND. MST1.GT.0) WRITE INDC3835 * (LOUT,99947) BLMX(MST1), BCRIT INDC3840 99947 FORMAT (///40H SINCE THE MAXIMUM BLOS ACROSS SUBSETS (, F10.8, INDC3845 * 53H) IS LESS THAN OR EQUAL TO THE STOPPING CRITERION OF , INDC3850 * 8HBLCRIT =, F12.8/14H THERE WILL BE, 24H NO MORE PRE-POLISHING I,INDC3855 * 10HTERATIONS.///) INDC3860 C INDC3865 C TAKING CARE OF NUMBER OF INNER ITERATIONS CONSTITUTING AN OUTER INDC3870 C ITERATION, DEPENDING ON THE NATURE OF CURRENT MAJOR ITERATION. INDC3875 C INDC3880 320 IF (MST-ITPRE1) 325, 330, 335 INDC3885 325 L15 = IUPPER INDC3890 K15 = IUPPER INDC3895 GO TO 345 INDC3900 330 L15 = JPOST + 30 INDC3905 WRITE (LOUT,99946) MST INDC3910 99946 FORMAT (/49H 30 EXTRA INNER ITERATIONS FOR THE FIRST MAJOR, INDC3915 * 12H ITERATION (, I3, 29H) WITH POLISHING (BETA = 1.0)//) INDC3920 GO TO 340 INDC3925 335 L15 = JPOST INDC3930 340 K15 = JPOST INDC3935 C INDC3940 C NO EFFECTIVE LIMIT ON ADJUST CALLS FOR POST-POLISHING ITERATIONS INDC3945 JLIM = 1000 INDC3950 C INDC3955 C INDC3960 C INDC3965 C THE FOLLOWING DO LOOP SETS UP THE OUTER ITERATIONS. INDC3970 C INDC3975 C INDC3980 345 DO 545 IIW=1,LD1 INDC3985 C INDC3990 C KC DETERMINES PRINTING AT THE END OF A MAJOR ITERATION INDC3995 C INDC4000 KC = 0 INDC4005 IF (IIW.EQ.LD1) KC = 1 INDC4010 C INDC4015 C PROCEDURE FOR REVERSING THE ORDER IN WHICH SUBSETS INDC4020 C (AND THEIR WEIGHTS) ARE COMPUTED, FOR EVEN-NUMBERED INDC4025 C MAJOR ITERATIONS. (ALSO USED BELOW IN COMBINATORIAL INDC4030 C OPTIMIZATION ITERATIONS). INDC4035 C INDC4040 IF (KMOD.EQ.1) GO TO 350 INDC4045 IW = LD2 - IIW INDC4050 GO TO 355 INDC4055 350 IW = IIW INDC4060 C INDC4065 C USER-SUPPLIED TRIAL STEP SIZE FOR FIRST INNER ITERATION INDC4070 355 ALF0 = DLAM INDC4075 IFLOP = 2 INDC4080 OBSLF(2) = 1.0E20 INDC4085 ASTEP(2) = 1.0E20 INDC4090 RCORR = 0.0 INDC4095 C INDC4100 C GET THE RESIDUALS AND CENTER THEM AT THE BEGINNING OF OUTER IT. INDC4105 C INDC4110 CALL RESID(IW) INDC4115 JADJ = 0 INDC4120 C INDC4125 C INDC4130 C INDC4135 C NOW SET UP THE INNER ITERATIONS. INDC4140 C INDC4145 C INDC4150 C INDC4155 DO 480 KL1=1,L15 INDC4160 KL = KL1 INDC4165 IFLIP = IFLOP INDC4170 IFLOP = 3 - IFLOP INDC4175 GRSUM = 0.0 INDC4180 GRMEN(IFLOP) = 0.0 INDC4185 C INDC4190 C ON THE VERY FIRST POLISHING ITERATION, CALL ADJUST INDC4195 C ON EVERY THIRD INNER ITERATION DURING THE EXTRA INDC4200 C INNER ITERATIONS INDC4205 C INDC4210 IF (KL.GT.K15 .AND. MOD(KL,3).EQ.0 .AND. ALFRAY(IW).GT.1.0E-6) INDC4215 * CALL ADJUST(MST, IW) INDC4220 C INDC4225 C NOW GET THE INITIAL CONFIGURATION FOR (ONLY) THE VERY FIRST INDC4230 C INNER ITERATION OF SUBSET IW (DURING THE FIRST MAJOR ITERATION). INDC4235 C INDC4240 IF (KL+MST.EQ.2 .AND. ICON.LE.0 .AND. ITPRE.GT.0) CALL INDC4245 * INICON(MST, IW, ICON) INDC4250 IF (MST.LE.ITPRE .OR. KL.NE.1) GO TO 360 INDC4255 C INDC4260 C INDC4265 C THE 'DE NOVO' APPROACH TO COMPUTING P(I) AFTER INITIAL POLISHING INDC4270 C INDC4275 CALL INICON(MST, IW, 0) INDC4280 ALFRAY(IW) = ALPHAI INDC4285 BETRAY(IW) = BETAI INDC4290 C INDC4295 C INDC4300 C INDC4305 360 BCONST = BSUM(IW) INDC4310 C INDC4315 C INDC4320 C NOW ESTIMATE THE VALUE OF RISE(M,IW) AND THE INDC4325 C REGRESSION CONSTANT, AND THEN TRANSLATE THE RESIDUALS. INDC4330 C INDC4335 CALL COMPW(IW, 0) INDC4340 C INDC4345 C INDC4350 C COMPUTE THE A PART OF THE LOSS FUNCTION (USED IN ESTIMATION INDC4355 C OF OPTIMAL LAMBDA) AND THE B PART (EARLIER) VIA FUNCTION CALLS. INDC4360 C (THE LATTER GIVES COEFFICIENT C FOR THE QUADRATIC INTERPOLATION INDC4365 C FUNCTION.) INDC4370 C INDC4375 IF (KL.GT.1) GO TO 365 INDC4380 CONST = ALFRAY(IW)*ASUM(IW) + BETRAY(IW)*BCONST INDC4385 GO TO 370 INDC4390 365 CONST = OBSLF(IFLIP) INDC4395 C INDC4400 C THEN EVALUATE THE PART OF THE EXPRESSIONS FOR EACH I = 1, N INDC4405 C INDC4410 370 CALL CDAPI(IW) INDC4415 C INDC4420 C INDC4425 C NOTE THAT DB IS AVAILABLE AS A FUNCTION CALL INDC4430 C INDC4435 DO 375 KLM=1,N INDC4440 GL(KLM,IFLOP) = BETRAY(IW)*DB(KLM,IW) + ALFRAY(IW)*DAPI(KLM) INDC4445 GLA = GL(KLM,IFLOP) INDC4450 IF (ABS(GLA).LE.1.0E-12) GO TO 375 INDC4455 GRMEN(IFLOP) = GRMEN(IFLOP) + GLA INDC4460 GRSUM = GRSUM + GLA**2 INDC4465 375 CONTINUE INDC4470 GRLEN(IFLOP) = GRSUM INDC4475 GRSUM = SQRT(GRSUM) INDC4480 IF (GRSUM.LE.1.0E-12) GO TO 485 INDC4485 IF (KL.LE.1) GO TO 385 INDC4490 RCORR = 0.0 INDC4495 C = 0. INDC4500 D = 0. INDC4505 DO 380 I=1,N INDC4510 C = C + GL(I,1) INDC4515 D = D + GL(I,2) INDC4520 RCORR = RCORR + GL(I,1)*GL(I,2) INDC4525 380 CONTINUE INDC4530 RCORR = (BN*RCORR-C*D)/(SQRT((BN*GRLEN(1)-GRMEN(1)**2)*(BN* INDC4535 * GRLEN(2)-GRMEN(2)**2))) INDC4540 385 CLAM = ALF0/GRSUM INDC4545 ILOOP = 0 INDC4550 DO 390 KLM=1,N INDC4555 P(KLM,IW) = P(KLM,IW) - CLAM*GL(KLM,IFLOP) INDC4560 390 CONTINUE INDC4565 IF (KP(MST).EQ.1 .AND. KL.EQ.1) WRITE (LOUT,99945) MST, KL, INDC4570 * ALFRAY(IW), BETRAY(IW), CLAM, ACRIT, GRSUM, RCORR, IW, INDC4575 * (RISE(M,IW),M=1,M3) INDC4580 99945 FORMAT (6H MAJOR, I3, 7H, INNER, I4, 7H ALPHA=, F6.4, 6H BETA=, INDC4585 * F6.4, 8H LAMBDA=, F11.4, 7H ACRIT=, F9.7, 7H GRSUM=, F10.5, INDC4590 * 5H COR=, F8.4, 11H; WTS. (IW=, I3, 9H) FOLLOW: /64(1X, 16F8.4/))INDC4595 C INDC4600 C EVALUATE THE LOSS FUNCTION NOW THAT A STEP OF LENGTH ALF0 HAS INDC4605 C BEEN TAKEN (USED FOR OPTIMAL LAMBDA). INDC4610 C INDC4615 FALPHA = ALFRAY(IW)*ASUM(IW) + BETRAY(IW)*BSUM(IW) INDC4620 C INDC4625 C THE FOLLOWING TWO STATEMENTS GIVE RISE TO BENIGN UNDERFLOW AND INDC4630 C EXP DIAGNOSTICS INDC4635 C INDC4640 395 IF (ABS(ALF0).LE.1.0E-12) GO TO 400 INDC4645 AA = (FALPHA-CONST+ALF0*GRSUM)/ALF0**2 INDC4650 ALFMIN = GRSUM*0.5/AA INDC4655 IF (AA.GT.0.001) GO TO 410 INDC4660 400 IF (ABS(RCORR).GT.1.0E-12) GO TO 405 INDC4665 ALFMIN = ALF0 INDC4670 GO TO 410 INDC4675 405 ALFMIN = ALF0*2.0**RCORR INDC4680 410 ASTEP(IFLOP) = (ALFMIN-ALF0)/GRSUM INDC4685 IF (AA.GT.0.001) GO TO 415 INDC4690 IF (KP(MST).EQ.1) WRITE (LOUT,99944) AA, ALF0 INDC4695 99944 FORMAT (27H QUADRATIC COEFFICIENT A =, F16.5, 15H; SUBSTITUTED V,INDC4700 * 14HALUE OF ALF0 =, F11.5) INDC4705 IF (KL.EQ.1) GO TO 425 INDC4710 IF (KP(MST).EQ.1) WRITE (LOUT,99943) ALF0, ASTEP(IFLOP), ALFMIN INDC4715 99943 FORMAT (34H USED COS MULTIPLIER OF OLD ALF0 =, F10.5, 9H TO GET N,INDC4720 * 17HORMALIZED ASTEP =, F10.5, 4X, 13H NEW ALFMIN =, F10.5) INDC4725 C INDC4730 C IN CASE OVERFLOW CAUSES HAVOC ... INDC4735 C INDC4740 415 IF (ASTEP(IFLOP).GT.1.0E10) GO TO 485 INDC4745 C INDC4750 C INDC4755 C NOW USE INTERMEDIATE STEP SIZE THAT SHOULD MINIMIZE THE LOSS INDC4760 C FUNCTION. INDC4765 C INDC4770 DO 420 KLM=1,N INDC4775 P(KLM,IW) = P(KLM,IW) - ASTEP(IFLOP)*GL(KLM,IFLOP) INDC4780 420 CONTINUE INDC4785 C INDC4790 C COMPUTE THE OBSERVED VALUE OF THE LOSS FUNCTION INDC4795 C INDC4800 425 ALOS = ASUM(IW) INDC4805 BCONST = BSUM(IW) INDC4810 OBSLF(IFLOP) = ALFRAY(IW)*ALOS + BETRAY(IW)*BCONST INDC4815 C INDC4820 C GET THE ESTIMATED VALUE OF THE LOSS FUNCTION AT ALPHA MIN INDC4825 C INDC4830 PRELF = ALFMIN*(AA*ALFMIN-GRSUM) + CONST INDC4835 C INDC4840 C STRATEGY WHEN OBSLF IS GETTING WORSE: INDC4845 C INDC4850 IF (OBSLF(IFLOP).LE.OBSLF(IFLIP)) GO TO 430 INDC4855 C CHECK TO SEE THAT PROCEDURE IS NOT CAUGHT IN A LOOP INDC4860 ILOOP = ILOOP + 1 INDC4865 IF (ILOOP.GT.3) GO TO 430 INDC4870 ALF0 = ALFMIN INDC4875 FALPHA = OBSLF(IFLOP) INDC4880 GO TO 395 INDC4885 C INDC4890 C INDC4895 C GET LENGTH OF 'RELATIVE GRADIENT', WHICH IS THE PRODUCT OF INDC4900 C LENGTH OF GRADIENT TIMES LENGTH OF P VECTOR. INDC4905 C INDC4910 430 RELGR = 0.0 INDC4915 DO 435 KLM=1,N INDC4920 RELGR = RELGR + P(KLM,IW)**2 INDC4925 435 CONTINUE INDC4930 RELGR = GRSUM*SQRT(RELGR) INDC4935 C INDC4940 IF (KP(MST).EQ.1 .AND. KL.EQ.1) WRITE (LOUT,99942) KL, INDC4945 * (P(KLM,IW),KLM=1,N) INDC4950 99942 FORMAT (21H P VALUES FOLLOW (KL=, I3, 2H):/64(1X, 16F8.4/)) INDC4955 IF (KP(MST).EQ.1 .AND. KL.EQ.1) WRITE (LOUT,99941) PRELF, INDC4960 * OBSLF(IFLOP), ASTEP(IFLOP), RELGR, DCON, ALOS, BCONST INDC4965 99941 FORMAT (14H F(ALPHA MIN)=, F11.5, 7H OBSLF=, F11.5, 7H ASTEP=, INDC4970 * F12.5, 8H RELGR =, F11.4, 7H DCON =, F11.4, 7H ALOS =, F10.4, INDC4975 * 6H BLOS=, F10.6) INDC4980 C INDC4985 C PUT IN NEW VALUE FOR LIMIT OF OPTIMAL STEP SIZE FOR INDC4990 C NEXT INNER ITERATION. INDC4995 C INDC5000 ALF0 = ALFMIN INDC5005 C INDC5010 C CHECK FOR CONVERGENCE, UNLESS THIS IS THE LAST INNER ITERATION. INDC5015 C INDC5020 C IF THIS IS THE LAST (L15-TH) INNER ITERATION, NO TESTS ARE INDC5025 C NECESSARY. INDC5030 C INDC5035 C INDC5040 C INDC5045 IF (KL.GT.K15) GO TO 455 INDC5050 C INDC5055 C INDC5060 C IF THE 'RELATIVE' GRADIENT IS NOT SHORT ENOUGH, DO NOT MAKE ANY INDC5065 C CHANGES. INDC5070 C INDC5075 IF (RELGR.GT.ACRIT) GO TO 440 INDC5080 C INDC5085 C IF BLOS IS SMALL ENOUGH, AND THIS IS A PRE-ITERATION, INDC5090 C TERMINATE THIS OUTER ITERATION. INDC5095 C INDC5100 IF (BCONST.GT.BLCRIT) GO TO 450 INDC5105 C INDC5110 IF (MST.GT.ITPRE) GO TO 445 INDC5115 IF (KP(MST).EQ.1) WRITE (LOUT,99940) KL INDC5120 99940 FORMAT (46H RELGR AND BCONST CRITERIA ARE BOTH SATISFIED., INDC5125 * 46H THIS OUTER ITERATION OF A MAJOR PRE-ITERATION, 10H TERMINATE,INDC5130 * 7HS AFTER, I3, 11H INNER ITS.) INDC5135 GO TO 485 INDC5140 C INDC5145 440 IF (ABS(ASTEP(1)).LE.1.0E-5 .AND. ABS(ASTEP(2)).LE.1.0E-5) GO TO INDC5150 * 450 INDC5155 GO TO 460 INDC5160 C INDC5165 C INDC5170 C SINCE THIS IS NOT A PRE-ITERATION, POSSIBLY CALL ADJUST AND INDC5175 C CONTINUE THE INNER ITERATIONS. INDC5180 C INDC5185 445 IF (KP(MST).EQ.1 .AND. ALFRAY(IW).GT.1.0E-6) WRITE (LOUT,99939) INDC5190 * RELGR, ACRIT, ALFRAY(IW), BETRAY(IW), BCONST, BLCRIT INDC5195 99939 FORMAT (7H RELGR=, F10.8, 10H.LE.ACRIT=, F10.8, 2X, 10H CALLED AD,INDC5200 * 12HJUST: ALPHA=, F9.6, 6H BETA=, F9.6, 8H BCONST=, F10.8, INDC5205 * 8H BLCRIT=, F10.8) INDC5210 450 IF (JADJ.EQ.JLIM) GO TO 485 INDC5215 IF (ALFRAY(IW).GT.1.0E-6) CALL ADJUST(MST, IW) INDC5220 C INDC5225 C INDC5230 C RETURN TO THE ORIGINAL (USER-SUPPLIED) ALPHA STEP SIZE AFTER INDC5235 C CONVERGENCE HAS OCCURRED. INDC5240 C INDC5245 455 ALF0 = DLAM INDC5250 C INDC5255 C RE-CENTER THE RESIDUALS. INDC5260 C INDC5265 460 DO 475 M=1,M3 INDC5270 DO 470 I=2,N INDC5275 I1 = I - 1 INDC5280 DO 465 J=1,I1 INDC5285 DATA(M,J,I) = DATA(M,J,I) + WCON(M) INDC5290 465 CONTINUE INDC5295 470 CONTINUE INDC5300 475 CONTINUE INDC5305 480 CONTINUE INDC5310 C INDC5315 485 IF (KP(MST).EQ.0) GO TO 490 INDC5320 WRITE (LOUT,99945) MST, KL, ALFRAY(IW), BETRAY(IW), CLAM, ACRIT, INDC5325 * GRSUM, RCORR, IW, (RISE(M,IW),M=1,M3) INDC5330 WRITE (LOUT,99941) PRELF, OBSLF(IFLOP), ASTEP(IFLOP), RELGR, INDC5335 * DCON, ALOS, BCONST INDC5340 C INDC5345 C GET THE VAF BEFORE AS WELL AS AFTER POLISHING AT END OF OUTER INDC5350 C ITERATION. INDC5355 C INDC5360 490 IF (KP(MST).EQ.1) CALL VARIAN(MST, 1, IW) INDC5365 IF (BCONST.GT.BLMX(MST)) BLMX(MST) = BCONST INDC5370 IF (BCONST.LT.BLMN(MST)) BLMN(MST) = BCONST INDC5375 C INDC5380 C INDC5385 C SINCE BLOS IS USING PAIRWISE PRODUCTS OF THE P(I), THE SIGNS INDC5390 C ARE POTENTIALLY REVERSED. SET THE SIGN OF THE LARGEST INDC5395 C OF THE ABS(P(I)) TO BE POSITIVE, AND MULTIPLY THE REST BY INDC5400 C THE SAME SIGN USED TO MAKE THE LARGEST POSITIVE. INDC5405 C INDC5410 I3 = 1 INDC5415 AM = ABS(P(1,IW)) INDC5420 DO 495 I=2,N INDC5425 AN = ABS(P(I,IW)) INDC5430 IF (AN.LE.AM) GO TO 495 INDC5435 AM = AN INDC5440 I3 = I INDC5445 495 CONTINUE INDC5450 IF (P(I3,IW).GE.0.0) GO TO 510 INDC5455 DO 500 I=1,N INDC5460 P(I,IW) = -P(I,IW) INDC5465 500 CONTINUE INDC5470 B3 = 0.0 INDC5475 DO 505 KLM=1,N INDC5480 B3 = B3 + P(KLM,IW) INDC5485 505 CONTINUE INDC5490 B3 = B3*BIN INDC5495 B4 = TM/B3 INDC5500 510 IF (MST.LE.ITPRE) GO TO 515 INDC5505 IF (KP(MST).EQ.1) WRITE (LOUT,99938) IW, B3, TM, B4, INDC5510 * (P(KLM,IW),KLM=1,N) INDC5515 99938 FORMAT (17H PRE-POLISHED P (, I2, 15H): (MEAN P(J) =, F12.6, 1H), INDC5520 * 17H MEAN(P(I)P(J)) =, F12.6, 8H RATIO =, F12.6, 11H; USED .707, INDC5525 * 17H P VALUES FOLLOW:/64(1X, 16F8.4/)) INDC5530 C INDC5535 C FORCE THE T VECTOR TO BE BINARY. INDC5540 C INDC5545 CALL POLISH(MST, IW) INDC5550 C INDC5555 C IF SOME OF THE FOLLOWING WRITE STATEMENTS WERE DELETED, INDC5560 C THE FOLLOWING CALL TO COMPW COULD BE ELIMINATED, SINCE INDC5565 C THE SUBSEQUENT CALL TO REGR TAKES CARE OF ESTIMATING ALL INDC5570 C THE WEIGHTS SIMULTANEOUSLY. INDC5575 C INDC5580 CALL COMPW(IW, 1) INDC5585 C INDC5590 C INDC5595 C INDC5600 515 IF (MST+IW.GT.2) CALL REGR(MST, IW, MST) INDC5605 IF (MST.GT.ITPRE) GO TO 520 INDC5610 IF (KP(MST).NE.1) GO TO 545 INDC5615 IF (KC.EQ.1) WRITE (LOUT,99937) TITL INDC5620 99937 FORMAT (/1X, 72A1) INDC5625 WRITE (LOUT,99936) MST, IW, B3, TM, B4, (P(KLM,IW),KLM=1,N) INDC5630 99936 FORMAT (/5H IT #, I3, 12H P MATRIX (, I2, 2H):, 6X, 9H( MEAN P(, INDC5635 * 4HJ) =, F12.6, 1H), 17H MEAN(P(I)P(J)) =, F12.6, 8H RATIO =, INDC5640 * F12.6/64(1X, 16F8.4/)) INDC5645 GO TO 545 INDC5650 520 IF (KC.NE.1) GO TO 525 INDC5655 99935 FORMAT (////) INDC5660 GO TO 530 INDC5665 525 IF (KP(MST).NE.1) GO TO 540 INDC5670 530 DENS(MST) = 0.0 INDC5675 DO 535 M=1,M3 INDC5680 CALL PRINTD(MST, IW, M) INDC5685 535 CONTINUE INDC5690 540 IF (IPUN.GE.1 .AND. KC.EQ.1) CALL PUN(MST) INDC5695 545 CONTINUE INDC5700 DO 555 M=1,M3 INDC5705 DO 550 IW=1,LD1 INDC5710 IF (RISE(M,IW).LT.0.0 .AND. NNEG.EQ.0) NW(MST) = NW(MST) + 1 INDC5715 IF (RISE(M,IW).LE.0.0 .AND. NNEG.GT.0) NW(MST) = NW(MST) + 1 INDC5720 550 CONTINUE INDC5725 555 CONTINUE INDC5730 WRITE (LOUT,99934) MST, IADJ(MST), NW(MST) INDC5735 99934 FORMAT (/24H DURING MAJOR ITERATION , I4, 19H ADJUST WAS CALLED , INDC5740 * I4, 22H TIMES AND THERE WERE , I4, 18H NEGATIVE WEIGHTS./////) INDC5745 IF (IERR.NE.4) SVAF(MST) = VAF INDC5750 C INDC5755 C OPTION TO TERMINATE PRE-ITERATIONS, BASED ON VAF FROM INDC5760 C ITERATIVE REGRESSION. INDC5765 C INDC5770 IF (BLMX(MST).LE.BCRIT .AND. MST.LE.ITPRE .AND. IPRE.NE.1) IPRE = INDC5775 * 1 INDC5780 C INDC5785 IF (MST.LT.ITPRE1 .OR. MST1.LE.1) GO TO 565 INDC5790 C INDC5795 C CHECK FOR OSCILLATION ON EVERY OTHER MAJOR ITERATION. INDC5800 C (NOTE THAT MST2 WILL BE POSITIVE BECAUSE OF THE PRECEDING INDC5805 C CHECK TO SEE THAT MST1 IS .GT. 1. ALSO, CHECKING FOR INDC5810 C DENSITY AS WELL AS VAF SHOULD MAKE IT UNLIKELY THAT THIS INDC5815 C TEST WOULD APPLY DURING PRE-ITERATIONS WHEN DENS(MST) = 0.0 INDC5820 C INDC5825 MST2 = MST1 - 1 INDC5830 C INDC5835 C NOTE THAT IF K, INSTEAD OF 1, WERE SUBTRACTED IN THE INDC5840 C PRECEDING STATEMENT, WE COULD TEST FOR A MORE GENERAL INDC5845 C PERIOD OF OSCILLATION. INDC5850 C INDC5855 IF (SVAF(MST)-SVAF(MST2).LT.0.001 .AND. DENS(MST)-DENS(MST2) INDC5860 * .LT.0.001 .AND. NW(MST).LE.NW(MST1)) GO TO 560 INDC5865 IF (SVAF(MST)-SVAF(MST1).GE.0.001 .OR. SVAF(MST)-SVAF(MST1) INDC5870 * .LT.0.0 .OR. ABS(DENS(MST)-DENS(MST1)).GE.0.005) GO TO 565 INDC5875 560 WRITE (LOUT,99933) INDC5880 99933 FORMAT (///43H TERMINATED THIS SERIES OF MAJOR ITERATIONS, INDC5885 * 34H BECAUSE OF NEGLIGIBLE IMPROVEMENT///) INDC5890 GO TO 570 INDC5895 565 CONTINUE INDC5900 C INDC5905 C INDC5910 C INDC5915 C INDC5920 570 IF (M24.EQ.M15) WRITE (LOUT,99932) M24 INDC5925 99932 FORMAT (//49H TERMINATED MAJOR ITERATIONS AFTER REACHING USER-, INDC5930 * 19HSPECIFIED LIMIT OF , I3//) INDC5935 WRITE (LOUT,99931) INDC5940 99931 FORMAT (1H1) INDC5945 WRITE (LOUT,99930) INDC5950 99930 FORMAT (//50H SUMMARY OF ALTERNATING LEAST SQUARES COMPUTATION , INDC5955 * 35HPRIOR TO COMBINATORIAL OPTIMIZATION//) INDC5960 WRITE (LOUT,99937) TITL INDC5965 M15 = MIN0(M24,M15) INDC5970 WRITE (LOUT,99929) INDC5975 99929 FORMAT (//20H VAF BY ITERATION, 19X, 13H DENSITY OF, 13X, INDC5980 * 57HB-PART OF LOSS FUNCTION NUMBER OF PREMATURE/43X,INDC5985 * 8H SUBSETS, 41H AT END OF OUTER ITS. , INDC5990 * 30H ADJUSTMENTS POLISHING/41X, 13H(=0. PRIOR TO, 58X, INDC5995 * 12H(=0 PRIOR TO/42X, 10HPOLISHING), 61X, 10HPOLISHING)) INDC6000 IN = 2 INDC6005 IF (NNEG.LE.0) IN = 1 INDC6010 WRITE (LOUT,99928) (MST,SVAF(MST),NW(MST),(ING(IN,I),I = 1, 4), INDC6015 * DENS(MST),BLMX(MST), INDC6020 * BLMN(MST),IADJ(MST),MP(MST), MST = 1, M15) INDC6025 99928 FORMAT (//(2X, 4HVAF(, I3, 2H)=, F12.6, 2X, I3, 1X, 4A1, 6H WTS ,INDC6030 * 6HDENS =, F10.5, 12H MAX BLOS =, F10.5, 12H MIN BLOS =, F10.7, INDC6035 * 6H ADJ =, I4, 11H MID POL =, I5/)) INDC6040 C INDC6045 IF (IPRE.EQ.1 .OR. M24.GE.ITPRE1) GO TO 575 INDC6050 WRITE (LOUT,99927) INDC6055 99927 FORMAT (//47H SINCE THE SUBSETS HAVE NOT YET BEEN POLISHED, , INDC6060 * 32HTHE P MATRIX IS NOT BINARY, AND /18H THERE WILL BE NO , INDC6065 * 28H COMBINATORIAL OPTIMIZATION./////) INDC6070 GO TO 5 INDC6075 C INDC6080 C TIME FOR COMBINATORIAL OPTIMIZATION INDC6085 C INDC6090 575 IF (IERR.EQ.4) GO TO 580 INDC6095 IF (M15.LE.0) GO TO 585 INDC6100 IF (SVAF(M15).GT.0.0 .AND. SVAF(M15).LE.1.0) GO TO 585 INDC6105 580 WRITE (LOUT,99926) INDC6110 99926 FORMAT (///49H ******** SINCE LINEAR DEPENDENCY WAS PRESENT IN, INDC6115 * 59H THE SOLUTION FROM THE ALTERNATING LEAST SQUARES PROCEDURE,/ INDC6120 * 12X, 48HCOMBINATORIAL OPTIMIZATION WILL NOT BE EXECUTED.///) INDC6125 GO TO 5 INDC6130 585 KC = 1 INDC6135 IXV = 1 INDC6140 IC1 = 0 INDC6145 DO 620 MST=1,15 INDC6150 NST = MST INDC6155 KMOD = MOD(NST,2) INDC6160 IC0 = 0 INDC6165 WRITE (LOUT,99925) NST INDC6170 99925 FORMAT (///22H THIS IS ITERATION, I3, 18H OF COMBINATORIAL , INDC6175 * 36HOPTIMIZATION, PHASE ONE (DOUBLETONS)///) INDC6180 DO 600 IIW=1,LD1 INDC6185 IF (KMOD.EQ.1) GO TO 590 INDC6190 IW = LD2 - IIW INDC6195 GO TO 595 INDC6200 590 IW = IIW INDC6205 595 CALL COMDUB(NST, IW) INDC6210 IC0 = IC0 + ICOUN INDC6215 600 CONTINUE INDC6220 DENS(NST) = 0.0 INDC6225 IF (IC0.EQ.LD1 .AND. IC1.EQ.LD1) GO TO 625 INDC6230 IC1 = 0 INDC6235 WRITE (LOUT,99924) NST INDC6240 99924 FORMAT (///22H THIS IS ITERATION, I3, 18H OF COMBINATORIAL , INDC6245 * 36HOPTIMIZATION, PHASE TWO (SINGLETONS)///) INDC6250 DO 615 IIW=1,LD1 INDC6255 IF (KMOD.EQ.1) GO TO 605 INDC6260 IW = LD2 - IIW INDC6265 GO TO 610 INDC6270 605 IW = IIW INDC6275 610 CALL COMSIN(NST, IW) INDC6280 IC1 = IC1 + ICOUN INDC6285 615 CONTINUE INDC6290 DENS(NST) = 0.0 INDC6295 IF (IC0.EQ.LD1 .AND. IC1.EQ.LD1) GO TO 625 INDC6300 IXV = 0 INDC6305 CALL REGR(0, 0, NST) INDC6310 IXV = 1 INDC6315 620 CONTINUE INDC6320 625 IXV = 0 INDC6325 CALL REGR(0, 0, NST) INDC6330 DO 630 M=1,M3 INDC6335 CALL PRINTD(NST, 0, M) INDC6340 630 CONTINUE INDC6345 IF (NST.EQ.15) WRITE (LOUT,99923) INDC6350 99923 FORMAT (///42H PROGRAM REACHED LIMIT OF 15 ITERATIONS OF, INDC6355 * 28H COMBINATORIAL OPTIMIZATION.///) INDC6360 IF (NST.LT.15) WRITE (LOUT,99922) INDC6365 99922 FORMAT (///42H SINCE THE PRECEDING ITERATION PRODUCED NO, INDC6370 * 50H IMPROVEMENT IN VAF, COMBINATORIAL OPTIMIZATION IS, 7H TERMIN,INDC6375 * 5HATED.///) INDC6380 IF (IPUN.GE.1) CALL PUN(NST) INDC6385 C INDC6390 C SORT THE SUBSETS BY WEIGHT AND PRINT THEM OUT. INDC6395 C INDC6400 635 DO 645 M=1,M3 INDC6405 DO 640 IW=1,LD1 INDC6410 IZ(IW) = IW INDC6415 A(IW) = RISE(M,IW) INDC6420 640 CONTINUE INDC6425 CALL SORT(LD1,A,IZ) INDC6430 C INDC6435 C NOTE THAT IF REGR WERE LATER CALLED, THE FOLLOWING DEFINITION OF INDC6440 C KC = 2 WOULD FOUL UP NNEG = 0 OR 1. INDC6445 C INDC6450 KC = 2 INDC6455 WRITE (LOUT,99921) INDC6460 99921 FORMAT (////45H SUBSETS RANKED ACCORDING TO NUMERICAL WEIGHT///) INDC6465 WRITE (LOUT,99937) TITL INDC6470 CALL PRINTD(NST, IW, M) INDC6475 645 CONTINUE INDC6480 C INDC6485 C THE FOLLOWING IPRN CONDITION ON CALLING FINIS DIFFERS FROM MAPCLUSINDC6490 C INDC6495 CALL FINIS(BMN, C3, C1, IPRN) INDC6500 GO TO 5 INDC6505 650 WRITE (LOUT,99920) INDC6510 99920 FORMAT (///34H ******** TROUBLE EXIT **********//9X, 9H CHECK F,INDC6515 * 21HORMAT OF INPUT DATA. ///) INDC6520 CX XX = ITIMEZ(0)/100000. INDC6525 655 CONTINUE INDC6528 CX WRITE (LOUT,829) XX INDC6530 CX829 FORMAT(//' AT THE END OF EXECUTION, THE TIME IS: ',F8.6///) INDC6535 STOP INDC6540 END INDC6545 CC$ FORTY FORM,XREF NRML 5 SUBROUTINE NRMLZE NRML 10 DOUBLE PRECISION DCON NRML 15 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, NRML 20 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, NRML 25 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, NRML 30 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON NRML 35 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), NRML 40 * ALOS, IXV NRML 45 COMMON /A9/ ILAB(030,6), MLAB(018,6), ALD1, SIGN NRML 50 C NRML 55 C FOR MATRIX CONDITIONAL DATA, EQUATE MEANS AND STANDARD NRML 60 C DEVIATIONS TO BE (0,1) FOR EACH SUBJECT. IF DATA WERE NRML 65 C DISSIMILARITIES, REVERSE THE SIGN. NRML 70 C NRML 75 C IF IT TURNS OUT THAT THE SUBJECTS' MEANS AND STANDARD DEVIATIONS NRML 80 C ARE NOT USED LATER, SDDAT AND DMEAN SHOULD BE CHANGED TO NRML 85 C SCALARS. NRML 90 C NRML 95 C NRML 100 IF (ISWIT.NE.1) GO TO 5 NRML 105 SIGN = -1.0 NRML 110 WRITE (LOUT,99999) NRML 115 99999 FORMAT (//36H THE INPUT DATA ARE DISSIMILARITIES.//) NRML 120 GO TO 10 NRML 125 5 SIGN = 1.0 NRML 130 WRITE (LOUT,99998) NRML 135 99998 FORMAT (//33H THE INPUT DATA ARE SIMILARITIES.//) NRML 140 10 DO 35 M=1,M3 NRML 145 AM1 = 0. NRML 150 AM2 = 0. NRML 155 DO 20 I=2,N NRML 160 I1 = I - 1 NRML 165 DO 15 J=1,I1 NRML 170 D = DATA(M,J,I) NRML 175 AM1 = AM1 + D NRML 180 AM2 = AM2 + D*D NRML 185 15 CONTINUE NRML 190 20 CONTINUE NRML 195 SDDAT(M) = SQRT((AM2-AM1*AM1*ABNV)*ABNV) NRML 200 DMEAN(M) = AM1*ABNV NRML 205 DO 30 I=2,N NRML 210 I1 = I - 1 NRML 215 DO 25 J=1,I1 NRML 220 DATA(M,I,J) = (DATA(M,J,I)-DMEAN(M))/SDDAT(M)*SIGN NRML 225 25 CONTINUE NRML 230 30 CONTINUE NRML 235 35 CONTINUE NRML 240 WRITE (LOUT,99997) NRML 245 99997 FORMAT (///50H USER HAS SPECIFIED THAT THE INPUT DATA ARE MATRIX, NRML 250 * 54H CONDITIONAL (MCON = 1). THUS, THE PROXIMITIES MATRIX, NRML 255 * 24H FOR EACH SOURCE OF DATA/29H WILL BE LINEARLY TRANSFORMED, NRML 260 * 45H (WHEN NECESSSARY) INTO A SIMILARITIES MATRIX, 11H AND NORMAL,NRML 265 * 23HIZED SEPARATELY TO HAVE, 19H ZERO MEAN AND UNIT/10H STANDARD ,NRML 270 * 10HDEVIATION.///1H1, 4X, 11HDATA SOURCE, 6X, 8HRAW MEAN, 6X, NRML 275 * 22HRAW STANDARD DEVIATION//) NRML 280 DO 40 M=1,M3 NRML 285 WRITE (LOUT,99996) (MLAB(M,I),I=1,6), DMEAN(M), SDDAT(M) NRML 290 40 CONTINUE NRML 295 99996 FORMAT (/8X, 6A1, 6X, F10.5, 11X, F10.5) NRML 300 RETURN NRML 305 END NRML 310 CC$ FORTY FORM,XREF FINI 5 SUBROUTINE FINIS(BMN, C3, C1, IPRN) FINI 10 DOUBLE PRECISION DCON FINI 15 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, FINI 20 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, FINI 25 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, FINI 30 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON FINI 35 COMMON /A2/ G(018,030), KT1(10), KT2(10), TK(10), LP(15), IPUN, FINI 40 * IERR, LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), FINI 45 * TRISE(018,017), E(018) FINI 50 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), FINI 55 * ALOS, IXV FINI 60 COMMON /A9/ ILAB(030,6), MLAB(018,6), ALD1, SIGN FINI 65 DO 20 M=1,M3 FINI 70 WK = RISE(M,LD2) FINI 75 DO 15 I=1,N FINI 80 DO 10 J=1,I FINI 85 DATA(M,J,I) = 0. FINI 90 DO 5 K=1,LD1 FINI 95 DATA(M,J,I) = DATA(M,J,I) + RISE(M,K)*P(I,K)*P(J,K) FINI 100 5 CONTINUE FINI 105 DATA(M,J,I) = DATA(M,J,I) + WK FINI 110 10 CONTINUE FINI 115 15 CONTINUE FINI 120 20 CONTINUE FINI 125 IF (IPRN.EQ.1) WRITE (LOUT,99999) FINI 130 99999 FORMAT (1H1, 49H TRANSFORMED TRANSFORMED RAW UNTRANSFORM,FINI 135 * 2HED, 48H RESIDUAL SUBSCRIPTS LABELS SOURCE/ FINI 140 * 26H OBSERVED ESTIMATED, 30H OBSERVED ESTIMATED ,FINI 145 * 4H , 29H OF STIMULI OF STIMULI/18H SIMILARITY SI ,FINI 150 * 32HMILARITY PROXIMITY PROXIMITY, 16X, 1HI, 5X, 1HJ, 6X, 1HI, FINI 155 * 8X, 1HJ//) FINI 160 IF (MCON.EQ.0) C4 = 1.0/(C3*C1) FINI 165 K = 0 FINI 170 DO 50 M=1,M3 FINI 175 DM = DMEAN(M) CONS = RISE(M,LD2) FINI 180 IF (MCON.EQ.0) GO TO 25 FINI 185 C4 = -SDDAT(M)*SIGN FINI 190 BMN = DMEAN(M) FINI 195 25 DO 45 I=1,N FINI 200 DO 40 J=1,I FINI 205 K = K + 1 FINI 210 IF (MOD(K,49).NE.0) GO TO 30 FINI 215 IF (IPRN.EQ.1) WRITE (LOUT,99996) FINI 220 IF (IPRN.EQ.1) WRITE (LOUT,99999) FINI 225 30 IF (MCON .EQ. 0) DATA(M,I,J) = DATA(M,I,J) + DM IF (MCON .EQ. 0) D1 = BMN - DATA(M,I,J)*C4 IF (MCON .EQ. 0) D2 = BMN - (DATA(M,J,I) + DM)*C4 IF (MCON. EQ. 1) D1 = BMN - C4*DATA(M,I,J) FINI 236 IF (MCON. EQ. 1) D2 = BMN - C4*DATA(M,J,I) FINI 238 D3 = D1 - D2 FINI 240 IF (I.NE.J) GO TO 35 FINI 250 D1 = 0.0 FINI 255 D3 = 0.0 FINI 260 35 IF (IPRN.EQ.1) WRITE (LOUT,99998) DATA(M,I,J), DATA(M,J,I), D1, FINI 265 * D2, D3, I, J, (ILAB(I,L),L=1,6), (ILAB(J,L),L=1,6), FINI 270 * (MLAB(M,L),L=1,6) FINI 275 40 CONTINUE FINI 280 45 CONTINUE FINI 285 99998 FORMAT (5(F12.4), 1X, 2I6, 3(3X, 6A1)) FINI 290 K = K + 2 FINI 295 IF (IPRN.EQ.1) WRITE (LOUT,99997) FINI 300 50 CONTINUE FINI 305 99997 FORMAT (//) FINI 310 IF (IPRN.EQ.1) WRITE (LOUT,99996) FINI 315 99996 FORMAT (//49H NOTE: SELF-SIMILARITY AND SELF-PROXIMITY VALUES, FINI 320 * 60H ARE ONLY PREDICTED BY, BUT NOT OBSERVED IN, THE INPUT DATA.) FINI 325 IF (IERR.NE.0) RETURN FINI 330 IF (MCON.EQ.1) GO TO 70 FINI 335 C FINI 340 C THE UNCONDITIONAL CASE FINI 345 C FINI 350 DO 65 M=1,M3 FINI 355 D1 = 0.0 FINI 360 D3 = 0.0 FINI 370 D4 = 0.0 FINI 375 DO 60 I=2,N FINI 380 I1 = I - 1 FINI 385 DO 55 J=1,I1 FINI 390 DZ1 = DATA(M,J,I) FINI 395 DZ3 = DATA(M,I,J) FINI 400 D1 = D1 + (DZ1 - DZ3)**2 D3 = D3 + DZ3 FINI 415 D4 = D4 + DZ3*DZ3 FINI 420 55 CONTINUE FINI 425 60 CONTINUE FINI 430 SDDAT(M) = 1.0 - D1/(D4-D3*D3*ABNV) 65 CONTINUE FINI 440 WRITE (LOUT,99995) TITL FINI 445 99995 FORMAT (1H1, 1X, 72A1, 5X, 22HUNCONDITIONAL ANALYSIS//) FINI 450 GO TO 90 FINI 455 C FINI 460 C THE CONDITIONAL CASE FINI 465 C FINI 470 70 DO 85 M=1,M3 FINI 475 C D1 = 0.0 DO 80 I=2,N FINI 490 I1 = I - 1 FINI 495 DO 75 J=1,I1 FINI 500 D1 = D1 + (DATA(M,J,I) - DATA(M,I,J))**2 75 CONTINUE FINI 520 80 CONTINUE FINI 525 SDDAT(M) = 1.0 - D1*ABNV 85 CONTINUE FINI 535 WRITE (LOUT,99994) TITL FINI 540 99994 FORMAT (1H1, 1X, 72A1, 5X, 20HCONDITIONAL ANALYSIS//) FINI 545 90 WRITE (LOUT,99993) FINI 550 99993 FORMAT (/////5X, 12H DATA SOURCE, 8X,12HSOURCE'S VAF) FINI 555 DO 95 M=1,M3 FINI 560 WRITE (LOUT,99992) M, (MLAB(M,I),I=1,6), SDDAT(M) FINI 565 95 CONTINUE FINI 570 99992 FORMAT (//1X, I3, 4X, 6A1, 10X, F10.5) FINI 575 RETURN FINI 580 END FINI 585 CC$ FORTY FORM HF00 5 FUNCTION HF(X) HF00 10 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, HF00 15 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, HF00 20 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, HF00 25 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON HF00 30 C HF00 35 C HF00 40 HF = X - ABNV*X*X HF00 45 RETURN HF00 50 END HF00 55 CC$ FORTY FORM,XREF COMS 5 SUBROUTINE COMSIN(MST, IW) COMS 10 DOUBLE PRECISION DCON COMS 15 INTEGER TITL COMS 20 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, COMS 25 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, COMS 30 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, COMS 35 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON COMS 40 COMMON /A2/ G(018,030), KT1(10), KT2(10), TK(10), LP(15), IPUN, COMS 45 * IERR, LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), COMS 50 * TRISE(018,017), E(018) COMS 55 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), COMS 60 * ALOS, IXV COMS 65 COMMON /A8/ W(018,017), IPRE, NNEG, NDOT, NLD1 COMS 70 COMMON /A9/ ILAB(030,6), MLAB(018,6), ALD1, SIGN COMS 75 ICOUN = 1 COMS 80 5 CALL RESID(IW) COMS 85 C COMS 90 C THAT GOT THE CENTERED RESIDUALS; NOW COMPUTE NUMERATOR (E) AND COMS 95 C DENOMINATOR (F) FOR V**2 COMS 100 C COMS 105 F = 0.0 COMS 110 DO 20 M=1,M3 COMS 115 E(M) = 0.0 COMS 120 DO 15 I=2,N COMS 125 I1 = I - 1 COMS 130 PI = P(I,IW) COMS 135 DO 10 J=1,I1 COMS 140 PIJ = P(J,IW)*PI COMS 145 IF (M.EQ.1) F = F + PIJ COMS 150 E(M) = E(M) + DATA(M,J,I)*PIJ COMS 155 10 CONTINUE COMS 160 15 CONTINUE COMS 165 20 CONTINUE COMS 170 VK2 = 0.0 COMS 175 DO 25 M=1,M3 COMS 180 VK2 = VK2 + E(M)**2 COMS 185 25 CONTINUE COMS 190 VK2 = VK2/HF(F) COMS 195 99999 FORMAT (//47H VARIANCE IN THE RESIDUALS CURRENTLY ACCOUNTED , COMS 200 * 15H FOR BY SUBSET , I3, 2H =, F10.5) COMS 205 HTOT = 0.0 COMS 210 DO 30 I=1,N COMS 215 HTOT = HTOT + P(I,IW) COMS 220 30 CONTINUE COMS 225 IF (LP(MST).EQ.1) WRITE (LOUT,99999) IW, VK2 COMS 230 DO 45 M=1,M3 COMS 235 DO 40 I=1,N COMS 240 G(M,I) = 0.0 COMS 245 DO 35 J=1,N COMS 250 IF (J.EQ.I) GO TO 35 COMS 255 K1 = MIN0(I,J) COMS 260 K2 = MAX0(I,J) COMS 265 G(M,I) = G(M,I) + P(J,IW)*DATA(M,K1,K2) COMS 270 35 CONTINUE COMS 275 40 CONTINUE COMS 280 45 CONTINUE COMS 285 DO 55 M=1,M3 COMS 290 DO 50 K=1,LD2 COMS 295 TRISE(M,K) = RISE(M,K) COMS 300 50 CONTINUE COMS 305 55 CONTINUE COMS 310 IC10 = 0 COMS 315 RMAX = -1.0E20 COMS 320 DO 70 I=1,N COMS 325 PI2 = 2.0*P(I,IW) - 1.0 COMS 330 HT2 = HTOT - PI2 COMS 335 C COMS 340 C DO NOT CONSIDER DROPPING EITHER MEMBER OF A DYAD SUBSET COMS 345 C OR FORMING THE COMPLETE SUBSET. COMS 350 C COMS 355 IF (HT2.LT.1.98 .OR. HT2.GT.AN1) GO TO 70 COMS 360 R2 = 0.0 COMS 365 DO 60 M=1,M3 COMS 370 ANUM = E(M) - PI2*G(M,I) COMS 375 IF (ABS(ANUM).LE.1.0E-12) GO TO 65 COMS 380 R2 = R2 + ANUM**2 COMS 385 60 CONTINUE COMS 390 R2 = R2/HF(F-PI2*(HTOT-P(I,IW))) COMS 395 IF (R2.LE.RMAX) GO TO 65 COMS 400 C COMS 405 C UPDATE BEST V**2 COMS 410 C COMS 415 MR2 = I COMS 420 RMAX = R2 COMS 425 C COMS 430 C PRINT OUT (10 AT A TIME) THE ALLOWABLE SINGLETON CHANGES COMS 435 C AND RESULTING V**2 COMS 440 C COMS 445 65 CONTINUE COMS 450 IF (LP(MST).NE.1) GO TO 70 COMS 455 IC10 = IC10 + 1 COMS 460 KT1(IC10) = -FLOAT(I)*PI2 COMS 465 TK(IC10) = R2 COMS 470 IF ((IC10.LT.10 .AND. I.LT.N) .OR. (I.EQ.N .AND. IC10.EQ.0)) GO COMS 475 * TO 70 COMS 480 WRITE (LOUT,99998) (KT1(K),TK(K),K=1,IC10) COMS 485 99998 FORMAT (1X, 10(1X, 1H(, I3, F7.2, 1H))) COMS 490 IC10 = 0 COMS 495 70 CONTINUE COMS 500 IF (IC10.NE.0) WRITE (LOUT,99998) (KT1(K),TK(K),K=1,IC10) COMS 505 IF (RMAX.LE.VK2) RETURN COMS 510 P(MR2,IW) = 1.0 - P(MR2,IW) COMS 515 UVAF = TVAF COMS 520 CALL REGR(0, 0, MST) COMS 525 IF (IERR.EQ.4) GO TO 95 COMS 530 IF (NNEG.GT.0) GO TO 90 COMS 535 C COMS 540 C FIRST CHECK TO SEE IF THE WEIGHT OF ANY SUBSET WOULD COMS 545 C BECOME NEGATIVE. COMS 550 C COMS 555 DO 80 M=1,M3 COMS 560 LM = M COMS 565 DO 75 IK=1,LD1 COMS 570 LK = IK COMS 575 C COMS 580 C COMS 585 C IN MAPCLUS, THE FOLLOWING STATEMENT USES .LE. INSTEAD COMS 590 C OF .LT. AND .GT. INSTEAD OF .GE. COMS 595 C COMS 600 IF (RISE(M,IK).LT.0.0 .AND. TRISE(M,IK).GE.0.0) GO TO 85 COMS 605 75 CONTINUE COMS 610 80 CONTINUE COMS 615 GO TO 110 COMS 620 85 WRITE (LOUT,99997) (ILAB(MR2,I),I=1,6), IW, TVAF, UVAF, LK, COMS 625 * (MLAB(LM,I),I=1,6), RISE(LM,LK) COMS 630 99997 FORMAT (//19H REVERSING ELEMENT , 1X, 6A1, 11H OF SUBSET , I3, COMS 635 * 31H WOULD INCREASE OVERALL VAF TO , F10.6, 12H, REPLACING , COMS 640 * F10.6/55H BUT WILL NOT BE DONE SINCE THE CURRENTLY POSITIVE WEIG,COMS 645 * 6HHT FOR, 7H SUBSET, I3, 13H FROM SOURCE , 6A1, 13H BECOMES NEG,COMS 650 * 7HATIVE (, F12.7, 2H).) COMS 655 GO TO 95 COMS 660 90 IF (TVAF.GT.UVAF) GO TO 110 COMS 665 C COMS 670 C UNSUCCESSFUL CHANGE COMS 675 C COMS 680 95 P(MR2,IW) = 1.0 - P(MR2,IW) COMS 685 TVAF = UVAF COMS 690 DO 105 M=1,M3 COMS 695 DO 100 K=1,LD2 COMS 700 RISE(M,K) = TRISE(M,K) COMS 705 100 CONTINUE COMS 710 105 CONTINUE COMS 715 RETURN COMS 720 C COMS 725 C FOR A SUCCESSFUL CHANGE COMS 730 C COMS 735 110 ICOUN = 2 COMS 740 WRITE (LOUT,99996) (ILAB(MR2,I),I=1,6), IW, TVAF, UVAF COMS 745 99996 FORMAT (//19H REVERSING ELEMENT , 1X, 6A1, 11H OF SUBSET , I3, COMS 750 * 33H GAVE AN IMPROVED OVERALL VAF OF , F10.6, 11H REPLACING , COMS 755 * F10.6//) COMS 760 UVAF = TVAF COMS 765 C COMS 770 C NOW GO BACK AND LOOK FOR ANY OTHER SUCCESSFUL CHANGES WITHIN COMS 775 C SUBSET IW. COMS 780 C COMS 785 GO TO 5 COMS 790 END COMS 795 CC$ FORTY FORM,XREF COMD 5 SUBROUTINE COMDUB(MST, IW) COMD 10 DOUBLE PRECISION DCON COMD 15 INTEGER TITL COMD 20 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, COMD 25 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, COMD 30 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, COMD 35 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON COMD 40 COMMON /A2/ G(018,030), KT1(10), KT2(10), TK(10), LP(15), IPUN, COMD 45 * IERR, LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), COMD 50 * TRISE(018,017), E(018) COMD 55 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), COMD 60 * ALOS, IXV COMD 65 COMMON /A8/ W(018,017), IPRE, NNEG, NDOT, NLD1 COMD 70 COMMON /A9/ ILAB(030,6), MLAB(018,6), ALD1, SIGN COMD 75 ICOUN = 1 COMD 80 5 CALL RESID(IW) COMD 85 C COMD 90 C THAT GOT THE CENTERED RESIDUALS; NOW COMPUTE NUMERATOR (E) AND COMD 95 C DENOMINATOR (F) FOR V**2 COMD 100 C COMD 105 F = 0.0 COMD 110 DO 20 M=1,M3 COMD 115 E(M) = 0.0 COMD 120 DO 15 I=2,N COMD 125 I1 = I - 1 COMD 130 PI = P(I,IW) COMD 135 DO 10 J=1,I1 COMD 140 PIJ = P(J,IW)*PI COMD 145 IF (M.EQ.1) F = F + PIJ COMD 150 E(M) = E(M) + DATA(M,J,I)*PIJ COMD 155 10 CONTINUE COMD 160 15 CONTINUE COMD 165 20 CONTINUE COMD 170 VK2 = 0.0 COMD 175 DO 25 M=1,M3 COMD 180 VK2 = VK2 + E(M)**2 COMD 185 25 CONTINUE COMD 190 VK2 = VK2/HF(F) COMD 195 99999 FORMAT (//47H VARIANCE IN THE RESIDUALS CURRENTLY ACCOUNTED , COMD 200 * 15H FOR BY SUBSET , I3, 2H =, F10.5) COMD 205 HTOT = 0.0 COMD 210 DO 30 I=1,N COMD 215 HTOT = HTOT + P(I,IW) COMD 220 30 CONTINUE COMD 225 IF (LP(MST).EQ.1) WRITE (LOUT,99999) IW, VK2 COMD 230 DO 45 M=1,M3 COMD 235 DO 40 I=1,N COMD 240 G(M,I) = 0.0 COMD 245 DO 35 J=1,N COMD 250 IF (J.EQ.I) GO TO 35 COMD 255 K1 = MIN0(I,J) COMD 260 K2 = MAX0(I,J) COMD 265 G(M,I) = G(M,I) + P(J,IW)*DATA(M,K1,K2) COMD 270 35 CONTINUE COMD 275 40 CONTINUE COMD 280 45 CONTINUE COMD 285 DO 55 M=1,M3 COMD 290 DO 50 K=1,LD2 COMD 295 TRISE(M,K) = RISE(M,K) COMD 300 50 CONTINUE COMD 305 55 CONTINUE COMD 310 IC10 = 0 COMD 315 RMAX = -1.0E20 COMD 320 DO 75 I=2,N COMD 325 I1 = I - 1 COMD 330 PI = P(I,IW) COMD 335 PI2 = 2.0*PI - 1.0 COMD 340 DO 70 J=1,I1 COMD 345 PJ = P(J,IW) COMD 350 PJ2 = 2.0*PJ - 1.0 COMD 355 C COMD 360 C DO NOT CONSIDER DROPPING EITHER MEMBER OF A DYAD SUBSET COMD 365 C OR FORMING THE COMPLETE SUBSET. COMD 370 C COMD 375 HT2 = HTOT - PI2 - PJ2 COMD 380 IF (HT2.LT.1.98 .OR. HT2.GT.AN1) GO TO 70 COMD 385 VIJ = PI2*PJ2 COMD 390 C COMD 395 C COMPUTE VARIANCE ACCOUNTED FOR BY SUBSET IW, WHEN A COMD 400 C DOUBLETON CHANGE IS MADE, REVERSING ELEMENTS I AND J OF SUBSET COMD 405 C IW. COMD 410 C COMD 415 C COMD 420 R3 = 0. COMD 425 DO 60 M=1,M3 COMD 430 ANUM = E(M) - PI2*G(M,I) - PJ2*G(M,J) + DATA(M,J,I)*VIJ COMD 435 IF (ABS(ANUM).LT.1.0E-12) GO TO 65 COMD 440 R3 = R3 + (ANUM**2) COMD 445 60 CONTINUE COMD 450 R3 = R3/HF(F-PI2*(HTOT-PI)-PJ2*(HTOT-PJ)+VIJ) COMD 455 IF (R3.LE.RMAX) GO TO 65 COMD 460 C COMD 465 C UPDATE BEST V**2 COMD 470 C COMD 475 MR2 = I COMD 480 MR3 = J COMD 485 RMAX = R3 COMD 490 C COMD 495 C PRINT OUT (8 AT A TIME) THE PERMISSIBLE DOUBLETON CHANGES AND COMD 500 C RESULTING V**2. COMD 505 C COMD 510 65 CONTINUE COMD 515 IF (LP(MST).NE.1) GO TO 70 COMD 520 IC10 = IC10 + 1 COMD 525 KT2(IC10) = -FLOAT(J)*PJ2 COMD 530 KT1(IC10) = -FLOAT(I)*PI2 COMD 535 TK(IC10) = R3 COMD 540 IF ((IC10.LT.8 .AND. (I.LT.N .OR. J.LT.N1)) .OR. (I.EQ.N .AND. COMD 545 * J.EQ.N1 .AND. IC10.EQ.0)) GO TO 70 COMD 550 WRITE (LOUT,99998) (KT1(K),KT2(K),TK(K),K=1,IC10) COMD 555 99998 FORMAT (1X, 8(1X, 1H(, 2I3, F7.2, 1H))) COMD 560 IC10 = 0 COMD 565 70 CONTINUE COMD 570 75 CONTINUE COMD 575 IF (IC10.NE.0) WRITE (LOUT,99998) (KT1(K),KT2(K),TK(K),K=1,IC10) COMD 580 IF (RMAX.LE.VK2) RETURN COMD 585 P(MR3,IW) = 1.0 - P(MR3,IW) COMD 590 P(MR2,IW) = 1.0 - P(MR2,IW) COMD 595 UVAF = TVAF COMD 600 CALL REGR(0, 0, MST) COMD 605 IF (IERR.EQ.4) GO TO 100 COMD 610 IF (NNEG.GT.0) GO TO 95 COMD 615 C COMD 620 C FIRST CHECK TO SEE IF THE WEIGHT OF ANY SUBSET WOULD COMD 625 C BECOME NEGATIVE. COMD 630 C COMD 635 DO 85 M=1,M3 COMD 640 LM = M COMD 645 DO 80 IK=1,LD1 COMD 650 LK = IK COMD 655 C COMD 660 C COMD 665 C IN MAPCLUS, THE FOLLOWING STATEMENT USES .LE. INSTEAD COMD 670 C OF .LT. AND .GT. INSTEAD OF .GE. COMD 675 C COMD 680 IF (RISE(M,IK).LT.0.0 .AND. TRISE(M,IK).GE.0.0) GO TO 90 COMD 685 80 CONTINUE COMD 690 85 CONTINUE COMD 695 GO TO 115 COMD 700 90 WRITE (LOUT,99997) (ILAB(MR2,I),I=1,6), (ILAB(MR3,I),I=1,6), IW, COMD 705 * TVAF, UVAF, LK, (MLAB(LM,I),I=1,6), RISE(LM,LK) COMD 710 99997 FORMAT (/20H REVERSING ELEMENTS , 1X, 6A1, 1X, 6A1, 10H OF SUBSET,COMD 715 * 1H , I3, 31H WOULD INCREASE OVERALL VAF TO , F10.6, 9H, REPLACI, COMD 720 * 3HNG , F10.6/48H BUT WILL NOT BE DONE SINCE THE CURRENTLY POSITI,COMD 725 * 13HVE WEIGHT FOR, 7H SUBSET, I3, 13H FROM SOURCE , 6A1, 6H BECOM,COMD 730 * 13HES NEGATIVE (, F12.7, 2H).) COMD 735 GO TO 100 COMD 740 95 IF (TVAF.GT.UVAF) GO TO 115 COMD 745 C COMD 750 C UNSUCCESSFUL CHANGE COMD 755 C COMD 760 100 P(MR2,IW) = 1.0 - P(MR2,IW) COMD 765 P(MR3,IW) = 1.0 - P(MR3,IW) COMD 770 TVAF = UVAF COMD 775 DO 110 M=1,M3 COMD 780 DO 105 K=1,LD2 COMD 785 RISE(M,K) = TRISE(M,K) COMD 790 105 CONTINUE COMD 795 110 CONTINUE COMD 800 RETURN COMD 805 C COMD 810 C FOR A SUCCESSFUL CHANGE COMD 815 C COMD 820 115 ICOUN = 2 COMD 825 WRITE (LOUT,99996) (ILAB(MR2,I),I=1,6), (ILAB(MR3,I),I=1,6), IW, COMD 830 * TVAF, UVAF COMD 835 99996 FORMAT (/20H REVERSING ELEMENTS , 6A1, 1X, 6A1, 11H OF SUBSET , COMD 840 * I3, 33H GAVE AN IMPROVED OVERALL VAF OF , F10.6, 11H REPLACING , COMD 845 * F10.6//) COMD 850 UVAF = TVAF COMD 855 C COMD 860 C NOW GO BACK AND LOOK FOR ANY OTHER SUCCESSFUL CHANGES WITHIN COMD 865 C SUBSET IW. COMD 870 C COMD 875 GO TO 5 COMD 880 END COMD 885 CC$ FORTY FORM,XREF PUN0 5 SUBROUTINE PUN(MST) PUN0 10 INTEGER FMT4, FMT5 PUN0 15 DIMENSION LIST(030), FMT4(12), FMT5(24) PUN0 20 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, PUN0 25 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, PUN0 30 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, PUN0 35 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON PUN0 40 COMMON /A2/ G(018,030), KT1(10), KT2(10), TK(10), LP(15), IPUN, PUN0 45 * IERR, LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), PUN0 50 * TRISE(018,017), E(018) PUN0 55 DATA FMT4(1), FMT4(2), FMT4(3), FMT4(4), FMT4(5), FMT4(6), PUN0 60 * FMT4(7), FMT4(8), FMT4(9), FMT4(10), FMT4(11), FMT4(12) /1H(,1H1,PUN0 65 * 1H4,1HX,1H,,1H5,1H8,1HF,1H1,1H.,1H0,1H)/ PUN0 70 DATA FMT5(1), FMT5(2), FMT5(3), FMT5(4), FMT5(5), FMT5(6), PUN0 75 * FMT5(7), FMT5(8), FMT5(9), FMT5(10), FMT5(11), FMT5(12), PUN0 80 * FMT5(13), FMT5(14), FMT5(15), FMT5(16), FMT5(17), FMT5(18), PUN0 85 * FMT5(19), FMT5(20), FMT5(21), FMT5(22), FMT5(23), FMT5(24) /1H(, PUN0 90 * 1H4,1HX,1H,,1H6,1HF,1H9,1H.,1H5,1H/,1H(,1HF,1H1,1H3,1H.,1H5,1H,, PUN0 95 * 1H5,1HF,1H9,1H.,1H5,1H),1H)/ PUN0 100 C PUN0 105 C SUBROUTINE TO PUNCH (TO DISK) POST-POLISHED SUBSETS PUN0 110 C PUN0 115 C PUN0 120 C NOTE THAT IF THE NUMBER OF STIMULI EVER EXCEEDS 58, OR IF PUN0 125 C CHANGES ARE MADE IN OUTPUT FORMAT STATEMENT 2, THEN FORMAT PUN0 130 C STATEMENT 1 AND ARRAY FMT4 MUST BE MODIFIED. PUN0 135 C PUN0 140 WRITE (LPUNCH,99999) MST, (TITL(I),I=1,61) PUN0 145 99999 FORMAT (4H IT., I3, 12H CONFIG FOR , 61A1) PUN0 150 WRITE (LPUNCH,99998) FMT4 PUN0 155 99998 FORMAT (27A1) PUN0 160 DO 10 I=1,LD1 PUN0 165 DO 5 J=1,N PUN0 170 LIST(J) = 0 PUN0 175 IF (P(J,I).GT.0.98 .AND. P(J,I).LT.1.02) LIST(J) = 1 PUN0 180 5 CONTINUE PUN0 185 C PUN0 190 C SUBSCRIPT M ARBITRARILY SET = 1 IN RISE(M,I) BELOW: PUN0 195 C PUN0 200 WRITE (LPUNCH,99997) I, RISE(1,I), (LIST(J),J=1,N) PUN0 205 10 CONTINUE PUN0 210 99997 FORMAT (1X, I3, 1X, F8.5, 1X, 58I1) PUN0 215 IF (IPUN.EQ.1) RETURN PUN0 220 WRITE (LPUNCH,99998) FMT5 PUN0 225 DO 15 I=1,M3 PUN0 230 WRITE (LPUNCH,99996) I, (RISE(I,J),J=1,LD2) PUN0 235 15 CONTINUE PUN0 240 99996 FORMAT (1X, I3, 6F9.5/(F13.5, 5F9.5)) PUN0 245 RETURN PUN0 250 END PUN0 255 CC$ FORTY FORM,XREF REGR 5 SUBROUTINE REGR(MST, IW, MRT) REGR 10 DOUBLE PRECISION DCON, Y REGR 15 DIMENSION Q(0435,016), S(0435), X(016), W2(016), ZZ(0435), REGR 20 * INDEX(016) REGR 25 DIMENSION Y(018,017), IPVT(017), B(017,017), R(017,017) REGR 30 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, REGR 35 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, REGR 40 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, REGR 45 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON REGR 50 COMMON /A2/ G(018,030), KT1(10), KT2(10), TK(10), LP(15), IPUN, REGR 55 * IERR, LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), REGR 60 * TRISE(018,017), E(018) REGR 65 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), REGR 70 * ALOS, IXV REGR 75 COMMON /A8/ W(018,017), IPRE, NNEG, NDOT, NLD1 REGR 80 COMMON /A9/ ILAB(030,6), MLAB(018,6), ALD1, SIGN REGR 85 C REGR 90 C REGR 95 C REGR 100 C NOTE THAT SUBROUTINE REGR USES THE TRANSFORMED DATA REGR 105 C (LOWERHALF OF ARRAY DATA), UNLIKE MOST OF THE OTHER REGR 110 C SUBROUTINES, WHICH USE THE RESIDUALS (UPPERHALF). REGR 115 C REGR 120 C REGR 125 IF (MST.NE.1 .OR. IW.EQ.LD1 .OR. ICON.GT.0) GO TO 10 REGR 130 C REGR 135 C REGR 140 C FOR GLOBAL REGRESSION DURING THE FIRST PRE-POLISHING REGR 145 C MAJOR ITERATION WHEN LESS THAN ALL LD1 SUBSETS HAVE BEEN REGR 150 C FORMED, THE ADDITIVE CONSTANT HAS TO BE TACKED ON AS SUBSET REGR 155 C (IW + 1), UNLESS A USER-SUPPLIED INITIAL CONFIGURATION WAS GIVEN. REGR 160 C REGR 165 C REGR 170 C REGR 175 LD0 = IW + 1 REGR 180 DO 5 J=1,N REGR 185 P(J,LD0) = 1.0 REGR 190 5 CONTINUE REGR 195 GO TO 15 REGR 200 10 LD0 = LD2 REGR 205 15 LDA = LD0 - 1 REGR 210 DO 35 M=1,M3 REGR 215 DO 30 I=1,LDA REGR 220 Y(M,I) = 0.0D0 REGR 225 DO 25 J=2,N REGR 230 PJI = P(J,I) REGR 235 J1 = J - 1 REGR 240 DO 20 L=1,J1 REGR 245 Y(M,I) = Y(M,I) + PJI*P(L,I)*DATA(M,J,L) REGR 250 20 CONTINUE REGR 255 25 CONTINUE REGR 260 30 CONTINUE REGR 265 35 CONTINUE REGR 270 DO 45 I=1,N2 REGR 275 DO 40 J=1,N2 REGR 280 R(I,J) = 0. REGR 285 40 CONTINUE REGR 290 45 CONTINUE REGR 295 DO 65 I=1,LD0 REGR 300 DO 60 J=1,I REGR 305 DO 55 L=2,N REGR 310 PLI = P(L,I) REGR 315 PLJ = P(L,J) REGR 320 L1 = L - 1 REGR 325 DO 50 M1=1,L1 REGR 330 R(I,J) = R(I,J) + PLI*PLJ*P(M1,I)*P(M1,J) REGR 335 50 CONTINUE REGR 340 55 CONTINUE REGR 345 60 CONTINUE REGR 350 65 CONTINUE REGR 355 DO 75 J=1,LDA REGR 360 DO 70 L=1,J REGR 365 R(J,L) = R(J,L) - ABNV*R(LD0,J)*R(LD0,L) REGR 370 70 CONTINUE REGR 375 75 CONTINUE REGR 380 CALL INVS2(LDA, R, N2, B, N2, IPVT, IERR) REGR 385 IF (IERR.EQ.0) GO TO 80 REGR 390 WRITE (LOUT,99999) REGR 395 99999 FORMAT (///51H ******** TROUBLE ENCOUNTERED WHILE INVERTING MATRI,REGR 400 * 1HX, 24H FOR REGRESSION ********/10X, 23HLINEAR DEPENDENCIES LIK,REGR 405 * 18HELY TO BE PRESENT.//) REGR 410 IF (LD0.NE.LD2) GO TO 190 REGR 415 RETURN REGR 420 80 DO 95 M=1,M3 REGR 425 DO 90 I=1,LDA REGR 430 W(M,I) = 0.0 REGR 435 DO 85 J=1,LDA REGR 440 W(M,I) = W(M,I) + Y(M,J)*B(J,I) REGR 445 85 CONTINUE REGR 450 90 CONTINUE REGR 455 95 CONTINUE REGR 460 DO 105 M=1,M3 REGR 465 DM = 0.0 REGR 466 IF (MCON.EQ.0) DM = DMEAN(M) REGR 468 SUMK = 0.0 REGR 470 DO 100 I=1,LDA REGR 475 SUMK = SUMK + W(M,I)*R(LD0,I) REGR 480 100 CONTINUE REGR 485 RISE(M,LD2) = DM - ABNV*SUMK REGR 490 105 CONTINUE REGR 495 C REGR 500 JB = LD1 REGR 505 IF (LD0.NE.LD2) JB = IW REGR 510 DO 170 M=1,M3 REGR 515 IF (NNEG+KC.LT.2) GO TO 160 REGR 520 C REGR 525 C CHECK ARRAY BOUNDS USED FOR NON-NEGATIVE REGRESSION REGR 530 C REGR 535 IF (NDO.LE.NDOT) GO TO 110 REGR 540 WRITE (LOUT,99998) NDOT, N, NDO REGR 545 99998 FORMAT (51H1IN THIS VERSION OF INDCLUS, THE UPPER BOUND ON THE, REGR 550 * 31H NUMBER OF PAIRS OF STIMULI IS , I3/19H HOWEVER, SINCE THE, REGR 555 * 10H USER HAS , I3, 37H STIMULI, THE NUMBER OF PAIRS MUST BE, REGR 560 * 13H INCREASED TO, I5, 1H./34H INCREASE THE VALUE OF NDOT IN THE, REGR 565 * 47H MAIN PROGRAM AND THE ARRAYS IN SUBROUTINE REGR, 9H ACCORDIN, REGR 570 * 4HGLY./22H EXECUTION TERMINATES.) REGR 575 WRITE (LOUT,99997) REGR 580 99997 FORMAT (//45H SEE DOCUMENTATION ON INDCLUS FOR DETAILS ON , REGR 585 * 32HINCREASING DIMENSION STATEMENTS.//) REGR 590 STOP REGR 595 110 DO 115 I=1,JB REGR 600 IF (W(M,I).LT.0.0) GO TO 120 REGR 605 115 CONTINUE REGR 610 GO TO 160 REGR 615 120 M1 = 0 REGR 620 DO 130 I=2,N REGR 625 I1 = I - 1 REGR 630 DO 125 J=1,I1 REGR 635 M1 = M1 + 1 REGR 640 C REGR 645 C PUT CENTERED DATA IN COLUMN VECTOR S REGR 650 C REGR 655 S(M1) = DATA(M,I,J) REGR 660 125 CONTINUE REGR 665 130 CONTINUE REGR 670 C REGR 675 C IN PRINCIPLE, Q SHOULD BE COMPUTED ONLY ONCE FOR A GIVEN REGR 680 C VALUE OF M. HOWEVER, SINCE SUBROUTINE NNLS OVERWRITES Q REGR 685 C IT SEEMS PREFERABLE TO RECOMPUTE Q THAN TO CREATE ANOTHER REGR 690 C LARGE ARRAY TO HOLD A COPY OF Q. REGR 695 C REGR 700 DO 150 K=1,LD1 REGR 705 QSUM = 0.0 REGR 710 M1 = 0 REGR 715 DO 140 I=2,N REGR 720 PIK = P(I,K) REGR 725 I1 = I - 1 REGR 730 DO 135 J=1,I1 REGR 735 Q1 = PIK*P(J,K) REGR 740 QSUM = QSUM + Q1 REGR 745 M1 = M1 + 1 REGR 750 Q(M1,K) = Q1 REGR 755 135 CONTINUE REGR 760 140 CONTINUE REGR 765 QSUM = QSUM*ABNV REGR 770 DO 145 M1=1,NDO REGR 775 Q(M1,K) = Q(M1,K) - QSUM REGR 780 145 CONTINUE REGR 785 150 CONTINUE REGR 790 C WRITE (LOUT,400) M,(W(M,I), I = 1, LD1) REGR 795 C 400 FORMAT(//' PRIOR TO CALLING NNLS, W(',I2,',I) = ' REGR 800 C 2 /(2X,13F10.5/)) REGR 805 CALL NNLS(Q, NDOT, NDO, LD1, S, X, RNORM, W2, ZZ, INDEX, MODE, REGR 810 * LOUT) REGR 815 C WRITE (LOUT,410) (X(I), I = 1, LD1) REGR 820 C 410 FORMAT(//' AFTER CALLING NNLS, X(I) = '/(2X,13F10.5/)) REGR 825 C WRITE (LOUT,324) M,MODE,(INDEX(I),I = 1, LD0) REGR 830 C 324 FORMAT (///' FOR SUBJECT ',I2,' MODE =',I3/1X,30I4//) REGR 835 DO 155 I=1,JB REGR 840 RISE(M,I) = X(I) REGR 845 155 CONTINUE REGR 850 GO TO 170 REGR 855 160 DO 165 I=1,JB REGR 860 RISE(M,I) = W(M,I) REGR 865 165 CONTINUE REGR 870 170 CONTINUE REGR 875 IF (LD0.EQ.LD2) GO TO 185 REGR 880 DO 180 M=1,M3 REGR 885 DO 175 J=LD0,LD1 REGR 890 RISE(M,J) = 0.0 REGR 895 175 CONTINUE REGR 900 180 CONTINUE REGR 905 185 MRT1 = MST REGR 910 IF (MRT1.EQ.0) MRT1 = MRT REGR 915 IF (LD0.NE.LD2) GO TO 190 REGR 920 CALL VARIAN(MRT1, 2, IW) REGR 925 RETURN REGR 930 190 DO 195 J=1,N REGR 935 P(J,LD0) = 0. REGR 940 195 CONTINUE REGR 945 RETURN REGR 950 END REGR 955 CC$ FORTY FORM,XREF NNLS 5 C NON-NEGATIVE LEAST SQUARES SUBROUTINE FROM LAWSON & HANSON NNLS 10 C SUBROUTINE NNLS (A,MDA,M,N,B,X,RNORM,W,ZZ,INDEX,MODE) NNLS 15 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY. 1973 JUNE 15NNLS 20 C TAKEN FROM "SOLVING LEAST SQUARES PROBLEMS". PRENTICE-HALL. 1974 NNLS 25 C NNLS 30 C ********** NONNEGATIVE LEAST SQUARES ********** NNLS 35 C NNLS 40 C GIVEN AN M BY N MATRIX, A, AND AN M-VECTOR, B, COMPUTE AN NNLS 45 C N-VECTOR, X, WHICH SOLVES THE LEAST SQUARES PROBLEM NNLS 50 C NNLS 55 C A * X = B SUBJECT TO X .GE. 0 NNLS 60 C NNLS 65 C A(),MDA,M,N MDA IS THE FIRST DIMENSIONING PARAMETER FOR THE NNLS 70 C ARRAY, A(). ON ENTRY A() CONTAINS THE M BY N NNLS 75 C MATRIX, A. ON EXIT A() CONTAINS NNLS 80 C THE PRODUCT MATRIX, Q*A WHERE Q IS AN NNLS 85 C M BY M ORTHOGONAL MATRIX GENERATED IMPLICITLY BY NNLS 90 C THIS SUBROUTINE. NNLS 95 C B() ON ENTRY B() CONTAINS THE M-VECTOR. B. ON EXIT B() CON-NNLS 100 C TAINS Q*B NNLS 105 C X() ON ENTRY X() NEED NOT BE INITIALIZED. ON EXIT X() WILL NNLS 110 C CONTAIN THE SOLUTION VECTOR. NNLS 115 C RNORM ON EXIT RNORM CONTAINS THE EUCLIDEAN NORM OF THE NNLS 120 C RESIDUAL VECTOR. NNLS 125 C W() AN N-ARRAY OF WORKING SPACE. ON EXIT W() WILL CONTAIN NNLS 130 C THE DUAL SOLUTION VECTOR. W WILL SATISFY W(I) = 0. NNLS 135 C FOR ALL I IN SET P AND W(I) .LE. 0. FOR ALL I IN SET Z NNLS 140 C ZZ() AN M-ARRAY OF WORKING SPACE. NNLS 145 C INDEX() AN INTEGER WORKING ARRAY OF LENGTH AT LEAST N. NNLS 150 C ON EXIT THE CONTENTS OF THIS ARRAY DEFINE THE SETS NNLS 155 C P AND Z AS FOLLOWS.. NNLS 160 C NNLS 165 C INDEX(1) THRU INDEX(NSETP) = P. NNLS 170 C INDEX(IZ1) THRU INDEX(IZ2) = Z. NNLS 175 C IZ1 = NSETP + = NPP1 NNLS 180 C IZ2 = N NNLS 185 C MODE THIS IS A SUCCESS-FAILURE FLAG WITH THE FOLLOWING NNLS 190 C MEANINGS NNLS 195 C 1 THE SOLUTION HAS BEEN COMPUTED SUCCESSFULLY. NNLS 200 C 2 THE DIMENSIONS OF THE PROBLEM ARE BAD. NNLS 205 C EITHER M .LE. 0 OR N .LE. 0. NNLS 210 C 3 ITERATION COUNT EXCEEDED. MORE THAN 3*N ITERATIONS NNLS 215 C NNLS 220 SUBROUTINE NNLS(A, MDA, M, N, B, X, RNORM, W, ZZ, INDEX, MODE, NNLS 225 * LOUT) NNLS 230 DIMENSION A(MDA,N), B(M), X(N), W(N), ZZ(M),DUMMY(1) NNLS 235 INTEGER INDEX(N) NNLS 240 ZERO = 0. NNLS 245 TWO = 2. NNLS 250 FACTOR = 0.01 NNLS 255 C NNLS 260 MODE = 1 NNLS 265 IF (M.GT.0 .AND. N.GT.0) GO TO 5 NNLS 270 MODE = 2 NNLS 275 RETURN NNLS 280 5 ITER = 0 NNLS 285 ITMAX = 3*N NNLS 290 C NNLS 295 C INITIALIZE THE ARRAYS INDEX() AND X(). NNLS 300 C NNLS 305 DO 10 I=1,N NNLS 310 X(I) = ZERO NNLS 315 INDEX(I) = I NNLS 320 10 CONTINUE NNLS 325 C NNLS 330 IZ2 = N NNLS 335 IZ1 = 1 NNLS 340 NSETP = 0 NNLS 345 NPP1 = 1 NNLS 350 C ****** MAIN LOOP BEGINS HERE ****** NNLS 355 15 CONTINUE NNLS 360 C QUIT IF ALL COEFFICIENTS ARE ALREADY IN THE SOLUTION.NNLS 365 C OR IF M COLS OF A HAVE BEEN TRIANGULARIZED. NNLS 370 C NNLS 375 IF (IZ1.GT.IZ2 .OR. NSETP.GE.M) GO TO 175 NNLS 380 C NNLS 385 C COMPUTE COMPONENTS OF THE DUAL (NEGATIVE GRADIENT) VECTOR W().NNLS 390 C NNLS 395 DO 25 IZ=IZ1,IZ2 NNLS 400 J = INDEX(IZ) NNLS 405 SM = ZERO NNLS 410 DO 20 L=NPP1,M NNLS 415 SM = SM + A(L,J)*B(L) NNLS 420 20 CONTINUE NNLS 425 W(J) = SM NNLS 430 25 CONTINUE NNLS 435 C FIND LARGEST POSITIVE W(J). NNLS 440 30 WMAX = ZERO NNLS 445 DO 35 IZ=IZ1,IZ2 NNLS 450 J = INDEX(IZ) NNLS 455 IF (W(J).LE.WMAX) GO TO 35 NNLS 460 WMAX = W(J) NNLS 465 IZMAX = IZ NNLS 470 35 CONTINUE NNLS 475 C NNLS 480 C IF WMAX .LE. 0. GO TO TERMINATION. NNLS 485 C THIS INDICATES SATISFACTION OF THE KUHN-TUCKER CONDITIONS.NNLS 490 C NNLS 495 IF (WMAX) 175, 175, 40 NNLS 500 40 IZ = IZMAX NNLS 505 J = INDEX(IZ) NNLS 510 C NNLS 515 C THE SIGN OF W(J) IS OK FOR J TO BE MOVED TO SET P. NNLS 520 C BEGIN THE TRANSFORMATION AND CHECK NEW DIAGONAL ELEMENT TO AVOID NNLS 525 C NEAR LINEAR DEPENDENCE. NNLS 530 C NNLS 535 ASAVE = A(NPP1,J) NNLS 540 CALL H12(1, NPP1, NPP1+1, M, A(1,J), 1, UP, DUMMY, 1, 1, 0) NNLS 545 UNORM = ZERO NNLS 550 IF (NSETP.EQ.0) GO TO 50 NNLS 555 DO 45 L=1,NSETP NNLS 560 UNORM = UNORM + A(L,J)**2 NNLS 565 45 CONTINUE NNLS 570 50 UNORM = SQRT(UNORM) NNLS 575 IF (DIFF(UNORM+ABS(A(NPP1,J))*FACTOR,UNORM)) 65, 65, 55 NNLS 580 C NNLS 585 C COL J IS SUFFICIENTLY INDEPENDENT. COPY B INTO ZZ, UPDATE ZZ AND NNLS 590 C SOLVE FOR ZTEST ( = PROPOSED NEW VALUE FOR X(J) ). NNLS 595 C NNLS 600 55 DO 60 L=1,M NNLS 605 ZZ(L) = B(L) NNLS 610 60 CONTINUE NNLS 615 CALL H12(2, NPP1, NPP1+1, M, A(1,J), 1, UP, ZZ, 1, 1, 1) NNLS 620 ZTEST = ZZ(NPP1)/A(NPP1,J) NNLS 625 C NNLS 630 C SEE IF ZTEST IS POSITIVE NNLS 635 C NNLS 640 IF (ZTEST) 65, 65, 70 NNLS 645 C NNLS 650 C REJECT J AS A CANDIDATE TO BE MOVED FROM SET Z TO SET P. NNLS 655 C RESTORE A(NPP1,J), SETW(J)=0., AND LOOP BACK TO TEST DUAL NNLS 660 C COEFFS AGAIN. NNLS 665 C NNLS 670 65 A(NPP1,J) = ASAVE NNLS 675 W(J) = ZERO NNLS 680 GO TO 30 NNLS 685 C NNLS 690 C THE INDEX J=INDEX(IZ) HAS BEEN SELECTED TO BE MOVED FROM NNLS 695 C SET Z TO SET P. UPDATE B, UPDATE INDICES. APPLY HOUSEHOLDER NNLS 700 C TRANSFORMATIONS TO COLS IN NEW SET Z. ZERO SUBDIAGONAL ELTS IN NNLS 705 C COL J. SET W(J)=0. NNLS 710 C NNLS 715 70 DO 75 L=1,M NNLS 720 B(L) = ZZ(L) NNLS 725 75 CONTINUE NNLS 730 C NNLS 735 INDEX(IZ) = INDEX(IZ1) NNLS 740 INDEX(IZ1) = J NNLS 745 IZ1 = IZ1 + 1 NNLS 750 NSETP = NPP1 NNLS 755 NPP1 = NPP1 + 1 NNLS 760 C NNLS 765 IF (IZ1.GT.IZ2) GO TO 85 NNLS 770 DO 80 JZ=IZ1,IZ2 NNLS 775 JJ = INDEX(JZ) NNLS 780 CALL H12(2, NSETP, NPP1, M, A(1,J), 1, UP, A(1,JJ), 1, MDA, 1) NNLS 785 80 CONTINUE NNLS 790 85 CONTINUE NNLS 795 C NNLS 800 IF (NSETP.EQ.M) GO TO 95 NNLS 805 DO 90 L=NPP1,M NNLS 810 A(L,J) = ZERO NNLS 815 90 CONTINUE NNLS 820 95 CONTINUE NNLS 825 C NNLS 830 W(J) = ZERO NNLS 835 C SOLVE THE TRIANGULAR SYSTEM. NNLS 840 C STORE THE SOLUTION TEMPORARILY IN ZZ() NNLS 845 ASSIGN 100 TO NEXT NNLS 850 GO TO 200 NNLS 855 100 CONTINUE NNLS 860 C NNLS 865 C ****** SECONDARY LOOP BEGINS HERE ****** NNLS 870 C NNLS 875 C INTERATION COUNTER. NNLS 880 C NNLS 885 105 ITER = ITER + 1 NNLS 890 IF (ITER.LE.ITMAX) GO TO 110 NNLS 895 MODE = 3 NNLS 900 WRITE (LOUT,99999) NNLS 905 GO TO 175 NNLS 910 110 CONTINUE NNLS 915 C NNLS 920 C SEE IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE. NNLS 925 C IF NOT COMPUTE ALPHA. NNLS 930 C NNLS 935 ALPHA = TWO NNLS 940 DO 120 IP=1,NSETP NNLS 945 L = INDEX(IP) NNLS 950 IF (ZZ(IP)) 115, 115, 120 NNLS 955 C NNLS 960 115 T = -X(L)/(ZZ(IP)-X(L)) NNLS 965 IF (ALPHA.LE.T) GO TO 120 NNLS 970 ALPHA = T NNLS 975 JJ = IP NNLS 980 120 CONTINUE NNLS 985 C NNLS 990 C IF ALL NEW CONSTRAINED COEFFS ARE FEASIBLE THEN ALPHA WILL NNLS 995 C STILL = 2. IF SO EXIT FROM SECONDARY LOOP. TO MAIN LOOP. NNLS1000 C NNLS1005 IF (ALPHA.EQ.TWO) GO TO 165 NNLS1010 C NNLS1015 C OTHERWISE USE ALPHA WHICH WILL BE BETWEEN 0. AND 1. TO NNLS1020 C INTERPOLATE BETWEEN THE OLD X AND THE NEW ZZ. NNLS1025 C NNLS1030 DO 125 IP=1,NSETP NNLS1035 L = INDEX(IP) NNLS1040 X(L) = X(L) + ALPHA*(ZZ(IP)-X(L)) NNLS1045 125 CONTINUE NNLS1050 C NNLS1055 C MODIFY A AND B AND THE INDEX ARRAYS TO MOVE COEFFICIENT I NNLS1060 C FROM SET P TO SET Z. NNLS1065 C NNLS1070 I = INDEX(JJ) NNLS1075 130 X(I) = ZERO NNLS1080 C NNLS1085 IF (JJ.EQ.NSETP) GO TO 145 NNLS1090 JJ = JJ + 1 NNLS1095 DO 140 J=JJ,NSETP NNLS1100 II = INDEX(J) NNLS1105 INDEX(J-1) = II NNLS1110 CALL G1(A(J-1,II), A(J,II), CC, SS, A(J-1,II)) NNLS1115 A(J,II) = ZERO NNLS1120 DO 135 L=1,N NNLS1125 IF (L.NE.II) CALL G2(CC, SS, A(J-1,L), A(J,L)) NNLS1130 135 CONTINUE NNLS1135 CALL G2(CC, SS, B(J-1), B(J)) NNLS1140 140 CONTINUE NNLS1145 145 NPP1 = NSETP NNLS1150 NSETP = NSETP - 1 NNLS1155 IZ1 = IZ1 - 1 NNLS1160 INDEX(IZ1) = I NNLS1165 C NNLS1170 C SEE IF THE REMAINING COEFFS IN SET P ARE FEASIBLE. THEY SHOULDNNLS1175 C BE BECAUSE OF THE WAY ALPHA WAS DETERMINED. NNLS1180 C IF ANY ARE INFEASIBLE IT IS DUE TO ROUND-OFF ERROR. ANY NNLS1185 C THAT ARE NONPOSITIVE WILL BE SET TO ZERO NNLS1190 C AND MOVED FROM SET P TO SET Z. NNLS1195 C NNLS1200 DO 150 JJ=1,NSETP NNLS1205 I = INDEX(JJ) NNLS1210 IF (X(I)) 130, 130, 150 NNLS1215 150 CONTINUE NNLS1220 C NNLS1225 C COPY B( ) INTO ZZ( ). THEN SOLVE AGAIN AND LOOP BACK. NNLS1230 C NNLS1235 DO 155 I=1,M NNLS1240 ZZ(I) = B(I) NNLS1245 155 CONTINUE NNLS1250 ASSIGN 160 TO NEXT NNLS1255 GO TO 200 NNLS1260 160 CONTINUE NNLS1265 GO TO 105 NNLS1270 C ****** END OF SECONDARY LOOP ****** NNLS1275 C NNLS1280 165 DO 170 IP=1,NSETP NNLS1285 I = INDEX(IP) NNLS1290 X(I) = ZZ(IP) NNLS1295 170 CONTINUE NNLS1300 C ALL NEW COEFFS ARE POSITIVE. LOOP BACK TO BEGINNING. NNLS1305 GO TO 15 NNLS1310 C NNLS1315 C ****** END OF MAIN LOOP ****** NNLS1320 C NNLS1325 C COME TO HERE FOR TERMINATION. NNLS1330 C COMPUTE THE NORM OF THE FINAL RESIDUAL VECTOR. NNLS1335 C NNLS1340 175 SM = ZERO NNLS1345 IF (NPP1.GT.M) GO TO 185 NNLS1350 DO 180 I=NPP1,M NNLS1355 SM = SM + B(I)**2 NNLS1360 180 CONTINUE NNLS1365 GO TO 195 NNLS1370 185 DO 190 J=1,N NNLS1375 W(J) = ZERO NNLS1380 190 CONTINUE NNLS1385 195 RNORM = SQRT(SM) NNLS1390 RETURN NNLS1395 C NNLS1400 C THE FOLLOWING BLOCK OF CODE IS USED AS AN INTERNAL SUBROUTINE NNLS1405 C TO SOLVE THE TRIANGULAR SYSTEM. PUTTING THE SOLUTION IN ZZ(). NNLS1410 C NNLS1415 200 DO 215 L=1,NSETP NNLS1420 IP = NSETP + 1 - L NNLS1425 IF (L.EQ.1) GO TO 210 NNLS1430 DO 205 II=1,IP NNLS1435 ZZ(II) = ZZ(II) - A(II,JJ)*ZZ(IP+1) NNLS1440 205 CONTINUE NNLS1445 210 JJ = INDEX(IP) NNLS1450 ZZ(IP) = ZZ(IP)/A(IP,JJ) NNLS1455 215 CONTINUE NNLS1460 GO TO NEXT, (100, 160) NNLS1465 99999 FORMAT (//46H DIAGNOSTIC: NNLS HAS EXHAUSTED THE SPECIFIED , NNLS1470 * 21HNUMBER OF ITERATIONS.//) NNLS1475 END NNLS1480 CC$ FORTY FORM DIFF 5 FUNCTION DIFF(X, Y) DIFF 10 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUNE 7 DIFF 15 C TO APPEAR IN "SOLVING LEAST SQUARES PROBLEMS", PRENTICE-HALL 1974 DIFF 20 DIFF = X - Y DIFF 25 RETURN DIFF 30 END DIFF 35 CC$ FORTY FORM,XREF H120 5 C SUBROUTINE H12 (MODE,LPIVOT,L1,M,U,IUE,UP,C,ICE,ICV,NCV) H120 10 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABOROTATORY,1973 JUN 12H120 15 C TAKEN FROM "SOLVING LEAST SQUARES PROBLEMS", PRENTICE-HALL, 1974 H120 20 C H120 25 C CONSTRUCTION AND/OR APPLICATION OF A SINGLE H120 30 C HOUSEHOLDER TRANSFORMATION.. Q = I + U*(U**T)/B H120 35 C H120 40 C MODE = 1 OR 2 TO SELECT ALGORITHM H1 OR H2 H120 45 C LPIVOT IS THE INDEX OF THE PIVOT ELEMENT. H120 50 C L1.M IF L1 .LE. M THE TRANSFORMATION WILL BE CONSTRUCTED TO H120 55 C ZERO ELEMENTS INDEXED FROM L1 THROUGH M. IF L1 GT. M H120 60 C THE SUBROUTINE DOES AN IDENTITY TRANSFORMATION. H120 65 C U(),IUE,UP ON ENTRY TO H1 U() CONTAINS THE PIVOT VECTOR. H120 70 C IUE IS THE STORAGE INCREMENT BETWEEN ELEMENTS. H120 75 C ON EXIT FROM H1 U() AND UP H120 80 C CONTAIN QUANTITIES DEFINING THE VECTOR U OF THE H120 85 C HOUSEHOLDER TRANSFORMATION. ON ENTRY TO H2 U() H120 90 C AND UP SHOULD CONTAIN QUANTITIES PREVIOUSLY COMPUTEDH120 95 C BY H1. THESE WILL NOT BE MODIFIED BY H2. H120 100 C C() ON ENTRY TO H1 OR H2 C() CONTAINS A MATRIX WHICH WILL BE H120 105 C REGARDED AS A SET OF VECTORS TO WHICH THE HOUSEHOLDER H120 110 C TRANSFORMATION IS TO BE APPLIED. ON EXIT C() CONTAINS THE H120 115 C SET OF TRANSFORMED VECTORS. H120 120 C ICE STORAGE INCREMENT BETWEEN ELEMENTS OF VECTORS IN C(). H120 125 C ICV STORAGE INCREMENT BETWEEN VECTORS. IN C(). H120 130 C NCV NUMBER OF VECTORS IN C() TO BE TRANSFORMED. IF NCV .LE. 0 H120 135 C NO OPERATIONS WILL BE DONE ON C(). H120 140 C H120 145 SUBROUTINE H12(MODE, LPIVOT, L1, M, U, IUE, UP, C, ICE, ICV, NCV) H120 150 DIMENSION U(IUE,M), C(1) H120 155 DOUBLE PRECISION SM, B H120 160 ONE = 1. H120 165 C H120 170 IF (0.GE.LPIVOT .OR. LPIVOT.GE.L1 .OR. L1.GT.M) RETURN H120 175 CL = ABS(U(1,LPIVOT)) H120 180 IF (MODE.EQ.2) GO TO 30 H120 185 C ****** CONSTRUCT THE TRANSFORMATION ****** H120 190 DO 5 J=L1,M H120 195 CL = AMAX1(ABS(U(1,J)),CL) H120 200 5 CONTINUE H120 205 IF (CL) 65, 65, 10 H120 210 10 CLINV = ONE/CL H120 215 SM = (DBLE(U(1,LPIVOT))*CLINV)**2 H120 220 DO 15 J=L1,M H120 225 SM = SM + (DBLE(U(1,J))*CLINV)**2 H120 230 15 CONTINUE H120 235 C CONVERT DBLE. PREC. SM TO SNGL. PREC. SM1H120 240 SM1 = SM H120 245 CL = CL*SQRT(SM1) H120 250 IF (U(1,LPIVOT)) 25, 25, 20 H120 255 20 CL = -CL H120 260 25 UP = U(1,LPIVOT) - CL H120 265 U(1,LPIVOT) = CL H120 270 GO TO 35 H120 275 C ****** APPLY THE TRANSFORMATION I+U*(U**T)/B TO C ******H120 280 C H120 285 30 IF (CL) 65, 65, 35 H120 290 35 IF (NCV.LE.0) RETURN H120 295 B = DBLE(UP)*U(1,LPIVOT) H120 300 C B MUST BE NONPOSITIVE HERE. IF B = 0., RETURN H120 305 C H120 310 IF (B) 40, 65, 65 H120 315 40 B = ONE/B H120 320 I2 = 1 - ICV + ICE*(LPIVOT-1) H120 325 INCR = ICE*(L1-LPIVOT) H120 330 DO 60 J=1,NCV H120 335 I2 = I2 + ICV H120 340 I3 = I2 + INCR H120 345 I4 = I3 H120 350 SM = C(I2)*DBLE(UP) H120 355 DO 45 I=L1,M H120 360 SM = SM + C(I3)*DBLE(U(1,I)) H120 365 I3 = I3 + ICE H120 370 45 CONTINUE H120 375 IF (SM) 50, 60, 50 H120 380 50 SM = SM*B H120 385 C(I2) = C(I2) + SM*DBLE(UP) H120 390 DO 55 I=L1,M H120 395 C(I4) = C(I4) + SM*DBLE(U(1,I)) H120 400 I4 = I4 + ICE H120 405 55 CONTINUE H120 410 60 CONTINUE H120 415 65 RETURN H120 420 END H120 425 CC$ FORTY FORM,XREF G100 5 SUBROUTINE G1(A, B, COS, SIN, SIG) G100 10 C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 G100 15 C TAKEN FROM "SOLVING LEAST SQUARES PROBLEMS". PRENTICE-HALL, 1974 G100 20 C G100 25 C G100 30 C COMPUTE ORTHOGONAL ROTATION MATRIX. G100 35 C COMPUTE.. MATRIX (C, S) SO THAT (C, S)(A) = (SQRT(A**2+B**2)) G100 40 C (-S,C) (-S,C)(B) ( 0 ) G100 45 C COMPUTE SIG = SQRT(A**2+B**2) G100 50 C SIG IS COMPUTED LAST TO ALLOW FOR THE POSSIBILITY THAT G100 55 C SIG MAY BE IN THE SAME LOCATION AS A OR B. G100 60 C G100 65 ZERO = 0. G100 70 ONE = 1. G100 75 IF (ABS(A).LE.ABS(B)) GO TO 5 G100 80 XR = B/A G100 85 YR = SQRT(ONE+XR**2) G100 90 COS = SIGN(ONE/YR,A) G100 95 SIN = COS*XR G100 100 SIG = ABS(A)*YR G100 105 RETURN G100 110 5 IF (B) 10, 15, 10 G100 115 10 XR = A/B G100 120 YR = SQRT(ONE+XR**2) G100 125 SIN = SIGN(ONE/YR,B) G100 130 COS = SIN*XR G100 135 SIG = ABS(B)*YR G100 140 RETURN G100 145 15 SIG = ZERO G100 150 COS = ZERO G100 155 SIN = ONE G100 160 RETURN G100 165 END G100 170 SUBROUTINE G2(COS, SIN, X, Y) G200 5 C C.L.LAWSON AND R.J.HANSON, JET POPULSION LABORATORY, 1972 DEC 15 G200 10 C TAKEN FROM "SOLVING LEAST SQUARES PROBLEMS", PRENTICE-HALL, 1974 G200 15 C APPLY THE ROTATION COMPUTED BY G1 TO (X,Y). G200 20 XR = COS*X + SIN*Y G200 25 Y = -SIN*X + COS*Y G200 30 X = XR G200 35 RETURN G200 40 END G200 45 CC$ 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 CC$ 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 CC$ 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 CC$ 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 EPS = 0.0 SSP2 145 DO 15 L=1,N SSP2 150 DO 10 M=1,N SSP2 155 EPS = EPS + A(L,M)**2 SSP2 160 10 CONTINUE SSP2 165 15 CONTINUE SSP2 170 EPS = SQRT(EPS)*1.0E-20 SSP2 175 99999 FORMAT (30H IN SSP2 JUST BEFORE 15, EPS =, E20.10, 9H A(I,I) =, SSP2 180 * E20.10) SSP2 185 IF (ABS(A(I,I)).GE.EPS) GO TO 20 SSP2 190 WRITE (6,99999) EPS, A(I,I) SSP2 195 FAIL = .TRUE. SSP2 200 RETURN SSP2 205 20 B(I) = SAVE/A(I,I) SSP2 210 IF (SAVE.EQ.0.0) GO TO 30 SSP2 215 DO 25 J=IP1,N SSP2 220 B(J) = B(J) + A(J,I)*SAVE SSP2 225 25 CONTINUE SSP2 230 30 I = IP1 SSP2 235 GO TO 5 SSP2 240 35 TEMP = B(I) SSP2 245 B(ICH) = B(IP1) SSP2 250 DENOM = A(I,I)*A(IP1,IP1)/A(IP1,I) - A(IP1,I) SSP2 255 B(IP1) = (SAVE*A(I,I)/A(IP1,I)-TEMP)/DENOM SSP2 260 B(I) = (SAVE-A(IP1,IP1)*B(IP1))/A(IP1,I) SSP2 265 IP2 = I + 2 SSP2 270 IF (IP2.GT.N) GO TO 45 SSP2 275 DO 40 J=IP2,N SSP2 280 B(J) = B(J) + A(J,I)*TEMP + A(J,IP1)*SAVE SSP2 285 40 CONTINUE SSP2 290 45 I = IP2 SSP2 295 GO TO 5 SSP2 300 C SSP2 305 C SSP2 310 50 IF (I.NE.N) GO TO 60 SSP2 315 99998 FORMAT (23H JUST BEFORE 105, EPS =, E20.10, 10H A(I,I) = , E20.10)SSP2 320 IF (ABS(A(I,I)).GE.EPS) GO TO 55 SSP2 325 WRITE (6,99998) EPS, A(I,I) SSP2 330 FAIL = .TRUE. SSP2 335 RETURN SSP2 340 55 B(I) = B(I)/A(I,I) SSP2 345 C NOW SOLVE M(TRANSPOSE) X = Y FOR X, WHERE Y IS STORED SSP2 350 C IN THE VECTOR B, AND STORE X IN B SSP2 355 C SSP2 360 60 I = N - 1 SSP2 365 IF (INTER(N).LT.0) I = N - 2 SSP2 370 65 IF (I.LE.0) RETURN SSP2 375 I1 = I SSP2 380 IF (INTER(I).LT.0) I1 = I - 1 SSP2 385 IP1 = I + 1 SSP2 390 DO 75 K=I1,I SSP2 395 SAVE = B(K) SSP2 400 DO 70 J=IP1,N SSP2 405 SAVE = SAVE + A(J,K)*B(J) SSP2 410 70 CONTINUE SSP2 415 B(K) = SAVE SSP2 420 75 CONTINUE SSP2 425 ICH = INTER(I1) SSP2 430 B(I) = B(ICH) SSP2 435 B(ICH) = SAVE SSP2 440 I = I1 - 1 SSP2 445 GO TO 65 SSP2 450 END SSP2 455 CC$ FORTY FORM ASUM 5 FUNCTION ASUM(IW) ASUM 10 DOUBLE PRECISION DCON ASUM 15 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, ASUM 20 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, ASUM 25 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, ASUM 30 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON ASUM 35 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), ASUM 40 * ALOS, IXV ASUM 45 C ASUM 50 C COMPUTE THE A PART OF THE LOSS FUNCTION, WHICH IS A CONSTANT ASUM 55 C FOR ALL I. ASUM 60 C ASUM 65 ASUM = 0.0 ASUM 70 DO 15 M=1,M3 ASUM 75 WIW = RISE(M,IW) ASUM 80 DO 10 I=2,N ASUM 85 PI = P(I,IW) ASUM 90 I1 = I - 1 ASUM 95 DO 5 J=1,I1 ASUM 100 TEST1 = DATA(M,J,I) - WIW*PI*P(J,IW) ASUM 105 IF (ABS(TEST1).LE.1.0E-12) GO TO 5 ASUM 110 ASUM = ASUM + TEST1**2 ASUM 115 5 CONTINUE ASUM 120 10 CONTINUE ASUM 125 15 CONTINUE ASUM 130 ASUM = ASUM*DCON*RM ASUM 135 RETURN ASUM 140 END ASUM 145 CC$ FORTY FORM,XREF BSUM 5 FUNCTION BSUM(IW) BSUM 10 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, BSUM 15 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, BSUM 20 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, BSUM 25 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON BSUM 30 COMMON /A4/ DAPI(030), BCONST, TM, U, V BSUM 35 C BSUM 40 C COMPUTE THE B PART OF THE LOSS FUNCTION, WHICH IS A CONSTANT BSUM 45 C FOR ALL I. BSUM 50 C BSUM 55 TM = 0.0 BSUM 60 DO 10 I=2,N BSUM 65 PI = P(I,IW) BSUM 70 I1 = I - 1 BSUM 75 DO 5 J=1,I1 BSUM 80 TM = TM + PI*P(J,IW) BSUM 85 5 CONTINUE BSUM 90 10 CONTINUE BSUM 95 TM = TM*ABNV BSUM 100 C BSUM 105 C TM IS THE MEAN (P(I,IW)*P(J,IW)) I .GT. J BSUM 110 C BSUM 115 U = 0.0 BSUM 120 V = 0.0 BSUM 125 DO 25 I=2,N BSUM 130 PI = P(I,IW) BSUM 135 I1 = I - 1 BSUM 140 DO 20 J=1,I1 BSUM 145 PIPJ = PI*P(J,IW) BSUM 150 TEST1 = (PIPJ-1.0)*PIPJ BSUM 155 IF (ABS(TEST1).LE.1.0E-12) GO TO 15 BSUM 160 U = U + TEST1**2 BSUM 165 15 CONTINUE BSUM 170 TEST2 = PIPJ - TM BSUM 175 IF (ABS(TEST2).LE.1.0E-12) GO TO 20 BSUM 180 V = V + TEST2**2 BSUM 185 20 CONTINUE BSUM 190 25 CONTINUE BSUM 195 C NOW ADD THE DIAGONAL ELEMENTS TO U BSUM 200 DO 30 I=1,N BSUM 205 PI2 = P(I,IW)**2 BSUM 210 TEST3 = (PI2-1.0)*PI2 BSUM 215 IF (ABS(TEST3).LE.1.0E-12) GO TO 30 BSUM 220 U = U + 0.5*TEST3**2 BSUM 225 30 CONTINUE BSUM 230 BSUM = U/V BSUM 235 RETURN BSUM 240 END BSUM 245 CC$ FORTY FORM,XREF VARI 5 SUBROUTINE VARIAN(MST, IND, IW) VARI 10 DOUBLE PRECISION DCON VARI 15 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, VARI 20 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, VARI 25 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, VARI 30 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON VARI 35 COMMON /A2/ G(018,030), KT1(10), KT2(10), TK(10), LP(15), IPUN, VARI 40 * IERR, LPUNCH, LPUN7, ICOUN, TVAF, ICON, TITL(72), VARI 45 * TRISE(018,017), E(018) VARI 50 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), VARI 55 * ALOS, IXV VARI 60 COMMON /A8/ W(018,017), IPRE, NNEG, NDOT, NLD1 VARI 65 C VARI 70 C VARI 75 C VARI 80 C VARI 85 IF (IND.NE.1) GO TO 5 VARI 90 VAF = 1.0 - 4.0*ALOS*AMV VARI 95 WRITE (LOUT,99999) MST, VAF, IW VARI 100 99999 FORMAT (/17H DURING MAJOR IT., I4, 25H VARIANCE ACCOUNTED FOR =, VARI 105 * F12.6, 29H IN THE RESIDUALS BY SUBSET (, I2, 16H), WITHOUT POLIS,VARI 110 * 4HHING) VARI 115 RETURN VARI 120 C VARI 125 C IF REGRESSION HAS JUST FAILED, THE FOLLOWING DEFINITION VARI 130 C IS INVALID. VARI 135 C VARI 140 5 SSQ = 0.0 VARI 145 DO 25 M=1,M3 VARI 150 DM = 0.0 VARI 152 IF (MCON.EQ.0) DM = DMEAN(M) VARI 154 CONS = RISE(M,LD2) VARI 155 DO 20 I=2,N VARI 160 I1 = I - 1 VARI 165 DO 15 J=1,I1 VARI 170 DATAJI = 0.0 VARI 175 DO 10 K=1,LD1 VARI 180 DATAJI = DATAJI + P(I,K)*P(J,K)*RISE(M,K) VARI 185 10 CONTINUE VARI 190 SSQ = SSQ + (DATA(M,I,J) + DM - CONS - DATAJI)**2 VARI 195 15 CONTINUE VARI 200 20 CONTINUE VARI 205 25 CONTINUE VARI 210 VAF = (VARSIM-SSQ)/VARSIM VARI 215 TVAF = VAF VARI 220 IF (IXV.EQ.1) GO TO 35 VARI 225 IF (MST.LE.0) GO TO 30 VARI 230 IF (KP(MST).NE.1 .AND. KC.NE.1) RETURN VARI 235 30 WRITE (LOUT,99998) MST, VAF VARI 240 99998 FORMAT (/24H DURING MAJOR ITERATION , I4, 19H, OVERALL VARIANCE , VARI 245 * 15HACCOUNTED FOR =, F12.6, 31H AFTER USING REGRESSION WEIGHTS) VARI 250 RETURN VARI 255 35 IF (NNEG.EQ.0) WRITE (LOUT,99997) MST, VAF VARI 260 99997 FORMAT (17H DURING ITERATION, I4, 19H, OVERALL VARIANCE , VARI 265 * 31H--TENTATIVELY-- ACCOUNTED FOR =, F12.6, 18H AFTER USING REGRE,VARI 270 * 13HSSION WEIGHTS) VARI 275 RETURN VARI 280 END VARI 285 CC$ FORTY FORM POLI 5 SUBROUTINE POLISH(MST, IW) POLI 10 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, POLI 15 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, POLI 20 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, POLI 25 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON POLI 30 COMMON /A4/ DAPI(030), BCONST, TM, U, V POLI 35 COMMON /A7/ A(030), IZ(030), ADUM(017), IZDUM(017) POLI 40 C POLI 45 C SUBROUTINE TO FORCE P(,IW) TO BE ZEROES AND ONES. POLI 50 C POLI 55 C POLI 60 M1 = 0 POLI 65 DO 10 KLM=1,N POLI 70 A(KLM) = P(KLM,IW) POLI 75 IF (A(KLM).GE.0.25 .AND. A(KLM).LE.0.75) MP(MST) = MP(MST) + 1 POLI 80 IF (A(KLM).GE.0.707) GO TO 5 POLI 85 P(KLM,IW) = 0.0 POLI 90 GO TO 10 POLI 95 5 M1 = M1 + 1 POLI 100 P(KLM,IW) = 1.0 POLI 105 10 CONTINUE POLI 110 IF (M1.LE.1 .OR. M1.EQ.N) GO TO 15 POLI 115 RETURN POLI 120 C POLI 125 C IF A VECTOR OF ALL ONES OR ZEROES HAS BEEN GENERATED, PRINT POLI 130 C A DIAGNOSTIC. POLI 135 C POLI 140 15 IF (KP(MST).EQ.1) WRITE (LOUT,99999) IW, M1 POLI 145 99999 FORMAT (/49H ******** WARNING FROM SUBROUTINE POLISH ********/ POLI 150 * 7H SUBSET, I3, 24H IS TRIVIAL SINCE IT HAS, I3, 7H ONE(S), POLI 155 * 39H ON THE PRESENT POLISHING; FIX-UP TAKEN//) POLI 160 DO 20 KLM=1,N POLI 165 IZ(KLM) = KLM POLI 170 20 CONTINUE POLI 175 CALL SORT(N,A,IZ) POLI 180 IF (M1.EQ.N) GO TO 25 POLI 185 C POLI 190 C FOR A NULL SUBSET, CHANGE THE TWO LARGEST ENTRIES TO UNITIES. POLI 195 C POLI 200 IZ1 = IZ(1) POLI 205 P(IZ1,IW) = 1.0 POLI 210 IZ2 = IZ(2) POLI 215 P(IZ2,IW) = 1.0 POLI 220 RETURN POLI 225 C POLI 230 C OR IF THE SUBSET CONTAINS THE WHOLE SET OF STIMULI, CHANGE THE POLI 235 C SMALLEST ENTRY IN THE VECTOR TO ZERO. POLI 240 C POLI 245 25 IZN = IZ(N) POLI 250 P(IZN,IW) = 0.0 POLI 255 RETURN POLI 260 END POLI 265 CC$ FORTY FORM SORT 5 SUBROUTINE SORT(N,A,IZ) SORT 10 DIMENSION A(N),IZ(N) 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 CC$ 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,017), DENS(050), NW(050), MP(050), M3, RM, AMV, INIC 25 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, INIC 30 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, INIC 35 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON INIC 40 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), INIC 45 * ALOS, IXV INIC 50 COMMON /A8/ W(018,017), IPRE, NNEG, NDOT, NLD1 INIC 55 C INIC 60 C INIC 65 IF (ICON) 90, 5, 65 INIC 70 C INIC 75 C FOR A RATIONAL INITIAL CONFIGURATION INIC 80 C INIC 85 5 NP = 0 INIC 90 NN = 0 INIC 95 SUM = 0.0 INIC 100 SUMP = 0.0 INIC 105 SUMN = 0.0 INIC 110 DO 20 M=1,M3 INIC 115 DO 15 I=1,N INIC 120 U(I) = 0.0 INIC 125 DO 10 J=1,N INIC 130 IF (I.EQ.J) GO TO 10 INIC 135 K1 = MIN0(I,J) INIC 140 K2 = MAX0(I,J) INIC 145 U(I) = U(I) + DATA(M,K1,K2) INIC 150 10 CONTINUE INIC 155 15 CONTINUE INIC 160 20 CONTINUE INIC 165 DO 25 J=1,N INIC 170 SUM = SUM + U(J) INIC 175 25 CONTINUE INIC 180 SUM = SUM*BIN INIC 185 DO 40 J=1,N INIC 190 U(J) = U(J) - SUM INIC 195 IF (U(J)) 30, 40, 35 INIC 200 30 NN = NN + 1 INIC 205 SUMN = SUMN + U(J) INIC 210 GO TO 40 INIC 215 35 NP = NP + 1 INIC 220 SUMP = SUMP + U(J) INIC 225 40 CONTINUE INIC 230 IF (NN*NP.EQ.0) GO TO 50 INIC 235 SUMP = SUMP/FLOAT(NP) INIC 240 SUMN = SUMN/FLOAT(NN) INIC 245 BB = 1.0/(SUMP-SUMN) INIC 250 A = -SUMN/(SUMP-SUMN) INIC 255 C WRITE (LOUT,200) NN,NP,SUMP,SUMN,A,BB,(U(K),K=1, N) INIC 260 C 200 FORMAT(' NN =',I3,' NP =',I3,' SUMP =',F12.6,' SUMN =',F12.6, INIC 265 C 2' A =',F12.6,' BB =',F12.6/(1X,10F12.6/)) INIC 270 DO 45 J=1,N INIC 275 P(J,IW) = A + BB*U(J) INIC 280 45 CONTINUE INIC 285 IF (MST.EQ.0) RETURN INIC 290 IF (MST.EQ.1 .OR. KP(MST).EQ.1) WRITE (LOUT,99999) IW, INIC 295 * (P(KLM,IW),KLM=1,N) INIC 300 99999 FORMAT (/53H RATIONAL INITIAL CONFIGURATION P VECTOR FOR SUBSET ,INIC 305 * 1H(, I2, 2H):/64(1X, 16F8.4/)) INIC 310 RETURN INIC 315 C INIC 320 C INIC 325 C ISSUE DIAGNOSTIC IF NN (NUMBER OF NEG. U(J)) OR NP (POS.) INIC 330 C IS ZERO. THAT IS, DATA(M,J,I) WAS A DOUBLE-CENTERED MATRIX. INIC 335 C INIC 340 50 WRITE (LOUT,99998) NN, NP INIC 345 99998 FORMAT (///50H ********DIAGNOSTIC FROM SUBROUTINE INICON********//INIC 350 * 19H NO. OF NEG U(J) =, I3, 15H AND POS U(J) =, I3//20X, INIC 355 * 30HLISTING OF MATRIX OF RESIDUALS//) INIC 360 DO 60 M=1,M3 INIC 365 DO 55 I=2,N INIC 370 I1 = I - 1 INIC 375 WRITE (LOUT,99997) (M,I,J,DATA(M,J,I),J=1,I1) INIC 380 55 CONTINUE INIC 385 60 CONTINUE INIC 390 99997 FORMAT (2H (, 3I3, 1H), F9.4, 2H (, 3I3, 1H), F9.4, 2H (, 3I3, INIC 395 * 1H), F9.4, 2H (, 3I3, 1H), F9.4, 2H (, 3I3, 1H), F9.4/) INIC 400 STOP INIC 405 C INIC 410 C INIC 415 C INIC 420 C IF (ICON .GT.0) READ IN A USER-SUPPLIED INITIAL CONFIGURATION FOR INIC 425 C P(,IW) VIA CHANNEL ICON. INIC 430 C INIC 435 65 WRITE (LOUT,99996) LD1, ICON INIC 440 99996 FORMAT (/////49H USER HAS SPECIFIED THAT INITIAL CONFIGURATION OF,INIC 445 * I4, 44H CLUSTERS IS TO BE READ FROM LOGICAL DEVICE , I3/6H WITH ,INIC 450 * 32HTHE TITLE AND FORMAT AS FOLLOWS:///) INIC 455 C INIC 460 C READ THE TITLE CARD AND PRINT IT; THEN DO THE SAME FOR THE FORMAT.INIC 465 C INIC 470 READ (ICON,99995) FMT2 INIC 475 WRITE (LOUT,99994) FMT2 INIC 480 99995 FORMAT (72A1) INIC 485 99994 FORMAT (//1X, 46HTITLE FOR USER-SUPPLIED INITIAL CONFIGURATION: , INIC 490 * 72A1//) INIC 495 99993 FORMAT (//1X, 33HUSER-SPECIFIED FORMAT OF INITIAL , 10HCONFIGURAT,INIC 500 * 8HION IS: , 72A1//) INIC 505 READ (ICON,99995) FMT2 INIC 510 WRITE (LOUT,99993) FMT2 INIC 515 WRITE (LOUT,99992) (I, I = 1, N) INIC 520 99992 FORMAT (///64(1X, I6, 15I8/)) INIC 525 DO 70 IX=1,LD1 INIC 530 READ (ICON,FMT2,END=100) (P(KLM,IX), KLM = 1, N) INIC 535 WRITE (LOUT,99991) IX, (P(KLM,IX),KLM=1,N) INIC 540 70 CONTINUE INIC 545 99991 FORMAT (//34H INITIAL CONFIGURATION P VECTOR (, I2, 2H):/64(1X, INIC 550 * 16F8.4/)) INIC 555 C INIC 560 C INIC 565 C THIS CHECK SPOTS SINGLETONS AND LUMPER, AND ALSO LOOKS INIC 570 C TO SEE IF USER-SUPPLIED P MATRIX WERE 0,1. INIC 575 C INIC 580 IDANG = 0 INIC 585 C INIC 590 DO 85 J=1,LD1 INIC 595 IK = 0 INIC 600 DO 80 I=1,N INIC 605 PIJ = P(I,J) INIC 610 IF (PIJ.LT.1.01 .AND. PIJ.GT.0.99) GO TO 75 INIC 615 IF (-1.0E-04.LE.PIJ .AND. PIJ.LT.1.0E-04) GO TO 80 INIC 620 WRITE (LOUT,99990) I, J, P(I,J) INIC 625 99990 FORMAT (/35H WARNING: IN USER-SUPPLIED INITIAL , 13HCONFIGURATION,INIC 630 * 4H, P(, I3, 1H,, I3, 4H) = , F10.5, 25H INSTEAD OF BEING A BINAR,INIC 635 * 14HY COEFFICIENT.) INIC 640 75 IK = IK + 1 INIC 645 80 CONTINUE INIC 650 IF (IK.GT.1 .AND. IK.LT.N) GO TO 85 INIC 655 IDANG = 1 INIC 660 WRITE (LOUT,99989) J, IK INIC 665 99989 FORMAT (///17H ******** SUBSET , I3, 24H HAS TOO FEW (0 OR 1) OR, INIC 670 * 30H TOO MANY (N) ELEMENTS, NAMELY, I3, 15H AND WILL CAUSE/10X, INIC 675 * 47HTHE MULTIPLE LINEAR REGRESSION TO FAIL ********//) INIC 680 85 CONTINUE INIC 685 IF (IDANG.NE.0) STOP INIC 690 RETURN INIC 695 100 WRITE (LOUT,99988) ICON INIC 700 99988 FORMAT (///49H ******** READ AN END OF FILE WHILE ATTEMPTING TO, INIC 705 * 26H READ INIT. CONF. FOR T(J)/10X, 26HIN SUBROUTINE INICON, USIN,INIC 710 * 1HG, 13H CHANNEL NO. , I3, 9H ********) INIC 715 RETURN INIC 720 C INIC 725 C FOR A RANDOM INITIAL CONFIGURATION INIC 730 C INIC 735 90 DO 95 KLM=1,N INIC 740 P(KLM,IW) = RUNIFV(0) INIC 745 95 CONTINUE INIC 750 WRITE (LOUT,99987) IW, (P(KLM,IW),KLM=1,N) INIC 755 99987 FORMAT (//33H RANDOM INITIAL CONFIGURATION P (, I2, 2H):/64(1X, INIC 760 * 16F8.4/)) INIC 765 RETURN INIC 770 END INIC 775 CC$ FORTY FORM RUNI 5 FUNCTION RUNIFV(L) RUNI 10 DATA MODULO, FLMOD, K /8192,8192.0,1/ RUNI 15 C RANDOM NUMBERS MODULO 2**13 RUNI 20 C BASED ON A MODIFICATION OF A ROUTINE SUGGESTED BY J. B. KRUSKAL RUNI 25 C (CACM, 1969, 12, 93-94), THIS SUBROUTINE IS SIMILAR TO THE ONE RUNI 30 C USED IN THE MULTIDIMENSIONAL SCALING PROGRAM KYST. RUNI 35 C RUNI 40 DO 5 I=1,3 RUNI 45 K = MOD(5*K,MODULO) RUNI 50 5 CONTINUE RUNI 55 RUNIFV = FLOAT(K)/FLMOD RUNI 60 RETURN RUNI 65 END RUNI 70 CC$ FORTY FORM,XREF RESI 5 SUBROUTINE RESID(IW) RESI 10 DOUBLE PRECISION DCON, SUMA RESI 15 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, RESI 20 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, RESI 25 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, RESI 30 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON RESI 35 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), RESI 40 * ALOS, IXV RESI 45 C RESI 50 C RESI 55 C IN THE ALTERNATING LEAST SQUARES PROCEDURE, COMPUTE RESI 60 C THE RESIDUALS FOR SUBSET IW, AND CENTER THEM RESI 65 C RESI 70 C RESI 75 DO 15 M=1,M3 RESI 80 DO 10 I=2,N RESI 85 I1 = I - 1 RESI 90 DO 5 J=1,I1 RESI 95 DATA(M,J,I) = 0.0 RESI 100 5 CONTINUE RESI 105 10 CONTINUE RESI 110 15 CONTINUE RESI 115 DO 35 K=1,LD1 RESI 120 IF (K.EQ.IW) GO TO 35 RESI 125 DO 30 M=1,M3 RESI 130 WIW = RISE(M,K) RESI 135 DO 25 I=2,N RESI 140 I1 = I - 1 RESI 145 PIK = P(I,K) RESI 150 DO 20 J=1,I1 RESI 155 DATA(M,J,I) = DATA(M,J,I) + PIK*P(J,K)*WIW RESI 160 20 CONTINUE RESI 165 25 CONTINUE RESI 170 30 CONTINUE RESI 175 35 CONTINUE RESI 180 SUMA = 0.0 RESI 185 DO 60 M=1,M3 RESI 190 SUM = 0.0 RESI 195 DO 45 I=2,N RESI 200 I1 = I - 1 RESI 205 DO 40 J=1,I1 RESI 210 DATA(M,J,I) = DATA(M,I,J) - DATA(M,J,I) RESI 215 SUM = SUM + DATA(M,J,I) RESI 220 40 CONTINUE RESI 225 45 CONTINUE RESI 230 SUM = SUM*ABNV RESI 235 DO 55 I=2,N RESI 240 I1 = I - 1 RESI 245 DO 50 J=1,I1 RESI 250 DATA(M,J,I) = DATA(M,J,I) - SUM RESI 255 C RESI 260 C SINCE THE UPPER TRIANGULAR HALF OF DATA CONSISTS OF CENTERED RESI 265 C RESIDUALS, RESI 270 C THE SUM OF SQUARES EQUALS THE NUMERATOR OF THE VARIANCE. RESI 275 C RESI 280 SUMA = SUMA + DATA(M,J,I)**2 RESI 285 50 CONTINUE RESI 290 55 CONTINUE RESI 295 60 CONTINUE RESI 300 C RESI 305 C COMPUTE THE CONSTANT THAT NORMALIZES THE SUM OF SQUARED ERROR. RESI 310 C RESI 315 DCON = 0.25*ABM/SUMA RESI 320 RETURN RESI 325 END RESI 330 CC$ FORTY FORM ADJU 5 SUBROUTINE ADJUST(MST, IW) ADJU 10 COMMON /A5/ AK, ALFRAY(017), BETRAY(017), 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 CC$ FORTY FORM,XREF CDAP 5 SUBROUTINE CDAPI(IW) CDAP 10 DOUBLE PRECISION DCON CDAP 15 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, CDAP 20 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, CDAP 25 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, CDAP 30 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON CDAP 35 COMMON /A4/ DAPI(030), BCONST, TM, U, V CDAP 40 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), CDAP 45 * ALOS, IXV CDAP 50 C CDAP 55 C COMPUTE THE PARTIAL OF A WITH RESPECT TO P(I,IW) CDAP 60 C CDAP 65 DO 15 I=1,N CDAP 70 PI = P(I,IW) CDAP 75 AS = 0. CDAP 80 DAPI(I) = 0.0 CDAP 85 DO 10 M=1,M3 CDAP 90 WIW = RISE(M,IW) CDAP 95 DO 5 J=1,N CDAP 100 IF (J.EQ.I) GO TO 5 CDAP 105 PJ = P(J,IW) CDAP 110 K1 = MIN0(I,J) CDAP 115 K2 = MAX0(I,J) CDAP 120 AS = AS + PJ*(DATA(M,K1,K2)-WIW*PJ*PI) CDAP 125 5 CONTINUE CDAP 130 DAPI(I) = WIW*AS + DAPI(I) CDAP 135 10 CONTINUE CDAP 140 15 CONTINUE CDAP 145 DO 20 I=1,N CDAP 150 DAPI(I) = -2.0*DAPI(I)*DCON*RM CDAP 155 20 CONTINUE CDAP 160 RETURN CDAP 165 END CDAP 170 CC$ FORTY FORM,XREF DB00 5 FUNCTION DB(KLM, IW) DB00 10 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, DB00 15 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, DB00 20 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, DB00 25 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON DB00 30 COMMON /A4/ DAPI(030), BCONST, TM, U, V DB00 35 C DB00 40 C FUNCTION TO ASSEMBLE ALL THE EXPRESSIONS FOR DB00 45 C EVALUATING DB/DP(I) DB00 50 C DB00 55 B1 = 0.0 DB00 60 B2 = 0.0 DB00 65 DO 5 I=1,N DB00 70 PI = P(I,IW) DB00 75 B1 = B1 + PI DB00 80 B2 = B2 + PI*PI DB00 85 5 CONTINUE DB00 90 EIK = ABNV*(B1-P(KLM,IW)) DB00 95 C DB00 100 C TM IS THE MEAN (P(I,IW)*P(J,IW)) I .GT. J DB00 105 C DB00 110 TM = ABNV*0.5*(B1*B1-B2) DB00 115 DUDPI = 0.0 DB00 120 DVDPI = 0.0 DB00 125 PI = P(KLM,IW) DB00 130 DO 15 J=1,N DB00 135 PJ = P(J,IW) DB00 140 PIPJ = PI*PJ DB00 145 IF (J.EQ.KLM) GO TO 10 DB00 150 DVDPI = DVDPI + (PJ-EIK)*(PIPJ-TM) DB00 155 10 DUDPI = DUDPI + PJ*PJ*(PIPJ*(2.0*PIPJ-3.0)+1.0) DB00 160 15 CONTINUE DB00 165 DUDPI = 2.0*PI*DUDPI DB00 170 DVDPI = 2.0*DVDPI DB00 175 DB = (V*DUDPI-U*DVDPI)/(V*V) DB00 180 RETURN DB00 185 END DB00 190 CC$ FORTY FORM,XREF COMP 5 SUBROUTINE COMPW(IW, IAV) COMP 10 DOUBLE PRECISION DCON, DELBAR COMP 15 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, COMP 20 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, COMP 25 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, COMP 30 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON COMP 35 COMMON /A4/ DAPI(030), BCONST, TM, U, V COMP 40 COMMON /A6/ DCON, DATA(018,030,030), SDDAT(018), DMEAN(018), COMP 45 * ALOS, IXV COMP 50 C COMP 55 C THIS SUBROUTINE COMPUTES THE ESTIMATE OF THE WEIGHT FOR COMP 60 C SUBSET IW AS WELL AS THE INNER ITERATION CONSTANT. THEN COMP 65 C TRANSLATE THE RESIDUALS BY THE CONSTANT FOR THE FIRST COMP 70 C CALL WHEN IAV = 0 COMP 75 C COMP 80 C COMP 85 C TM IS THE MEAN (P(I,IW)*P(J,IW)) I .GT. J COMP 90 C COMP 95 C COMP 100 C COMPUTE THE MEAN AND SQUARIANCE OF PAIRWISE PRODUCTS COMP 105 C IN STRAIGHTFORWARD MANNER TO AVOID TRUNCATION. COMP 110 C COMP 115 TM = 0.0 COMP 120 DO 10 I=2,N COMP 125 TA = P(I,IW) COMP 130 I1 = I - 1 COMP 135 DO 5 J=1,I1 COMP 140 TM = TM + TA*P(J,IW) COMP 145 5 CONTINUE COMP 150 10 CONTINUE COMP 155 TM = TM*ABNV COMP 160 C COMP 165 C DENOMINATOR OF RISE(M,IW) IS THE SUM OF SQUARES ABOUT TM COMP 170 C AND IS CONSTANT ACROSS SUBJECTS. COMP 175 C 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 C COMP 225 C COMP 230 C GET THE NUMERATOR (AW) AND THE WEIGHT -RISE- FOR THE M-TH COMP 235 C SUBJECT ON SUBSET IW. COMP 240 C COMP 245 DO 35 M=1,M3 COMP 250 AW = 0.0 COMP 255 DO 30 I=2,N COMP 260 I1 = I - 1 COMP 265 PI = P(I,IW) COMP 270 DO 25 J=1,I1 COMP 275 AW = AW + DATA(M,J,I)*PI*P(J,IW) COMP 280 25 CONTINUE COMP 285 30 CONTINUE COMP 290 RISE(M,IW) = AW/SQA COMP 295 35 CONTINUE COMP 300 C COMP 305 C THE REST OF THE COMPUTATION IS NOT NECESSARY FOR THE COMP 310 C CLEANUP CALL TO COMPW. COMP 315 C COMP 320 IF (IAV.EQ.1) RETURN COMP 325 C COMP 330 C COMP 335 C NOW FOR COMPUTING THE REGRESSION CONSTANT FOR SUBSET IW. COMP 340 C COMP 345 DO 60 M=1,M3 COMP 350 DELBAR = 0.0D0 COMP 355 DO 45 I=2,N COMP 360 I1 = I - 1 COMP 365 DO 40 J=1,I1 COMP 370 DELBAR = DELBAR + DATA(M,J,I) COMP 375 40 CONTINUE COMP 380 45 CONTINUE COMP 385 DELBAR = ABNV*DELBAR COMP 390 WCON(M) = -RISE(M,IW)*TM + DELBAR COMP 395 C COMP 400 C TRANSLATE THE RESIDUALS. COMP 405 C COMP 410 DO 55 I=2,N COMP 415 I1 = I - 1 COMP 420 DO 50 J=1,I1 COMP 425 DATA(M,J,I) = DATA(M,J,I) - WCON(M) COMP 430 50 CONTINUE COMP 435 55 CONTINUE COMP 440 60 CONTINUE COMP 445 RETURN COMP 450 END COMP 455 CC$ FORTY FORM,XREF PRIN 5 SUBROUTINE PRINTD(MST, IW, M) PRIN 10 COMMON /A1/ P(030,017), DENS(050), NW(050), MP(050), M3, RM, AMV, PRIN 15 * BLMX(050), BLMN(050), KP(050), AB, WCON(018), ABIN, ABM, BIN, PRIN 20 * N2, N, LD2, N1, AN1, LOUT, VARSIM, VAF, ABNV, NDO, PRIN 25 * RISE(018,017), LD1, IFLOP, KMOD, KC, ISWIT, MCON PRIN 30 COMMON /A7/ AZ(030), IZA(030), A(017), IZ(017) PRIN 35 COMMON /A9/ ILAB(030,6), MLAB(018,6), ALD1, SIGN PRIN 40 DIMENSION LSTIM(030) PRIN 45 C PRIN 50 C THIS SUBROUTINE PRINTS OUT THE STIMULI CONTAINED IN SUBSETS. PRIN 55 C PRIN 60 C IF KC .EQ. 0, SUBSET IW IS PRINTED. IF KC .EQ. 1, SUBSETS 1,...NPRIN 65 C ARE PRINTED. IF KC .EQ. 2, THE SUBSETS ARE PRINTED AFTER PRIN 70 C BEING RANKED ACCORDING TO WEIGHT. PRIN 75 C PRIN 80 WRITE (LOUT,99999) M, (MLAB(M,I),I=1,6) PRIN 85 99999 FORMAT (/1X, 6HSUBSET, 4X, 6HWEIGHT, 3X, 21HSTIMULI CONTAINED IN ,PRIN 90 * 6HSUBSET, 15X, 13HFROM SOURCE #, I3, 2H: , 3X, 6A1) PRIN 95 IF (KC.NE.0) GO TO 5 PRIN 100 IA = IW PRIN 105 IB = IW PRIN 110 IZ(IW) = IW PRIN 115 GO TO 15 PRIN 120 5 IA = 1 PRIN 125 IB = LD1 PRIN 130 IF (KC.EQ.2) GO TO 15 PRIN 135 DO 10 KX=1,LD1 PRIN 140 IZ(KX) = KX PRIN 145 10 CONTINUE PRIN 150 15 DO 25 KX=IA,IB PRIN 155 LX = IZ(KX) PRIN 160 LCNT = 0 PRIN 165 DO 20 I=1,N PRIN 170 IF (P(I,LX).LE.0.00001) GO TO 20 PRIN 175 LCNT = LCNT + 1 PRIN 180 LSTIM(LCNT) = I PRIN 185 20 CONTINUE PRIN 190 C PRIN 195 C IF YOUR COMPILER CANNOT HANDLE THE ILAB PART OF THE FOLLOWING PRIN 200 C WRITE STATEMENT, THEN WRITING PRIN 205 C PRIN 210 C (LSTIM(KK), KK = 1, LCNT) PRIN 215 C PRIN 220 C IN I FORMAT WILL LIST THE ORDINAL NUMBERS (1, ..., N) OF THE PRIN 225 C STIMULI (INSTEAD OF THEIR USER-SUPPLIED LABELS) FOR EACH PRIN 230 C SUBSET. ALTERNATIVELY, A SCRATCH ARRAY COULD BE CREATED FOR PRIN 235 C THE LINE OF LABELS BEING PRINTED IN THE PRESENT ARRANGEMENT. PRIN 240 C PRIN 245 WRITE (LOUT,99998) LX,RISE(M,LX), PRIN 250 * ((ILAB(LSTIM(KK),L), L = 1, 6), KK = 1, LCNT) PRIN 255 99998 FORMAT (/1X, 1H(, I4, 1H), F10.5, 2X, 13(6A1, 1X)/(19X, 6A1, 1X, PRIN 260 * 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1, 1X, PRIN 265 * 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1, 1X, 6A1/)) PRIN 270 IF (KC.EQ.1 .AND. MST.GT.0 .AND. M.EQ.1) DENS(MST) = DENS(MST) + PRIN 275 * FLOAT(LCNT) PRIN 280 25 CONTINUE PRIN 285 IF (KC.EQ.0) RETURN PRIN 290 IF (MST.GT.0) GO TO 30 PRIN 295 WRITE (LOUT,99997) RISE(M,LD2) PRIN 300 99997 FORMAT (/31H PLUS AN ADDITIVE CONSTANT OF , F10.5) PRIN 305 RETURN PRIN 310 30 IF (KC.EQ.1 .AND. M.EQ.1) DENS(MST) = DENS(MST)*ALD1 PRIN 315 WRITE (LOUT,99996) RISE(M,LD2), DENS(MST) PRIN 320 99996 FORMAT (/30H PLUS AN ADDITIVE CONSTANT OF , F10.5, 15X, 7HDENSITY,PRIN 325 * 2H =, F10.5) PRIN 330 RETURN PRIN 335 END PRIN 34 CUT HERE------------ indclus.f echo processing file indclus.data 1>&2 cat > indclus.data <<'CUT HERE------------ indclus.data' # indclus.data Sample Data. C The authors of this software are J.Douglas Carroll and Phipps Arabie. 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----+----@----+----@----+----@----+----@----+----@----+----@----+----@ 005 0 5 0 15 0 1 2 17 000000000000000000000000000000000000000000000000000000000000000000000000000000 0 0 0 0 ROSENBERG-KIM 1975 S-MEAS FOR FEMALE,MALE MULT SORTS1,2,