SUBROUTINE DECLAR(SUBR,DUMMY,RSUM) ************************************************************************ * * * Function : Checks for enough memory. * * Version : 1.0 * * Date : August 1996 * * Author : Patrick Groenen * * Calls : SUBR * * Used by : Adapted for Prefmap3 * * * ************************************************************************ C----------------------------------------------------------------------- C Parameters. C INTEGER RSUM EXTERNAL SUBR REAL DUMMY C----------------------------------------------------------------------- C Local variables. C INTEGER RMAXME PARAMETER (RMAXME = 50000) REAL A(RMAXME) C----------------------------------------------------------------------- C Program start. C DO 10 I=1,RMAXME A(I) = 0 10 CONTINUE IF (RSUM .LE. RMAXME) THEN CALL SUBR(A,RSUM) ELSE PRINT *,' ' PRINT *,' Not enough memory available for this task.' END IF RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C P R E F M A P - 3 VERSION 1.0 FEBRUARY 1985 C C C C E X T E R N A L U N F O L D I N G C C C C C C C C INPUT DESCRIPTION AND OPTIONS C C C C C C -CARD 1: TITLE (20A4) C C C C -CARD 2: NROW,NCOL,IOPT,(IDIF(K),K=1,4) (7I5) C C C C NROW = NUMBER OF ROW POINTS; C C NCOL = NUMBER OF COLUMN POINTS C C IOPT = 0: ALL ROWS SAME OPTIONS (DEFAULT) C C = 1: ALL ROWS DIFFERENT OPTIONS C C = 2: INDICATED ROWS SAME OPTIONS, AS SPECIFIED IN C C IDIF = STARTING POINTS FOR DIFFERENT OPTIONS C C C C -CARD 3: NDIM,NANA,IPRE,ISTA,MOD4,MOD3,MOD2,MOD1,MAX,CRIT C C (9I5,F10.8) C C C C NDIM = NUMBER OF DIMENSIONS C C NANA = NUMBER OF SEPARATE ANALYSES C C IPRE = 0: NO PRELIMINARY TRANSFORMATION (DEFAULT) C C = 1: CONFIGURATION WILL BE WEIGHTED C C = 2: CONFIGURATION WILL BE ROTATED AND WEIGHTED C C ISTA = 0: STANDARDIZATION EXTERNAL DATA (DEFAULT) C C = 1: CENTERING ONLY C C = 2: NORMALIZATION ONLY C C = 3: NONE OF THE ABOVE C C MOD4 = 1: VECTOR MODEL WILL BE APPLIED C C MOD3 = 1: UNFOLDING MODEL WILL BE APPLIED C C MOD2 = 1: WEIGHTED UNFOLDING MODEL WILL BE APPLIED C C MOD1 = 1: GENERAL UNFOLDING MODEL WILL BE APPLIED C C MAX = MAXIMUM NR OF NON-METRIC ITERATIONS (DEFAULT=50)C C CRIT = CONVERGENCE CRITERION (DEFAULT=0.00001) C C C C -CARD 4: IPRI0,IPRI1,IPRI2,ISHE,IPLO,IHIS,ISTAT,IWRI1,IWRI2 (9I5) C C C C IPRI0 = 0: NO PRINT OF INPUT C C = 1: PRINT CONFIGURATION C C = 2: PRINT EXTERNAL DATA C C = 3: PRINT BOTH C C IPRI1 = 0: NO RESULTS ACROSS OPTIONS C C = 1: PRINT COMPLETE RESULTS ACROSS OPTIONS C C = 2: PRINT FIT ACROSS OPTIONS C C = 3: PRINT FIT, CRITERION AND PREDICTED VALUES C C IPRI2 = 0: NO RESULTS ACROSS ANALYSES C C = 1: PRINT RESULTS ACROSS ANALYSES C C ISHE = 0: NO SCATTER PLOTS C C = N: FOR N ROWS PLOT OF CRITERION C C AND DISTANCES VERSUS THE DATA C C IPLO = 0: NO PLOT OF FINAL CONFIGURATION C C = 1: PLOT FIRST DIMENSION C C = 2: PLOT FIRST TWO DIMENSIONS C C = K: PLOT (PAIRWISE) K DIMENSIONS C C IHIS = 0: NO PRINT OF HISTORY C C = 1: PRINT HISTORY OF ITERATION C C IST = 0: NO STATISTICS C C = 1: PRINT STATISTICS FOR THE METRIC MODELS C C IWRI1 = 0: NO STORAGE C C = 1: STORE COORDINATES C C = 2: STORE 1 + WEIGHTS C C = 3: STORE 2 + ROTATION MATRICES C C IWRI2 = 0: NO STORAGE C C = 1: PREDICETD VALUES C C = 2: STORE 1 + CRITERION C C C C -CARD 5: INCON,INDATA,INOPT,ISCR1,ISCR2,(LUN(I),I=1,4) (9I5) C C C C INCON = INPUT UNIT CONFIGURATION C C INDATA = INPUT UNIT EXTERNAL DATA C C INOPT = INPUT UNIT OPTIONS C C ISCR1 = UNIT NUMBER SCRATCH FILE 1 C C ISCR2 = UNIT NUMBER SCRATCH FILE 2 C C LUN(1) = OUTPUT UNIT NUMBER COORDINATES C C LUN(2) = OUTPUT UNIT NUMBER WEIGHTS C C LUN(3) = OUTPUT UNIT NUMBER ROTATION MATRICES C C LUN(4) = OUTPUT UNIT NUMBER PREDICTED VALUES C C C C C C -CARD 6: FORMAT OF TARGET CONFIGURATION (20A4) C C C C -FOLLOWING CARDS: TARGET DATA C C C C -FOLLOWING CARD : FORMAT OF EXTERNAL DATA (20A4) C C C C -FOLLOWING CARDS: EXTERNAL DATA C C C C CARDS 1 - 5 AND THE TWO FORMAT CARDS ARE READ FROM INPUT MEDIUM C C INPARA (IN THIS VERSION FIXED TO BE 5). C C C C TO ESTIMATE THE TOTAL NUMBER OF WORDS NEEDED FOR THE ARRAY AREA, C C USE THE FORMULA NWORDS=NCOL*(2*NPRT+NDIM+5)+NPRT+2*NPRM+NCR1*NCR1 C C +NDIM*NANA*(2*NPL1+NPL2+NDIM+2)+NDIM+NW2+NW3+NW4 C C +NPR1**2+NPR2**2+NPR3**2+NPR4**2+NANA C C IN THE STATIC ALLOCATION VERSION, NWORDS IS FIXED AT 25100. C C C C AUTHORS : JACQUELINE MEULMAN & WILLEM HEISER C C BELL LABORATORIES C C 600 MOUNTAIN AVENUE - MURRAY HILL - NEW JERSEY 07974 C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC EXTERNAL MAIN2 COMMON/CDECLA/NDIM,MAX,CRIT,IPRI1,IPRI2,IPLO,IHIS,ISHE,IWRI1,IST, X IWRI2,IPRE,ISTA,NV,NPR1,NPR2,NPR3,NPR4,NPRM,NPL1,NPL2, X NCR1,NW2,NW3,NW4,ND,NANA,IPRI0,MOD1,MOD2,MOD3,MOD4 COMMON/CFORM /FORM1(20) COMMON/CHEAD /TITLE(20) COMMON/CPARAN/NROW,NCOL,IOPT,IDIF(4) COMMON/CPARAT/INCON,INDATA,INOPT,INPARA,IWRITE,ISCR1,ISCR2,LUN(4) C C READER AND PRINTER SPECIFICATION; CHANGE HERE AND IN SUBROUTINE C DECLAR IF OTHER CONVENTIONS ARE USED ON YOUR MACHINE C INPARA=5 IWRITE=6 C C C C READ AND WRITE PARAMETER CARDS C READ(INPARA,500) TITLE READ(INPARA,501) NROW,NCOL,IOPT,(IDIF(K),K=1,4) READ(INPARA,501) NDIM,NANA,IPRE,ISTA,MOD4,MOD3,MOD2,MOD1,MAX,CRIT READ(INPARA,501) IPRI0,IPRI1,IPRI2,ISHE,IPLO,IHIS,IST,IWRI1,IWRI2 READ(INPARA,501) INCON,INDATA,INOPT,ISCR1,ISCR2,(LUN(I),I=1,4) WRITE(IWRITE,600) TITLE WRITE(IWRITE,601) WRITE(IWRITE,602) NROW,NCOL,IOPT,(IDIF(K),K=1,4), 1 NDIM,NANA,IPRE,ISTA,MOD4,MOD3,MOD2,MOD1,MAX,CRIT, 2 IPRI0,IPRI1,IPRI2,ISHE,IPLO,IHIS,IST,IWRI1,IWRI2, 3 INCON,INDATA,INOPT,ISCR1,ISCR2,(LUN(I),I=1,4) C C SET DEFAULT VALUES AND EXTRA PARAMETERS C IF (IOPT.GT.2) IOPT=0 IDIF(1)=1 IF (IOPT.EQ.1.AND.IDIF(2).EQ.0) IDIF(2)=NROW IF (IOPT.EQ.1.AND.IDIF(3).EQ.0) IDIF(3)=NROW IF (IOPT.EQ.1.AND.IDIF(4).EQ.0) IDIF(4)=NROW IF (IOPT.EQ.2.AND.IDIF(2).LE.IDIF(1)) IDIF(2)=NROW+1 IF (IOPT.EQ.2.AND.IDIF(3).LE.IDIF(2)) IDIF(3)=NROW+1 IF (IOPT.EQ.2.AND.IDIF(4).LE.IDIF(3)) IDIF(4)=NROW+1 IF (INDATA.LE.0) INDATA=INPARA IF (INCON.LE.0) INCON=INPARA IF (INOPT.LE.0) INOPT=INPARA IF (LUN(1).LE.0.AND.IWRI1.GT.0) LUN(1)=IWRITE IF (LUN(2).LE.0.AND.IWRI1.GT.1) LUN(2)=IWRITE IF (LUN(3).LE.0.AND.IWRI1.GT.2) LUN(3)=IWRITE IF (LUN(4).LE.0.AND.IWRI2.GT.0) LUN(4)=IWRITE IF (NDIM.LE.0) NDIM=2 IF (NANA.LE.0.OR.NANA.GT.4) NANA=4 IF (MAX.LE.0) MAX=50 IF (CRIT.LT.1.0E-10) CRIT=1.0E-5 IF (IPRE.GT.2) IPRE=0 IF (MOD4.NE.1) MOD4=0 IF (MOD3.NE.1) MOD3=0 IF (MOD2.NE.1) MOD2=0 IF (MOD1.NE.1) MOD1=0 C C GIVE SUMMARY OF OPTIONS C WRITE(IWRITE,603) NROW,NCOL,IOPT IN=1 IR=NROW+1 IF (IDIF(2).NE.IR) IN=IN+1 IF (IDIF(3).NE.IR) IN=IN+1 IF (IDIF(4).NE.IR) IN=IN+1 IF (IOPT.EQ.2) WRITE(IWRITE,1603) (IDIF(K),K=1,IN) WRITE(IWRITE,2603) IF (NROW.EQ.0.OR.NCOL.EQ.0) WRITE(IWRITE,700) IF (NROW.EQ.0.OR.NCOL.EQ.0) STOP WRITE(IWRITE,604) NDIM,NANA,IPRE,ISTA WRITE(IWRITE,1604) WRITE(IWRITE,605) MOD4,MOD3,MOD2,MOD1,MAX,CRIT WRITE(IWRITE,606) IPRI0,IPRI1,IPRI2 WRITE(IWRITE,1606) ISHE,IPLO WRITE(IWRITE,607) IHIS,IST,IWRI1,IWRI2 WRITE(IWRITE,608) INCON,INDATA,INOPT,ISCR1,ISCR2,(LUN(I),I=1,4) C IF (ISCR1.EQ.0.OR.ISCR2.EQ.0) WRITE(IWRITE,700) IF (ISCR1.EQ.0.OR.ISCR2.EQ.0) STOP C C COMPUTE THE NUMBER OF PREDICTOR VARIABLES FOR THE VARIOUS C MODELS. IF A MODEL WILL NOT BE APPLIED, THE NUMBER OF COLUMNS C OF ITS PREDICTOR MATRIX WILL BE DIMINISHED TO 1, UNLESS A PRE- C LIMINARY TRANSFORMATION IS ASKED FOR (NPR1 AND NPR4) C WHEN THE NUMBER OF COLUMNS IS NOT LARGE ENOUGH TO ESTIMATE C THE PARAMETERS OF A MODEL, THE MODEL WILL NOT BE APPLIED. C WHEN THE NUMBER OF COLUMNS IS NOT LARGE ENOUGH THE PRELIMI- C NARY TRANSFORMATION WILL NOT BE PERFORMED. C WHEN THE NUMBER OF DIMENSIONS EQUALS 1, THE GENERAL MODEL AND C THE PRELIMINARY ROTATION WILL NOT BE APPLIED C NPRE=IPRE KPT1=NDIM+2 KPT2=(NDIM*(NDIM+1))/2+1 IF (KPT1.GT.NCOL) NPRE=0 IF (KPT2.GT.NCOL.AND.IPRE.EQ.2) NPRE=0 IF (NDIM.EQ.1) MOD1=0 IF (NDIM.EQ.1.AND.IPRE.EQ.2) NPRE=0 IF (NPRE.NE.IPRE) WRITE (IWRITE,701) IPRE=NPRE C NPR1=(NDIM*(NDIM+3))/2+1 IF (NPR1.GT.NCOL) MOD1=0 KPR1=MAX0(MOD1,IPRE) IF (KPR1.GT.0) KPR1=1 NPR1=KPR1*(NPR1-1)+1 NPR2=2*NDIM+1 IF (NPR2.GT.NCOL) MOD2=0 NPR2=MOD2*(NPR2-1)+1 NPR3=NDIM+2 IF (NPR3.GT.NCOL) MOD3=0 NPR3=MOD3*(NPR3-1)+1 NPR4=NDIM+1 IF (NPR4.GT.NCOL) MOD4=0 KPR4=MAX0(MOD4,IPRE) IF (KPR4.GT.0) KPR4=1 NPR4=KPR4*(NPR4-1)+1 NPRT=NPR1+NPR2+NPR3+NPR4 NPRM=MAX0(NPR1,NPR2,NPR3,NPR4) N=MAX0(IPLO,IPRI2) IF (N.NE.0) N=1 NW2=MAX0((NROW+NCOL)*N,NCOL+NCOL,NDIM*NDIM) NW3=NW2 NW4=NW2 ND =NDIM*NANA C C WHEN PLOTTING OF IDEAL-POINTS AND/OR VECTORS IS NOT DESIRED, C THE DIMENSIONS OF TWEIGH,TCOOR AND TROTAT WILL BE DIMINISHED C N=MAX0(IPLO,IPRI2) IF (N.NE.0) N=1 NPL1=(NROW-1)*N+1 NPL2=(NROW-1)*NDIM*MOD1*N+NDIM C C WHEN PREPROCESSING IS NOT DESIRED, C THE DIMENSION OF ARRAY CROSS WILL BE DIMINISHED C N=0 IF (IPRE.NE.0) N=1 NCR1=(NCOL-1)*N+1 C C ALLOCATE STORAGE AND START COMPUTATIONS C C C C CON = NCOL*NDIM C PR1 = NCOL*NPR1 C PR2 = NCOL*NPR2 C PR3 = NCOL*NPR3 C PR4 = NCOL*NPR4 C SIVA1 = NPR1 C SIVA2 = NPR2 C SIVA3 = NPR3 C SIVA4 = NPR4 C LSV1 = NCOL*NPR1 C LSV2 = NCOL*NPR2 C LSV3 = NCOL*NPR3 C LSV4 = NCOL*NPR4 C RSV1 = NPR1*NPR1 C RSV2 = NPR2*NPR2 C RSV3 = NPR3*NPR3 C RSV4 = NPR4*NPR4 C TCOOR = NPL1*NDIM*NANA C TWEIGH = NPL1*NDIM*NANA C TROTAT = NPL2*NDIM*NANA C COEF = NPRM C COOR = NDIM*NANA C WEIGH = NDIM*NANA C ROTAT = NDIM*NANA*NDIM C CROSS = NCR1*NCR1 C WSIGN = NDIM C DATA = NCOL C DISP = NCOL C DSTAR = NCOL C IORD = NCOL C LBK = NCOL C WORK1 = NPRM C WORK2 = NW2 C WORK3 = NW3 C WORK4 = NW4 C OPTION = NANA C DANA = NCOL C C NWORDS=NCOL*(2*NPRT+NDIM+6)+NPRT+2*NPRM+NCR1*NCR1 X +NDIM*NANA*(2*NPL1+NPL2+NDIM+2)+NDIM+NW2+NW3+NW4 X +NPR1**2+NPR2**2+NPR3**2+NPR4**2+NANA CALL DECLAR(MAIN2,A,NWORDS) STOP C C FORMAT SPECIFICATIONS C 500 FORMAT(20A4) 501 FORMAT(9I5,F10.8) 600 FORMAT(1H1//,5X,17HP R E F M A P - 3/ 1 1H ,47X, 31H *** EXTERNAL UNFOLDING *** ,25X, 1 18HJACQUELINE MEULMAN / 2 5X,15H VERSION - 1.0 ,28X,31H *** & *** ,25X, 3 13HWILLEM HEISER/ 3 5X,14H FEBRUARY 1985, 29X,31H *** PROPERTY FITTING *** ,25X, 3 18HJ. DOUGLAS CARROLL/ 4 1H ,103X,17HBELL LABORATORIES/ 5 1H ,103X,15HMURRAY HILL, NJ//// 6 1H ,4X,13H1 JOB TITLE: ,20A4//) 601 FORMAT(1H ,4X,24HECHO OF PARAMETER CARDS://) 602 FORMAT(1H ,4X,2H2 ,7I5/ 1 1H ,4X,2H3 ,9I5,F10.8/ 2 1H ,4X,2H4 ,9I5/ 3 1H ,4X,2H5 ,9I5///) 603 FORMAT(1H ,4X,22H2 DATA SPECIFICATIONS:// 1 1H ,8X,44HNUMBER OF ROW POINTS (THESE ARE FITTED) IS ,I5/ 2 1H ,8X,44HNUMBER OF COLUMN POINTS (THESE ARE GIVEN) IS,I5/ 3 1H ,8X,44HOPTION SET SELECTION: ,I5/ 4 1H ,8X,44H 0 = ALL ROWS SAME OPTION SET / 5 1H ,8X,44H 1 = ALL ROWS DIFFERENT OPTION SETS / 6 1H ,8X,44H 2 = SPECIFIED ROWS SAME OPTION SET ) 1603 FORMAT(1H ,8X,27HOPTION SETS START WITH ROWS,17X,4I5) 2603 FORMAT(1H ,//) 604 FORMAT(1H ,4X,26H3 ANALYSIS SPECIFICATIONS:// 1 1H ,8X,44HTHE NUMBER OF DIMENSIONS IS ,I5/ 2 1H ,8X,44HMAXIMUM NUMBER OF ANALYSES IN ANY OPTION SET,I5/ 3 1H ,8X,44HPRELIMINARY TRANSFORMATION OF CONFIGURATION ,I5/ 4 1H ,8X,44H 0 = REMAINS UNCHANGED / 5 1H ,8X,44H 1 = WEIGHTED / 6 1H ,8X,44H 2 = ROTATED AND WEIGHTED / 7 1H ,8X,44HSTANDARDIZE EXTERNAL DATA ,I5/ 8 1H ,8X,44H 0 = YES (=1+2) / 9 1H ,8X,44H 1 = CENTER ONLY ) 1604 FORMAT(1H ,8X,44H 2 = NORMALIZE ONLY / 2 1H ,8X,44H 3 = NONE OF THE ABOVE ) 605 FORMAT(1H ,8X,44HMODELS: 0 = NOT APPLIED, 1 = APPLIED / 1 1H ,8X,44H VECTOR MODEL ,I5/ 2 1H ,8X,44H UNFOLDING MODEL ,I5/ 3 1H ,8X,44H WEIGHTED UNFOLDING MODEL ,I5/ 4 1H ,8X,44H GENERAL UNFOLDING MODEL ,I5/ 5 1H ,8X,44HTHE NUMBER OF NON-METRIC ITERATIONS IS ,I5/ 6 1H ,8X,39HTHE CONVERGENCE CRITERION IS ,E10.2///) 606 FORMAT(1H ,4X,44H4 PRINT/PLOT OPTIONS: 0 = NONE // 1 1H ,8X,44HPRINT INPUT: ,I5/ 2 1H ,8X,44H 1 = TARGET CONFIGURATION / 3 1H ,8X,44H 2 = EXTERNAL DATA / 4 1H ,8X,44H 3 = 1 & 2 / 5 1H ,8X,44HPRINT RESULTS ACROSS OPTION SETS ,I5/ 6 1H ,8X,44H 1 = COMPLETE RESULTS / 7 1H ,8X,44H 2 = FIT / 8 1H ,8X,44H 3 = 2 & CRITERION-PREDICTED VALUES / 9 1H ,8X,44HSELECTED RESULTS ACROSS ANALYSES ,I5) 1606 FORMAT(1H ,8X,44HSCATTER&TRANSFORMATION PLOT FOR N ROWS, N =,I5/ 2 1H ,8X,44HPAIRWISE PLOTS IDEAL POINTS/VECTORS, DIM =,I5) 607 FORMAT(1H ,8X,44HHISTORY OF NON-METRIC REGRESSION ,I5/ 1 1H ,8X,44HSTATISTICS FOR METRIC ANALYSES ,I5/ 2 1H ,8X,44HSTORE OUTPUT: ,I5/ 3 1H ,8X,44H 1 = COORDINATES / 4 1H ,8X,44H 2 = LIKE 1 + WEIGHTS / 5 1H ,8X,44H 3 = LIKE 2 + ROTATION MATRICES / 6 1H ,8X,44HSTORE OUTPUT: ,I5/ 8 1H ,8X,44H 1 = PREDICTED VALUES / 9 1H ,8X,44H 2 = LIKE 1 + CRITERION ///) 608 FORMAT(1H ,4X,31H5 UNIT NUMBERS FOR INPUT/OUTPUT// 1 1H ,8X,44HTHE CONFIGURATION WILL BE READ FROM UNIT ,I5/ 2 1H ,8X,44HTHE EXTERNAL DATA WILL BE READ FROM UNIT ,I5/ 3 1H ,8X,44HTHE OPTIONS WILL BE READ FROM UNIT ,I5/ 4 1H ,8X,44HUNIT NUMBER FOR SCRATCH FILE 1 ,I5/ 5 1H ,8X,44HUNIT NUMBER FOR SCRATCH FILE 2 ,I5/ 6 1H ,8X,44HOUTPUT UNIT FOR THE COORDINATES ,I5/ 7 1H ,8X,44HOUTPUT UNIT FOR THE WEIGHTS ,I5/ 8 1H ,8X,44HOUTPUT UNIT FOR THE ROTATION MATRICES ,I5/ 9 1H ,8X,44HOUTPUT UNIT FOR PREDICTED&CRITERION VALUES ,I5/) 700 FORMAT(1H0///,9X,43HTHE PROGRAM STOPS BECAUSE OF A FATAL ERROR , X 43HIN THE PARAMETER CARD JUST PRINTED. / X 1H ,8X, 41HCHECK NUMBER OF ROWS OR COLUMNS, OR UNIT , X 43HNUMBERS FOR 2 OBLIGATORY SCRATCH FILES. ) 701 FORMAT(1H0///,5X,43H*** WARNING *** THE PRELIMINARY TRANSFORMAT, X 43HION WILL NOT BE PERFORMED, BECAUSE THE / X 1H0,20X, 42HNUMBER OF TARGET POINTS IS NOT SUFFICIENT.) C END SUBROUTINE ESTIMA(DATA,DSTAR,PRM,PSI,NCOL,NPR,COEF,RR,KK) C ****************************************************************** C * * C * E S T I M A * C * * C * PURPOSE: COMPUTES THE REGRESSION COEFFICIENTS, THE PREDICTED * C * VALUES AND THE MULTIPLE CORRELATION * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION DATA(NCOL),DSTAR(NCOL),PRM(NCOL,NPR),PSI(NCOL,NPR), X COEF(NPR) C C FOR THE VECTOR MODEL THE DISSIMILARITIES WILL BE REVERSED C DSYM=1. IF (KK.EQ.4) DSYM=-1.0 C DO 10 J=1,NPR SUM=0.0 DO 20 I=1,NCOL SUM=SUM+PSI(I,J)*DATA(I)*DSYM 20 CONTINUE COEF(J)=SUM 10 CONTINUE C RMEAN=0.0 DO 30 I=1,NCOL SUM=0.0 DO 40 K=1,NPR SUM=SUM+PRM(I,K)*COEF(K)*DSYM 40 CONTINUE DSTAR(I)=SUM RMEAN=RMEAN+SUM 30 CONTINUE RMEAN=RMEAN/FLOAT(NCOL) STANDP=0.0 DO 60 I=1,NCOL A=DSTAR(I)-RMEAN STANDP=STANDP+A*A 60 CONTINUE STANDP=SQRT(STANDP) C RMEAN=0.0 STANDD=0.0 DO 70 I=1,NCOL RMEAN=RMEAN+DATA(I) 70 CONTINUE RMEAN=RMEAN/FLOAT(NCOL) DO 80 I=1,NCOL A=DATA(I)-RMEAN STANDD=STANDD+A*A 80 CONTINUE STANDD=SQRT(STANDD) C RR=STANDP/STANDD RETURN END SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR) C ****************************************************************** C * * C * I M T Q L 2 * C * * C * PURPOSE: DETERMINES THE EIGENVALUES AND EIGENVECTORS OF A * C * SYMMETRIC TRIDIAGONAL MATRIX USING THE IMPLICIT * C * QL METHOD. * C * * C * SUBROUTINES CALLED: NONE * C * * C * ADAPTED FROM EISPACK * C * * C ****************************************************************** DIMENSION D(N),E(N),Z(NM,N) REAL MACHEP C C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPCIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C C ********** C MACHEP = 2.0**(-20) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C E(N) = 0.0 C DO 240 L = 1, N J = 0 C ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT ********** 105 DO 110 M = L, N IF (M .EQ. N) GO TO 120 IF (ABS(E(M)) .LE. MACHEP* (ABS(D(M)) + ABS(D(M+1)))) X GO TO 120 110 CONTINUE C 120 P = D(L) IF (M .EQ. L) GO TO 240 IF (J .EQ.30) GO TO 1000 J = J + 1 C ********** FORM SHIFT ********** G = (D(L+1) - P) / (2.0 * E(L)) R = SQRT(G*G+1.0) G = D(M) - P + E(L) / (G + SIGN(R,G)) S = 1.0 C = 1.0 P = 0.0 NML = M - L C ********** FOR I=M-1 STEP -1 UNTIL L DO -- ********** DO 200 II = 1, NML I = M - II F = S * E(I) B = C * E(I) IF (ABS(F) .LT. ABS(G)) GO TO 150 C = G / F R = SQRT(C*C+1.0) E(I+1) = F * R S = 1.0 / R C = C * S GO TO 160 150 S = F / G R = SQRT(S*S+1.0) E(I+1) = G* R C = 1.0 / R S = S * C 160 G = D(I+1) - P R = (D(I) - G) * S + 2.0 * C * B P = S * R D(I+1) = G + P G = C * R - B C ********** FORM VECTOR ********** DO 180 K = 1, N F = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * F Z(K,I) = C * Z(K,I) - S * F 180 CONTINUE C 200 CONTINUE C D(L) = D(L) - P E(L) = G E(M) = 0.0 GO TO 105 240 CONTINUE C ********** ORDER EIGENVALUES AND EIGENVECTORS ********** DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P ) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C ********** SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS ********** 1000 IERR = L 1001 RETURN END SUBROUTINE MAIN2(A,NWORDS) C ****************************************************************** C * * C * M A I N 2 * C * * C * PURPOSE: COMPUTES THE RELATIVE ADRESSES FOR ALL ARRAYS IN * C * THE REST OF THE PROGRAM PREFMAP-III * C * * C * SUBROUTINES CALLED: MAIN3 * C * * C ****************************************************************** DIMENSION A(NWORDS) COMMON/CDECLA/NDIM,MAX,CRIT,IPRI1,IPRI2,IPLO,IHIS,ISHE,IWRI1,IST, X IWRI2,IPRE,ISTA,NV,NPR1,NPR2,NPR3,NPR4,NPRM,NPL1,NPL2, X NCR1,NW2,NW3,NW4,ND,NANA,IPRI0,MOD1,MOD2,MOD3,MOD4 COMMON/CPARAN/NROW,NCOL,IOPT,IDIF(4) C ICON = 1 IPR1 = 1 + NCOL*NDIM IPR2 = IPR1 + NCOL*NPR1 IPR3 = IPR2 + NCOL*NPR2 IPR4 = IPR3 + NCOL*NPR3 ISV1 = IPR4 + NCOL*NPR4 ISV2 = ISV1 + NPR1 ISV3 = ISV2 + NPR2 ISV4 = ISV3 + NPR3 ILSV1 = ISV4 + NPR4 ILSV2 = ILSV1 + NCOL*NPR1 ILSV3 = ILSV2 + NCOL*NPR2 ILSV4 = ILSV3 + NCOL*NPR3 IRSV1 = ILSV4 + NCOL*NPR4 IRSV2 = IRSV1 + NPR1*NPR1 IRSV3 = IRSV2 + NPR2*NPR2 IRSV4 = IRSV3 + NPR3*NPR3 ITCOOR = IRSV4+ NPR4*NPR4 ITWEIG = ITCOOR+ NPL1*NDIM*NANA ITROTA = ITWEIG+ NPL1*NDIM*NANA ICOEF = ITROTA+ NPL2*NDIM*NANA ICOOR = ICOEF + NPRM IWEIGH = ICOOR + NDIM*NANA IROTAT = IWEIGH+ NDIM*NANA ICROSS = IROTAT+ NDIM*NANA*NDIM IWSIGN = ICROSS+ NCR1*NCR1 IDATA = IWSIGN+ NDIM IDISP = IDATA + NCOL IDSTAR = IDISP + NCOL IIORD = IDSTAR+ NCOL ILBK = IIORD + NCOL IW1 = ILBK + NCOL IW2 = IW1 + NPRM IW3 = IW2 + NW2 IW4 = IW3 + NW3 IOO = IW4 + NW4 IDANA = IOO + NANA C CALL MAIN3(A(ICON),A(IPR1),A(IPR2),A(IPR3),A(IPR4), X A(ISV1),A(ISV2),A(ISV3),A(ISV4), X A(ILSV1),A(ILSV2),A(ILSV3),A(ILSV4), X A(IRSV1),A(IRSV2),A(IRSV3),A(IRSV4),A(ITCOOR), X A(ITWEIG),A(ITROTA),A(ICOEF),A(ICOOR),A(IWEIGH), X A(IROTAT),A(ICROSS),A(IWSIGN),A(IDATA),A(IDISP), X A(IDSTAR),A(IIORD),A(ILBK),A(IW1),A(IW2),A(IW3),A(IW4), X A(IOO),A(IDANA),NDIM,MAX,CRIT, X IPRI0,IPRI1,IPRI2,IPLO,IHIS,IST,ISHE,IWRI1,IWRI2,IPRE,ISTA, X NROW,NCOL,NPR1,NPR2,NPR3,NPR4,NPRM,NPL1,NPL2,NCR1, X NW2,NW3,NW4,ND,NANA,IOPT,IDIF,MOD1,MOD2,MOD3,MOD4) C RETURN END SUBROUTINE MAIN3(CON,PR1,PR2,PR3,PR4,SIVA1,SIVA2,SIVA3,SIVA4, X LSV1,LSV2,LSV3,LSV4,RSV1,RSV2,RSV3,RSV4, X TCOOR,TWEIGH,TROTAT,COEF,COOR,WEIGH,ROTAT, X CROSS,WSIGN,DATA,DISP,DSTAR,IORD,LBK,WORK1, X WORK2,WORK3,WORK4,OPTION,DANA, X NDIM,MAX,CRIT,IPRI0,IPRI1,IPRI2, X IPLO,IHIS,IST,ISHE,IWRI1,IWRI2,IPRE,ISTA, X NROW,NCOL,NPR1,NPR2,NPR3,NPR4,NPRM,NPL1,NPL2, X NCR1,NW2,NW3,NW4,ND,NANA,IOPT,IDIF,MOD1,MOD2, X MOD3,MOD4) C ****************************************************************** C * * C * M A I N 3 * C * * C * PURPOSE: CONTROLS THE FLOW OF THE PREFMAP-III PROGRAM * C * * C * SUBROUTINES CALLED: PREREC,ZERO,PREPRO,PRED1(-4),PSINV,ZEROI, * C * ESTIMA,MONO,UNRAV,RESUL2,STATIS,OUTP3 * C * * C ****************************************************************** COMMON/CFORM /FORM1(20) COMMON/CHEAD /TITLE(20) COMMON/CPARAT/INCON,INDATA,INOPT,INPARA,IWRITE,ISCR1,ISCR2,LUN(4) C DIMENSION CON(NCOL,NDIM),PR1(NCOL,NPR1),PR2(NCOL,NPR2), X PR3(NCOL,NPR3),PR4(NCOL,NPR4) DIMENSION SIVA1(NPR1),SIVA2(NPR2),SIVA3(NPR3),SIVA4(NPR4), X LSV1(NCOL,NPR1),LSV2(NCOL,NPR2),LSV3(NCOL,NPR3), X LSV4(NCOL,NPR4),RSV1(NPR1,NPR1),RSV2(NPR2,NPR2), X RSV3(NPR3,NPR3),RSV4(NPR4,NPR4),TCOOR(NPL1,ND), X TWEIGH(NPL1,ND),TROTAT(NPL2,ND),COEF(NPRM),DANA(NCOL), X COOR(ND),WEIGH(ND),ROTAT(NDIM,ND),CROSS(NCR1,NCR1), X WSIGN(NDIM),OPTION(NANA),OPTLAB(12),OPTLBS(12),IDIF(4), X DATA(NCOL),DISP(NCOL),DSTAR(NCOL),IORD(NCOL),LBK(NCOL), X WORK1(NPRM),WORK2(NW2),WORK3(NW3),WORK4(NW4),FITT(4), X HULP(4),HULP2(4),FITM(12),KCOUNT(12),VAF(12),PVAF(12), X SFITM(12) REAL LSV1,LSV2,LSV3,LSV4 INTEGER OPTION,OPTLAB,OPTLBS,HULP,HULP2 DATA OPTLAB/'GM ','WM ','UM ','VM ','GP ','WP ','UP ', X 'VP ','GS ','WS ','US ','VS '/ DATA OPTLBS/'gm ','wm ','um ','vm ','gp ','wp ','up ', X 'vp ','gs ','ws ','us ','vs '/ C C READ FORMAT AND TARGET CONFIGURATION C READ(INPARA,500) FORM1 WRITE(IWRITE,635) FORM1 READ(INCON,FORM1) ((CON(J,M),M=1,NDIM),J=1,NCOL) C C CENTER TARGET CONFIGURATION C DO 61 M=1,NDIM RMEAN=0.0 DO 62 J=1,NCOL RMEAN=RMEAN+CON(J,M) 62 CONTINUE RMEAN=RMEAN/FLOAT(NCOL) DO 63 J=1,NCOL CON(J,M)=CON(J,M)-RMEAN 63 CONTINUE WRITE (IWRITE,629) M,RMEAN 61 CONTINUE C IF (IPRI0.EQ.1.OR.IPRI0.EQ.3) WRITE (IWRITE,600) IF (IPRI0.EQ.1.OR.IPRI0.EQ.3) X CALL PREREC(CON,NCOL,NDIM,NCOL,IWRITE,0) C C READ THE EXTERNAL DATA AND STORE ON SCRATCH FILE. IF PREPROCESSING C OF THE TARGET CONFIGURATION IS DESIRED, ACCUMULATE THE MATRIX OF C SUMS OF SQUARES AND CROSS PRODUCTS. C REWIND ISCR1 REWIND ISCR2 C C DUMMY WRITING TO SCRATCH FILES (NECESSARY FOR VAX-UNIX) C WRITE(ISCR1) (K,K=1,NCOL) WRITE(ISCR2) (K,K=1,NANA) REWIND ISCR1 REWIND ISCR2 C CALL ZERO(CROSS,NCR1*NCR1) CALL ZERO(TWEIGH,NPL1*ND) CALL ZERO(TCOOR,NPL1*ND) CALL ZERO(TROTAT,NPL2*ND) READ (INPARA,500) FORM1 WRITE(IWRITE,636) FORM1 IF (IPRI0.EQ.0.OR.IPRI0.EQ.1) KN=MIN0(NROW,10) IF (IPRI0.EQ.2.OR.IPRI0.EQ.3) KN=NROW KNN=MIN0(NCOL,10) WRITE(IWRITE,604) KN,(K,K=1,KNN) NN=NROW DO 5 I=1,NROW READ (INDATA,FORM1) (DATA(J),J=1,NCOL) IF (IPRI0.EQ.2.OR.IPRI0.EQ.3.OR.I.LE.10) X WRITE (IWRITE,605) I,(DATA(J),J=1,NCOL) C C NORMALIZE DATA AND/OR COMPUTE SSQ C C THE PRELIMINARY TRANSFORMATION AND THE METRIC ANALYSES ARE PERFORMED C ACCORDING TO THE DESIRED NORMALIZATION OF THE EXTERNAL DATA. C THE NONMETRIC ANALYSES ARE PERFORMED ON THE BASIS OF THE EXTERNAL C DATA, NORMALIZED TO HAVE ZERO MEAN AND VARIANCE EQUAL TO ONE. C THEREFORE THE RAW DATA WILL BE WRITTEN TO THE SCRATCH FILE (ARRAY C DATA) WHILE ARRAY DANA WILL FIRST CONTAIN THE DATA, NORMALIZED C ACCORDING TO OPTION ISTA. FOR NONMETRIC ANALYSES THE RAW DATA WILL C BE NORMALIZED AS DESCRIBED ABOVE, NO MATTER WHAT VALUE THE PARAMETER C ISTA TAKES. ARRAY DISP WILL THEN BE USED TO OBTAIN THE OPTIMALLY C TRANSFORMED NORMALIZED DATA. IN THIS WAY THE NORMALIZATION OF THE ORI- C GINAL DATA IS CONTROLLED (E.G. PRESERVING THE ORIGINAL VARIANCE), WHILE C FOR NONMETRIC ANALYSES WE OBTAIN AN APPROPRIATE SCALING TO MAKE TRANS- C FORMATION PLOTS. C RMEAN=0.0 STAND=0.0 VAR=0.0 DO 19 J=1,NCOL RMEAN=RMEAN+DATA(J) 19 CONTINUE RMEAN=RMEAN/FLOAT(NCOL) DO 29 J=1,NCOL A=DATA(J)-RMEAN IF (ISTA.EQ.0.OR.ISTA.EQ.1) DANA(J)=A IF (ISTA.EQ.2.OR.ISTA.EQ.3) DANA(J)=DATA(J) STAND=STAND+DANA(J)*DANA(J) VAR=VAR+A*A 29 CONTINUE IF (VAR.EQ.0.0) WRITE(IWRITE,900) I IF (VAR.EQ.0.0) GO TO 6 IF (ISTA.EQ.1) GO TO 18 IF (ISTA.EQ.3) GO TO 18 STAND=SQRT(STAND/FLOAT(NCOL)) DO 39 J=1,NCOL DANA(J)=DANA(J)/STAND 39 CONTINUE 18 CONTINUE C IF (IPRE.LE.0) GO TO 6 DO 7 J=1,NCOL DO 8 K=1,NCOL CROSS(J,K)=CROSS(J,K)+DANA(J)*DANA(K) 8 CONTINUE 7 CONTINUE 6 CONTINUE WRITE (ISCR1) (DATA(J),J=1,NCOL) 5 CONTINUE C NROW=NN REWIND ISCR1 C DO 9 I=1,NDIM WSIGN(I)=1.0 9 CONTINUE C IF (IPRE.GT.0) CALL PREPRO(CON,CROSS,PR1,LSV1,PR4,LSV4,RSV1,RSV4, X WORK1,WSIGN,IPRE,IWRITE,NCOL,NDIM,NPR1, X NPR4,NPRM) IF (IPRE.GT.0) WRITE (IWRITE,601) IF (IPRE.GT.0) CALL PREREC(CON,NCOL,NDIM,NCOL, X IWRITE,0) C C COMPUTE MAXIMAL DISTANCE FROM ORIGIN C CMAX=0.0 DO 52 I=1,NCOL A=0.0 DO 53 J=1,NDIM A=A+CON(I,J)*CON(I,J) 53 CONTINUE IF (A.GT.CMAX) CMAX=A 52 CONTINUE CMAX=SQRT(CMAX) C C CHECK MODELS THAT WILL BE USED C AND ACCORDINGLY CONSTRUCT PREDICTOR MATRICES C IF (MOD1.EQ.1) CALL PRED1(CON,PR1,LSV1,NCOL,NDIM,NPR1) IF (MOD2.EQ.1) CALL PRED2(CON,PR2,LSV2,NCOL,NDIM,NPR2) IF (MOD3.EQ.1) CALL PRED3(CON,PR3,LSV3,WSIGN,NCOL,NDIM,NPR3) IF (MOD4.EQ.1) CALL PRED4(CON,PR4,LSV4,NCOL,NDIM,NPR4) C C C COMPUTE THE PSEUDO INVERSES OF THE PREDICTOR MATRICES C ON BASIS OF THEIR SINGULAR VALUE DECOMPOSITIONS C IF (MOD1.EQ.1) CALL PSINV(LSV1,RSV1,SIVA1,NPR1,NCOL,WORK1) IF (MOD2.EQ.1) CALL PSINV(LSV2,RSV2,SIVA2,NPR2,NCOL,WORK1) IF (MOD3.EQ.1) CALL PSINV(LSV3,RSV3,SIVA3,NPR3,NCOL,WORK1) IF (MOD4.EQ.1) CALL PSINV(LSV4,RSV4,SIVA4,NPR4,NCOL,WORK1) C CALL ZEROI(HULP2,4) CALL ZEROI(KCOUNT,12) CALL ZERO(FITM,12) CALL ZERO(SFITM,12) CALL ZERO(VAF,12) CALL ZERO(PVAF,12) C IO=1 II=1 WRITE (IWRITE,632) WRITE (IWRITE,631) (J,J=1,NANA) DO 1 I=1,NROW IF (IOPT.LE.0.AND.I.GT.1) GO TO 17 IF (IOPT.EQ.2.AND.I.NE.IDIF(IO)) GO TO 17 IF (IOPT.GT.1.AND.IO.LT.4) IO=IO+1 READ (INOPT,501) (HULP(J),J=1,NANA) WRITE (IWRITE,633) II,(HULP(J),J=1,NANA) II=II+1 DO 2 K=1,NANA DO 3 J=1,12 IF (HULP(K).EQ.OPTLAB(J)) GO TO 4 IF (HULP(K).EQ.OPTLBS(J)) GO TO 4 3 CONTINUE HULP(K)=0 GO TO 2 4 HULP(K)=J IF (J.GT.8) HULP(K)=HULP(K)+12 IF (J.GT.4.AND.J.LT.9) HULP(K)=HULP(K)+6 IF (MOD(HULP(K),10).EQ.1.AND.MOD1.EQ.0) X WRITE(IWRITE,630) K IF (MOD(HULP(K),10).EQ.1.AND.MOD1.EQ.0) HULP(K)=0 IF (MOD(HULP(K),10).EQ.2.AND.MOD2.EQ.0) X WRITE(IWRITE,630) K IF (MOD(HULP(K),10).EQ.2.AND.MOD2.EQ.0) HULP(K)=0 IF (MOD(HULP(K),10).EQ.3.AND.MOD3.EQ.0) X WRITE(IWRITE,630) K IF (MOD(HULP(K),10).EQ.3.AND.MOD3.EQ.0) HULP(K)=0 IF (MOD(HULP(K),10).EQ.4.AND.MOD4.EQ.0) X WRITE(IWRITE,630) K IF (MOD(HULP(K),10).EQ.4.AND.MOD4.EQ.0) HULP(K)=0 2 CONTINUE 17 CONTINUE WRITE(ISCR2) (HULP(K),K=1,NANA) 1 CONTINUE REWIND ISCR2 C C MAIN LOOP ACROSS ROWS C C FOR EACH ROW SEPARATELY : C REREAD THE EXTERNAL DATA AND THE MODEL DESIRED C DO 10 I=1,NROW READ(ISCR1) (DATA(J),J=1,NCOL) READ(ISCR2) (OPTION(J),J=1,NANA) C CALL ZERO(COOR,NANA*NDIM) CALL ZERO(WEIGH,NANA*NDIM) CALL ZERO(ROTAT,NANA*NDIM*NDIM) CALL ZERO(FITT,NANA) C C C INNER LOOP ACROSS OPTIONS C DO 20 J=1,NANA C C NORMALIZE EXTERNAL DATA C RMEAN=0.0 STAND=0.0 VAR=0.0 DO 119 M=1,NCOL RMEAN=RMEAN+DATA(M) 119 CONTINUE RMEAN=RMEAN/FLOAT(NCOL) DO 129 M=1,NCOL A=DATA(M)-RMEAN IF (ISTA.EQ.0.OR.ISTA.EQ.1) DANA(M)=A IF (ISTA.EQ.2.OR.ISTA.EQ.3) DANA(M)=DATA(M) STAND=STAND+DANA(M)*DANA(M) VAR=VAR+A*A 129 CONTINUE IF (VAR.EQ.0.0) GO TO 10 STAND=SQRT(STAND/FLOAT(NCOL)) VAR=VAR/FLOAT(NCOL) IF (ISTA.EQ.1) GO TO 118 IF (ISTA.EQ.3) GO TO 118 RMEAN=0.0 VAR=0.0 DO 139 M=1,NCOL DANA(M)=DANA(M)/STAND RMEAN=RMEAN+DANA(M) 139 CONTINUE RMEAN=RMEAN/FLOAT(NCOL) DO 140 M=1,NCOL A=DANA(M)-RMEAN VAR=VAR+A*A 140 CONTINUE VAR=VAR/FLOAT(NCOL) 118 CONTINUE C C COMPUTE REGRESSION COEFFICIENTS, PREDICTED VALUES AND MULTIPLE C CORRELATION FOR THE METRIC CASE C KK=MOD(OPTION(J),10) IF (KK.EQ.0) GO TO 20 IF (KK.EQ.1.OR.KK.EQ.2) HULP2(J)=HULP2(J)+1 IF (IPRI1.NE.1) GO TO 38 IF (KK.EQ.1) WRITE (IWRITE,611) I,J IF (KK.EQ.2) WRITE (IWRITE,612) I,J IF (KK.EQ.3) WRITE (IWRITE,613) I,J IF (KK.EQ.4) WRITE (IWRITE,614) I,J IF (OPTION(J).LT.11.AND.OPTION(J).GT.0) WRITE (IWRITE,608) IF (OPTION(J).GT.10) WRITE (IWRITE,609) IF (OPTION(J).GT.10.AND.OPTION(J).LT.20) WRITE (IWRITE,620) IF (OPTION(J).GT.20) WRITE (IWRITE,621) C 38 IF(KK.EQ.1) CALL ESTIMA(DANA,DSTAR,PR1,LSV1,NCOL,NPR1,COEF,FIT,KK) IF(KK.EQ.2) CALL ESTIMA(DANA,DSTAR,PR2,LSV2,NCOL,NPR2,COEF,FIT,KK) IF(KK.EQ.3) CALL ESTIMA(DANA,DSTAR,PR3,LSV3,NCOL,NPR3,COEF,FIT,KK) IF(KK.EQ.4) CALL ESTIMA(DANA,DSTAR,PR4,LSV4,NCOL,NPR4,COEF,FIT,KK) C WF=(FIT**2)*VAR IK=OPTION(J) IF (IK.GT.10.AND.IK.LT.15) IK=IK-6 IF (IK.GT.20) IK=IK-12 VAF(IK)=VAF(IK)+WF PVAF(IK)=PVAF(IK)+VAR KCOUNT(IK)=KCOUNT(IK)+1 IF (IPRI1.EQ.1) WRITE (IWRITE,615) FIT,VAR,WF IF (IPRI1.GT.1) WRITE (IWRITE,617) I,J,FIT,VAR,WF FITT(J)=FIT C C COPY THE DATA, APPROPRIATELY NORMALIZED IN ARRAY DANA, C INTO ARRAY DISP (FOR TRANSFORMATION PLOTS) C DO 37 M=1,NCOL DISP(M)=DANA(M) 37 CONTINUE C C NON-METRIC REGRESSION; STANDARDIZE ORIGINAL DATA C IF (OPTION(J).LE.4) GO TO 36 RMEAN=0.0 VAR=0.0 DO 219 M=1,NCOL RMEAN=RMEAN+DATA(M) 219 CONTINUE RMEAN=RMEAN/FLOAT(NCOL) DO 229 M=1,NCOL DANA(M)=DATA(M)-RMEAN VAR=VAR+DANA(M)*DANA(M) 229 CONTINUE STAND=SQRT(VAR/FLOAT(NCOL)) SSQ=FLOAT(NCOL) DO 239 M=1,NCOL DANA(M)=DANA(M)/STAND DISP(M)=DANA(M) 239 CONTINUE IF (OPTION(J).LE.14) IAPP=1 IF (OPTION(J).GE.21) IAPP=2 IF (OPTION(J).EQ.11.OR.OPTION(J).EQ.21) X CALL MONO (DANA,DISP,PR1,LSV1,IORD,LBK,NCOL,NPR1,MAX,CRIT, X IAPP,IWRITE,IHIS,WORK2,FIT,SSQ) IF (OPTION(J).EQ.12.OR.OPTION(J).EQ.22) X CALL MONO (DANA,DISP,PR2,LSV2,IORD,LBK,NCOL,NPR2,MAX,CRIT, X IAPP,IWRITE,IHIS,WORK2,FIT,SSQ) IF (OPTION(J).EQ.13.OR.OPTION(J).EQ.23) X CALL MONO (DANA,DISP,PR3,LSV3,IORD,LBK,NCOL,NPR3,MAX,CRIT, X IAPP,IWRITE,IHIS,WORK2,FIT,SSQ) IF (OPTION(J).EQ.14.OR.OPTION(J).EQ.24) X CALL MONO (DANA,DISP,PR4,LSV4,IORD,LBK,NCOL,NPR4,MAX,CRIT, X IAPP,IWRITE,IHIS,WORK2,FIT,SSQ) C IF (OPTION(J).EQ.11.OR.OPTION(J).EQ.21) X CALL ESTIMA(DISP,DSTAR,PR1,LSV1,NCOL,NPR1,COEF,FIT,KK) IF (OPTION(J).EQ.12.OR.OPTION(J).EQ.22) X CALL ESTIMA(DISP,DSTAR,PR2,LSV2,NCOL,NPR2,COEF,FIT,KK) IF (OPTION(J).EQ.13.OR.OPTION(J).EQ.23) X CALL ESTIMA(DISP,DSTAR,PR3,LSV3,NCOL,NPR3,COEF,FIT,KK) IF (OPTION(J).EQ.14.OR.OPTION(J).EQ.24) X CALL ESTIMA(DISP,DSTAR,PR4,LSV4,NCOL,NPR4,COEF,FIT,KK) C IF (IPRI1.EQ.1) WRITE (IWRITE,616) FIT IF (IPRI1.GT.1) WRITE (IWRITE,618) FIT C 36 CONTINUE C IK=OPTION(J) IF (IK.GT.10.AND.IK.LT.15) IK=IK-6 IF (IK.GT.20) IK=IK-12 FITM(IK)=FITM(IK)+FIT SFITM(IK)=SFITM(IK)+FIT**2 C C UNRAVEL THE REGRESSION COEFFICIENTS INTO THE PARAMETERS C OF THE MODEL DESIRED C JJ=J CALL UNRAV(COEF,NPRM,COOR,WEIGH,ROTAT,OPTION(J),NDIM,ND,JJ,I, X IWRITE,WSIGN,IPRE,WORK1,WORK3,WORK2,WORK4,CMAX,IPRI1,SLOP,CEPT, X FIT) CALL RESUL2 (DANA,DISP,DSTAR,WORK2,WORK3,I,J,OPTION(J), X WORK4,NCOL,NCOL*2,IWRITE,IWRI2,LUN(4),ISHE,IPRI1,SLOP,CEPT) 20 CONTINUE C C STORE COORDINATES,WEIGHTS,ROTATION MATRICES C IF (IWRI1.EQ.0) GO TO 21 LU=LUN(1) DO 23 J=1,NCOL IF (I.EQ.1) WRITE(LU,679) J,(CON(J,K),K=1,NDIM) 23 CONTINUE KE=0 DO 24 J=1,NANA KB=KE+1 KE=KB+NDIM-1 WRITE(LU,680) I,J,(COOR(K),K=KB,KE) 24 CONTINUE IF (IWRI1.EQ.1) GO TO 21 LU=LUN(2) KE=0 DO 25 J=1,NANA KB=KE+1 KE=KB+NDIM-1 KK=MOD(OPTION(J),10) IF (KK.LT.4.AND.KK.GT.0) WRITE(LU,680) I,J,(WEIGH(K),K=KB,KE) 25 CONTINUE IF (IWRI1.EQ.2) GO TO 21 LU=LUN(3) KE=0 DO 26 J=1,NANA KB=KE+1 KE=KB+NDIM-1 KK=MOD(OPTION(J),10) DO 27 M=1,NDIM IF (KK.EQ.1) WRITE(LU,685) I,J,M,(ROTAT(M,K),K=KB,KE) 27 CONTINUE 26 CONTINUE C C COMPUTE STATISTICS FOR METRIC ANALYSES C 21 IF(IST.EQ.1) CALL STATIS(OPTION,NROW,NCOL,NDIM,FITT,NANA,I,IWRITE) C C IF (NPL1.EQ.1.AND.NROW.GT.1) GO TO 10 DO 30 N=1,ND TCOOR(I,N)=COOR(N) TWEIGH(I,N)=WEIGH(N) 30 CONTINUE IF (NPL2.EQ.NDIM.AND.NROW.GT.1) GO TO 10 DO 35 M=1,NDIM MM=M+(I-1)*NDIM DO 40 N=1,ND TROTAT(MM,N)=ROTAT(M,N) 40 CONTINUE 35 CONTINUE C C 10 CONTINUE C WRITE(IWRITE,634) DO 54 K=1,12 IF (KCOUNT(K).EQ.0) GO TO 54 SFITM(K)=SQRT(SFITM(K)/FLOAT(KCOUNT(K))) FITM(K)=FITM(K)/FLOAT(KCOUNT(K)) WRITE(IWRITE,6343) OPTLAB(K),KCOUNT(K),FITM(K),SFITM(K) IF (K.GT.4) GO TO 54 WRITE(IWRITE,6342) PVAF(K),VAF(K) PVAF(K)=VAF(K)/PVAF(K) WRITE(IWRITE,6341) PVAF(K) 54 CONTINUE NTOT=NROW+NCOL N=MAX0(IPLO,IPRI2) IF (N.GT.0) CALL OUTP3(TCOOR,CON,WORK2,WORK3,WORK4,IPLO,NROW, X NCOL,NDIM,NTOT,ND,NPL2,TWEIGH,TROTAT,NANA,HULP2,IPRI2, X OPTION,CMAX) C 500 FORMAT(20A4) 501 FORMAT(4A4) 629 FORMAT(1H ,4X,42HTHE TARGET CONFIGURATION WILL BE CENTERED:, X 1X,23HMEAN ORIGINAL DIMENSION,I5,F11.3) 600 FORMAT(1H0//,5X,31HTARGET CONFIGURATION (CENTERED)/, X 1H , 4X,31H-------------------------------) 601 FORMAT(1H0//,5X,41HTARGET CONFIGURATION AFTER TRANSFORMATION/ X 1H , 4X,41H-----------------------------------------) 604 FORMAT(1H0//,5X,I5,26H ROWS OF THE EXTERNAL DATA/ X 1H ,4X,31H-------------------------------/1H ,9X,10I9) 605 FORMAT(1H ,5X,I5,2X,10F9.3/(1H ,12X,10F9.3)) 608 FORMAT(1H+,51X,18H METRIC REGRESSION/1H ,4X, X 52H====================================================, X 13H=============) 609 FORMAT(1H+,51X,22H NON-METRIC REGRESSION) 611 FORMAT(1H0//,5X,4HROW=,I5,10H ANALYSIS=,I2, X 25H GENERAL UNFOLDING MODEL ) 612 FORMAT(1H0//,5X,4HROW=,I5,10H ANALYSIS=,I2, X 25H WEIGHTED UNFOLDING MODEL) 613 FORMAT(1H0//,5X,4HROW=,I5,10H ANALYSIS=,I2, X 25H UNFOLDING MODEL ) 614 FORMAT(1H0//,5X,4HROW=,I5,10H ANALYSIS=,I2, X 25H VECTOR MODEL ) 615 FORMAT(1H0/,16X,14HMETRIC FIT ,F11.3,10H VARIANCE,F11.3, X 10H V.A.F.,F11.3) 617 FORMAT(1H0,3X,3HROW,I5,7H OPTION,I2,14H METRIC FIT ,F6.3, X 10H VARIANCE,F11.3,10H V.A.F.,F11.3) 618 FORMAT(1H0,16X,15H NON-METRIC FIT,3X,F6.3) 616 FORMAT(1H0/,16X,14HNON-METRIC FIT,F11.3) 620 FORMAT(1H+,74X,27H PRIMARY APPROACH TO TIES/1H ,4X, X 52H====================================================, X 45H=============================================) 621 FORMAT(1H+,74X,27H SECONDARY APPROACH TO TIES/1H ,4X, X 52H====================================================, X 45H=============================================) 630 FORMAT(1H0/,5X,22H*** WARNING *** OPTION,I2,12H WILL NOT BE, X 51H APPLIED. EITHER THE MODEL HAS NOT BEEN REFERRED TO/ X 1H0,20X,9HON CARD 3, X 49H OR THE NUMBER OF TARGET POINTS IS NOT SUFFICIENT, X 19H TO FIT THIS MODEL.) 631 FORMAT(1H0//,19X,8HANALYSIS/1H ,4X,11HOPTION SET:,4I4) 632 FORMAT(1H0/,5X,35HMODEL: V = VECTOR / X 1H ,4X, 35H U = UNFOLDING / X 1H ,4X, 35H W = WEIGHTED UNFOLDING/ X 1H ,4X, 35H G = GENERAL UNFOLDING / X 1H ,4X, 35HREGRESSION: M = METRIC / X 1H ,4X, 29H P = NON-METRIC, , X 24HPRIMARY APPROACH TO TIES/ X 1H ,4X, 29H S = NON-METRIC, , X 26HSECONDARY APPROACH TO TIES) 633 FORMAT(1H ,10X,I5,2X,4A4) 634 FORMAT(1H0//,5X,41HSUMMARY OF RESULTS: (PROPORTION OF TOTAL , X 53HVARIANCE ACCOUNTED FOR ONLY GIVEN FOR METRIC OPTIONS)/) 6343 FORMAT(1H0,4X,9HMODEL ' ,1A4,5H' N =,I5,15H AVERAGE FIT =,F6.3/ X 18X,25H ROOT MEAN SQUARED FIT =,F6.3) 6342 FORMAT(1H ,4X,16HTOTAL VARIANCE =,F9.3,3X,14HTOTAL VARIANCE, X 16H ACCOUNTED FOR =,F9.3) 6341 FORMAT(1H ,4X,44HPROPORTION OF TOTAL VARIANCE ACCOUNTED FOR =, X F6.3) 635 FORMAT(1H0//,5X,30HTHE CONFIGURATION WILL BE READ, X 13H WITH FORMAT ,20A4/(1H ,4X,20A4)) 636 FORMAT(1H0//,5X,30HTHE EXTERNAL DATA WILL BE READ, X 13H WITH FORMAT ,20A4/(1H ,4X,20A4)) 679 FORMAT(I10,5F12.5/(10X,5F12.5)) 680 FORMAT(2I5,5F12.5/(10X,5F12.5)) 685 FORMAT(3I5,4F12.5/(15X,4F12.5)) 900 FORMAT(1H0//,5X,3HROW,I5,35H WILL BE OMITTED FROM THE ANALYSIS , X 26HSINCE IT HAS ZERO VARIANCE//) C C RETURN END SUBROUTINE MONO(DATA,DISP,PRM,LSV,IORD,LBK,NCOL,NPR,IMAX,CRIT, X IAPP,IWRITE,IHIS,HELP,FIT,DSTAND) C ****************************************************************** C * * C * M O N O * C * * C * PURPOSE: CONTROLS THE MONOTONE REGRESSION * C * * C * SUBROUTINES CALLED: RANK1,PROJEC,MONOR * C * * C ****************************************************************** DIMENSION DATA(NCOL),DISP(NCOL),IORD(NCOL),LBK(NCOL),HELP(NCOL) DIMENSION PRM(NCOL,NPR),LSV(NCOL,NPR) REAL LSV C OLDFIT=FIT NITER=0 IF (IHIS.EQ.1) WRITE (IWRITE,605) CALL RANK1(DATA,IORD,LBK,HELP,IAPP,NCOL) DO 10 NITER=1,IMAX CALL PROJEC(PRM,LSV,DISP,HELP,NCOL,NPR,STAND1) CALL MONOR(DISP,HELP,IORD,LBK,IAPP,NCOL,STAND2,DSTAND) FIT=STAND2/STAND1 DIFF=FIT-OLDFIT IF (IHIS.EQ.1) WRITE(IWRITE,610) NITER,FIT,DIFF IF (ABS(DIFF).LT.CRIT) GO TO 60 OLDFIT=FIT 10 CONTINUE 60 RETURN C 605 FORMAT (1H0/,20X,32HHISTORY OF NON-METRIC REGRESSION/, X 1H ,19X,32H--------------------------------/, X 1H ,50X,15HDIFFERENCE WITH/, X 1H ,19X,9HITERATION,7X,3HFIT,12X,19HPRECEDING ITERATION) 610 FORMAT (1H0,25X,I4,1X,F12.5,3X,F12.5) END SUBROUTINE MONOR(DISP,HELP,IORD,LBK,IAPP,NCOL,FIT,DSTAND) C ****************************************************************** C * * C * M O N O R * C * * C * PURPOSE: MONITOR FOR MONOTONE REGRESSION * C * * C * SUBROUTINES CALLED: PRAR,SCAR,MRMNH * C * * C ****************************************************************** DIMENSION DISP(NCOL),HELP(NCOL),IORD(NCOL),LBK(NCOL) REAL MEAN DO 10 I=1,NCOL HELP(I)=DISP(IORD(I)) 10 CONTINUE IF (IAPP.EQ.1) CALL PRAR(HELP,IORD,LBK,NCOL) IF (IAPP.EQ.2) CALL SCAR(HELP,LBK,NCOL) CALL MRMNH(HELP,NCOL) STAND=0.0 MEAN=0.0 FIT=0.0 DO 15 I=1,NCOL STAND=STAND+HELP(I)*HELP(I) MEAN=MEAN+HELP(I) 15 CONTINUE STAND=SQRT(DSTAND/STAND) MEAN=MEAN/FLOAT(NCOL) DO 20 I=1,NCOL DISP(IORD(I))=HELP(I)*STAND HELP(I)=HELP(I)-MEAN FIT=FIT+HELP(I)*HELP(I) 20 CONTINUE FIT=SQRT(FIT) RETURN END SUBROUTINE MRMNH (DISP,N) C ****************************************************************** C * * C * M R M N H * C * * C * PURPOSE: MONOTONE REGRESSION. * C * * C * SUBROUTINES CALLED: NONE. * C * * C * AUTHOR: ERNST VAN WANING. RELEASED: JULY 1976 * C * * C ****************************************************************** C C DIMENSION DISP(N) C LOVBKH = 0 LUPH = 0 DO 1 IUP = 2, N IF (DISP(IUP) .GE. DISP(IUP-1)) GOTO 1 SDS = DISP(IUP) SW = 1. IDOWN = IUP 5 IDOWN = IDOWN - 1 IF (LUPH .EQ. IDOWN) GOTO 2 SDS = SDS + DISP(IDOWN) SW = SW + 1. GOTO 3 2 SDS = SDS + (LOVBKH * DISP(LUPH)) SW = SW + LOVBKH IDOWN = IDOWN - LOVBKH + 1 3 TRIALV = SDS / SW IF (IDOWN .EQ. 1) GOTO 4 IF (DISP(IDOWN - 1) .GT. TRIALV) GOTO 5 4 DO 6 J = IDOWN, IUP 6 DISP(J) = TRIALV LOVBKH = IUP - IDOWN + 1 LUPH = IUP 1 CONTINUE RETURN END SUBROUTINE OUTP3(TCOOR,CON,H1,H2,IND,IPLO,NROW,NCOL,NDIM,NTOT,ND, X NPL2,TWEIGH,TROTAT,NANA,IWE,IPRI2,OPTION,CMAX) C ****************************************************************** C * * C * O U T P 3 * C * * C * PURPOSE: PROVIDES FOR THE PRINTING, AND PLOTTING OF RESULTS * C * * C * SUBROUTINES CALLED: PLOT * C * * C ****************************************************************** COMMON/CPARAT/INCON,INDATA,INOPT,INPARA,IWRITE,ISCR1,ISCR2,LUN(4) DIMENSION TCOOR(NROW,ND),CON(NCOL,NDIM),H1(NTOT),H2(NTOT),IWE(4), X IND(NTOT),TWEIGH(NROW,ND),TROTAT(NPL2,ND), X OPTION(NANA),OPTLAB(13) INTEGER OPTION,OPTLAB DATA OPTLAB/'GM ','WM ','UM ','VM ','GP ','WP ','UP ', X 'VP ','GS ','WS ','US ','VS ','N.A.'/ REWIND ISCR2 IF (IPRI2.EQ.0) GO TO 21 KE=0 DO 1 J=1,NANA KB=KE+1 KE=KB+NDIM-1 DO 5 I=1,NROW READ (ISCR2) (OPTION(JN),JN=1,NANA) IF (I.EQ.1) WRITE(IWRITE,600) J IO=OPTION(J) IF (IO.GT.10.AND.IO.LT.15) IO=IO-6 IF (IO.GT.20) IO=IO-12 IF (IO.EQ.0) IO=13 WRITE (IWRITE,603) I,OPTLAB(IO),(TCOOR(I,K),K=KB,KE) 5 CONTINUE REWIND ISCR2 IC=0 DO 10 I=1,NROW READ (ISCR2) (OPTION(JN),JN=1,NANA) IO=OPTION(J) IF (IO.GT.10.AND.IO.LT.15) IO=IO-6 IF (IO.GT.20) IO=IO-12 IF (IO.EQ.0) IO=13 IK=MOD(OPTION(J),10) IF (IK.LT.4.AND.IK.GT.0) IC=IC+1 IF (IC.EQ.1) WRITE(IWRITE,601) J IF (IC.EQ.1) IC=IC+1 IF (IK.LT.4.AND.IK.GT.0) WRITE(IWRITE,603) X I,OPTLAB(IO),(TWEIGH(I,K),K=KB,KE) 10 CONTINUE REWIND ISCR2 IC=0 MM=1 DO 15 I=1,NROW READ (ISCR2) (OPTION(JN),JN=1,NANA) IO=OPTION(J) IF (IO.GT.10.AND.IO.LT.15) IO=IO-6 IF (IO.GT.20) IO=IO-12 IF (IO.EQ.0) IO=13 IK=MOD(OPTION(J),10) IF (IK.LT.2.AND.IK.GT.0) IC=IC+1 IF (IC.EQ.1) WRITE(IWRITE,602) J IF (IC.EQ.1) IC=IC+1 DO 19 M=1,NDIM IF (IK.LT.2.AND.IK.GT.0.AND.M.EQ.1) WRITE(IWRITE,604) X I,OPTLAB(IO),M,(TROTAT(MM,K),K=KB,KE) IF (IK.LT.2.AND.IK.GT.0.AND.M.GT.1) WRITE(IWRITE,614) X M,(TROTAT(MM,K),K=KB,KE) MM=MM+1 19 CONTINUE 15 CONTINUE REWIND ISCR2 1 CONTINUE 21 IF (IPLO.EQ.0) RETURN DO 25 M=1,NANA NLB=(M-1)*NDIM NDPL=MIN0(NDIM,IPLO) IF (NDPL.EQ.1) GO TO 3 DO 20 I=2,NDPL IMIN=I-1 DO 30 J=1,IMIN JB=J+NLB IB=I+NLB NA=0 35 L=0 NC=0 NC2=0 DO 40 K=1,NROW READ (ISCR2) (OPTION(JN),JN=1,NANA) IO=OPTION(M) IF (IO.GT.10.AND.IO.LT.15) IO=IO-6 IF (IO.GT.20) IO=IO-12 IF (IO.EQ.0) IO=13 NOUT=0 NT=1-NA IF (TWEIGH(K,JB).LE.0.0.OR.TWEIGH(K,IB).LE.0.0) X NT=NA NC=NC+NT NC2=NC2+NT L=L+1 H1(L)=TCOOR(K,JB)*FLOAT(NT) H2(L)=TCOOR(K,IB)*FLOAT(NT) TDIS=SQRT(H1(L)**2+H2(L)**2) IF (TDIS.GT.2*CMAX) NOUT=1 IF (NOUT.EQ.1) H1(L)=0.0 IF (NOUT.EQ.1) H2(L)=0.0 IF (NOUT.EQ.1) NC2=NC2-1 IF (NOUT.EQ.1) WRITE (IWRITE,605) M,K, X OPTLAB(IO),TCOOR(K,JB),TCOOR(K,IB) 40 CONTINUE REWIND ISCR2 DO 50 K=1,NCOL L=L+1 H1(L)=CON(K,J) H2(L)=CON(K,I) 50 CONTINUE IF (NC.EQ.0) GO TO 51 IF (NC2.EQ.0) GO TO 51 WRITE (IWRITE,110) M,J,I CALL PRPLOT(H1,H2,IND,NTOT,NROW,1,IWRITE,1) IF (NC.EQ.NROW) GO TO 30 51 NA=1 L=0 NC=0 NC2=0 DO 52 K=1,NROW READ (ISCR2) (OPTION(JN),JN=1,NANA) IO=OPTION(M) IF (IO.GT.10.AND.IO.LT.15) IO=IO-6 IF (IO.GT.20) IO=IO-12 IF (IO.EQ.0) IO=13 NOUT=0 NT=1-NA IF (TWEIGH(K,JB).LT.0.0.OR.TWEIGH(K,IB).LT.0.0) X NT=NA NC=NC+NT NC2=NC2+NT L=L+1 H1(L)=TCOOR(K,JB)*FLOAT(NT) H2(L)=TCOOR(K,IB)*FLOAT(NT) TDIS=SQRT(H1(L)**2+H2(L)**2) IF (TDIS.GT.2*CMAX) NOUT=1 IF (NOUT.EQ.1) H1(L)=0.0 IF (NOUT.EQ.1) H2(L)=0.0 IF (NOUT.EQ.1) NC2=NC2-1 IF (NOUT.EQ.1) WRITE (IWRITE,605) M,K, X OPTLAB(IO),TCOOR(K,JB),TCOOR(K,IB) 52 CONTINUE REWIND ISCR2 IF (NC.EQ.0) GO TO 30 IF (NC2.EQ.0) GO TO 30 WRITE (IWRITE,111) M,J,I CALL PRPLOT(H1,H2,IND,NTOT,NROW,1,IWRITE,1) 30 CONTINUE 20 CONTINUE IF (IWE(M).LE.1) GO TO 25 DO 61 I=2,NDPL IMIN=I-1 DO 62 J=1,IMIN JB=J+NLB IB=I+NLB WRITE (IWRITE,210) M,J,I CALL PRPLOT(TWEIGH(1,JB),TWEIGH(1,IB),IND,NROW,NROW, X 1,IWRITE,1) 62 CONTINUE 61 CONTINUE GO TO 25 C C C ONE-DIMENSIONAL PLOT C C C 3 L=0 NA=0 NC=0 NC2=0 JB=NLB+1 DO 63 K=1,NROW READ (ISCR2) (OPTION(JN),JN=1,NANA) IO=OPTION(M) IF (IO.GT.10.AND.IO.LT.15) IO=IO-6 IF (IO.GT.20) IO=IO-12 IF (IO.EQ.0) IO=13 NOUT=0 NT=1-NA IF (TWEIGH(K,JB).LE.0.0) NT=NA L=L+1 NC=NC+NT NC2=NC2+NT H1(L)=TCOOR(K,JB)*FLOAT(NT) H2(L)=0.0 TDIS=ABS(H1(L)) IF (TDIS.GT.2*CMAX) NOUT=1 IF (NOUT.EQ.1) H1(L)=0.0 IF (NOUT.EQ.1) NC2=NC2-1 IF (NOUT.EQ.1) WRITE (IWRITE,606) M,K, X OPTLAB(IO),TCOOR(K,JB) 63 CONTINUE REWIND ISCR2 DO 64 K=1,NCOL L=L+1 H1(L)=CON(K,1) H2(L)=0.0 64 CONTINUE IF (NC.EQ.0) GO TO 26 IF (NC2.EQ.0) GO TO 26 WRITE (IWRITE,120) M CALL PRPLOT(H1,H2,IND,NTOT,NROW,1,IWRITE,1) IF (NC.EQ.NROW) GO TO 69 26 NA=1 L=0 NC=0 NC2=0 DO 65 K=1,NROW READ (ISCR2) (OPTION(JN),JN=1,NANA) IO=OPTION(M) IF (IO.GT.10.AND.IO.LT.15) IO=IO-6 IF (IO.GT.20) IO=IO-12 IF (IO.EQ.0) IO=13 NOUT=0 NT=1-NA IF (TWEIGH(K,JB).LT.0.0) NT=NA L=L+1 NC=NC+NT NC2=NC2+NT H1(L)=TCOOR(K,JB)*FLOAT(NT) H2(L)=0.0 TDIS=ABS(H1(L)) IF (TDIS.GT.2*CMAX) NOUT=1 IF (NOUT.EQ.1) H1(L)=0.0 IF (NOUT.EQ.1) NC2=NC2-1 IF (NOUT.EQ.1) WRITE (IWRITE,606) M,K, X OPTLAB(IO),TCOOR(K,JB) 65 CONTINUE REWIND ISCR2 IF (NC.EQ.0) GO TO 69 IF (NC2.EQ.0) GO TO 69 WRITE(IWRITE,121) M CALL PRPLOT(H1,H2,IND,NTOT,NROW,1,IWRITE,1) 69 IF (IWE(NANA).LE.1) GO TO 25 DO 70 K=1,NROW H2(K)=0.0 70 CONTINUE WRITE(IWRITE,220) M CALL PRPLOT(TWEIGH(1,JB),H2,IND,NROW,NROW,1,IWRITE,1) 25 CONTINUE C RETURN C 600 FORMAT(1H0///,12X,8HANALYSIS,I2,12H COORDINATES/1H ,11X, X 22H----------------------) 601 FORMAT(1H0//,12X,8HANALYSIS,I2,8H WEIGHTS/1H ,11X, X 18H------------------) 602 FORMAT(1H0//,12X,8HANALYSIS,I2,18H ROTATION MATRICES/1H ,11X, X 28H----------------------------) 603 FORMAT(1H ,25X,I6,5X,A4,10F9.3/(1H ,40X,10F9.3)) 604 FORMAT(1H ,25X,I6,5X,A4,I2,10F7.3/(1H ,42X,10F7.3)) 614 FORMAT(1H ,40X,I2,10F7.3/(1H ,42X,10F7.3)) 605 FORMAT(1H0//,6X,8HANALYSIS,I2,20H *** NOT PLOTTED ***,I6,5X,A4, X 2F9.3,30H *** SEVERELY OUT OF RANGE ***) 606 FORMAT(1H0//,6X,8HANALYSIS,I2,20H *** NOT PLOTTED ***,I6,5X,A4, X F9.3,30H *** SEVERELY OUT OF RANGE ***) 110 FORMAT(1H1,3X,38HIDEAL POINTS AND/OR VECTORS (INTEGERS), X 26H IN TARGET CONFIGURATION. , X 8HANALYSIS,I2,4H DIM,I2,16H (X-AXIS) VS DIM,I2,9H (Y-AXIS)) 111 FORMAT(1H1,3X,42HANTI-IDEAL AND/OR SADDLE POINTS (INTEGERS), X 26H IN TARGET CONFIGURATION. , X 8HANALYSIS,I2,4H DIM,I2,16H (X-AXIS) VS DIM,I2,9H (Y-AXIS)) 120 FORMAT(1H1,3X,43HONE-DIMENSIONAL PLOT OF IDEAL POINTS AND/OR, X 32H VECTORS (INTEGERS) FOR ANALYSIS,I2) 121 FORMAT(1H1,3X,48HONE-DIMENSIONAL PLOT OF ANTI-IDEAL AND/OR SADDLE, X 20H POINTS FOR ANALYSIS,I2) 210 FORMAT(1H1,3X,28HPLOT OF WEIGHTS FOR ANALYSIS, 1 I2,10H DIMENSION,I2,26H (X-AXIS) VERSUS DIMENSION,I2, 2 9H (Y-AXIS)) 220 FORMAT(1H1,3X,31HONE-DIMENSIONAL PLOT OF WEIGHTS, X 13H FOR ANALYSIS,I2) END SUBROUTINE PRAR(DIST,IORD,LBK,N) C ****************************************************************** C * * C * P R A R * C * * C * PURPOSE: TO SORT THE DISTANCES WITHIN TIEBLOCKS AND REORDER * C * THEIR INDICES IN IORD ACCORDINGLY * C * * C * SUBROUTINES CALLED: SHEL9 * C * * C ****************************************************************** DIMENSION DIST(N),IORD(N),LBK(N) C IDLBK=0 JB=1 1 IDLBK=IDLBK+1 NLBK=LBK(IDLBK) IF (NLBK.GT.1) CALL SHEL9(DIST(JB),IORD(JB),NLBK) JB=JB+NLBK IF (JB.LT.N) GO TO 1 RETURN END SUBROUTINE PRED1(CON,PRM,LSV1,NCOL,NDIM,NPR1) C ****************************************************************** C * * C * P R E D 1 * C * * C * PURPOSE: CONSTRUCTS THE PREDICTOR MATRIX FOR THE * C * GENERAL UNFOLDING MODEL * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION CON(NCOL,NDIM),PRM(NCOL,NPR1),LSV1(NCOL,NPR1) REAL LSV1 C DO 10 J=1,NCOL PRM(J,1)=1.0 LSV1(J,1)=PRM(J,1) DO 20 I=1,NDIM L=I+1 LL=L+NDIM PRM(J,L)=CON(J,I) LSV1(J,L)=PRM(J,L) PRM(J,LL)=CON(J,I)**2 LSV1(J,LL)=PRM(J,LL) 20 CONTINUE LB=2 LE=NDIM-1 L=2*NDIM+2 DO 30 M=1,LE DO 40 I=LB,NDIM PRM(J,L)=CON(J,M)*CON(J,I) LSV1(J,L)=PRM(J,L) L=L+1 40 CONTINUE LB=LB+1 30 CONTINUE 10 CONTINUE RETURN END SUBROUTINE PRED2(CON,PRM,LSV2,NCOL,NDIM,NPR2) C ****************************************************************** C * * C * P R E D 2 * C * * C * PURPOSE: CONSTRUCTS THE PREDICTOR MATRIX FOR THE * C * WEIGHTED UNFOLDING MODEL * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION CON(NCOL,NDIM),PRM(NCOL,NPR2),LSV2(NCOL,NPR2) REAL LSV2 C DO 10 J=1,NCOL PRM(J,1)=1.0 LSV2(J,1)=1.0 DO 20 I=1,NDIM L=I+1 LL=L+NDIM PRM(J,L)=CON(J,I) LSV2(J,L)=PRM(J,L) PRM(J,LL)=CON(J,I)**2 LSV2(J,LL)=PRM(J,LL) 20 CONTINUE 10 CONTINUE RETURN END SUBROUTINE PRED3(CON,PRM,LSV3,WSIGN,NCOL,NDIM,NPR3) C ****************************************************************** C * * C * P R E D 3 * C * * C * PURPOSE: CONSTRUCTS THE PREDICTOR MATRIX FOR THE * C * UNFOLDING MODEL * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION CON(NCOL,NDIM),PRM(NCOL,NPR3),LSV3(NCOL,NPR3) DIMENSION WSIGN(NDIM) REAL LSV3 C DO 10 J=1,NCOL PRM(J,1)=1.0 LSV3(J,1)=1.0 SSQ=0. DO 20 I=1,NDIM L=I+1 PRM(J,L)=CON(J,I) LSV3(J,L)=PRM(J,L) SSQ=SSQ+WSIGN(I)*CON(J,I)**2 20 CONTINUE L=NDIM+2 PRM(J,L)=SSQ LSV3(J,L)=SSQ 10 CONTINUE RETURN END SUBROUTINE PRED4(CON,PRM,LSV4,NCOL,NDIM,NPR4) C ****************************************************************** C * * C * P R E D 4 * C * * C * PURPOSE: CONSTRUCTS THE PREDICTOR MATRIX FOR THE * C * VECTOR MODEL * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION CON(NCOL,NDIM),PRM(NCOL,NPR4),LSV4(NCOL,NPR4) REAL LSV4 C DO 10 J=1,NCOL PRM(J,1)=1.0 LSV4(J,1)=1.0 DO 20 I=1,NDIM L=I+1 PRM(J,L)=CON(J,I) LSV4(J,L)=PRM(J,L) 20 CONTINUE 10 CONTINUE RETURN END SUBROUTINE PREPRO(CON,CROSS,U,UH,V,PSEU,WORK1,WORK4,W,WNEG, X IPRE,IWRITE,NCOL,NDIM,NPR1,NPR4,NPRM) C ****************************************************************** C * * C * P R E P R O * C * * C * PURPOSE: PREPROCESSING THE TARGET CONFIGURATION * C * IF IPRE EQUALS 1, ONLY WEIGHTING OF AXES * C * IF IPRE EQUALS 2, CANONICAL REFERENCE FRAME * C * * C * SUBROUTINES CALLED: PRED4,PSINV,TRED2,IMTQL2,PREREC * C * * C ****************************************************************** DIMENSION CON(NCOL,NDIM),CROSS(NCOL,NCOL),U(NCOL,NPR1), X UH(NCOL,NPR1),V(NCOL,NPR4),PSEU(NCOL,NPR4), X WORK1(NPR1,NPR1),WORK4(NPR4,NPR4),W(NPRM),WNEG(NDIM) C C CHECK OPTION AND SET UP U,V AND PSEU C IF (IPRE.EQ.1) NRAV=NDIM IF (IPRE.EQ.2) NRAV=NPR1-(NDIM+1) CALL PRED4(CON,V,PSEU,NCOL,NDIM,NPR4) CALL PSINV(PSEU,WORK4,W,NPR4,NCOL,WORK1) DO 10 J=1,NCOL DO 20 L=1,NDIM U(J,L)=CON(J,L)*CON(J,L) 20 CONTINUE IF (IPRE.EQ.1) GO TO 10 ME=NDIM-1 KB=2 L=NDIM+1 DO 30 M=1,ME DO 40 K=KB,NDIM U(J,L)=2*CON(J,M)*CON(J,K) L=L+1 40 CONTINUE KB=KB+1 30 CONTINUE 10 CONTINUE C C FORM UH C DO 50 I=1,NCOL DO 60 L=1,NRAV UH(I,L)=0.0 DO 70 J=1,NCOL DO 80 K=1,NPR4 UH(I,L)=UH(I,L)+PSEU(I,K)*V(J,K)*U(J,L) 80 CONTINUE 70 CONTINUE UH(I,L)=U(I,L)-UH(I,L) 60 CONTINUE 50 CONTINUE C C C USE U AS WORKING ARRAY FOR INVERSE SQUARE ROOT OF UHU C DO 90 L=1,NRAV DO 100 M=1,NRAV WORK1(L,M)=0.0 DO 110 I=1,NCOL WORK1(L,M)=WORK1(L,M)+UH(I,L)*UH(I,M) 110 CONTINUE 100 CONTINUE 90 CONTINUE CALL TRED2(NPR1,NRAV,W,WORK4,WORK1) CALL IMTQL2(NPR1,NRAV,W,WORK4,WORK1,IERR) IF (IERR.GT.0) GO TO 900 SMALL=0.0000001 DO 120 L=1,NRAV IF (W(L).LE.SMALL) W(L)=0.0 IF (W(L).GT.SMALL) W(L)=1.0/SQRT(W(L)) 120 CONTINUE DO 130 L=1,NRAV DO 140 M=1,NRAV U(L,M)=0.0 DO 150 N=1,NRAV U(L,M)=U(L,M)+WORK1(L,N)*W(N)*WORK1(M,N) 150 CONTINUE 140 CONTINUE 130 CONTINUE C C C FORM MATRIX FOR GENERALIZED EIGENVECTOR PROBLEM C DO 160 L=1,NRAV DO 170 M=1,NRAV WORK1(L,M)=0.0 DO 180 K1=1,NRAV DO 190 I=1,NCOL DO 200 J=1,NCOL DO 210 K2=1,NRAV P=U(L,K1)*UH(I,K1)*CROSS(I,J)*UH(J,K2)*U(K2,M) WORK1(L,M)=WORK1(L,M)+P 210 CONTINUE 200 CONTINUE 190 CONTINUE 180 CONTINUE 170 CONTINUE 160 CONTINUE C C C FIND SOLUTION FOR W C CALL TRED2(NPR1,NRAV,W,WORK4,WORK1) CALL IMTQL2(NPR1,NRAV,W,WORK4,WORK1,IERR) IF (IERR.GT.0) GO TO 900 SUM=0.0 DO 220 L=1,NRAV W(L)=0.0 DO 230 M=1,NRAV W(L)=W(L)+U(L,M)*WORK1(M,1) 230 CONTINUE IF (L.LE.NDIM) SUM=SUM+W(L) 220 CONTINUE DO 235 L=1,NRAV W(L)=W(L)*SIGN(1.0,SUM) 235 CONTINUE C IF (IPRE.EQ.2) GO TO 260 C C FINISH HERE IF WEIGHTING ONLY C STAND=0.0 DO 236 L=1,NRAV IF (L.LE.NDIM) STAND=STAND+ABS(W(L)) 236 CONTINUE STAND=FLOAT(NDIM)/STAND DO 240 L=1,NDIM WNEG(L)=SIGN(1.,W(L)) W(L)=SQRT(WNEG(L)*W(L)*STAND) DO 250 I=1,NCOL CON(I,L)=CON(I,L)*W(L) 250 CONTINUE 240 CONTINUE WRITE (IWRITE,601) WRITE (IWRITE,602) (W(L),L=1,NDIM) RETURN C C FINISH HERE FOR FULL TRANSFORMATION C 260 CONTINUE DO 270 L=1,NDIM WORK4(L,L)=W(L) 270 CONTINUE ME=NDIM-1 KB=2 L=NDIM+1 DO 280 M=1,ME DO 290 K=KB,NDIM WORK4(M,K)=W(L) WORK4(K,M)=W(L) L=L+1 290 CONTINUE KB=KB+1 280 CONTINUE CALL TRED2(NPR4,NDIM,W,WORK1,WORK4) CALL IMTQL2(NPR4,NDIM,W,WORK1,WORK4,IERR) IF (IERR.GT.0) GO TO 900 C C TRANSMIT SIGNS OF AXES THROUGH WNEG AND PERFORM REGULAR C TRANSFORMATION HERE. C STAND=0.0 DO 295 L=1,NDIM STAND=STAND+ABS(W(L)) 295 CONTINUE STAND=FLOAT(NDIM)/STAND DO 300 L=1,NDIM WNEG(L)=SIGN(1.,W(L)) W(L)=SQRT(WNEG(L)*W(L)*STAND) 300 CONTINUE WRITE (IWRITE,601) WRITE (IWRITE,602) (W(L),L=1,NDIM) WRITE (IWRITE,603) CALL PREREC(WORK4,NDIM,NDIM,NPR4,IWRITE,0) DO 310 I=1,NCOL DO 320 L=1,NDIM WORK1(1,L)=0.0 DO 330 M=1,NDIM WORK1(1,L)=WORK1(1,L)+CON(I,M)*WORK4(M,L) 330 CONTINUE 320 CONTINUE DO 340 L=1,NDIM CON(I,L)=WORK1(1,L)*W(L) 340 CONTINUE 310 CONTINUE RETURN C C TROUBLE REPORT C 900 WRITE(IWRITE,901) 901 FORMAT(1H1,33HT R O U B L E , X /1H0,34HONE OF THE EIGENVALUE ROUTINES , X /1H0,34HIN TRANSFORMATION STEP BREAKS DOWN, X /1H0,34HCONFIGURATION LEFT UNCHANGED ) 601 FORMAT(1H0//,5X,32HNORMALIZED WEIGHTS PER DIMENSION, X /1H , 4X,32H--------------------------------) 602 FORMAT(1H ,12X,12F9.3) 603 FORMAT(1H0///5X,15HROTATION MATRIX/1H ,4X,15H---------------) RETURN END SUBROUTINE PREREC(X,M,N,MM,IWRITE,NR) C ****************************************************************** C * * C * P R E F R E C * C * * C * PURPOSE: PRINTS A M*N ARRAY WITH ROW- AND COLUMN IDENTIFI- * C * CATION, EACH LINE CONTAINING AT MOST 10 VALUES (F9.3)* C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION X(MM,N) C C C CHECK HOW MANY TIMES THE NUMBER OF COLUMNS EXCEEDS 10 C C C NTRU=N/10 NREM=N-NTRU*10 IF (NREM.GT.0) GO TO 1 NREM=10 NTRU=NTRU-1 1 NREP=NTRU+1 C C C PRINT NREP TIMES A BLOCK OF ELEMENTS C C C NK=10 KB=1 DO 2 I=1,NREP IF (I.EQ.NREP) NK=NREM KE=KB+NK-1 WRITE(IWRITE,60) (K,K=KB,KE) WRITE(IWRITE,61) DO 3 L=1,M IF (NR.EQ.0) WRITE(IWRITE,62) L,(X(L,K),K=KB,KE) IF (NR.NE.0) WRITE(IWRITE,62) NR,(X(L,K),K=KB,KE) 3 CONTINUE KB=KB+NK 2 CONTINUE RETURN 60 FORMAT(1H0,9X,10I9) 61 FORMAT(1H ) 62 FORMAT(1H ,I10,2X,10F9.3) END SUBROUTINE PROJEC(PRM,LSV,DISP,WORK,NCOL,NPR,STAND) C ****************************************************************** C * * C * P R O J E C T * C * * C * PURPOSE: COMPUTES THE PROJECTION OF THE PREDICTAND ON THE * C * PREDICTOR * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION DISP(NCOL),PRM(NCOL,NPR),LSV(NCOL,NPR),WORK(NCOL) REAL LSV,MEAN C DO 5 I=1,NCOL WORK(I)=0.0 5 CONTINUE DO 10 K=1,NCOL DO 20 I=1,NCOL SUM=0.0 DO 30 J=1,NPR SUM=SUM+PRM(K,J)*LSV(I,J)*DISP(K) 30 CONTINUE WORK(I)=WORK(I)+SUM 20 CONTINUE 10 CONTINUE MEAN=0.0 STAND=0.0 DO 40 I=1,NCOL DISP(I)=WORK(I) MEAN=MEAN+DISP(I) 40 CONTINUE MEAN=MEAN/FLOAT(NCOL) DO 50 I=1,NCOL A=WORK(I)-MEAN STAND=STAND+A*A 50 CONTINUE STAND=SQRT(STAND) RETURN END SUBROUTINE PRPLOT(X,Y,IND,NITEM,NRC,IDENT,IWRITE,NOR) C ****************************************************************** C * * C * P R P L O T * C * * C * PURPOSE: PRINTPLOT OF Y VERSUS X * C * IF IDENT.GT.0, THE FIRST NRC POINTS ARE IDENTIFIED BY* C * INTEGER NUMBERS 1-9,0,1-9... , THE REST OF THE POINTS* C * BY CHARACTERS A-I,J,A-I... * C * IF IDENT.EQ.0, ALL POINTS ARE PRINTED AS STARS * C * IF IDENT.LT.0, THE FIRST NRC POINTS ARE IDENTIFIED BY* C * STARS, THE OTHERS BY A 'D' * C * IF NOR.EQ.1, THE ORIGIN WILL BE PLOTTED * C * * C * SUBROUTINES CALLED: NONE * C * * C * AUTHOR: WILLEM HEISER RELEASED: MARCH 1978 * C * ADAPTED : NOVEMBER 1982 * C * ADAPTED : DECEMBER 1984 * C * * C ****************************************************************** C C TO ADJUST THE PLOT CHANGE THE CONSTANTS UNDER WHICH C ## IS PLACED ACCORDING TO THE COMMENTS C DIMENSION LINE(71),XPR(11),X(NITEM),Y(NITEM),LABEL(20),IND(NITEM) C ##=NN+1 (NN=INT(NC/7)) C ##=NUMBER OF COLUMNS PER LINE=NC INTEGER BLANK,STAR DATA BLANK,STAR,MORE/1H ,1H*,1HM/,NL/56/,NC/71/,NN/10/ C ##=NC##=NN C ##=NUMBER OF LINES PER PAGE DATA LABEL/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9, X 1HJ,1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI/ C 1 FORMAT(16X,3H.--,10(7H+------),5H+---.) C ##=NN 2 FORMAT(8X,F7.3,1X,1HI,2X,71A1,3X,1HI) C ##=NC C ##=NC+5 4 FORMAT(15X,11F7.3) C ##=NN+1 C DO 10 I=1,NITEM 10 IND(I)=I C C THE INDICES ARE ORDERED SUCH THAT THEY WOULD ORDER Y IN DESCENDING C ORDER C M=NITEM 20 M=M/2 IF (M.LT.1) GO TO 40 K=NITEM-M J=1 41 I=J 49 L=I+M IF (Y(IND(I)).GE.Y(IND(L))) GO TO 60 II=IND(I) IND(I)=IND(L) IND(L)=II I=I-M IF (I.GE.1) GO TO 49 60 J=J+1 IF (J-K) 41,41,20 40 CONTINUE C C THE PLOT IS ADJUSTED TO THE RANGE OF THE 'LONGEST' AXIS C XMAX=X(1) XMIN=X(1) DO 70 I=2,NITEM IF (X(I).GT.XMAX) XMAX=X(I) IF (X(I).LT.XMIN) XMIN=X(I) 70 CONTINUE IF (XMIN.GT. 0.0.AND.NOR.EQ.1) XMIN=0.0 IF (XMAX.LT. 0.0.AND.NOR.EQ.1) XMAX=0.0 SPANX=XMAX-XMIN YMAX=Y(IND(1)) IF (YMAX.LT. 0.0.AND.NOR.EQ.1) YMAX=0.0 YMIN=Y(IND(NITEM)) IF (YMIN.GT. 0.0.AND.NOR.EQ.1) YMIN=0.0 SPANY=YMAX-YMIN IF (SPANY.LE.SPANX) GO TO 75 XMIN=XMIN-(SPANY-SPANX)/2.0 SPANX=SPANY 75 STEPX=SPANX/(7.0*NN) STEPY=(STEPX*NC)/(NL+0.) HSTEP=STEPY/2.0 TOP=(YMIN+YMAX+SPANX)/2.0 C C THE POINTS ARE PLOTTED LINE AFTER LINE C WRITE(IWRITE,1) L=1 I=1 YPR=TOP+STEPY IF (NOR.EQ.1) IP=(0.0-XMIN)/STEPX+1.0 80 YPR=YPR-STEPY 90 DO 95 J=1,NC 95 LINE(J)=BLANK IF ((YPR-HSTEP.LE.0.0.AND.YPR+HSTEP.GT.0.0).AND.NOR.EQ.1) X LINE(IP)=STAR IF (I.GT.NITEM) GO TO 110 IF ((YPR-Y(IND(I))).GT.HSTEP) GO TO 110 100 JP=(X(IND(I))-XMIN)/STEPX+1.0 IF (LINE(JP).EQ.BLANK) GO TO 101 LINE(JP)=MORE GO TO 103 101 IF (IDENT.EQ.0) GO TO 102 IF (IDENT.GT.0) GO TO 104 IF (IND(I).LE.NRC) LINE(JP)=STAR IF (IND(I).GT.NRC) LINE(JP)=LABEL(15) GO TO 103 104 IF (IND(I).LE.NRC) LINE(JP)=LABEL(MOD(IND(I),10)+1) IF (IND(I).GT.NRC) X LINE(JP)=LABEL(MOD((IND(I)-NRC),10)+11) GO TO 103 102 LINE(JP)=STAR 103 I=I+1 IF (I.GT.NITEM) GO TO 105 IF ((YPR-Y(IND(I))).LT.HSTEP) GO TO 100 105 WRITE(IWRITE,2) YPR,LINE GO TO 115 110 WRITE(IWRITE,2) YPR,LINE 115 L=L+1 IF (L.LE.NL) GO TO 80 WRITE(IWRITE,1) C C VALUES OF X-AXIS ARE PRINTED C XPR(1)=XMIN DO 130 J=1,10 130 XPR(J+1)=XPR(J)+STEPX*7.0 WRITE(IWRITE,4) XPR RETURN END SUBROUTINE PSINV(LSV,RSV,SIVA,NPR,NCOL,WORK) C ****************************************************************** C * * C * P S I N V * C * * C * PURPOSE: COMPUTES THE PSEUDO INVERSE * C * * C * -1 * C * (A'A) A' ON THE BASIS OF THE * C * * C * SINGULAR VALUE DECOMPOSITION A = KSL' BY TAKING * C * * C * ******************** * C * -1 -1 * C * (A'A) A' = LS K' * C * ************ * C * * C * SUBROUTINES CALLED: SIVAD * C * * C ****************************************************************** DIMENSION LSV(NCOL,NPR),RSV(NPR,NPR),SIVA(NPR),WORK(NPR) REAL LSV C IF (NPR.GT.1) CALL SIVAD X (NCOL,NPR,SIVA,LSV,NCOL*NPR,RSV,NPR*NPR,WORK,ICV) C DO 5 I=1,NPR SMALL=.000001 IF(SIVA(I).GE.SMALL) SIVA(I)=1.0/SIVA(I) IF(SIVA(I).LT.SMALL) SIVA(I)=0. 5 CONTINUE C C DO 10 K=1,NCOL DO 20 I=1,NPR SUM=0.0 DO 30 J=1,NPR SUM=SUM+RSV(I,J)*SIVA(J)*LSV(K,J) 30 CONTINUE WORK(I)=SUM 20 CONTINUE DO 40 I=1,NPR LSV(K,I)=WORK(I) 40 CONTINUE 10 CONTINUE RETURN END SUBROUTINE RANK1(DATA,IORD,LBK,HELP,IAPP,NCOL) C ****************************************************************** C * * C * R A N K 1 * C * * C * PURPOSE: CONSTRUCT IORD AND LBK FROM DATA. * C * * C * SUBROUTINES CALLED: SHEL9, TIEBL * C * * C ****************************************************************** DIMENSION DATA(NCOL),IORD(NCOL),LBK(NCOL),HELP(NCOL) DO 10 I=1,NCOL HELP(I)=DATA(I) IORD(I)=I LBK(I)=0 10 CONTINUE CALL SHEL9(HELP,IORD,NCOL) CALL TIEBL(HELP,LBK,NCOL,IDLBK) IF (IDLBK.EQ.NCOL) IAPP=0 RETURN END SUBROUTINE RESUL1(COOR,WEIGH,ROTAT,NDIM,ND,NOP,IWRITE,OPT,NR,IP) C ****************************************************************** C * * C * R E S U L 1 * C * * C * PURPOSE: PRINTS RESULTS PER SUBJECT * C * * C * SUBROUTINES CALLED: PREREC * C * * C ****************************************************************** DIMENSION COOR(ND),WEIGH(ND),ROTAT(NDIM,ND) INTEGER OPT C NOPT=NOP-1 C WRITE(IWRITE,600) ILB=NOPT*NDIM+1 IF (OPT.NE.3) WRITE(IWRITE,606) IF (OPT.NE.3) GO TO 5 IF (IP.EQ.0.AND.WEIGH(ILB).GE.0.0) WRITE(IWRITE,603) IF (IP.EQ.0.AND.WEIGH(ILB).LT.0.0) WRITE(IWRITE,604) IF (IP.EQ.0) GO TO 5 WRITE (IWRITE,605) K=ILB KE=ILB+NDIM-1 WRITE(IWRITE,607) (WEIGH(K),K=ILB,KE) WRITE(IWRITE,608) 5 CALL PREREC(COOR(ILB),1,NDIM,1,IWRITE,NR) IF (OPT.LT.3) WRITE(IWRITE,601) IF (OPT.LT.3) CALL PREREC(WEIGH(ILB),1,NDIM,1,IWRITE,NR) IF (OPT.LT.2) WRITE(IWRITE,602) IF (OPT.LT.2) CALL PREREC(ROTAT(1,ILB),NDIM,NDIM,NDIM,IWRITE,0) C 600 FORMAT(1H0/,16X,11HCOORDINATES) 606 FORMAT(1H ,15X,11H-----------) 601 FORMAT(1H0/,16X,18HNORMALIZED WEIGHTS/ X 1H , 15X,18H------------------) 602 FORMAT(1H0/,16X,8HROTATION/1H ,15X,8H--------) 603 FORMAT(1H+,27X,11HIDEAL POINT/1H ,15X,23H-----------------------) 604 FORMAT(1H+,27X,16HANTI-IDEAL POINT/1H ,15X, X 28H----------------------------) C 605 FORMAT(1H+,27X,13HWITH WEIGHTS:) 607 FORMAT(1H+,40X,20F4.0) 608 FORMAT(1H ,15X,25H-------------------------) RETURN END SUBROUTINE RESUL2(DATA,DISP,DSTAR,H1,H2,NI,NJ,OPTION,IND,NCOL, X NC,IWRITE,IW,LU,ISHE,IPRI1,SLOPE,CEPT) C ****************************************************************** C * * C * R E S U L 2 * C * * C * PURPOSE: PLOTS DATA VERSUS DISP AND DSTAR (OPTIONAL) * C * STORES DSTAR AND DISP (OPTIONAL) * C * * C * SUBROUTINES CALLED: PLOT * C * * C ****************************************************************** DIMENSION DATA(NCOL),DISP(NCOL),DSTAR(NCOL),H1(NC),H2(NC),IND(NC) INTEGER OPTION,OPT C OPT=MOD(OPTION,10) IPRI=0 IF (IPRI1.EQ.1.OR.IPRI1.EQ.3) IPRI=1 IF (IPRI.EQ.1) WRITE (IWRITE,601) IF (IPRI.EQ.1) CALL PREREC(DISP,1,NCOL,1,IWRITE,NI) IF (IPRI.EQ.1) WRITE (IWRITE,603) IF (IPRI.EQ.1) CALL PREREC(DSTAR,1,NCOL,1,IWRITE,NI) IF (IPRI.EQ.1) WRITE (IWRITE,604) SLOPE,CEPT IF (ISHE.LT.NI) GO TO 30 L=0 DO 10 I=1,NCOL L=L+1 H1(L)=DATA(I) H2(L)=DISP(I) 10 CONTINUE DO 20 I=1,NCOL L=L+1 H1(L)=DATA(I) H2(L)=DSTAR(I) 20 CONTINUE WRITE(IWRITE,600) NI IF (OPT.EQ.1) WRITE(IWRITE,611) IF (OPT.EQ.2) WRITE(IWRITE,612) IF (OPT.EQ.3) WRITE(IWRITE,613) IF (OPT.EQ.4) WRITE(IWRITE,614) CALL PRPLOT(H1,H2,IND,NC,NCOL,-1,IWRITE,0) 30 IF (IW.EQ.0) RETURN DO 31 J=1,NCOL IF (IW.EQ.1) WRITE (LU,620) NI,NJ,J,DSTAR(J) IF (IW.EQ.2) WRITE (LU,620) NI,NJ,J,DSTAR(J),DATA(J),DISP(J) 31 CONTINUE C 600 FORMAT(1H1,3X,44HSCATTER PLOT. DATA (X-AXIS) VS CRITERION (*), X 38H AND PREDICTED VALUES (D) (M=*+D), ROW,I5) 611 FORMAT(1H+,90X,25H GENERAL UNFOLDING MODEL ) 612 FORMAT(1H+,90X,25H WEIGHTED UNFOLDING MODEL) 613 FORMAT(1H+,90X,25H UNFOLDING MODEL ) 614 FORMAT(1H+,90X,25H VECTOR MODEL ) 601 FORMAT(1H0/,16X,16HCRITERION VALUES, X 59H (EXTERNAL DATA, POSSIBLY TRANSFORMED BY MONOTONE FUNCTION)/ X 1H ,15X,16H----------------) 603 FORMAT(1H0/,16X,16HPREDICTED VALUES/1H ,15X,16H----------------) 604 FORMAT(1H0/,16X,5HSLOPE,F11.5,3X,9HINTERCEPT,F11.5) 620 FORMAT(3I5,3F13.6) RETURN END SUBROUTINE SCAR(DIST,LBK,N) C ****************************************************************** C * * C * S C A R * C * * C * PURPOSE: POOL DISTANCES WITHIN TIEBLOCKS * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION DIST(N),LBK(N) C IDLBK=0 JB=0 3 IDLBK=IDLBK+1 NLBK=LBK(IDLBK) MAX=JB+NLBK MIN=JB+1 DMEAN=0.0 DO 1 K=MIN,MAX DMEAN=DMEAN+DIST(K) 1 CONTINUE DMEAN=DMEAN/FLOAT(NLBK) DO 2 K=MIN,MAX DIST(K)=DMEAN 2 CONTINUE JB=MAX IF (JB.LT.N) GO TO 3 RETURN END SUBROUTINE SHEL9(DIS,IND,N) C ****************************************************************** C DIMENSION DIS(N) ,IND(N) C M=N 4 M=M/2 IF (M .LT. 1) RETURN K=N-M J=1 3 I=J 2 L=I+M IF ( DIS(L) .GE. DIS(I)) GOTO 1 B = DIS(I) DIS(I) = DIS(L) DIS(L) = B KK = IND(I) IND(I) = IND(L) IND(L) = KK I=I-M IF (I .GE. 1) GOTO 2 1 J=J+1 IF (J-K) 3,3,4 END SUBROUTINE SIVAD(M,N,Q,U,IU,V,IV,E,ICV) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C**** C C**** NAME: SIVAD(M,N,Q,U,IU,V,IV,E,ICV) C C**** C C**** C C**** PURPOSE : COMPUTES AND ORDERS THE SINGULAR VALUES OF A REAL C C**** MATRIX A AND GIVES A COMPLETE ORTHOGONAL DECOMPO- C C**** SITION: C C**** A = U * Q * V(T) C C**** C C**** C C**** PARAMETERS : M NUMBER OF ROWS OF THE MATRIX A C C**** C C**** N NUMBER OF COLUMNS OF THE MATRIX A C C**** C C**** IU LENGTH OF THE INPUT- AND OUTPUTVECTOR C C**** U,SO IU EQUALS M*N C C**** C C**** IV LENGTH OF THE OUTPUTVECTOR V,SO IV C C**** EQUALS N*N C C**** C C**** ICV OUTPUT PARAMETER,CONTAINS THE MAXIMUM C C**** NUMBER OF ITERATIONS USED,IF ICV = 0 C C**** THERE WAS NO CONVERGENCE FOR ONE OF THE C C**** SINGULAR VALUES AND THE SUBROUTINE C C**** HAS STOPPED C C**** C C**** C C**** DIMENSIONS: U(IU) CONTAINS ON INPUT THE M*N MATRIX TO BE C C**** COMPOSED,ON OUTPUT IT CONTAINS THE M*N C C**** ORTHONORMALIZED MATRIX U C C**** C C**** V(IV) REPRESENTS THE ORTHOGONAL MATRIX V C C**** C C**** Q(N) A VECTOR HOLDING THE NONNEGATIVE C C**** SINGULAR VALUES OF A,IN DECREASING C C**** SEQUENCE C C**** C C**** E(N) A VECTOR USED IN THE SUBROUTINE,IT C C**** SHOULD BE DECLARED IN THE MAIN PROGRAM C C**** C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C**** C**** * * * * * * START OF PROGRAM * * * * * * * * * * * * C**** C**** DIMENSION U(IU),V(IV),Q(N),E(N) DATA ZERO,ONE,TWO/0.0,1.0,2.0/ EPS=2.0**(-24) TOL=10.0**(-38) C**** C**** C**** * * * * * * HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM C**** C**** G=ZERO X=ZERO DO 130 I=1,N E(I)=G S=ZERO L=I+1 DO 20 J=I,M II=J+(I-1)*M S=S+U(II)*U(II) 20 CONTINUE IF(S.GE.TOL) GOTO 30 G=ZERO GOTO 60 30 CONTINUE JJ=I+(I-1)*M F=U(JJ) G=-SQRT(S) IF(F.LT.ZERO) G=-G H=F*G-S U(JJ)=F-G IF(L.GT.N) GOTO 60 DO 50 J=L,N S=ZERO DO 40 K=I,M S=S+U(K+(I-1)*M)*U(K+(J-1)*M) 40 CONTINUE F=S/H DO 50 K=I,M KJ=K+(J-1)*M U(KJ)=U(KJ)+F*U(K+(I-1)*M) 50 CONTINUE 60 CONTINUE Q(I)=G S=ZERO IF(L.GT.N) GOTO 75 DO 70 J=L,N IJ=I+(J-1)*M S=S+U(IJ)*U(IJ) 70 CONTINUE 75 CONTINUE IF(S.GE.TOL) GOTO 80 G=ZERO GOTO 120 80 CONTINUE I1=(M+1)*I F=U(I1) G=-SQRT(S) IF(F.LT.ZERO) G=-G H=F*G-S U(I1)=F-G DO 90 J=L,N E(J)=U(I+(J-1)*M)/H 90 CONTINUE DO 110 J=L,M S=ZERO DO 100 K=L,N S=S+U(J+(K-1)*M)*U(I+(K-1)*M) 100 CONTINUE DO 110 K=L,N JK=J+(K-1)*M U(JK)=U(JK)+S*E(K) 110 CONTINUE 120 CONTINUE Y=ABS(Q(I))+ABS(E(I)) IF(Y.GT.X) X=Y 130 CONTINUE C**** C**** C**** * * * * * * ACCUMULATION OF RIGHT HAND TRANSFORMATIONS C**** C**** DO 210 II=1,N I=(N+1)-II IF(G.NE.ZERO) GOTO 150 GOTO 190 150 CONTINUE IPLUS=(M+1)*I H=U(IPLUS)*G DO 160 J=L,N V(J+(I-1)*N)=U(I+(J-1)*M)/H 160 CONTINUE DO 180 J=L,N S=ZERO DO 170 K=L,N S=S+U(I+(K-1)*M)*V(K+(J-1)*N) 170 CONTINUE DO 180 K=L,N KJ=K+(J-1)*N V(KJ)=V(KJ)+S*V(K+(I-1)*N) 180 CONTINUE 190 CONTINUE IF(L.GT.N) GOTO 205 DO 200 J=L,N V(J+(I-1)*N)=ZERO V(I+(J-1)*N)=ZERO 200 CONTINUE 205 CONTINUE V(I+(I-1)*N)=ONE G=E(I) L=I 210 CONTINUE C**** C**** C**** * * * * * * ACCUMULATION OF LEFT HAND TRANSFORMATIONS C**** C**** DO 320 II=1,N I=(N+1)-II JJ=I+(I-1)*M L=I+1 G=Q(I) IF(L.GT.N) GOTO 245 DO 240 J=L,N U(I+(J-1)*M)=ZERO 240 CONTINUE 245 CONTINUE IF(G.NE.ZERO) GOTO 250 GOTO 290 250 CONTINUE H=U(JJ)*G IF(L.GT.N) GOTO 275 DO 270 J=L,N S=ZERO DO 260 K=L,M S=S+U(K+(I-1)*M)*U(K+(J-1)*M) 260 CONTINUE F=S/H DO 270 K=I,M KJ=K+(J-1)*M U(KJ)=U(KJ)+F*U(K+(I-1)*M) 270 CONTINUE 275 CONTINUE DO 280 J=I,M JI=J+(I-1)*M U(JI)=U(JI)/G 280 CONTINUE GOTO 310 290 CONTINUE DO 300 J=I,M U(J+(I-1)*M)=ZERO 300 CONTINUE 310 CONTINUE U(JJ)=U(JJ)+ONE 320 CONTINUE C**** C**** C**** * * * * * * DIAGONALIZATION OF THE BIDIAGONAL FORM * * C**** C**** EPS=EPS*X ICV=0 ICONV=1 DO 470 KK=1,N K=(N+1)-KK IF(ICV.LT.ICONV) ICV=ICONV ICONV=0 340 CONTINUE ICONV=ICONV+1 IF(ICONV.GT.50) GOTO 550 DO 350 LL=1,K L=(K+1)-LL IF(ABS(E(L)).LE.EPS) GOTO 390 LMIN=L-1 IF(ABS(Q(LMIN)).LE.EPS) GOTO 360 350 CONTINUE 360 CONTINUE C**** C**** C**** * * * * * * CANCELLATION OF E(L) IF L>1 * * * * * * * C**** C**** C=ZERO S=ONE L1=L-1 DO 380 I=L,K F=S*E(I) E(I)=C*E(I) IF(ABS(F).LE.EPS) GOTO 390 G=Q(I) IF(ABS(F).LT.ABS(G)) GOTO 362 IF (F.EQ.ZERO) GOTO 365 H=G/F H=ABS(F)*SQRT(ONE+H*H) GOTO 368 362 CONTINUE H=F/G H=ABS(G)*SQRT(ONE+H*H) GOTO 368 365 CONTINUE H=ZERO GOTO 369 368 C=G/H S=-F/H 369 Q(I)=H DO 370 J=1,M JL=J+(L1-1)*M JI=J+(I-1)*M Y=U(JL) Z=U(JI) U(JL)=Y*C+Z*S U(JI)=-Y*S+Z*C 370 CONTINUE 380 CONTINUE 390 CONTINUE C**** C**** C**** * * * * * * TEST CONVERGENCE * * * * * * * * * * * * C**** C**** Z=Q(K) IF(L.EQ.K) GOTO 450 C**** C**** C**** * * * * * * SHIFT FROM BOTTOM 2 * 2 MINOR * * * * * * C**** C**** X=Q(L) KMIN=K-1 Y=Q(KMIN) G=E(KMIN) H=E(K) F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(TWO*H*Y) G=SQRT(F*F+ONE) R=F+G IF(F.LT.ZERO) R=F-G F=((X-Z)*(X+Z)+H*(Y/R-H))/X C**** C**** C**** * * * * * * NEXT QR TRANSFORMATION * * * * * * * * * C**** C**** C=ONE S=ONE L2=L+1 DO 430 I=L2,K G=E(I) Y=Q(I) H=S*G G=C*G IMIN=I-1 IF(ABS(F).LT.ABS(H)) GOTO 392 IF(F.EQ.ZERO) GOTO 395 Z=H/F Z=ABS(F)*SQRT(ONE+Z*Z) GOTO 398 392 CONTINUE Z=F/H Z=ABS(H)*SQRT(ONE+Z*Z) GOTO 398 395 CONTINUE Z=ZERO GOTO 399 398 C=F/Z S=H/Z 399 E(IMIN)=Z F=X*C+G*S G=-X*S+G*C H=Y*S Y=Y*C DO 400 J=1,N JM=J+(IMIN-1)*N JI=J+(I-1)*N X=V(JM) Z=V(JI) V(JM)=X*C+Z*S V(JI)=-X*S+Z*C 400 CONTINUE IF(ABS(F).LT.ABS(H)) GOTO 412 IF(F.EQ.ZERO) GOTO 415 Z=H/F Z=ABS(F)*SQRT(ONE+Z*Z) GOTO 418 412 CONTINUE Z=F/H Z=ABS(H)*SQRT(ONE+Z*Z) GOTO 418 415 CONTINUE Z=ZERO GOTO 419 418 C=F/Z S=H/Z 419 Q(IMIN)=Z F=C*G+S*Y X=-S*G+C*Y DO 420 J=1,M JM=J+(IMIN-1)*M JI=J+(I-1)*M Y=U(JM) Z=U(JI) U(JM)=Y*C+Z*S U(JI)=-Y*S+Z*C 420 CONTINUE 430 CONTINUE E(L)=ZERO E(K)=F Q(K)=X GOTO 340 C**** C**** C**** * * * * * * CONVERGENCE * * * * * * * * * * * * * * * C**** C**** 450 CONTINUE IF(Z.GE.ZERO) GOTO 470 C**** C**** C**** * * * * * * Q(K) IS MADE NON-NEGATIVE * * * * * * * * C**** C**** Q(K)=-Z DO 460 J=1,N JK=J+(K-1)*N V(JK)=-V(JK) 460 CONTINUE 470 CONTINUE C**** C**** C**** * * * * * * ORDERING OF SINGULAR VALUES * * * * * * * C**** C**** NMIN=N-1 DO 540 I=1,NMIN IPLUS=I+1 480 CONTINUE DO 530 J=IPLUS,N IF(Q(I).GE.Q(J)) GOTO 530 R=Q(I) Q(I)=Q(J) Q(J)=R DO 490 K=1,N KI=K+(I-1)*N KJ=K+(J-1)*N R=V(KI) V(KI)=V(KJ) V(KJ)=R 490 CONTINUE DO 510 K=1,M KI=K+(I-1)*M KJ=K+(J-1)*M R=U(KI) U(KI)=U(KJ) U(KJ)=R 510 CONTINUE JPLUS=J+1 IF(JPLUS.GT.N) GOTO 540 GOTO 480 530 CONTINUE 540 CONTINUE GOTO 570 550 CONTINUE ICV=0 570 CONTINUE RETURN END SUBROUTINE STATIS(OPTION,NROW,NCOL,NDIM,FITT,NANA,K,IWRITE) C ****************************************************************** C * * C * S T A T I S * C * * C * PURPOSE: COMPUTE STATISTICS * C * * C * SUBROUTINES CALLED: SHEL9 * C * * C ****************************************************************** DIMENSION OPTION(NANA),IORD(4),HELP(4),FITT(NANA) INTEGER OPTION WRITE(IWRITE,640) K DO 5 J=1,NANA C C WARNING C IF (OPTION(J).GT.4) WRITE (IWRITE,660) IF (OPTION(J).GT.4) GO TO 6 5 CONTINUE 6 DO 10 J=1,NANA FITT(J)=FITT(J)**2 HELP(J)=FLOAT(MOD(OPTION(J),10)) IORD(J)=J 10 CONTINUE CALL SHEL9(HELP,IORD,NANA) DO 15 J=1,NANA R=HELP(J) IF (R.EQ.0.) GO TO 15 IF (FITT(IORD(J)).GT. .9995) WRITE(IWRITE,650) IORD(J) IF (FITT(IORD(J)).GT. .9995) GO TO 15 F1=FITT(IORD(J))/(1.-FITT(IORD(J))) ID=NCOL-1 IF (J.EQ.NANA) RR=5 IF (J.EQ.NANA) GO TO 16 RR=HELP(J+1) IF (RR.EQ.R) GO TO 15 F2=(FITT(IORD(J))-FITT(IORD(J+1)))/(1.-FITT(IORD(J))) 16 IF (R.EQ.4.) GO TO 40 IF (R.EQ.3.) GO TO 30 IF (R.EQ.2.) GO TO 20 IDF1=IFIX((FLOAT(NDIM)/2)*(NDIM+3)) IF (RR.EQ.2.) IDF12=.5*NDIM*(NDIM-1) IF (RR.EQ.3.) IDF12=.5*NDIM*(NDIM+1)-1 IF (RR.EQ.4.) IDF12=.5*NDIM*(NDIM+1) GO TO 50 20 IDF1=2*NDIM IF (RR.EQ.3.) IDF12=NDIM-1 IF (RR.EQ.4.) IDF12=NDIM GO TO 50 30 IDF1=NDIM+1 IF (RR.EQ.4.) IDF12=1 GO TO 50 40 IDF1=NDIM IDF2=ID-IDF1 GO TO 60 50 IDF2=ID-IDF1 IF (RR.EQ.0.) GO TO 60 IF (RR.EQ.5.) GO TO 60 FF2=(FLOAT(IDF2)/FLOAT(IDF12))*F2 60 FF1=(FLOAT(IDF2)/FLOAT(IDF1))*F1 FIT=FITT(IORD(J)) IF (R.EQ.1. .AND.RR.NE.1.) WRITE (IWRITE,600) FIT,FF1,IDF1,IDF2 IF (R.EQ.1. .AND.RR.EQ.2.) WRITE (IWRITE,601) FF2,IDF12,IDF2 IF (R.EQ.1. .AND.RR.EQ.3.) WRITE (IWRITE,602) FF2,IDF12,IDF2 IF (R.EQ.1. .AND.RR.EQ.4.) WRITE (IWRITE,603) FF2,IDF12,IDF2 IF (R.EQ.2. .AND.RR.NE.2.) WRITE (IWRITE,604) FIT,FF1,IDF1,IDF2 IF (R.EQ.2. .AND.RR.EQ.3.) WRITE (IWRITE,602) FF2,IDF12,IDF2 IF (R.EQ.2. .AND.RR.EQ.4.) WRITE (IWRITE,603) FF2,IDF12,IDF2 IF (R.EQ.3. .AND.RR.NE.3.) WRITE (IWRITE,605) FIT,FF1,IDF1,IDF2 IF (R.EQ.3. .AND.RR.EQ.4.) WRITE (IWRITE,603) FF2,IDF12,IDF2 IF (R.EQ.4. .AND.RR.NE.4.) WRITE (IWRITE,606) FIT,FF1,IDF1,IDF2 15 CONTINUE WRITE(IWRITE,607) C 600 FORMAT(1H0,12X,34HGENERAL UNFOLDING MODEL, P.VAF. =,F6.3,5H, F =, X F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM) 601 FORMAT(1H ,15X,42HVERSUS WEIGHTED UNFOLDING MODEL F =, X F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM) 602 FORMAT(1H ,15X,42HVERSUS UNFOLDING MODEL F =, X F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM) 603 FORMAT(1H ,15X,42HVERSUS VECTOR MODEL F =, X F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM) 604 FORMAT(1H0,12X,34HWEIGHTED UNFOLDING MODEL, P.VAF. =,F6.3,5H, F =, X F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM) 605 FORMAT(1H0,12X,34HUNFOLDING MODEL, P.VAF. =,F6.3,5H, F =, X F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM) 606 FORMAT(1H0,12X,34HVECTOR MODEL, P.VAF. =,F6.3,5H, F =, X F8.3,5H WITH,I4,4H AND, I4,19H DEGREES OF FREEDOM) 607 FORMAT(1H0//) 640 FORMAT(1H0//,13X,33HSTATISTICS ACROSS OPTIONS FOR ROW,I5/ X 1H ,12X,38H--------------------------------------) 650 FORMAT(1H0/,13X,22HPERFECT FIT FOR OPTION,I2,16H. NO STATISTICS.) 660 FORMAT(1H ,12X,38HTHE STATISTICS GIVEN BELOW DEPEND ONLY, X 46H ON THE METRIC FIT, NOT ON THE NON-METRIC FIT.) RETURN END SUBROUTINE TEST(P1,N1,N2) C ****************************************************************** C * * C * T E S T * C * * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION P1(N1,N2) C DO 10 I=1,N1 WRITE(6,600) (P1(I,J),J=1,N2) 10 CONTINUE WRITE(6,601) N1,N2 600 FORMAT(10F11.4) 601 FORMAT(10X,2I10) RETURN END SUBROUTINE TIEBL (DISM, LBK, N, IDLBK ) C ****************************************************************** C * * C * T I E B L * C * * C * PURPOSE:- FIND A VECTOR LBK SO THAT EACH ELEMENT OF IT * C * REPRESENTS THE LENGTH OF A TIEBLOCK FOUND IN DISM. * C * - IDLBK GIVES THE NUMBER OF TIEBLOCKS * C * * C * SUBROUTINES CALLED: NONE. * C * * C * SEE ACCOMPANYING PAPER SECTIONS: 2.2. * C * * C * AUTHOR: ERNST VAN WANING. RELEASED: JULY 1976 * C * * C ****************************************************************** C C DIMENSION DISM(N),LBK(N) C IDLBK = 0 JBASE = 1 3 K = 1 2 IDDISM = JBASE + K IF (IDDISM .GT. N) GOTO 1 IF (DISM(JBASE) .NE. DISM(IDDISM)) GOTO 1 K = K + 1 GOTO 2 1 IDLBK = IDLBK + 1 LBK(IDLBK) = K JBASE = IDDISM IF (JBASE .LE. N) GOTO 3 RETURN END SUBROUTINE TRED2(NM,N,D,E,Z) C ****************************************************************** C * * C * T R E D 2 * C * * C * PURPOSE: REDUCES A REAL SYMMETRIC MATRIX TO A SYMMETRIC TRI- * C * DIAGONAL MATRIX USING ORTHOGONAL SIMILARITY TRANS- * C * FORMATIONS; THE ACCUMULATED ORTHOGONAL TRANSFORMA- * C * TIONS ARE RETAINED IN Z. * C * * C * SUBROUTINES CALLED: NONE * C * * C * ADAPTED FROM EISPACK * C * * C ****************************************************************** DIMENSION Z(NM,N),D(N),E(N) C IF (N .EQ. 1) GO TO 320 C ********** FOR I=N STEP -1 UNTIL 2 DO -- ********** DO 300 II =2, N I = N + 2 - II L = I - 1 H = 0.0 SCALE = 0.0 IF (L .LT.2) GO TO 130 C ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) ********** DO 120 K =1, L 120 SCALE = SCALE + ABS(Z(I,K)) C IF (SCALE .NE. 0.0) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L Z(I,K) = Z(I,K) / SCALE H = H + Z(I,K) * Z(I,K) 150 CONTINUE C F = Z(I,L) G = -SIGN(SQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0 C DO 240 J = 1, L Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0 C ********** FORM ELEMENT OF A*U ********** DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C ********** FORM ELEMENT OF P ********** 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C ********** FORM REDUCED A ********** DO 260 J =1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) -G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0 E(1) = 0.0 C ********** ACCUMULATION OF TRANSFORMATION MATRICES ********** DO 500 I = 1, N L = I - 1 IF (D(I) .EQ. 0.0) GO TO 380 C C DO 360 J = 1, L G = 0.0 DO 340 K =1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0 IF (L .LT. 1) GO TO 500 C DO 400 J =1, L Z(I,J) = 0.0 Z(J,I) = 0.0 400 CONTINUE C 500 CONTINUE C RETURN END SUBROUTINE UNRAV(COEF,NPRM,COOR,WEIGH,ROTAT,OPTION,NDIM,ND,JJ,NR, X IWRITE,WSIGN,IPRE,WORK1,WORK3,WORK2,WORK4,CMAX,IPRI1,SLOPE,CEPT, X FIT) C ****************************************************************** C * * C * U N R A V E L * C * * C * PURPOSE: UNRAVELS THE REGRESSION COEFFICIENTS AND COMPUTES * C * THE COORDINATES FOR THE IDEAL POINTS AND/OR VECTORS * C * * C * SUBROUTINES CALLED: UNRA1,UNRA2,UNRA3,UNRA4,RESUL1 * C * * C ****************************************************************** DIMENSION COEF(NPRM),COOR(ND),WEIGH(ND),ROTAT(NDIM,ND), X WSIGN(NDIM),WORK1(NDIM),WORK4(NDIM),WORK2(NDIM,NDIM), X WORK3(NDIM,NDIM) INTEGER OPTION,OPT C OPT=MOD(OPTION,10) IF (OPT.EQ.1) CALL UNRA1(COEF,NPRM,COOR,WEIGH,ROTAT,NDIM,ND,JJ, X WORK1,WORK3,WORK2,WORK4,SLOPE,CEPT) IF (OPT.EQ.2) CALL UNRA2(COEF,NPRM,COOR,WEIGH,NDIM,ND,JJ, X SLOPE,CEPT) IF (OPT.EQ.3) CALL UNRA3(COEF,NPRM,COOR,WEIGH,WSIGN,NDIM,ND,JJ, X SLOPE,CEPT) IF (OPT.EQ.4) CALL UNRA4(COEF,NPRM,COOR,WEIGH,NDIM,ND,JJ,CMAX, X SLOPE,CEPT,FIT) IF (OPT.NE.0.AND.IPRI1.EQ.1) X CALL RESUL1(COOR,WEIGH,ROTAT,NDIM,ND,JJ,IWRITE,OPT,NR,IPRE) C RETURN END SUBROUTINE UNRA1(COEF,NPR,COOR,WEIGH,ROTAT,NDIM,ND,NOP, X WORK1,LSV,RSV,SIVA,STAND,CEPT) C ****************************************************************** C * * C * U N R A 1 * C * * C * PURPOSE: UNRAVELS FOR THE GENERAL UNFOLDING MODEL * C * * C * SUBROUTINES CALLED: TRED2,IMTQL2 * C * * C ****************************************************************** DIMENSION COEF(NPR),COOR(ND),LSV(NDIM,NDIM),RSV(NDIM,NDIM), X SIVA(NDIM),WORK1(NDIM),WEIGH(ND),ROTAT(NDIM,ND) REAL LSV C C NOPT=NOP-1 C C NOPT WILL INDICATE ON WHAT POSITIONS ARRAY COOR HAS TO BE FILLED C C N=NDIM+2 K=2*NDIM+2 L=2 M=NDIM-1 DO 1 I=1,M DO 3 J=L,NDIM RSV(I,J)=.5*COEF(K) RSV(J,I)=RSV(I,J) K=K+1 3 CONTINUE L=L+1 1 CONTINUE DO 5 I=1,NDIM RSV(I,I)=COEF(N) N=N+1 5 CONTINUE C C C COMPUTE THE EIGENVALUE DECOMPOSITION OF THE TRANSFORMATION MATRIX C RSV. C ON THE BASIS THEREOF COMPUTE (IN LSV) THE INVERSE (FOR THE IDEAL C POINTS); THE EIGENVALUES (IN SIVA) WILL GIVE THE WEIGHTS; C THE EIGENVECTORS (IN RSV) WILL GIVE THE ROTATION. C CALL TRED2(NDIM,NDIM,SIVA,WORK1,RSV) CALL IMTQL2(NDIM,NDIM,SIVA,WORK1,RSV,IERR) C C COLLECT WEIGHTS FROM ARRAY SIVA AND ROTATION FROM ARRAY RSV C STANDARDIZE WEIGHTS AND COMPUTE SLOPE C STAND=0.0 DO 64 I=1,NDIM STAND=STAND+SIVA(I)*SIVA(I) 64 CONTINUE STAND=SQRT(STAND/NDIM) DO 65 I=1,NDIM IL=NOPT*NDIM+I WEIGH(IL)=SIVA(I)/STAND 65 CONTINUE DO 70 I=1,NDIM DO 80 J=1,NDIM IL=NOPT*NDIM+J ROTAT(I,IL)=RSV(I,J) 80 CONTINUE 70 CONTINUE C DO 6 I=1,NDIM SMALL=.0000001 IF (ABS(SIVA(I)).GT.SMALL) SIVA(I)=1.0/SIVA(I) IF (ABS(SIVA(I)).LE.SMALL) SIVA(I)=0.0 6 CONTINUE C DO 10 K=1,NDIM DO 20 I=1,NDIM SUM=0.0 DO 30 J=1,NDIM SUM=SUM+RSV(I,J)*RSV(K,J)*SIVA(J) 30 CONTINUE WORK1(I)=SUM 20 CONTINUE DO 40 I=1,NDIM LSV(K,I)=WORK1(I) 40 CONTINUE 10 CONTINUE C C C COMPUTE THE IDEAL POINTS AND INTERCEPT C CEPT=0.0 DO 50 I=1,NDIM IL=NOPT*NDIM+I SUM=0.0 DO 60 J=1,NDIM L=J+1 SUM=SUM+COEF(L)*LSV(J,I) 60 CONTINUE L=I+1 COOR(IL)=-.5*SUM ASIGN=SIGN(1.,COOR(IL)) IF (ABS(COOR(IL)).GT.9999.) COOR(IL)=9999.*ASIGN CEPT=CEPT+COOR(IL)*COEF(L) 50 CONTINUE CEPT=COEF(1)+.5*CEPT C RETURN END SUBROUTINE UNRA2(COEF,NPR,COOR,WEIGH,NDIM,ND,NOP,STAND,CEPT) C ****************************************************************** C * * C * U N R A 2 * C * * C * PURPOSE: UNRAVELS FOR THE WEIGHTED UNFOLDING MODEL * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION COEF(NPR),COOR(ND),WEIGH(ND) C NOPT=NOP-1 C C FOR COMMENT SEE UNRA1 C C COMPUTE SLOPE AND STANDARDIZE WEIGHTS C STAND=0.0 DO 5 I=1,NDIM K=I+NDIM+1 STAND=STAND+COEF(K)*COEF(K) 5 CONTINUE STAND=SQRT(STAND/NDIM) CEPT=0.0 DO 10 I=1,NDIM IL=NOPT*NDIM+I L=I+1 K=L+NDIM WEIGH(IL)=COEF(K)/STAND COOR(IL)=-.5*COEF(L)/COEF(K) ASIGN=SIGN(1.,COOR(IL)) IF (ABS(COOR(IL)).GT.9999.) COOR(IL)=9999.*ASIGN CEPT=CEPT+COOR(IL)*COEF(L) 10 CONTINUE CEPT=COEF(1)+.5*CEPT RETURN END SUBROUTINE UNRA3(COEF,NPR,COOR,WEIGH,WSIGN,NDIM,ND,NOP,STAND,CEPT) C ****************************************************************** C * * C * U N R A 3 * C * * C * PURPOSE: UNRAVELS FOR THE UNFOLDING MODEL * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION COEF(NPR),COOR(ND),WEIGH(ND),WSIGN(NDIM) C C NOPT=NOP-1 C C SEE UNRA1 FOR COMMENT C K=NDIM+2 STAND=ABS(COEF(K)) CEPT=0.0 ASIGN=1. DO 10 I=1,NDIM L=I+1 IL=NOPT*NDIM+I COOR(IL)=-.5*COEF(L)/COEF(K)*WSIGN(I) CEPT=CEPT+COOR(IL)*COEF(L) IF (COEF(K).LT.0.) ASIGN=-1. WEIGH(IL)=WSIGN(I)*ASIGN BSIGN=SIGN(1.,COOR(IL)) IF (ABS(COOR(IL)).GT.9999.) COOR(IL)=9999.*BSIGN 10 CONTINUE CEPT=COEF(1)+.5*CEPT RETURN END SUBROUTINE UNRA4(COEF,NPR,COOR,WEIGH,NDIM,ND,NOP,CMAX,STAND,CEPT, X FIT) C ****************************************************************** C * * C * U N R A 4 * C * * C * PURPOSE: UNRAVELS FOR THE VECTOR MODEL * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION COEF(NPR),COOR(ND),WEIGH(ND) C NOPT=NOP-1 C C C SEE UNRA1 FOR COMMENT C SUM=0.0 K=NDIM+1 DO 10 J=2,K SUM=SUM+COEF(J)*COEF(J) 10 CONTINUE STAND=SQRT(SUM) DO 20 I=1,NDIM IL=NOPT*NDIM+I L=I+1 COOR(IL)=FIT*CMAX*COEF(L)/STAND WEIGH(IL)=1. 20 CONTINUE CEPT=COEF(1) STAND=STAND/(FIT*CMAX) RETURN END SUBROUTINE ZERO(A,L) C ****************************************************************** C * * C * Z E R O * C * * C * PURPOSE: FILLS ARRAY A WITH ZERO'S * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION A(L) C DO 10 J=1,L A(J)=0.0 10 CONTINUE RETURN END SUBROUTINE ZEROI(N,L) C ****************************************************************** C * * C * Z E R O I * C * * C * PURPOSE: FILLS ARRAY N WITH ZERO'S * C * * C * SUBROUTINES CALLED: NONE * C * * C ****************************************************************** DIMENSION N(L) C DO 10 J=1,L N(J)=0 10 CONTINUE RETURN END