C sindscal.f C The author of this software is Sandra Pruzansky. 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 Arabie,P., Carroll,J.D., & DeSarbo,W.S. (1987). Three-Way Scaling and C Clustering. Newbury Park, CA: Sage, and C Carroll,J.D. & Chang,J.J. (1970), "Analysis of individual differences C in multidimensional scaling via an N-way generalization of C 'Eckart-Young' decomposition" in Psychometrika, 35,283-319, and C Carroll,J.D. (1972). Individual differences and multidimensional C scaling, in R.N.Shepard, A.K.Romney & S.B. Nerlove (Eds.), C Multidimensional Scaling: Theory and Applications in the Behavioral C Sciences (Vol.1, pp.105-155). New York and London: Seminar Press. C The latter two have been reprinted in C P. Davies & A.P.M. Coxon (Eds.), (1984) C "Key Texts in Multidimensional Scaling" Exeter, NH: Heinemann. C First one: pp. 229-252; second: pp. 267-301. C----+----@----+----@----+----@----+----@----+----@----+----@----+----@ C MACHINE DEPENDENT CONSTANT, XMAXN, AT LINE XA00 12 C CURRENTLY SET FOR THE VAX. C REPRESENTS MINIMUM OF LARGEST POSITIVE MACHINE NUMBER C AND MINUS LARGEST NEGATIVE MACHINE NUMBER C MODIFIED MAY, 1987; CHANGE IN RDAT C RDAT 44 COMMENTED OUT TO AVOID LOADER ERROR; C IF A USER-WRITTEN MYREAD SUBROUTINE IS C TO BE USED, REMOVE C FROM FIRST COLUMN C AND RECOMPILE C MODIFIED MAY, 1982; CHANGE IN CANDE C ADDED RELAXATION FACTOR TO SPEED UP CONVERGENCE; C CHANGED CONVERGENCE CRIT TO 1.E-6 FROM 1.E-5 C THIS WILL INCREASE CONSIDERABLY, THE NUMBER OF C ITERATIONS FOR CONVERGENCE BUT WILL REDUCE THE C POSSIBILITY OF LOCAL MINIMA C MODIFIED OCT., 1977; CHANGE IN SERCH X000 3 C MODIFIED AUG., 1977; CHANGES IN SORT, RUNIFV, SMTINV X000 4 C WRITTEN BY SANDRA PRUZANSKY X000 5 C BELL TELEPHONE LABORATORIES X000 6 C MURRAY HILL, NEW JERSEY X000 7 C X000 8 C BASED ON ALGORITHM BY CARROLL AND CHANG PUBLISHED IN X000 9 C PSYCHOMETRIKA,35:1970,283-319. X000 10 C X000 11 C X000 12 C PARAMETER DEFINITIONS X000 13 C*************************************** X000 14 C MAXDIM- MAXIMUM NUMBER OF DIMENSIONS (MAX. =10) X000 15 C MINDIM- MINIMUM NUMBER OF DIMENSIONS X000 16 C NWT(1) - NUMBER OF MATRICES X000 17 C NWT(2)- NUMBER OF STIMULI X000 18 C IRDATA- 1, READ IN LOWERHALF OF SIMILARITY MATRIX WITHOUT DIAGS. X000 19 C IRDATA- 2, READ IN DISSIMILARITIES, LOWERHALF WITHOUT DIAGS. X000 20 C IRDATA- 3, READ IN EUCLIDEAN DISTANCES, LOWERHALF WITHOUT DIAGS. X000 21 C IRDATA- 4, READ IN LOWERHALF CORRELATION MATRIX WITHOUT DIAGS. X000 22 C -4, SAME AS ABOVE BUT PROGRAM WILL NORMALIZE SO SSQ=1. X000 23 C IRDATA- 5, READ IN LOWERHALF COVARIANCE MATRIX WITH DIAGS. X000 24 C -5, SAME AS ABOVE BUT PROGRAM WILL NORMALIZE SO SSQ=1. X000 25 C IRDATA- 6, READ IN FULL SYMMETRIC SIMILARITY MATRIX, DIAG. IGNOREDX000 26 C IRDATA- 7, READ IN FULL SYMMETRIC DISSIM. MATRIX,DIAG. IGNORED X000 27 C ITMAX - MAXIMUM NUMBER OF ITERATIONS. X000 28 C IF ITMAX=0 PROGRAM WILL SOLVE FOR WEIGHTS ONLY. X000 29 C PUNCH OPTIONS - ALL PUNCH OPTIONS INCLUDE PUNCHING NORMALIZED X000 30 C SOLUTION EXCEPT -1. X000 31 C IPUNSP -1, NO PUNCHED OUTPUT X000 32 C IPUNSP - 0, PUNCH NORMALIZED SOLUTION X000 33 C IPUNSP - 1, PUNCH SCALAR PRODUCT MATRICES. X000 34 C IPUNSP - 2, PUNCH UNNORMALIZED MATRICES - A1,A2,A3 X000 35 C IPUNSP - 3, PUNCH ALL OF THE ABOVE X000 36 C IPLOT - -1, NO PLOT X000 37 C IPLOT - 0, PLOT WITH POINTS X000 38 C IRN - 0, READ IN INITIAL CONFIGURATION. X000 39 C IRN - AN INTEGER FOR GENERATING A RANDOM STARTING X000 40 C CONFIGURATION. (MAXIMUM 4 DIGITS) X000 41 C ARRAYS X000 42 C Y - DATA. Y IS A 2 SUBSCRIPTED VARIABLE. X000 43 C Y(N1,LSD), WHERE N1 IS NO. OF MATRICES X000 44 C AND LSD IS (N2(N2+1))/2;(N2=# OF STIM) X000 45 C A - COORDINATE MATRIX. X000 46 C A(I,J,K), WHERE I - ELEMENTS WITHIN MATRICES X000 47 C J - DIMENSIONS. X000 48 C K - MATRIX. X000 49 C AA - SCRATCH MATRIX X000 50 C AA(I,J) - WHERE I - MAXIMUM OF # STIM AND # OF MATRICES X000 51 C J - MAXIMUM # OF DIMENSIONS X000 52 C TIT - TITLE, USE NO MORE THAN ONE CARD. X000 53 C X000 54 C ----- DATA CARD INPUT FORMAT: ----- X000 55 C X000 56 C FIRST CARD (4I4): MAXDIM,MINDIM,NWT(1),NWT(2) X000 57 C SECOND CARD (5I4): ITMAX,IRDATA,IPUNSP,IPLOT,IRN X000 58 C THIRD CARD (80A1) : TITLE CARD FOR THIS DATA SET X000 59 C FOURTH CARD (80A1) : (DATA FORMAT ENCLOSED IN PARENTHESES) X000 60 C OR X000 61 C MYREAD (IN COLUMNS 1-6) X000 62 C FIFTH TO NTH CARD : DATA IN FORMAT GIVEN ON FOURTH CARD X000 63 C OPTIONAL CONFIGURATION CARDS - (IF IRN=0) X000 64 C FIRST CONFIG. CARD - (80A1) :FORMAT FOR STARTING CONFIGURATION X000 65 C SECOND THRU LAST CONFIG. CARD : STARTING STIMULUS X000 66 C CONFIGURATION IN FORMAT SPECIFIED BY X000 67 C PRECEDING CARD. X000 68 C LAST CARD: MUST BE BLANK IN COLUMNS 1-20 X000 69 C************************* X000 70 C X000 71 C X000 72 C MAIN PROGRAM FOR SINDSCAL - (FOR PORTABLE USE) X000 73 C WRITTEN BY S. PRUZANSKY X000 74 C LATEST CHANGE 9-16-75 X000 75 DIMENSION ARRAY(5000) X000 76 COMMON /PARAM/ MAXDIM, MINDIM, ITMAX, IRDATA, N, IPUNSP, X000 77 * IPLOT, IRN, CRIT, NWT(3) X000 78 COMMON /DEVICE/ LREAD, LPRINT, LPUNCH X000 79 DATA MAXCOR /5000/ X000 80 C DESIGNATE DEVICE NUMBERS FOR INPUT, OUTPUT X000 81 LREAD = 5 X000 82 LPRINT = 6 X000 83 LPUNCH = 7 X000 84 C READ INPUT PARMETERS NEEDED FOR STORAGE ALLOCATION X000 85 10 READ (LREAD,9999) MAXDIM, MINDIM, NWT(1), NWT(2) X000 86 IF (MAXDIM.EQ.0) STOP X000 87 C CHECK VALIDITY OF PARAMETER CARD X000 88 IF (MAXDIM.GT.10) GO TO 20 X000 89 IF (MINDIM.GT.10) GO TO 20 X000 90 9999 FORMAT (4I4) X000 91 READ (LREAD,9998) ITMAX, IRDATA, IPUNSP, IPLOT, IRN X000 92 9998 FORMAT (5I4) X000 93 C COMPUTE AMOUNT OF CORE NEEDED FOR ARRAYS X000 94 LY = NWT(1)*((NWT(2)*(NWT(2)+1))/2) X000 95 N = MAX0(NWT(1),NWT(2)) X000 96 LA = MAX0(N*MAXDIM*3,NWT(2)**2) X000 97 LAA = N*MAXDIM X000 98 MEMORY = LY + LA + LAA X000 99 IF (MEMORY.GT.MAXCOR) GO TO 30 X000 100 C COMPUTE STARTING LOCATIONS FOR THREE STORAGE ARRAYS X000 101 LOCY = 1 X000 102 LOCA = LOCY + LY X000 103 LOCAA = LOCA + LA X000 104 CRIT = 1.E-6 X000 105 CALL SINSCL(ARRAY(LOCY), LY, ARRAY(LOCA), LA, X000 106 * ARRAY(LOCAA), LAA) X000 107 GO TO 10 X000 108 20 WRITE (LPRINT,9997) X000 109 9997 FORMAT (45H0REQUEST EXCEEDS MAXIMUM NUMBER OF DIMENSIONS/ X000 110 * 23H ***CHECK DECK SETUP***) X000 111 STOP X000 112 30 WRITE (LPRINT,9996) X000 113 9996 FORMAT (46H0CORE REQUESTED EXCEEDS LIMITS OF DIMENSIONED , X000 114 * 5HARRAY, 32H ***CHECK VARIABLES LY,LA,LAA***) X000 115 STOP X000 116 END X000 117 C SINS 1 SUBROUTINE SINSCL(Y, LY, A, LA, AA, LAA) SINS 2 C SINSCL FOR SINDSCAL SEPTEMBER, 1975 SINS 3 C WRITTEN BY S. PRUZANSKY SINS 4 C CONTROL PROGRAM FOR INDSCAL ANALYSIS SINS 5 COMMON /PARAM/ MAXDIM, MINDIM, ITMAX, IRDATA, MNWT, SINS 6 * IPUNSP, IPLOT, IRN, CRIT, NWT(3) SINS 7 COMMON /DEVICE/ LREAD, LPRINT, LPUNCH SINS 8 DIMENSION Y(LY), A(LA), AA(LAA) SINS 9 INTEGER TIT(80) SINS 10 C READ TITLE SINS 11 READ (LREAD,9999) (TIT(I),I=1,80) SINS 12 C SET PARAMETERS SINS 13 NA = 2 SINS 14 NWT(3) = NWT(2) SINS 15 N = 3 SINS 16 N1 = NWT(1) SINS 17 N2 = NWT(2) SINS 18 LSD = (N2*(N2+1))/2 SINS 19 NDIM = 0 SINS 20 C DO INDSCAL ANALYSES IN MAXDIM TO MINDIM DIMENSIONS SINS 21 MMDIM = MAXDIM + MINDIM SINS 22 DO 60 IDIM=MINDIM,MAXDIM SINS 23 NF = MMDIM - IDIM SINS 24 MAXIT = ITMAX SINS 25 IF (ITMAX.EQ.0) GO TO 10 SINS 26 ISTI = 0 SINS 27 GO TO 20 SINS 28 10 ISTI = 1 SINS 29 20 IF (IDIM.GT.MINDIM) WRITE (LPRINT,9992) SINS 30 C PRINT TITLE AND PARAMETER INFO. SINS 31 WRITE (LPRINT,9997) SINS 32 WRITE (LPRINT,9998) (TIT(I),I=1,80) SINS 33 WRITE (LPRINT,9995) SINS 34 WRITE (LPRINT,9994) SINS 35 WRITE (LPRINT,9996) NF, IRDATA, MAXIT, IPUNSP, IPLOT, SINS 36 * IRN SINS 37 WRITE (LPRINT,9993) (NWT(I),I=1,NA) SINS 38 WRITE (LPRINT,9995) SINS 39 IF (ISTI.NE.0) MAXIT = 0 SINS 40 JRDATA = IABS(IRDATA) SINS 41 IF (NDIM.EQ.0) GO TO 30 SINS 42 C REDUCE DIMENSIONALITY OF STARTING CONFIG. SINS 43 NFO = NF + 1 SINS 44 NQ = NFO*N SINS 45 CALL DIMA(A, MNWT, NQ, NFO, NF) SINS 46 GO TO 40 SINS 47 C READ AND PREPROCESS INPUT DATA SINS 48 30 CALL RDATA(Y, N1, N2, LSD, JRDATA, A, IPUNSP) SINS 49 C PERFORM CANONICAL DECOMPOSITION SINS 50 40 CALL CANDE(Y, N1, N2, LSD, NA, A, MNWT, NF, N, NWT, SINS 51 * IRDATA, ISTI, MAXIT, CRIT, IRN, NDIM, FVAF, AA, SINS 52 * IPUNSP) SINS 53 C NORMALIZE FINAL WEIGHTS AND STIMULUS MATRICES SINS 54 CALL NORMA(A, AA, MNWT, NF, N1, N2, NA, FVAF, AA, SINS 55 * IPUNSP) SINS 56 C COMPUTE CORRELATIONS BETWEEN SCALAR PROD. AND SOLUTION FOR EACH SINS 57 C DATA MATRIX. SINS 58 CALL SUBR(A, MNWT, NF, N, Y, N1, N2, LSD) SINS 59 IF (IPLOT.LT.0) GO TO 50 SINS 60 C PLOT FINAL WEIGHTS AND STIMULUS SPACES SINS 61 CALL PPLOT(A, MNWT, NF, TIT, NWT) SINS 62 50 NDIM = 1 SINS 63 60 CONTINUE SINS 64 9999 FORMAT (80A1) SINS 65 9998 FORMAT (1H0, 80A1) SINS 66 9997 FORMAT (1H0, 23X, 17HSYMMETRIC INDSCAL) SINS 67 9996 FORMAT (1H , I3, 4I7, 2X, I7) SINS 68 9995 FORMAT (46H0*********************************************, SINS 69 * 5H*****) SINS 70 9994 FORMAT (11H PARAMETERS/20H DIM IRDATA ITMAX, 7H IPUNCH, SINS 71 * 14H IPLOT IRN) SINS 72 9993 FORMAT (18H NO. OF MATRICES =, I4, 16H NO. OF STIM. =, SINS 73 * I4) SINS 74 9992 FORMAT (1H1) SINS 75 RETURN SINS 76 END SINS 77 C RDAT 1 SUBROUTINE RDATA(Y, N1, N2, LSD, IRDATA, D, IPUNSP) RDAT 2 C RDATA FOR SINDSCAL SEPTEMBER, 1975 RDAT 3 C WRITTEN BY S. PRUZANSKY RDAT 4 COMMON /DEVICE/ LREAD, LPRINT, LPUNCH RDAT 5 DIMENSION Y(N1,LSD), D(N2,N2), MY(2) RDAT 6 INTEGER FMT(80) RDAT 7 DATA IC /1HI/, LP /1H)/, MY(1) /1HM/, MY(2) /1HY/ RDAT 8 C IRDATA- 1, READ IN LOWERHALF OF SIMILARITY MATRIX WITHOUT DIAGS. RDAT 9 C IRDATA- 2, READ IN DISSIMILARITIES, LOWERHALF WITHOUT DIAGS. RDAT 10 C IRDATA- 3, READ IN EUCLIDEAN DISTANCES, LOWERHALF WITHOUT DIAGS. RDAT 11 C IRDATA- 4, READ IN LOWERHALF CORRELATION MATRIX WITHOUT DIAGS. RDAT 12 C IRDATA- 5, READ IN LOWERHALF COVARIANCE MATRIX WITH DIAGS. RDAT 13 C IRDATA- 6, READ IN FULL SYMMETRIC SIMILARITY MATRIX, DIAG. IGNOREDRDAT 14 C IRDATA- 7, READ IN FULL SYMMETRIC DISSIMILARITY MATRIX,DIAG. IGNORRDAT 15 C IADDC=1, COMPUTE ADDITIVE CONSTANT. RDAT 16 C IADDC=0, DO NOT COMPUTE ADDITIVE CONSTANT. RDAT 17 C LSD - LENGTH OF SYM. DATA VECTOR RDAT 18 C********************************************* RDAT 19 MIP2 = MOD(IPUNSP,2) RDAT 20 IF (MIP2.LE.0) GO TO 10 RDAT 21 WRITE (LPUNCH,9997) RDAT 22 10 IADDC = 1 RDAT 23 IFIRST = 2 RDAT 24 II = N2 RDAT 25 IF (IRDATA.EQ.3) IADDC = 0 RDAT 26 IF (IRDATA.GE.5) IFIRST = 1 RDAT 27 READ (LREAD,9999) (FMT(I),I=1,80) RDAT 28 C BYPASS FORMAT CHECKING RDAT 29 9999 FORMAT (80A1) RDAT 30 IF (FMT(1).EQ.MY(1) .AND. FMT(2).EQ.MY(2)) GO TO 30 RDAT 31 C CHECK THAT DATA FORMAT IS PROPER RDAT 32 DO 20 I=1,80 RDAT 33 IF (FMT(I).EQ.LP) GO TO 30 RDAT 34 IF (FMT(I).NE.IC) GO TO 20 RDAT 35 WRITE (LPRINT,9998) RDAT 36 STOP RDAT 37 20 CONTINUE RDAT 38 9998 FORMAT (38H*** FORMAT FOR INPUT DATA IS IMPROPER, RDAT 39 * 24H - MUST BE E OR F FORMAT) RDAT 40 30 DO 190 I1=1,N1 RDAT 41 C READ USERS DATA READ SUBROUTINE RDAT 42 IF (FMT(1).NE.MY(1) .OR. FMT(2).NE.MY(2)) GO TO 40 RDAT 43 C CALL MYREAD(D, N2, I1) RDAT 44 GO TO 90 RDAT 45 40 IN = 0 RDAT 46 DO 80 J=IFIRST,N2 RDAT 47 IF (IRDATA-5) 50, 60, 70 RDAT 48 50 II = J - 1 RDAT 49 GO TO 70 RDAT 50 60 II = J RDAT 51 70 READ (LREAD,FMT) (D(J,M),M=1,II) RDAT 52 80 CONTINUE RDAT 53 90 DO 110 J=2,N2 RDAT 54 JJ = J - 1 RDAT 55 DO 100 M=1,JJ RDAT 56 IF (IRDATA.EQ.1 .OR. IRDATA.EQ.6) D(J,M) = -D(J,M) RDAT 57 D(M,J) = D(J,M) RDAT 58 100 CONTINUE RDAT 59 110 CONTINUE RDAT 60 GO TO (120, 120, 120, 140, 160, 120, 120), IRDATA RDAT 61 C COMPUTE ADDITIVE CONSTANT AND SCALAR PRODUCT MATRICES. RDAT 62 120 CALL TORGB(D, N2, IADDC) RDAT 63 IF (MIP2.LE.0) GO TO 160 RDAT 64 C PUNCH LOWER-HALF SCALAR PROD. MATRIX WITH DIAG RDAT 65 DO 130 I=1,N2 RDAT 66 WRITE (LPUNCH,9996) I1, I, (D(I,J),J=1,I) RDAT 67 130 CONTINUE RDAT 68 GO TO 160 RDAT 69 140 DO 150 I=1,N2 RDAT 70 D(I,I) = 1. RDAT 71 150 CONTINUE RDAT 72 C STORE LOWER HALF SCALAR PROD. MATRIX WITH DIAG IN Y RDAT 73 160 K = 0 RDAT 74 DO 180 I=1,N2 RDAT 75 DO 170 J=1,I RDAT 76 K = K + 1 RDAT 77 Y(I1,K) = D(I,J) RDAT 78 170 CONTINUE RDAT 79 180 CONTINUE RDAT 80 190 CONTINUE RDAT 81 9997 FORMAT (38H(6X,8F8.4) SCALAR PRODUCT MATRICES) RDAT 82 9996 FORMAT (2I3, 8F8.4/(F14.4, 7F8.4)) RDAT 83 RETURN RDAT 84 END RDAT 85 C TORG 1 SUBROUTINE TORGB(A, N, IADDC) TORG 2 C TORGB FOR SINDSCAL SEPTEMBER, 1975 TORG 3 C WRITTEN BY J.J. CHANG FOR INDSCAL TORG 4 C MODIFIED FOR SINDSCAL BY S. PRUZANSKY TORG 5 C TORGERSON METHOD TO COMPUTE SCALAR PRODUCT MATRICES TORG 6 DIMENSION DR(200), A(N,N) TORG 7 XN = N TORG 8 DO 10 I=1,N TORG 9 DR(I) = 0. TORG 10 A(I,I) = 0. TORG 11 10 CONTINUE TORG 12 SUMD = 0. TORG 13 IF (IADDC) 110, 110, 20 TORG 14 C COMPUTE ADDITIVE CONSTANT TORG 15 20 ADDC = -100. TORG 16 N1 = N - 1 TORG 17 N2 = N - 2 TORG 18 DO 100 I=1,N2 TORG 19 JJ = I + 1 TORG 20 DO 90 J=JJ,N1 TORG 21 KK = J + 1 TORG 22 DO 80 K=KK,N TORG 23 A1 = A(I,K) TORG 24 A2 = A(J,K) TORG 25 A3 = A(I,J) TORG 26 DO 70 M=1,3 TORG 27 GO TO (30, 40, 50), M TORG 28 30 CIJK = A1 - A2 - A3 TORG 29 GO TO 60 TORG 30 40 CIJK = A2 - A1 - A3 TORG 31 GO TO 60 TORG 32 50 CIJK = A3 - A2 - A1 TORG 33 60 IF (CIJK.GT.ADDC) ADDC = CIJK TORG 34 70 CONTINUE TORG 35 80 CONTINUE TORG 36 90 CONTINUE TORG 37 100 CONTINUE TORG 38 110 DO 140 I=2,N TORG 39 II = I - 1 TORG 40 DO 130 J=1,II TORG 41 IF (IADDC.EQ.0) GO TO 120 TORG 42 A(I,J) = A(I,J) + ADDC TORG 43 A(J,I) = A(I,J) TORG 44 120 SUMD = SUMD + A(I,J)**2 TORG 45 130 CONTINUE TORG 46 140 CONTINUE TORG 47 ADIJ = 2.*SUMD/XN**2 TORG 48 DO 160 I=1,N TORG 49 DO 150 J=1,N TORG 50 DR(I) = DR(I) + A(I,J)**2 TORG 51 150 CONTINUE TORG 52 160 CONTINUE TORG 53 C COMPUTE SCALAR PRODUCT MATRIX TORG 54 SUM = 0. TORG 55 SUMD = 0. TORG 56 DO 190 I=1,N TORG 57 DO 180 J=1,I TORG 58 QQ = (DR(I)+DR(J))/XN TORG 59 A(I,J) = (QQ-A(I,J)**2-ADIJ)*(.5) TORG 60 IF (I.EQ.J) GO TO 170 TORG 61 SUM = SUM + A(I,J)**2 TORG 62 GO TO 180 TORG 63 170 SUMD = SUMD + A(I,J)**2 TORG 64 180 CONTINUE TORG 65 190 CONTINUE TORG 66 C NORMALIZE SCALAR PRODUCT MATRIX TORG 67 Q = 1./SQRT(SUMD+2.*SUM) TORG 68 DO 210 I=1,N TORG 69 DO 200 J=1,I TORG 70 A(I,J) = A(I,J)*Q TORG 71 200 CONTINUE TORG 72 210 CONTINUE TORG 73 RETURN TORG 74 END TORG 75 C CAND 1 SUBROUTINE CANDE(Y, N1, N2, LSD, NA, A, MNWT, ND, N, NWT, CAND 2 * IRDATA, JSTI, NAXIT, CRIT, IRN, NDIM, FVAF, S, IPUN) CAND 3 C CANDE FOR SINDSCAL SEPTEMBER, 1975 CAND 4 C WRITTEN BY S. PRUZANSKY CAND 5 C LATEST CHANGE 5-20-82 BY SP CAND 6 C RELAXATION FACTOR INTRODUCED TO SPEED UP CONVERGENCE COMMON /DEVICE/ LREAD, LPRINT, LPUNCH CAND 7 DIMENSION Y(N1,LSD), A(MNWT,ND,3), S(MNWT,ND), R(55) CAND 8 DIMENSION NWT(3), SCR(55), SUMM(2) CAND 9 INTEGER FMT(80) CAND 10 ITCOM = 0 CAND 11 C THE FOLLOWING STATEMENT WAS INSERTED 5-20-82 ALPHA=1.5 MAXIT = NAXIT CAND 12 ISTI = JSTI CAND 13 NF = ND CAND 14 IF (NDIM.EQ.1) GO TO 40 CAND 15 SSQY = N1 CAND 16 C COMPUTE SSQ FOR COVAR. OR CORR. DATA (WITH OPTION TO NORMALIZE) CAND 17 IAD = IABS(IRDATA) CAND 18 IF (IAD.EQ.4 .OR. IAD.EQ.5) CALL SSQC(Y, N1, LSD, S, N2, CAND 19 * SSQY, IRDATA) CAND 20 IF (IRN) 10, 20, 10 CAND 21 C GENERATE RANDOM STARTING CONFIG. CAND 22 10 CALL RCON(A, MNWT, NF, N2, IRN) CAND 23 GO TO 70 CAND 24 C READ IN STARTING STIM. CONFIG. CAND 25 20 READ (LREAD,9997) (FMT(I),I=1,80) CAND 26 DO 30 K=1,NF CAND 27 READ (LREAD,FMT) (A(J,K,2),J=1,N2) CAND 28 30 CONTINUE CAND 29 C SET TABLE 2 EQ. TABLE 3 TO START CAND 30 40 DO 60 K=1,NF CAND 31 DO 50 J=1,N2 CAND 32 A(J,K,3) = A(J,K,2) CAND 33 50 CONTINUE CAND 34 60 CONTINUE CAND 35 70 WRITE (LPRINT,9991) CAND 36 DO 80 J=1,NF CAND 37 WRITE (LPRINT,9996) J, (A(M,J,2),M=1,N2) CAND 38 80 CONTINUE CAND 39 C BEGIN FITTING CYCLE CAND 40 90 IT = 0 CAND 41 100 DO 350 JZ=1,N CAND 42 C THE FOLLOWING STATEMENT MODIFIED 5-20-82 JB = JZ CAND 43 IF (IT.NE.0 .AND. ISTI.EQ.0) GO TO 110 CAND 44 IF (JB.GE.2) GO TO 350 CAND 45 C COMPUTE INVERSE CAND 46 110 KL = 0 CAND 47 DO 150 K1=1,NF CAND 48 DO 140 K2=1,K1 CAND 49 KL = KL + 1 CAND 50 SUM1 = 1. CAND 51 DO 130 I=1,N CAND 52 IF (I.EQ.JB) GO TO 130 CAND 53 MM = NWT(I) CAND 54 SUM = 0. CAND 55 DO 120 J=1,MM CAND 56 SUM = SUM + A(J,K1,I)*A(J,K2,I) CAND 57 120 CONTINUE CAND 58 SUM1 = SUM*SUM1 CAND 59 130 CONTINUE CAND 60 SCR(KL) = SUM1 CAND 61 140 CONTINUE CAND 62 150 CONTINUE CAND 63 CALL SMTINV(SCR, SCR, NF, R, IER) CAND 64 IF (IER.EQ.0) GO TO 160 CAND 65 C PRINT ERROR MESSAGE CAND 66 WRITE (LPRINT,9995) CAND 67 STOP CAND 68 160 IF (JB.GT.1) GO TO 230 CAND 69 C COMPUTE S FOR WEIGHTS MATRIX (JB=1) CAND 70 DO 220 I1=1,N1 CAND 71 DO 170 L=1,NF CAND 72 S(I1,L) = 0 CAND 73 170 CONTINUE CAND 74 LC = 0 CAND 75 DO 210 J=1,N2 CAND 76 DO 200 K=1,J CAND 77 LC = LC + 1 CAND 78 DO 190 L1=1,NF CAND 79 IF (J.EQ.K) GO TO 180 CAND 80 S(I1,L1) = Y(I1,LC)*A(J,L1,2)*A(K,L1,3) + CAND 81 * S(I1,L1) CAND 82 180 S(I1,L1) = Y(I1,LC)*A(K,L1,2)*A(J,L1,3) + CAND 83 * S(I1,L1) CAND 84 190 CONTINUE CAND 85 200 CONTINUE CAND 86 210 CONTINUE CAND 87 220 CONTINUE CAND 88 GO TO 300 CAND 89 C COMPUTE S FOR STIM MATRICES (JB>1) CAND 90 230 JC = 5 - JB CAND 91 JRB = 1 CAND 92 DO 290 J=1,N2 CAND 93 LC = JRB CAND 94 DO 240 L=1,NF CAND 95 S(J,L) = 0 CAND 96 240 CONTINUE CAND 97 DO 280 K=1,N2 CAND 98 DO 260 I1=1,N1 CAND 99 DO 250 L1=1,NF CAND 100 S(J,L1) = Y(I1,LC)*A(I1,L1,1)*A(K,L1,JC) + CAND 101 * S(J,L1) CAND 102 250 CONTINUE CAND 103 260 CONTINUE CAND 104 IF (K.LT.J) GO TO 270 CAND 105 LC = LC + K CAND 106 GO TO 280 CAND 107 270 LC = LC + 1 CAND 108 280 CONTINUE CAND 109 JRB = JRB + J CAND 110 290 CONTINUE CAND 111 C COMPUTE A CAND 112 300 NN = NWT(JB) CAND 113 DO 340 I=1,NN CAND 114 JRI = 1 CAND 115 DO 330 J=1,NF CAND 116 KL = JRI CAND 117 SUM = 0. CAND 118 DO 320 K=1,NF CAND 119 SUM = SUM + S(I,K)*R(KL) CAND 120 IF (K.LT.J) GO TO 310 CAND 121 KL = KL + K CAND 122 GO TO 320 CAND 123 310 KL = KL + 1 CAND 124 320 CONTINUE CAND 125 JRI = JRI + J CAND 126 C THE FOLLOWING STATEMENT WAS INSERTED 5-20-82 IF(IT.NE.0) A(I,J,JB)=ALPHA*SUM+A(I,J,JB)*(1.-ALPHA) C THE FOLLOWING STATEMENT WAS MODIFIED 5-20-82 IF(IT.EQ.0) A(I,J,JB) = SUM CAND 127 330 CONTINUE CAND 128 340 CONTINUE CAND 129 350 CONTINUE CAND 130 C COMPUTE Y HAT AND ERROR CAND 131 ERR = 0. CAND 132 SSQYH = 0. CAND 133 DO 410 I1=1,N1 CAND 134 LC = 0 CAND 135 DO 400 I2=1,N2 CAND 136 DO 390 I3=1,I2 CAND 137 KE = 2 CAND 138 LC = LC + 1 CAND 139 SUMM(1) = 0 CAND 140 SUMM(2) = 0. CAND 141 DO 370 I=1,NF CAND 142 IF (I3.EQ.I2) GO TO 360 CAND 143 SUMM(2) = SUMM(2) + A(I1,I,1)*A(I3,I,2)*A(I2,I,3) CAND 144 360 SUMM(1) = SUMM(1) + A(I1,I,1)*A(I2,I,2)*A(I3,I,3) CAND 145 370 CONTINUE CAND 146 YS = Y(I1,LC) CAND 147 IF (I2.EQ.I3) KE = 1 CAND 148 DO 380 I=1,KE CAND 149 DIF = YS - SUMM(I) CAND 150 SSQYH = SSQYH + SUMM(I)**2 CAND 151 ERR = ERR + DIF**2 CAND 152 380 CONTINUE CAND 153 390 CONTINUE CAND 154 400 CONTINUE CAND 155 410 CONTINUE CAND 156 YY = (ERR-SSQY-SSQYH)*(-.5) CAND 157 CORYY = YY/SQRT(SSQY*SSQYH) CAND 158 RSQ = CORYY**2 CAND 159 ERR = ERR/FLOAT(N1) CAND 160 IF (ISTI.LT.2 .AND. IT.EQ.0) WRITE (LPRINT,9992) CAND 161 IF (ISTI.GT.0) GO TO 420 CAND 162 WRITE (LPRINT,9990) IT, CORYY, RSQ, ERR CAND 163 GO TO 430 CAND 164 420 WRITE (LPRINT,9989) CORYY, RSQ, ERR CAND 165 FVAF = RSQ CAND 166 GO TO 520 CAND 167 430 CONTINUE CAND 168 IF (IT.LE.1) GO TO 440 CAND 169 ETEST = PERR - ERR CAND 170 IF (ETEST.GT.CRIT) GO TO 440 CAND 171 ITCOM = 1 CAND 172 WRITE (LPRINT,9994) IT CAND 173 GO TO 450 CAND 174 440 PERR = ERR CAND 175 IF (IT.LT.MAXIT) GO TO 460 CAND 176 WRITE (LPRINT,9993) CAND 177 ITCOM = 2 CAND 178 450 FVAF = RSQ CAND 179 460 IF ((IPUN-1).LE.0) GO TO 500 CAND 180 IF (ITCOM.GT.0) GO TO 470 CAND 181 C PRINT UNNORMALIZED A MATRICES EVERY 10 TH ITERATION CAND 182 IF (IT.EQ.0) GO TO 510 CAND 183 IF (MOD(IT,10).NE.0) GO TO 510 CAND 184 470 WRITE (LPUNCH,9999) IT CAND 185 DO 490 I=1,N CAND 186 NN = NWT(I) CAND 187 DO 480 J=1,NF CAND 188 WRITE (LPUNCH,9998) I, J, (A(M,J,I),M=1,NN) CAND 189 480 CONTINUE CAND 190 490 CONTINUE CAND 191 9999 FORMAT (38H UNNORMALIZED MATRICES AFTER ITERATION, I4) CAND 192 9998 FORMAT (I3, 1X, I3, 10F7.3/(7X, 10F7.3)) CAND 193 500 IF (ITCOM.GT.0) GO TO 520 CAND 194 510 IT = IT + 1 CAND 195 GO TO 100 CAND 196 C SET TABLE 2=TABLE 3; ITERATE AGAIN TO COMPUTE SUBJ. WEIGHTS CAND 197 520 IF (ISTI.NE.0) GO TO 550 CAND 198 DO 540 I=1,NF CAND 199 DO 530 J=1,N2 CAND 200 A(J,I,2) = A(J,I,3) CAND 201 530 CONTINUE CAND 202 540 CONTINUE CAND 203 MAXIT = 0 CAND 204 ISTI = 2 CAND 205 GO TO 90 CAND 206 550 RETURN CAND 207 9997 FORMAT (80A1) CAND 208 9996 FORMAT (1H , I3, 10F7.3/(1H 3X, 10F7.3)) CAND 209 9995 FORMAT (21H0R CANNOT BE INVERTED) CAND 210 9994 FORMAT (31H0REACHED CRITERION ON ITERATION, I3) CAND 211 9993 FORMAT (26HREACHED MAXIMUM ITERATIONS) CAND 212 9992 FORMAT (1H0, 23X, 22HHISTORY OF COMPUTATION/10H ITERATION, CAND 213 * 8X, 20HCORRELATIONS BETWEEN, 7X, 3HVAF, 17X, 4HLOSS/ CAND 214 * 19X, 16HY(DATA) AND YHAT, 8X, 6H(R**2), 12X, 6H(Y-YHA, CAND 215 * 5HT)**2) CAND 216 9991 FORMAT (24H0INITIAL STIMULUS MATRIX) CAND 217 9990 FORMAT (I7, 3X, 3F20.6) CAND 218 9989 FORMAT (4X, 5HFINAL, 1X, 3F20.6) CAND 219 END CAND 220 C SSQC 1 SUBROUTINE SSQC(Y, N1, LSD, USSQ, N2, TSSQ, IRDATA) SSQC 2 C SSQC FOR SINDSCAL SEPTEMBER, 1975 SSQC 3 C WRITTEN BY S. PRUZANSKY SSQC 4 C SUBROUTINE TO COMPUTE SSQ (WITH OPTION TO NORMALIZE) SSQC 5 C FOR CORR. AND COVAR. DATA SSQC 6 COMMON /DEVICE/ LREAD, LPRINT, LPUNCH SSQC 7 DIMENSION Y(N1,LSD), USSQ(N2) SSQC 8 TSSQ = 0 SSQC 9 DO 70 I1=1,N1 SSQC 10 SSQD = 0 SSQC 11 SSQ = 0 SSQC 12 K = 0 SSQC 13 DO 30 J=1,N2 SSQC 14 DO 20 L=1,J SSQC 15 K = K + 1 SSQC 16 IF (J.EQ.L) GO TO 10 SSQC 17 SSQ = Y(I1,K)**2 + SSQ SSQC 18 GO TO 20 SSQC 19 10 SSQD = Y(I1,K)**2 + SSQD SSQC 20 20 CONTINUE SSQC 21 30 CONTINUE SSQC 22 USSQ(I1) = SSQ*2. + SSQD SSQC 23 IF (IRDATA.GT.0) GO TO 60 SSQC 24 C NORMALIZE DATA SO SSQ=1 SSQC 25 TSSQ = N1 SSQC 26 Q = 1./SQRT(USSQ(I1)) SSQC 27 K1 = 0 SSQC 28 DO 50 I=1,N2 SSQC 29 DO 40 J=1,I SSQC 30 K1 = K1 + 1 SSQC 31 Y(I1,K1) = Y(I1,K1)*Q SSQC 32 40 CONTINUE SSQC 33 50 CONTINUE SSQC 34 GO TO 70 SSQC 35 60 TSSQ = TSSQ + USSQ(I1) SSQC 36 70 CONTINUE SSQC 37 IF (IRDATA.LT.0) RETURN SSQC 38 C PRINT SSQ FOR UNNORMALIZED DATA SSQC 39 WRITE (LPRINT,9999) SSQC 40 DO 80 I=1,N1,10 SSQC 41 IE = I + 9 SSQC 42 IF (IE.GT.N1) IE = N1 SSQC 43 WRITE (LPRINT,9998) I, (USSQ(L),L=I,IE) SSQC 44 80 CONTINUE SSQC 45 9999 FORMAT (46H0SUM OF SQUARES OF UNNORMALIZED CORR. OR COVAR, SSQC 46 * 2H. , 8HMATRICES) SSQC 47 9998 FORMAT (1H I3, 10F7.2) SSQC 48 RETURN SSQC 49 END SSQC 50 C NORM 1 SUBROUTINE NORMA(A, AA, MNWT, K, N1, N2, NA, FVAF, S, NORM 2 * IPUN) NORM 3 C NORMA FOR SINDSCAL SEPTEMBER, 1975 NORM 4 C WRITTEN BY J.J. CHANG NORM 5 C MODIFIED FOR SINDSCAL BY S. PRUZANSKY NORM 6 C PROGRAM TO NORMALIZE A MATRICES FROM THE CANNONICAL DECOMPOSITION NORM 7 C OF 3 WAY TABLES - TABLE 2 MUST EQ. TABLE 3 NORM 8 COMMON /DEVICE/ LREAD, LPRINT, LPUNCH NORM 9 DIMENSION NWT(3), A(MNWT,K,3), S(K,K) NORM 10 DIMENSION AA(MNWT,K), DIAG(10), LPAS(10), V(10) NORM 11 NWT(1) = N1 NORM 12 NWT(2) = N2 NORM 13 C NORMALIZE STIM VECTORS SO SSQ OF EACH DIM EQ 1 NORM 14 XN1 = NWT(1) NORM 15 I = 2 NORM 16 DO 30 J=1,K NORM 17 SUM = 0. NORM 18 DO 10 M=1,N2 NORM 19 SUM = SUM + A(M,J,I)**2 NORM 20 10 CONTINUE NORM 21 V(J) = SUM NORM 22 SUM = SQRT(SUM) NORM 23 DO 20 M=1,N2 NORM 24 A(M,J,I) = A(M,J,I)/SUM NORM 25 20 CONTINUE NORM 26 30 CONTINUE NORM 27 C MULT. ELEMENTS OF WEIGHTS VECTOR BY SSQ OF CORRES. DIM NORM 28 DO 50 J=1,K NORM 29 DO 40 M=1,N1 NORM 30 A(M,J,1) = A(M,J,1)*V(J) NORM 31 40 CONTINUE NORM 32 50 CONTINUE NORM 33 C COMPUTE SUM SQUARES OF WEIGHTS FOR EACH DIM. NORM 34 DO 70 I=1,K NORM 35 LPAS(I) = I NORM 36 DIAG(I) = 0. NORM 37 DO 60 J=1,N1 NORM 38 DIAG(I) = DIAG(I) + A(J,I,1)**2 NORM 39 60 CONTINUE NORM 40 70 CONTINUE NORM 41 CALL SORT(DIAG, K, LPAS, -1) NORM 42 C PERMUTE DIMENSIONS ACCORDING TO SUM SQUARES IN DESCENDING ORDER NORM 43 DO 120 I=1,NA NORM 44 NN = NWT(I) NORM 45 DO 90 J=1,K NORM 46 DO 80 M=1,NN NORM 47 AA(M,J) = A(M,J,I) NORM 48 80 CONTINUE NORM 49 90 CONTINUE NORM 50 DO 110 J=1,K NORM 51 JJ = LPAS(J) NORM 52 DO 100 M=1,NN NORM 53 A(M,J,I) = AA(M,JJ) NORM 54 100 CONTINUE NORM 55 110 CONTINUE NORM 56 120 CONTINUE NORM 57 C PRINT AND PUNCH NORMALIZED SOLUTION NORM 58 WRITE (LPRINT,9999) NORM 59 IF (IPUN.GE.0) WRITE (LPUNCH,9994) NORM 60 DO 190 I=1,NA NORM 61 NN = NWT(I) NORM 62 IF (IPUN.LT.0) GO TO 140 NORM 63 DO 130 J=1,K NORM 64 WRITE (LPUNCH,9993) J, I, (A(M,J,I),M=1,NN) NORM 65 130 CONTINUE NORM 66 140 GO TO (150, 160), I NORM 67 150 WRITE (LPRINT,9996) NORM 68 GO TO 170 NORM 69 160 WRITE (LPRINT,9992) NORM 70 GO TO 170 NORM 71 170 DO 180 J=1,K NORM 72 WRITE (LPRINT,9997) J, (A(M,J,I),M=1,NN) NORM 73 180 CONTINUE NORM 74 190 CONTINUE NORM 75 C COMPUTE SUMS OF PRODUCTS AND SUMS OF SQUARES OF A MATRIX NORM 76 DO 290 I=1,NA NORM 77 IF (I.GT.1) WRITE (LPRINT,9998) NORM 78 NN = NWT(I) NORM 79 DO 220 L=1,K NORM 80 DO 210 J=1,K NORM 81 S(L,J) = 0. NORM 82 DO 200 M=1,NN NORM 83 S(L,J) = S(L,J) + A(M,L,I)*A(M,J,I) NORM 84 200 CONTINUE NORM 85 210 CONTINUE NORM 86 220 CONTINUE NORM 87 SUM = 0. NORM 88 DO 240 L=1,K NORM 89 IF (I.EQ.1) GO TO 230 NORM 90 WRITE (LPRINT,9997) L, (S(L,J),J=1,L) NORM 91 GO TO 240 NORM 92 230 V(L) = S(L,L) NORM 93 SUM = SUM + V(L) NORM 94 240 CONTINUE NORM 95 IF (I.GT.1) GO TO 290 NORM 96 DO 260 II=1,K NORM 97 DO 250 JJ=1,II NORM 98 S(II,JJ) = S(II,JJ)/SQRT(V(II)*V(JJ)) NORM 99 250 CONTINUE NORM 100 260 CONTINUE NORM 101 WRITE (LPRINT,9990) NORM 102 DO 270 L=1,K NORM 103 WRITE (LPRINT,9997) L, (S(L,J),J=1,L) NORM 104 270 CONTINUE NORM 105 DO 280 II=1,K NORM 106 V(II) = (V(II)*FVAF)/SUM NORM 107 280 CONTINUE NORM 108 290 CONTINUE NORM 109 DO 300 II=1,K NORM 110 LPAS(II) = II NORM 111 300 CONTINUE NORM 112 WRITE (LPRINT,9991) (LPAS(II),II=1,K) NORM 113 WRITE (LPRINT,9995) (V(II),II=1,K) NORM 114 9999 FORMAT (1H0, 23X, 19HNORMALIZED SOLUTION) NORM 115 9998 FORMAT (26H0SUM OF PRODUCTS (STIMULI)) NORM 116 9997 FORMAT (I3, 3X, 10F7.3/(6X, 10F7.3)) NORM 117 9996 FORMAT (23H0SUBJECTS WEIGHT MATRIX) NORM 118 9995 FORMAT (6X, 10F7.3) NORM 119 9994 FORMAT (35H(2X,5E13.5) NORMALIZED SOLUTION) NORM 120 9993 FORMAT (2I1, 5E13.5/(2X, 5E13.5)) NORM 121 9992 FORMAT (16H0STIMULUS MATRIX) NORM 122 9991 FORMAT (46H0APPROXIMATE PROPORTION OF TOTAL VARIANCE ACCO, NORM 123 * 9HUNTED FOR, 18H BY EACH DIMENSION/6X, 10(I4, 3X)) NORM 124 9990 FORMAT (38H0NORMALIZED SUM OF PRODUCTS (SUBJECTS)) NORM 125 RETURN NORM 126 END NORM 127 C SUBR 1 SUBROUTINE SUBR(A, MNWT, NF, N, Y, N1, N2, LSD) SUBR 2 C SUBR FOR SINDSCAL SEPTEMBER, 1975 SUBR 3 C WRITTEN BY J.J. CHANG SUBR 4 C MODIFIED FOR SINDSCAL BY S. PRUZANSKY SUBR 5 COMMON /DEVICE/ LREAD, LPRINT, LPUNCH SUBR 6 DIMENSION A(MNWT,NF,N), Y(N1,LSD) SUBR 7 XN = N2*N2 SUBR 8 WRITE (LPRINT,9999) SUBR 9 DO 50 I1=1,N1 SUBR 10 SUMX = 0. SUBR 11 SUMY = 0. SUBR 12 SSQX = 0. SUBR 13 SSQY = 0. SUBR 14 SXY = 0. SUBR 15 LC = 0 SUBR 16 DO 40 I2=1,N2 SUBR 17 DO 30 I3=1,I2 SUBR 18 KE = 2 SUBR 19 LC = LC + 1 SUBR 20 SUM = 0 SUBR 21 YY = Y(I1,LC) SUBR 22 DO 10 I=1,NF SUBR 23 SUM = SUM + A(I1,I,1)*A(I2,I,2)*A(I3,I,2) SUBR 24 10 CONTINUE SUBR 25 IF (I2.EQ.I3) KE = 1 SUBR 26 DO 20 I=1,KE SUBR 27 SUMX = SUMX + SUM SUBR 28 SUMY = SUMY + YY SUBR 29 SSQX = SSQX + SUM**2 SUBR 30 SSQY = SSQY + YY**2 SUBR 31 SXY = SXY + SUM*YY SUBR 32 20 CONTINUE SUBR 33 30 CONTINUE SUBR 34 40 CONTINUE SUBR 35 R = (XN*SXY-SUMX*SUMY)/SQRT((XN*SSQX-SUMX**2)*(XN* SUBR 36 * SSQY-SUMY**2)) SUBR 37 WRITE (LPRINT,9998) I1, R SUBR 38 50 CONTINUE SUBR 39 9999 FORMAT (46H0CORRELATION BETWEEN COMPUTED SCORES AND SCALA, SUBR 40 * 8HR PROD. , 12HFOR SUBJECTS/1H0) SUBR 41 9998 FORMAT (I4, 3X, F10.6) SUBR 42 RETURN SUBR 43 END SUBR 44 C DIMA 1 SUBROUTINE DIMA(A, MNWT, NQ, NFO, NF) DIMA 2 C DIMA FOR SINDSCAL SEPTEMBER, 1975 DIMA 3 C WRITTEN BY S. PRUZANSKY DIMA 4 DIMENSION A(MNWT,NQ) DIMA 5 NN = NFO DIMA 6 L = NF DIMA 7 DO 20 I=1,NF DIMA 8 L = L + 1 DIMA 9 NN = NN + 1 DIMA 10 DO 10 J=1,MNWT DIMA 11 A(J,L) = A(J,NN) DIMA 12 10 CONTINUE DIMA 13 20 CONTINUE DIMA 14 RETURN DIMA 15 END DIMA 16 C SMTI 1 SUBROUTINE SMTINV(A, UL, N, AINV, IER) SMTI 2 C-SMTINV-------ADAPTED FROM IMSL LIBRARY 2--------------------------- SMTI 3 C SMTINV MERGES THE TWO IMSL ROUTINES LINV1P AND LUDECP. SMTI 4 C THESE IMSL ROUTINES ARE PROPRIETARY SOFTWARE, OWNED BY IMSL. SMTI 5 C THEY MAY BE USED ONLY IN THE CODE IN WHICH THEY ARE EMBEDDED. SMTI 6 C NO OTHER USE OF THESE ROUTINES IS PERMITTED. SMTI 7 C SMTI 8 C FUNCTION - INVERSION OF A POSITIVE DEFINITE SYMMETRIC SMTI 9 C MATRIX - SYMMETRIC STORAGE MODE - SPACE SMTI 10 C ECONOMIZER SOLUTION SMTI 11 C USAGE - CALL SMTINV(A,UL,N,AINV,IER) SMTI 12 C PARAMETERS A - INPUT VECTOR OF LENGTH N(N+1)/2 CONTAINING SMTI 13 C THE N BY N POSITIVE DEFINITE SYMMETRIC SMTI 14 C MATRIX TO BE INVERTED. A IS STORED IN SMTI 15 C SYMMETRIC STORAGE MODE. SMTI 16 C UL ON OUTPUT, A IS REPLACED BY THE LOWER SMTI 17 C TRIANGULAR MATRIX L WHERE A = L*L-TRANSPOSE.SMTI 18 C L IS STORED IN SYMMETRIC STORAGE MODE WITH SMTI 19 C THE DIAGONAL ELEMENTS OF L STORED IN SMTI 20 C RECIPROCAL FORM. SMTI 21 C N - ORDER OF A.(INPUT) SMTI 22 C AINV - OUTPUT VECTOR OF LENGTH N(N+1)/2 CONTAINING SMTI 23 C THE N BY N INVERSE OF A. AINV IS STORED IN SMTI 24 C SYMMETRIC STORAGE MODE. SMTI 25 C IER - ERROR PARAMETER. SMTI 26 C TERMINAL ERROR = 128+N. SMTI 27 C PRECISION - SINGLE SMTI 28 C REQ'D IMSL ROUTINES - LUELMP SMTI 29 C-----------------------------------------------------------------------SMTI 30 DIMENSION A(1), AINV(1), UL(1) SMTI 31 DOUBLE PRECISION X SMTI 32 DATA ZERO, ONE, FOUR, SIXTN, SIXTH /0.0,1.,4.,16.,.0625/ SMTI 33 IER = 0 SMTI 34 L = 1 SMTI 35 N1 = N - 1 SMTI 36 LN = N SMTI 37 C DECOMPOSE A SMTI 38 D1 = ONE SMTI 39 D2 = ZERO SMTI 40 IP = 1 SMTI 41 IER = 0 SMTI 42 DO 90 I=1,N SMTI 43 IQ = IP SMTI 44 IR = 1 SMTI 45 DO 80 J=1,I SMTI 46 X = A(IP) SMTI 47 IF (J.EQ.1) GO TO 20 SMTI 48 DO 10 K=IQ,IP1 SMTI 49 X = X - DBLE(UL(K))*DBLE(UL(IR)) SMTI 50 IR = IR + 1 SMTI 51 10 CONTINUE SMTI 52 20 IF (I.NE.J) GO TO 60 SMTI 53 D1 = D1*X SMTI 54 IF (X.LE.0.D0) GO TO 100 SMTI 55 30 IF (ABS(D1).LT.ONE) GO TO 40 SMTI 56 D1 = D1*SIXTH SMTI 57 D2 = D2 + FOUR SMTI 58 GO TO 30 SMTI 59 40 IF (ABS(D1).GE.SIXTH) GO TO 50 SMTI 60 D1 = D1*SIXTN SMTI 61 D2 = D2 - FOUR SMTI 62 GO TO 40 SMTI 63 50 UL(IP) = 1.D0/DSQRT(X) SMTI 64 GO TO 70 SMTI 65 60 UL(IP) = X*UL(IR) SMTI 66 70 IP1 = IP SMTI 67 IP = IP + 1 SMTI 68 IR = IR + 1 SMTI 69 80 CONTINUE SMTI 70 90 CONTINUE SMTI 71 GO TO 110 SMTI 72 100 IER = 129 SMTI 73 110 CONTINUE SMTI 74 IF (IER.NE.0) GO TO 140 SMTI 75 DO 130 I=1,N SMTI 76 DO 120 J=L,LN SMTI 77 AINV(J) = ZERO SMTI 78 120 CONTINUE SMTI 79 LI = L + I - 1 SMTI 80 AINV(LI) = ONE SMTI 81 C OBTAIN THE SOLUTION SMTI 82 CALL LUELMP(UL, AINV(L), N, AINV(L)) SMTI 83 L = L + I SMTI 84 LN = L + N1 SMTI 85 130 CONTINUE SMTI 86 140 RETURN SMTI 87 END SMTI 88 C LUEL 1 SUBROUTINE LUELMP(A, B, N, X) LUEL 2 C-LUELMP--------FROM IMSL LIBRARY 2---------------------------------- LUEL 3 C THE IMSL ROUTINE LUELMP IS PROPRIETARY SOFTWARE, OWNED BY IMSL. LUEL 4 C IT MAY BE USED ONLY IN THE CODE IN WHICH IT IS EMBEDDED. LUEL 5 C NO OTHER USE OF THIS ROUTINE IS PERMITTED. LUEL 6 C LUEL 7 C FUNCTION - ELIMINATION PART OF THE SOLUTION OF AX=B - LUEL 8 C SYMMETRIC STORAGE MODE LUEL 9 C USAGE - CALL LUELMP (A,B,N,X) LUEL 10 C PARAMETERS A - INPUT VECTOR OF LENGTH N(N+1)/2 CONTAINING LUEL 11 C THE N BY N MATRIX L WHERE A = L*L-TRANSPOSE.LUEL 12 C L IS A LOWER TRIANGULAR MATRIX STORED IN LUEL 13 C SYMMETRIC STORAGE MODE. THE MAIN DIAGONAL LUEL 14 C ELEMENTS OF L ARE STORED IN RECIPROCAL LUEL 15 C FORM. MATRIX L MAY BE OBTAINED FROM IMSL LUEL 16 C ROUTINE LUDECP. LUEL 17 C B - VECTOR OF LENGTH N CONTAINING THE RIGHT HAND LUEL 18 C SIDE OF THE EQUATION AX = B. (INPUT) LUEL 19 C N - ORDER OF A AND THE LENGTH OF B AND X. (INPUT) LUEL 20 C X - VECTOR OF LENGTH N CONTAINING THE SOLUTION TO LUEL 21 C THE EQUATION AX = B. (OUTPUT) LUEL 22 C PRECISION - SINGLE LUEL 23 C LANGUAGE - FORTRAN LUEL 24 C-----------------------------------------------------------------------LUEL 25 C LATEST REVISION - FEBRUARY 12,1973 LUEL 26 C LUEL 27 DIMENSION A(1), B(1), X(1) LUEL 28 DOUBLE PRECISION T LUEL 29 C SOLUTION OF LY = B LUEL 30 IP = 1 LUEL 31 DO 30 I=1,N LUEL 32 IQ = I - 1 LUEL 33 T = B(I) LUEL 34 IF (IQ.LT.1) GO TO 20 LUEL 35 DO 10 K=1,IQ LUEL 36 T = T - DBLE(A(IP))*DBLE(X(K)) LUEL 37 IP = IP + 1 LUEL 38 10 CONTINUE LUEL 39 20 X(I) = T*A(IP) LUEL 40 IP = IP + 1 LUEL 41 30 CONTINUE LUEL 42 C SOLUTION OF UX = Y LUEL 43 N1 = N + 1 LUEL 44 DO 60 I=1,N LUEL 45 II = N1 - I LUEL 46 IP = IP - 1 LUEL 47 IS = IP LUEL 48 IQ = II + 1 LUEL 49 T = X(II) LUEL 50 IF (N.LT.IQ) GO TO 50 LUEL 51 KK = N LUEL 52 DO 40 K=IQ,N LUEL 53 T = T - DBLE(A(IS))*DBLE(X(KK)) LUEL 54 KK = KK - 1 LUEL 55 IS = IS - KK LUEL 56 40 CONTINUE LUEL 57 50 X(II) = T*A(IS) LUEL 58 60 CONTINUE LUEL 59 RETURN LUEL 60 END LUEL 61 C SORT 1 CSORT SORT 2 SUBROUTINE SORT(A, N, LPAS, SWITCH) SORT 3 C SORT FOR KYST JANUARY,1973 SORT 4 C WRITTEN BY J.KRUSKAL MODIFIED FOR SINDSCAL 8-18-77 BY SP SORT 5 C SORT 6 C THIS ROUTINE SORTS INPUT ARRAY 'A' AND REARRANGES, OPTIONALLY, SORT 7 C ARRAY LPAS, IN ORDER CORRESPONDING TO 'A'. SORT 8 C N = NUMBER OF ITEMS IN 'A' (AND LPAS , IF USED) SORT 9 C IF 'SWITCH' IS POSITIVE, SORT WILL BE IN ASCENDING ORDER, SORT 10 C IF ZERO OR NEGATIVE, IN DESCENDING ORDER. SORT 11 C ALGORITHM FROM CACM, JULY 1959, PAGE 30 BY D. L. SHELL SORT 12 C SORT 13 DIMENSION A(1), LPAS(1) SORT 14 INTEGER SWITCH SORT 15 IF (N.LE.1) GO TO 90 SORT 16 M = 1 SORT 17 10 M = M + M SORT 18 IF (M.LE.N) GO TO 10 SORT 19 M = M - 1 SORT 20 20 M = M/2 SORT 21 IF (M.EQ.0) GO TO 90 SORT 22 KK = N - M SORT 23 J = 1 SORT 24 30 I = J SORT 25 40 IM = I + M SORT 26 IF (SWITCH) 60, 60, 50 SORT 27 50 IF (A(I).GT.A(IM)) GO TO 80 SORT 28 GO TO 70 SORT 29 60 IF (A(I).LT.A(IM)) GO TO 80 SORT 30 70 J = J + 1 SORT 31 IF (J.GT.KK) GO TO 20 SORT 32 GO TO 30 SORT 33 80 TEMP = A(I) SORT 34 A(I) = A(IM) SORT 35 A(IM) = TEMP SORT 36 ITEMP = LPAS(I) SORT 37 LPAS(I) = LPAS(IM) SORT 38 LPAS(IM) = ITEMP SORT 39 I = I - M SORT 40 IF (I.LT.1) GO TO 70 SORT 41 GO TO 40 SORT 42 90 RETURN SORT 43 END SORT 44 C RCON 1 SUBROUTINE RCON(A, MNWT, NF, N2, IRN) RCON 2 C RCON FOR SINDSCAL SEPTEMBER, 1975 RCON 3 C WRITTEN BY S. PRUZANSKY RCON 4 C SUBROUTINE TO GENERATE RANDOM STARTING CONFIGURATION RCON 5 C USES PORTABLE RANDOM NUMBER GENERATOR RCON 6 DIMENSION A(MNWT,NF,3) RCON 7 IF (IRN.LT.100) IRN = IRN*11 RCON 8 DO 10 I=1,IRN RCON 9 TEMP = RUNIFV(DUM) RCON 10 10 CONTINUE RCON 11 DO 40 L=2,3 RCON 12 DO 30 K=1,NF RCON 13 DO 20 J=1,N2 RCON 14 TEMP = RUNIFV(DUM) RCON 15 A(J,K,L) = TEMP - .5 RCON 16 20 CONTINUE RCON 17 30 CONTINUE RCON 18 40 CONTINUE RCON 19 RETURN RCON 20 END RCON 21 C RUNI 1 FUNCTION RUNIFV(A) RUNI 2 C RUNIFV FOR KYST JANUARY,1973 RUNI 3 C WRITTEN BY J.KRUSKAL RUNI 4 C THIS IS A SIMPLE ROUTINE FOR GENERATING RANDOM NUMBERS. RUNI 5 C THEY ARE UNIFORMLY DISTRIBUTED BETWEEN 0 AND 1. RUNI 6 C IT IS NOT FAST, NOR DOES IT PRODUCE VERY HIGH QUALITY NUMBERS. RUNI 7 C ITS MAIN VIRTUE IS THAT IT IS IN FORTRAN AND WILL WORK ON RUNI 8 C ALMOST ANY MACHINE ON WHICH FORTRAN IS AVAILABLE. RUNI 9 C RUNI 10 DATA MODULO, FLMOD, K /8192,8192.0,1/ RUNI 11 C MODULO 2**13 RUNI 12 C RUNI 13 DO 10 I=1,15 RUNI 14 K = MOD(5*K,MODULO) RUNI 15 10 CONTINUE RUNI 16 Z = FLOAT(K)/FLMOD RUNI 17 RUNIFV = Z RUNI 18 RETURN RUNI 19 END RUNI 20 C PPLO 1 SUBROUTINE PPLOT(X, MNWT, NF, TIT, N) PPLO 2 C PPLOT FOR SINDSCAL SEPTEMBER, 1975 PPLO 3 C WRITTEN BY S. PRUZANSKY PPLO 4 C SUBROUTINE TO PLOT CONFIGURATION ON PRINTER PPLO 5 C REQUIRES 120 CHARACTERS PER LINE PPLO 6 C WRITTEN 2-6-75 BY SP PPLO 7 DIMENSION X(MNWT,NF,3), N(1) PPLO 8 INTEGER TIT(1) PPLO 9 COMMON /DEVICE/ LREAD, LPRINT, LPUNCH PPLO 10 DO 30 L1=1,2 PPLO 11 DO 20 J=2,NF PPLO 12 J1 = J - 1 PPLO 13 DO 10 K=1,J1 PPLO 14 C DETERMINE MAX. ABSOLUTE VALUE FOR PLANE TO BE PLOTTED PPLO 15 CALL MAXI(X(1,K,L1), X(1,J,L1), XMAX, N(L1)) PPLO 16 XMIN = -XMAX PPLO 17 C PRINT TITLE INFORMATION AT TOP OF PAGE PPLO 18 IF (L1.EQ.1) WRITE (LPRINT,9999) J, K, PPLO 19 * (TIT(II),II=1,80) PPLO 20 IF (L1.EQ.2) WRITE (LPRINT,9998) J, K, PPLO 21 * (TIT(II),II=1,80) PPLO 22 C PLOT A PLANE AT A TIME PPLO 23 CALL PLOT(X(1,K,L1), X(1,J,L1), XMAX, XMAX, XMIN, PPLO 24 * XMIN, N(L1), 1, 2, 66, 1) PPLO 25 10 CONTINUE PPLO 26 20 CONTINUE PPLO 27 30 CONTINUE PPLO 28 9999 FORMAT (26H1WEIGHTS SPACE DIMENSION, I2, 11H (Y-AXIS) V, PPLO 29 * 12HS. DIMENSION, I3, 9H (X-AXIS), /, 30X, 80A1) PPLO 30 9998 FORMAT (27H1STIMULUS SPACE DIMENSION, I2, 10H (Y-AXIS) , PPLO 31 * 13HVS. DIMENSION, I3, 9H (X-AXIS), /, 30X, 80A1) PPLO 32 RETURN PPLO 33 END PPLO 34 C MAXI 1 SUBROUTINE MAXI(X, Y, Z, N) MAXI 2 C MAXI FOR SINDSCAL SEPTEMBER, 1975 MAXI 3 C WRITTEN BY J.J. CHANG MAXI 4 C MAXI FINDS THE MAGNITUDE OF ARRAYS X AND Y EACH OF LENGTH N. MAXI 5 C ON RETURN ABSOLUTE MAXIMUM IS STORED IN Z. MAXI 6 DIMENSION X(1), Y(1) MAXI 7 Z = ABS(X(1)) MAXI 8 DO 10 I=1,N MAXI 9 AXO = ABS(X(I)) MAXI 10 IF (AXO.GT.Z) Z = AXO MAXI 11 AYO = ABS(Y(I)) MAXI 12 IF (AYO.GT.Z) Z = AYO MAXI 13 10 CONTINUE MAXI 14 RETURN MAXI 15 END MAXI 16 C XA00 1 BLOCK DATA XA00 2 C BLDA FOR KYST JANUARY,1973 XA00 3 C WRITTEN BY J.KRUSKAL XA00 4 C MODIFIED FOR SINDSCAL BY S PRUZANSKY JUNE,1975 XA00 5 C XA00 6 REAL PTID(100), ITEM(101) XA00 7 C XA00 8 COMMON /ACCUR/ XMAXN XA00 9 COMMON /PLTCHR/ PTID, ITEM XA00 10 C XA00 11 C XMAXN = MINIMUM OF LARGEST POSITIVE MACHINE NUMBER AND MINUS C LARGEST NEGATIVE MACHINE NUMBER C DATA XMAXN /1.0E38/ XA00 12 DATA PTID /1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK, XA00 13 * 1HL,1HM,1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX, XA00 14 * 1HY,1HZ,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HA,1HB, XA00 15 * 1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO, XA00 16 * 1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H1,1H2, XA00 17 * 1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HA,1HB,1HC,1HD,1HE,1HF, XA00 18 * 1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO,1HP,1HQ,1HR,1HS, XA00 19 * 1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H1,1H2,1H3,1H4/ XA00 20 DATA ITEM /101*1H / XA00 21 END XA00 22 C PLOT 1 SUBROUTINE PLOT(X, Y, XA, YA, XI, YI, NPOI, NY, ID, LONG, PLOT 2 * MX) PLOT 3 C PLOT FOR KYST JANUARY,1973 PLOT 4 C WRITTEN BY FORREST W. YOUNG, NOVEMBER, 1965 PLOT 5 C UPDATED OCTOBER 11,1971 JAY R LEVINSOHN PLOT 6 C ALTERED FOR KYST BY J.KRUSKAL AND J.SEERY JANUARY,1973 PLOT 7 C PLOT 8 C ROUTINE TO GENERATE A ONE PAGE PLOT OF A MATRIX Y VS SOME VECTOR X PLOT 9 C PLOT 10 C Y MAY BE A MATRIX AND IT IS ASSUMED THAT X IS A VECTOR PLOT 11 C PLOT 12 C THE PARAMETERS XA,YA,XI,YI INDICATE THE UPPER PLOT 13 C AND LOWER BOUNDS FOR EACH AXIS. PLOT 14 C THE USER MAY PASS THE MIN AND MAX FOR X AND Y OR HE MAY LET THE PLOT 15 C PROGRAM FIND THEM. IF XA=XI THEN THE PROGRAM WILL FIND THE MAX AND MIPLOT 16 C FOR THE X VECTOR, IF YA=YI THE PROGRAM WILL FIND THE MAX AND MIN PLOT 17 C FOR THE Y MATRIX. PLOT 18 C PLOT 19 C IF -ID- IS POSITIVE, AXES WILL BE INCLUDED ON THE PLOT, PLOT 20 C IF -ID- IS NEGATIVE, NO AXES WILL APPEAR. PLOT 21 C IF ID IS PLUS OR MINUS 2, THE POINTS WILL BE IDENTIFIED BY ROWS PLOT 22 C IF ID IS NOT EQUAL TO 2 THEN POINTS WILL BE LABELED BY COLS. PLOT 23 C PLOT 24 C IF ID IS PLUS OR MINUS THREE THE PLOT IS ASSUMED TO BE 1 DIMENSIONAL. PLOT 25 C ONE DIMENSIONAL PLOTS ARE PLOTTED WITH NO AXES, AND ALL Y VALUES PLOT 26 C FIXED AT MAXY-MINY/2 PLOT 27 C POINTS IN THE 1 DIM. PLOTS ARE IDENDTIFIED BY ROW. PLOT 28 C PLOT 29 C DEFINITIONS PLOT 30 C PARAMETERS PLOT 31 C Y--A MATRIX (LONG,NX) IT CONTAINS VALUES TO BE PLOTTED AGAINST X PLOT 32 C X--A VECTOR (LONG); IT CONTAINS VALUES TO BE PLOTTED AGAINST Y PLOT 33 C XA-- X AXIS UPPER BOUND (MAY BE PASSED AS ZERO) PLOT 34 C YA-- Y AXIS " " " " " PLOT 35 C XI-- X AXIS LOWER BOUND (MAY BE PASSED AS ZERO) PLOT 36 C YI-- Y AXIS LOWER BOUND(MAY BE PASSED AS ZERO ) PLOT 37 C NPOI-- NO. OF ELEMENTS IN X, NO. OF ROWS IN Y PLOT 38 C NY--NUMBER OF COLUMNS IN Y MATRIX PLOT 39 C ID-- PRINT CONTROL PARAMETER, CONTROLS LABELING, AND AXIS GEN. PLOT 40 C LONG-- MAXIMUM NUMBER OF ELEMENTS IN X, MAX. NO. ROWS IN Y PLOT 41 C MX -- WORK WITH ELEMENTS MX THROUGH MX+NPOI-1 PLOT 42 C VARIABLES PLOT 43 C DELX-- X AXIS INCREMENT PLOT 44 C DELY-- Y AXIS INCREMENT PLOT 45 C FINC-- FLOATING INCREMENT FOR XLABEL PLOT 46 C IDAT-- COUNTER OVER ROWS OF X AND Y PLOT 47 C ILAB-- INTEGER LABEL FLAG PLOT 48 C ITEM(101)--ITEMS TO BE PRINTED INITIALLY BLANK PLOT 49 C LINE-- COUNTER OF CURRENT LINE PLOT 50 C NLPP-- NUMBER OF LINES PLOTTED PER PAGE PLOT 51 C SPANX-- RANGE OF X PLOT 52 C SPANY-- RANGE OF Y PLOT 53 C VALUE-- DATUM PRINTED ON Y AXIS PLOT 54 C XLABEL(11)--X-AXIS LABELS PLOT 55 C XNOW-- CURRENT VALUE OF X VECTOR NOW BEING WORKED ON PLOT 56 C YNOW-- CURRENT VALUE OF Y MATRIX NOW BEING WORKED ON PLOT 57 C PLOT 58 REAL ITEM(101), MINUS, PTID(100), YMINV(2), YMAXV(2) PLOT 59 DIMENSION XLABL(11), X(1), Y(LONG,NY), XS(10), YS(10), PLOT 60 * IND(2) PLOT 61 COMMON /DEVICE/ LREAD, LPRINT, LPUNCH PLOT 62 COMMON /ACCUR/ XMAXN PLOT 63 COMMON /PLTCHR/ PTID, ITEM PLOT 64 DATA MINUS, DD, STAR, BLANK, AIE /1H-,1H.,1H*,1H ,1HI/ PLOT 65 DATA XS /0.,0.,0.,0.,0.,4.,2.,4.,1.,0./ PLOT 66 DATA YS /0.,0.,0.,0.,2.,4.,3.,4.,7.,2./ PLOT 67 C PLOT 68 JY = 1 PLOT 69 KODE = 1 PLOT 70 INDEX = 0 PLOT 71 LINE = 1 PLOT 72 IY0L = 0 PLOT 73 IX0L = 0 PLOT 74 IDAT = 1 PLOT 75 NLPP = 53 PLOT 76 XMAX = XA PLOT 77 YMAX = YA PLOT 78 YMIN = YI PLOT 79 XMIN = XI PLOT 80 NPOI2 = NY*NPOI PLOT 81 NPOIRL = MX + NPOI - 1 PLOT 82 LXFLG = 0 PLOT 83 LYFLG = 0 PLOT 84 C PLOT 85 C DETERMINE MINIMUM AND MAXIMUM OF X AXIS, IF NECESSARY PLOT 86 C PLOT 87 IF (XMAX-XMIN) 60, 10, 60 PLOT 88 10 XMAX = -XMAXN PLOT 89 XMIN = XMAXN PLOT 90 DO 50 I=MX,NPOIRL PLOT 91 IF (X(I)-XMAX) 30, 20, 20 PLOT 92 20 XMAX = X(I) PLOT 93 30 IF (X(I)-XMIN) 40, 40, 50 PLOT 94 40 XMIN = X(I) PLOT 95 50 CONTINUE PLOT 96 C PLOT 97 C CALL SERCH TO INITIALIZE; IF DESIRED SET YMAX AND YMIN PLOT 98 C PLOT 99 C PLOT 100 C IF WE HAVE A 1D PLOT NO NEED FOR SERCH PLOT 101 C PLOT 102 IF (IABS(ID)-3) 60, 100, 60 PLOT 103 60 DO 70 I=1,NY PLOT 104 IY = I PLOT 105 CALL SERCH(Y(1,IY), NPOIRL, INDEX, KODE, IY, MX) PLOT 106 YMAXV(I) = Y(INDEX,I) PLOT 107 YMINV(I) = Y(KODE,I) PLOT 108 IND(I) = INDEX PLOT 109 IF (NY.EQ.1) GO TO 80 PLOT 110 KODE = 1 PLOT 111 70 CONTINUE PLOT 112 YMAX = AMAX1(YMAXV(1),YMAXV(2)) PLOT 113 YMIN = AMIN1(YMINV(1),YMINV(2)) PLOT 114 IF (YMAX.NE.YMAXV(1)) JY = 2 PLOT 115 YNOW = YMAX PLOT 116 YNOW1 = AMIN1(YMAXV(1),YMAXV(2)) PLOT 117 GO TO 100 PLOT 118 80 YNOW = Y(INDEX,1) PLOT 119 YNOW1 = -XMAXN PLOT 120 IF (YMAX-YMIN) 100, 90, 100 PLOT 121 90 YMAX = Y(INDEX,1) PLOT 122 YMIN = Y(KODE,1) PLOT 123 C PLOT 124 C DETERMINE RANGE AND INCREMENT OF BOTH AXES PLOT 125 100 P = ALOG10(XMAX-XMIN+.0001) PLOT 126 IP = IFIX(P) PLOT 127 IF (P.LT.0.) IP = IP - 1 PLOT 128 SPAN = (XMAX-XMIN)/(10.**IP) PLOT 129 IF (SPAN.LT.5.0) GO TO 110 PLOT 130 DEL = 10.**IP PLOT 131 GO TO 130 PLOT 132 110 IF (SPAN.LT.2.0) GO TO 120 PLOT 133 DEL = (10.**IP)/2.0 PLOT 134 GO TO 130 PLOT 135 120 DEL = (10.**IP)/5.0 PLOT 136 130 IF (XMAX.LT.0.) GO TO 160 PLOT 137 XMAX = FLOAT(IFIX(XMAX/DEL+.9999))*DEL PLOT 138 IF (XMIN) 150, 140, 140 PLOT 139 140 XMIN = FLOAT(IFIX(XMIN/DEL))*DEL PLOT 140 GO TO 170 PLOT 141 150 XMIN = FLOAT(IFIX(XMIN/DEL-.9999))*DEL PLOT 142 GO TO 170 PLOT 143 160 XMAX = FLOAT(IFIX(XMAX/DEL))*DEL PLOT 144 XMIN = FLOAT(IFIX(XMIN/DEL+.9999))*DEL PLOT 145 170 SPANX = XMAX - XMIN PLOT 146 NXUNIT = SPANX/DEL + .00001 PLOT 147 IF (NXUNIT.LE.10) GO TO 180 PLOT 148 LXFLG = LXFLG + 1 PLOT 149 IF (LXFLG-1) 830, 100, 830 PLOT 150 180 NXP1 = NXUNIT + 1 PLOT 151 FINC = DEL PLOT 152 DELX = SPANX/(100.-XS(NXUNIT)) PLOT 153 IXS = IFIX(XS(NXUNIT)/2.+.1) PLOT 154 XSF = IXS PLOT 155 IF (ID.LT.0) GO TO 190 PLOT 156 IF ((XMIN.GE.0.) .OR. (XMAX.LE.0.)) GO TO 190 PLOT 157 IX0L = -XMIN/DELX + 1.5 PLOT 158 IX0L = IX0L + IXS PLOT 159 C PLOT 160 190 IF (ID.EQ.3) GO TO 280 PLOT 161 P = ALOG10(YMAX-YMIN+.0001) PLOT 162 IP = IFIX(P) PLOT 163 IF (P.LT.0.) IP = IP - 1 PLOT 164 SPAN = (YMAX-YMIN)/(10.**IP) PLOT 165 IF (SPAN.LT.5.0) GO TO 200 PLOT 166 DEL = 10.**IP PLOT 167 GO TO 220 PLOT 168 200 IF (SPAN.LT.2.0) GO TO 210 PLOT 169 DEL = (10.**IP)/2.0 PLOT 170 GO TO 220 PLOT 171 210 DEL = (10.**IP)/5.0 PLOT 172 220 IF (YMAX.LT.0.) GO TO 250 PLOT 173 YMAX = FLOAT(IFIX(YMAX/DEL+.9999))*DEL PLOT 174 IF (YMIN) 240, 230, 230 PLOT 175 230 YMIN = FLOAT(IFIX(YMIN/DEL))*DEL PLOT 176 GO TO 260 PLOT 177 240 YMIN = FLOAT(IFIX(YMIN/DEL-.9999))*DEL PLOT 178 GO TO 260 PLOT 179 250 YMAX = FLOAT(IFIX(YMAX/DEL))*DEL PLOT 180 YMIN = FLOAT(IFIX(YMIN/DEL+.9999))*DEL PLOT 181 260 SPANY = YMAX - YMIN PLOT 182 NYUNIT = SPANY/DEL + .00001 PLOT 183 IF (NYUNIT.LE.10) GO TO 270 PLOT 184 LYFLG = LYFLG + 1 PLOT 185 IF (LYFLG-1) 830, 190, 830 PLOT 186 270 DELY = SPANY/(FLOAT(NLPP-1)-YS(NYUNIT)) PLOT 187 DELTY = DEL PLOT 188 VALUE = YMAX PLOT 189 YININC = (NLPP-1)/NYUNIT PLOT 190 LININC = IFIX(YININC) PLOT 191 IYLAB = LININC PLOT 192 IF (ID.LT.0) GO TO 290 PLOT 193 IF ((YMAX.LE.0.) .OR. (YMIN.GE.0.)) GO TO 290 PLOT 194 IY0L = YMAX/DELY + 1.5 PLOT 195 GO TO 290 PLOT 196 280 VALUE = 0.0 PLOT 197 DELTY = 0.0 PLOT 198 LININC = 27 PLOT 199 IYLAB = 1 PLOT 200 C PLOT 201 C LABEL PLOT AT TOP PLOT 202 C PLOT 203 290 XLABL(1) = XMIN PLOT 204 DO 300 I=1,NXUNIT PLOT 205 XLABL(I+1) = XLABL(I) + FINC PLOT 206 300 CONTINUE PLOT 207 GO TO (830, 830, 830, 310, 320, 330, 340, 350, 360, 370), PLOT 208 * NXUNIT PLOT 209 310 WRITE (LPRINT,9999) (XLABL(I),I=1,NXP1) PLOT 210 WRITE (LPRINT,9998) PLOT 211 9999 FORMAT (1H0, F20.4, 4F25.4) PLOT 212 9998 FORMAT (1H , 14X, 36H*.************************.*********, PLOT 213 * 50H***************.************************.*********, PLOT 214 * 17H***************.*) PLOT 215 GO TO 380 PLOT 216 320 WRITE (LPRINT,9997) (XLABL(I),I=1,NXP1) PLOT 217 WRITE (LPRINT,9996) PLOT 218 9997 FORMAT (1H0, 6F20.4) PLOT 219 9996 FORMAT (1H , 14X, 36H*.*******************.**************, PLOT 220 * 50H*****.*******************.*******************.****, PLOT 221 * 17H***************.*) PLOT 222 GO TO 380 PLOT 223 330 WRITE (LPRINT,9995) (XLABL(I),I=1,NXP1) PLOT 224 WRITE (LPRINT,9994) PLOT 225 9995 FORMAT (1H0, 6X, 7F16.4) PLOT 226 9994 FORMAT (1H , 14X, 36H***.***************.***************., PLOT 227 * 50H***************.***************.***************.**, PLOT 228 * 17H*************.***) PLOT 229 GO TO 380 PLOT 230 340 WRITE (LPRINT,9993) (XLABL(I),I=1,NXP1) PLOT 231 WRITE (LPRINT,9992) PLOT 232 9993 FORMAT (1H0, 7X, 8F14.4) PLOT 233 9992 FORMAT (1H , 14X, 36H**.*************.*************.*****, PLOT 234 * 50H********.*************.*************.*************, PLOT 235 * 17H.*************.**) PLOT 236 GO TO 380 PLOT 237 350 WRITE (LPRINT,9991) (XLABL(I),I=1,NXP1) PLOT 238 WRITE (LPRINT,9990) PLOT 239 9991 FORMAT (1H0, 10X, 9F12.4) PLOT 240 9990 FORMAT (1H , 14X, 36H***.***********.***********.********, PLOT 241 * 50H***.***********.***********.***********.**********, PLOT 242 * 17H*.***********.***) PLOT 243 GO TO 380 PLOT 244 360 WRITE (LPRINT,9989) (XLABL(I),I=1,NXP1) PLOT 245 WRITE (LPRINT,9988) PLOT 246 9989 FORMAT (1H0, 8X, 10F11.3) PLOT 247 9988 FORMAT (1H , 14X, 36H*.**********.**********.**********.*, PLOT 248 * 50H*********.**********.**********.**********.*******, PLOT 249 * 17H***.**********.**) PLOT 250 GO TO 380 PLOT 251 370 WRITE (LPRINT,9987) (XLABL(I),I=1,NXP1) PLOT 252 WRITE (LPRINT,9986) PLOT 253 9987 FORMAT (1H0, 9X, 11F10.3) PLOT 254 9986 FORMAT (1H , 14X, 36H*.*********.*********.*********.****, PLOT 255 * 50H*****.*********.*********.*********.*********.****, PLOT 256 * 17H*****.*********.*) PLOT 257 380 CONTINUE PLOT 258 C PLOT 259 C PREPARE PLOT PLOT 260 C PLOT 261 KODE = 0 PLOT 262 390 IF (LINE-NLPP) 400, 400, 740 PLOT 263 400 IF (IDAT-1) 410, 470, 410 PLOT 264 410 IF (IABS(ID)-3) 420, 470, 420 PLOT 265 420 IF (LFLAG) 430, 430, 470 PLOT 266 430 IF (IDAT-NPOI2) 440, 440, 470 PLOT 267 440 CALL SERCH(Y(1,JY), NPOIRL, INDEX, KODE, JY, MX) PLOT 268 IF (INDEX.EQ.(-1)) GO TO 450 PLOT 269 IND(JY) = INDEX PLOT 270 YNOW = Y(INDEX,JY) PLOT 271 IF (YNOW-YNOW1) 460, 470, 470 PLOT 272 450 JY = NY - JY + 1 PLOT 273 YNOW = YNOW1 PLOT 274 YNOW1 = -XMAXN PLOT 275 GO TO 470 PLOT 276 460 JY = NY - JY + 1 PLOT 277 TEMP = YNOW PLOT 278 YNOW = YNOW1 PLOT 279 YNOW1 = TEMP PLOT 280 470 CONTINUE PLOT 281 C PLOT 282 C CHECK TO SEE IF ANY POINTS IN Y REMAIN PLOT 283 C PLOT 284 IF (IDAT-NPOI2) 480, 480, 600 PLOT 285 480 IF (IABS(ID)-3) 490, 520, 490 PLOT 286 490 IF (YMAX-YNOW) 570, 500, 500 PLOT 287 500 IF (YNOW-YMIN) 570, 510, 510 PLOT 288 510 I = (YMAX-YNOW)/DELY + 1.50001 PLOT 289 GO TO 550 PLOT 290 520 IF (LINE-27) 550, 530, 550 PLOT 291 530 DO 540 II=MX,NPOIRL PLOT 292 I = II PLOT 293 CALL XITEM(I, XMAX, XMIN, ID, DELX, X, JY, XSF) PLOT 294 540 CONTINUE PLOT 295 GO TO 590 PLOT 296 550 IF (I-LINE) 580, 580, 560 PLOT 297 560 LFLAG = 1 PLOT 298 GO TO 600 PLOT 299 570 LFLAG = 0 PLOT 300 GO TO 590 PLOT 301 580 LFLAG = 0 PLOT 302 INDEX = IND(JY) PLOT 303 CALL XITEM(INDEX, XMAX, XMIN, ID, DELX, X, JY, XSF) PLOT 304 590 IDAT = IDAT + 1 PLOT 305 GO TO 390 PLOT 306 C PLOT 307 C PRINT ONE LINE OF PLOT PLOT 308 C PLOT 309 C CHECK ID -- IF MINUS NO AXES,IF + AXES PLOT 310 C IF ID=3 1D PLOT PLOT 311 C PLOT 312 600 IF (ID) 700, 610, 610 PLOT 313 610 IF (IABS(ID)-3) 620, 700, 620 PLOT 314 620 IF (IX0L) 650, 650, 630 PLOT 315 630 IF (ITEM(IX0L)-BLANK) 650, 640, 650 PLOT 316 640 ITEM(IX0L) = AIE PLOT 317 650 IF (IY0L) 700, 700, 660 PLOT 318 660 IF (LINE-IY0L) 700, 670, 700 PLOT 319 670 DO 690 I=1,101 PLOT 320 IF (ITEM(I)-BLANK) 690, 680, 690 PLOT 321 680 ITEM(I) = MINUS PLOT 322 690 CONTINUE PLOT 323 700 CONTINUE PLOT 324 C PLOT 325 IF (IYLAB.NE.LININC) GO TO 710 PLOT 326 WRITE (LPRINT,9985) VALUE, DD, (ITEM(J),J=1,101), DD, PLOT 327 * VALUE PLOT 328 9985 FORMAT (1H , F12.4, 2X, 103A1, F13.4) PLOT 329 VALUE = VALUE - DELTY PLOT 330 IYLAB = 0 PLOT 331 GO TO 720 PLOT 332 710 WRITE (LPRINT,9984) STAR, (ITEM(J),J=1,101), STAR PLOT 333 9984 FORMAT (1H , 14X, 103A1) PLOT 334 720 IYLAB = IYLAB + 1 PLOT 335 DO 730 I=1,101 PLOT 336 ITEM(I) = BLANK PLOT 337 730 CONTINUE PLOT 338 LINE = LINE + 1 PLOT 339 GO TO 390 PLOT 340 C PLOT 341 C LABEL PLOT ALONG THE BOTTOM PLOT 342 C PLOT 343 740 GO TO (830, 830, 830, 750, 760, 770, 780, 790, 800, 810), PLOT 344 * NXUNIT PLOT 345 750 WRITE (LPRINT,9998) PLOT 346 WRITE (LPRINT,9999) (XLABL(I),I=1,NXP1) PLOT 347 GO TO 820 PLOT 348 760 WRITE (LPRINT,9996) PLOT 349 WRITE (LPRINT,9997) (XLABL(I),I=1,NXP1) PLOT 350 GO TO 820 PLOT 351 770 WRITE (LPRINT,9994) PLOT 352 WRITE (LPRINT,9995) (XLABL(I),I=1,NXP1) PLOT 353 GO TO 820 PLOT 354 780 WRITE (LPRINT,9992) PLOT 355 WRITE (LPRINT,9993) (XLABL(I),I=1,NXP1) PLOT 356 GO TO 820 PLOT 357 790 WRITE (LPRINT,9990) PLOT 358 WRITE (LPRINT,9991) (XLABL(I),I=1,NXP1) PLOT 359 GO TO 820 PLOT 360 800 WRITE (LPRINT,9988) PLOT 361 WRITE (LPRINT,9989) (XLABL(I),I=1,NXP1) PLOT 362 GO TO 820 PLOT 363 810 WRITE (LPRINT,9986) PLOT 364 WRITE (LPRINT,9987) (XLABL(I),I=1,NXP1) PLOT 365 820 CONTINUE PLOT 366 RETURN PLOT 367 830 WRITE (LPRINT,9983) PLOT 368 9983 FORMAT (22H0ERROR EXIT FROM PLOT.) PLOT 369 STOP PLOT 370 END PLOT 371 C SERC 1 SUBROUTINE SERCH(VEC, LEN, INDEX, KODE, JY, MX) SERC 2 C SERCH FOR KYST JANUARY,1973 SERC 3 C OCTOBER 10,1971 JAY R LEVINSOHN SERC 4 C MODIFIED FOR KYST BY J.SEERY JANUARY,1973 SERC 5 C SERC 6 C SEARCH SERVES TWO FUNCTIONS SERC 7 C 1 FIND THE MIN AND MAX OF SOME VECTOR VEC SERC 8 C 2 FIND THE NEXT BIGGEST ITEM IN SOME VECTOR VEC SERC 9 C CAN KEEP TRACK OF TWO VECTORS (OF THE SAME LENGTH) AT A TIME SERC 10 C PARAMETERS SERC 11 C VEC-- VECTOR TO BE SEARCHED SERC 12 C LEN-- NUMBER OF ELEMENTS IN VEC SERC 13 C MX-- WORK WITH ELEMENTS MX THROUGH LEN SERC 14 C INDEX-- POINTER TO NEXT BIGGEST ITEM IN VEC SERC 15 C KODE-- OP. CODE 1=FIND MIN AND MAX,0= FIND NEXT BIGGEST ITEM SERC 16 C IF KODE=1 THEN POINTER TO MIN IS RETURNED IN KODE WHILE POINTER SERC 17 C TO MAX IS RETURNED IN INDEX SERC 18 C ROUTINE ASSUMES IT WILL BE CALLED TO FIRST FIND MIN AND MAX OF SERC 19 C VECTOR AND THEN CALLED ITERATIVELY TO FIND THE SECOND LARGEST SERC 20 C ELEMENT, THIRD LARGEST ELEMENT,, ETC.... UNTIL FINAL ELEMENT HAS SERC 21 C HAS BEEN EXTRACTED. ONE CAN TEST THE FINAL ELEMENT AGAINST THE SERC 22 C MIN TO CHECK FOR ACCURATE COMPLETION. SERC 23 C SERC 24 DIMENSION VEC(1) SERC 25 DIMENSION FMIN(2), IMIN(2), INDB4(2), KNT(2) SERC 26 COMMON /DEVICE/ LREAD, LPRINT, LPUNCH SERC 27 COMMON /ACCUR/ XMAXN SERC 28 COMMON /SCOM/ FMIN, IMIN, INDB4, KNT SERC 29 C SERC 30 IF (KODE) 170, 60, 10 SERC 31 C SERC 32 C OP CODE =1 THEN FIND MIN AND MAX SERC 33 C SERC 34 10 FMAX = -XMAXN SERC 35 FMIN(JY) = XMAXN SERC 36 KNT(JY) = MX SERC 37 DO 50 I=MX,LEN SERC 38 IF (VEC(I)-FMAX) 30, 20, 20 SERC 39 20 FMAX = VEC(I) SERC 40 INDEX = I SERC 41 30 IF (VEC(I)-FMIN(JY)) 40, 40, 50 SERC 42 40 FMIN(JY) = VEC(I) SERC 43 KODE = I SERC 44 50 CONTINUE SERC 45 C SERC 46 C SAVE MIN IN SAVE,MIN POINTER IN IMIN,RESET FMIN, SAVE MAX SERC 47 C SERC 48 IMIN(JY) = KODE SERC 49 FMIN(JY) = -XMAXN SERC 50 INDB4(JY) = INDEX SERC 51 RETURN SERC 52 C SERC 53 C OP. CODE NOT EQUAL 1 FIND NEXT LARGEST ELE IN VEC SERC 54 C SERC 55 60 IF (KNT(JY).LT.LEN) GO TO 70 SERC 56 INDEX = -1 SERC 57 GO TO 160 SERC 58 70 INDEX = INDB4(JY) SERC 59 FMB4 = VEC(INDEX) SERC 60 VEC(INDEX) = FMB4 + 1.0 SERC 61 KNT(JY) = KNT(JY) + 1 SERC 62 DO 110 I=MX,LEN SERC 63 IF (VEC(I)-FMIN(JY)) 110, 80, 80 SERC 64 80 IF (VEC(I)-FMB4) 90, 90, 110 SERC 65 90 IF (I-INDB4(JY)) 100, 110, 100 SERC 66 100 FMIN(JY) = VEC(I) SERC 67 INDEX = I SERC 68 110 CONTINUE SERC 69 C SERC 70 C SUBTRACT 1 FROM EACH ELEMENT OF VEC TO GET IT BACK AS BEFORE SERC 71 C DON'T DO SUBTRACTION UNLESS AT LAST ITEM OF VEC SERC 72 C SERC 73 IF (KNT(JY)-LEN) 150, 120, 150 SERC 74 120 DO 140 I=MX,LEN SERC 75 IF (I-INDEX) 130, 140, 130 SERC 76 130 VEC(I) = VEC(I) - 1.0 SERC 77 140 CONTINUE SERC 78 C SERC 79 C SAVE INDEX AND RESET FMIN SERC 80 C SERC 81 150 INDB4(JY) = INDEX SERC 82 FMIN(JY) = -XMAXN SERC 83 160 RETURN SERC 84 C SERC 85 C ERROR ROUTINE ENTER HERE IF KODE<0 SERC 86 C SERC 87 170 WRITE (LPRINT,9999) KODE SERC 88 9999 FORMAT (24H ERROR IN SERCH, KODE= ,I4) SERC 89 RETURN SERC 90 END SERC 91 C XITE 1 SUBROUTINE XITEM(INDEX, XMAX, XMIN, ID, DELX, X, JX, XSF) XITE 2 C XITEM FOR KYST JANUARY,1973 XITE 3 C JAY R LEVINSOHN OCTOBER 31,1971 XITE 4 C ADAPTED FOR KYST BY J.SEERY JANUARY,1973 XITE 5 C XITE 6 C ITS FUNCTION IS TO PLACE ELEMENTS OF THE X VECTOR IN THE XITE 7 C PRINT LINE. XITE 8 C XITE 9 C XITE 10 REAL ITEM(101), PTID(100), X(1) XITE 11 COMMON /PLTCHR/ PTID, ITEM XITE 12 C XITE 13 XNOW = X(INDEX) XITE 14 IF (XNOW-XMAX) 10, 10, 60 XITE 15 10 IF (XNOW-XMIN) 60, 20, 20 XITE 16 20 J = (XNOW-XMIN)/DELX + 1.5 + XSF XITE 17 C XITE 18 C DETERMINE LABEL FOR OBSERVATION XITE 19 C XITE 20 IF (IABS(ID)-3) 30, 50, 30 XITE 21 30 IF (IABS(ID)-2) 40, 50, 40 XITE 22 40 IF (ITEM(J).NE.PTID(2)) ITEM(J) = PTID(JX) XITE 23 GO TO 60 XITE 24 50 ITEM(J) = PTID(INDEX) XITE 25 60 CONTINUE XITE 26 RETURN XITE 27 END XITE 28