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
.