!C IMSL ROUTINE NAME - ODFISH !C !C----------------------------------------------------------------------- !C !C COMPUTER - HP1000/SINGLE !C !C LATEST REVISION - NOVEMBER 1, 1979 !C !C PURPOSE - LINEAR DISCRIMINANT ANALYSIS METHOD OF FISHER !C FOR REDUCING THE NUMBER OF VARIABLES !C !C USAGE - CALL ODFISH (X,IX,NG,NV,ND,XM,IXM,NNV,E,C,IC, !C SW,SB,ISB,EX,CX,ICX,IS,IER) !C !C ARGUMENTS X - INPUT/OUTPUT MATRIX OF DIMENSION NV BY THE !C SUM OF THE ND(I). ON INPUT, X CONTAINS THE !C GROUP DATA. !C ON OUTPUT, EACH COLUMN DATA VECTOR IS REPLACED !C BY THE TRANSFORMED NNV BY 1 DATA VECTOR. !C (SEE REMARK 1). !C IX - INPUT ROW DIMENSION OF MATRIX X EXACTLY AS !C SPECIFIED IN THE DIMENSION STATEMENT IN THE !C CALLING PROGRAM. !C NG - INPUT NUMBER OF DISTINCT GROUPS. !C NV - INPUT NUMBER OF DATA VARIABLES. !C ND - INPUT VECTOR OF LENGTH NG. ND(I) CONTAINS !C THE NUMBER OF DATA VECTORS FROM THE I-TH !C GROUP IN THE DATA MATRIX X. !C XM - OUTPUT MATRIX OF DIMENSION NV BY NG+1. !C THE I-TH COLUMN OF XM CONTAINS THE !C MEAN VECTOR FOR THE I-TH GROUP. COLUMN !C NG+1 CONTAINS THE OVERALL MEAN VECTOR. !C IXM - INPUT ROW DIMENSION OF MATRIX XM EXACTLY AS !C SPECIFIED IN THE DIMENSION STATEMENT IN THE !C CALLING PROGRAM. !C NNV - OUTPUT NUMBER OF NEW VARIABLES COMPUTED. !C NNV IS GIVEN BY THE MINIMUM OF NV AND NG-1. !C E - OUTPUT VECTOR OF LENGTH NV CONTAINING THE !C EIGENVALUES IN THE FIRST NNV LOCATIONS IN !C DESCENDING ORDER. THE REMAINDER OF E IS !C USED AS WORK AREA. (SEE REMARK 2). !C C - OUTPUT MATRIX OF DIMENSION NV BY NNV. THE !C I-TH COLUMN OF C CONTAINS THE VECTOR FOR !C COMPUTING THE I-TH NEW LINEAR VARIABLE. !C (SEE REMARK 1). !C IC - INPUT ROW DIMENSION OF MATRIX C EXACTLY AS !C SPECIFIED IN THE DIMENSION STATEMENT IN THE !C CALLING PROGRAM. !C SW - WORK VECTOR OF LENGTH ((NV+1)*NV)/2 OR 4, !C WHICHEVER IS GREATER. (SEE REMARK 3). !C SB - WORK MATRIX OF DIMENSION NV BY NV. !C (SEE REMARK 3). !C ISB - INPUT ROW DIMENSION OF MATRIX SB EXACTLY AS !C SPECIFIED IN THE DIMENSION STATEMENT IN THE !C CALLING PROGRAM. !C EX - COMPLEX WORK VECTOR OF LENGTH NV CONTAINING !C THE UNSORTED EIGENVALUES. (SEE REMARK 4). !C NOTE - THE ROUTINE TREATS EX AS A REAL VECTOR !C OF LENGTH 2*NV. AN APPROPRIATE EQUIVALENCE !C STATEMENT MAY BE REQUIRED. SEE DOCUMENT !C EXAMPLE. !C CX - COMPLEX WORK MATRIX OF DIMENSION NV BY NV !C CONTAINING THE EIGENVECTORS ASSOCIATED WITH !C EX. (SEE REMARK 4). !C NOTE - THE ROUTINE TREATS CX AS A REAL VECTOR !C OF LENGTH 2*NV*NV. AN APPROPRIATE EQUI- !C VALENCE STATEMENT MAY BE REQUIRED. SEE !C DOCUMENT EXAMPLE. !C ICX - INPUT ROW DIMENSION OF MATRIX CX EXACTLY AS !C SPECIFIED IN THE DIMENSION STATEMENT IN THE !C CALLING PROGRAM. !C IS - INTEGER WORK VECTOR OF LENGTH NV. !C IER - OUTPUT ERROR PARAMETER !C TERMINAL ERROR !C IER = 129 INDICATES THAT THE COVARIANCE !C MATRIX IS SINGULAR. CHECK FOR !C REDUNDANCIES IN THE ORIGINAL VARIABLES. !C IER = 130 INDICATES THAT THE EIGENVALUES !C COULD NOT BE COMPUTED. REDUCE THE !C NUMBER OF VARIABLES AND TRY AGAIN. !C IER = 131 INDICATES ONE OF IX,IXM,IC,ISB, OR !C ICX IS TOO SMALL. REDIMENSION MATRICES. !C !C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 !C - SINGLE/H36,H48,H60 !C !C REQD. IMSL ROUTINES - SINGLE/EBALAF,EBBCKF,EHBCKF,EHESSF,EIGRF, !C EQRH3F,LEQT1P,LUDECP,LUELMP,UERTST,UGETIO, !C VSRTR !C - DOUBLE/EBALAF,EBBCKF,EHBCKF,EHESSF,EIGRF, !C EQRH3F,LEQT1P,LUDECP,LUELMP,UERTST,UGETIO, !C VSRTRD !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C REMARKS 1. ON INPUT, THE MATRIX X CONTAINS THE ND(1)+ND(2)+... !C ...+ND(NG) SAMPLE (COLUMN) VECTORS. THE FIRST ND(1) !C COLUMNS OF MATRIX X SHOULD CONTAIN THE DATA FOR GROUP !C ONE. THE NEXT ND(2) COLUMNS OF MATRIX X SHOULD CONTAIN !C THE DATA FOR GROUP TWO, AND SO ON. !C ON OUTPUT, EACH ORIGINAL OBSERVATION X(I) (WHERE X(I) !C REPRESENTS A COLUMN VECTOR OF INPUT MATRIX X) IS !C REPLACED BY ITS TRANSFORMED OBSERVATION ( A NNV BY 1 !C DATA VECTOR). THE TRANSFORMED OBSERVATION IS GIVEN BY !C THE FORMULA (C TRANSPOSE) * X(I) WHERE C IS THE OUTPUT !C MATRIX C. !C 2. THE COLUMN VECTORS IN THE TRANSFORMATION MATRIX C ARE !C IN ORDER OF DECREASING IMPORTANCE, THAT IS, THE EIGEN- !C VECTOR ASSOCIATED WITH THE LARGEST EIGENVALUE IS IN !C THE FIRST COLUMN OF C, AND SO ON. !C 3. IF THE WITHIN AND BETWEEN SUMS OF SQUARES MATRICES, SW !C AND SB, ARE DESIRED, THEY MAY BE PRINTED BY MODIFYING !C THE SOURCE CODE AT THE DESIGNATED COMMENT CARD. !C 4. ONE SHOULD PRUDENTLY INSPECT THE COMPLEX EIGENVALUE !C VECTOR AND EIGENVECTOR MATRIX AND CHECK THAT THERE ARE !C ONLY NG-1 OR NV NONZERO EIGENVALUES AND THAT THE IMAG- !C INARY PARTS OF THE EIGENVECTORS ARE TRULY ZERO. !C !C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- !C SUBROUTINE ODFISH (X,IX,NG,NV,ND,XM,IXM,NNV,E,C,IC,SW,SB,ISB,EX,CX,ICX,IS,IER) !C SPECIFICATIONS FOR ARGUMENTS INTEGER IX,NG,NV,ND(1),IXM,NNV,IC,ISB,ICX,IS(1),IER REAL X(IX,1),XM(IXM,1),E(1),C(IC,1),SW(1),SB(ISB,1) REAL EX(2,1),CX(2,ICX,1) !C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER I,IDGT,IJOB,J,K,KHI,KLO,L,N,NGP1 REAL D1,D2,XMI,XMJ DOUBLE PRECISION DSUM,DTOT,DZERO DATA DZERO/0.0D0/ !C ERROR CHECKS !C FIRST EXECUTABLE STATEMENT IER = 0 IF (IX.GE.NV.AND.IXM.GE.NV.AND.IC.GE.NV.AND.ISB.GE.NV.AND.ICX.GE.NV) GO TO 5 IER = 131 GO TO 9000 !C COMPUTE MEAN VECTOR MATRIX 5 NGP1 = NG+1 DO 20 I=1,NV DTOT = DZERO KHI = 0 DO 15 J=1,NG KLO = KHI+1 KHI = KLO+ND(J)-1 DSUM = DZERO DO 10 K=KLO,KHI DSUM = DSUM+X(I,K) 10 CONTINUE XM(I,J) = DSUM/ND(J) DTOT = DTOT+DSUM 15 CONTINUE XM(I,NGP1) = DTOT/KHI 20 CONTINUE N = KHI !C COMPUTE POOLED SUMS OF SQUARES MATRIX IR = 0 DO 40 I=1,NV DO 35 J=1,I DSUM = DZERO KHI = 0 DO 30 L=1,NG KLO = KHI+1 KHI = KLO+ND(L)-1 XMI = XM(I,L) XMJ = XM(J,L) DO 25 K=KLO,KHI DSUM = DSUM+(X(I,K)-XMI)*(X(J,K)-XMJ) 25 CONTINUE 30 CONTINUE IR = IR+1 SW(IR) = DSUM 35 CONTINUE 40 CONTINUE !C PRINT WITHIN SUMS OF SQUARES MATRIX HERE !C COMPUTE BETWEEN SUMS OF SQUARES MATRIX DO 55 I=1,NV DO 50 J=1,I XMI = XM(I,NGP1) XMJ = XM(J,NGP1) DSUM = DZERO DO 45 K=1,NG DSUM = DSUM+ND(K)*(XM(I,K)-XMI)*(XM(J,K)-XMJ) 45 CONTINUE SB(I,J) = DSUM SB(J,I) = DSUM 50 CONTINUE 55 CONTINUE !C PRINT BETWEEN SUMS OF SQUARES MATRIX HERE !C FORM MATRIX TO GET EIGENVALUES AND VECTORS IDGT = 0 CALL LEQT1P (SW,NV,NV,SB,ISB,IDGT,D1,D2,IER) IF (IER.NE.0) GO TO 9000 !C COMPUTE EIGEN INFORMATION OF SB IJOB = 1 CALL EIGRF (SB,NV,ISB,IJOB,EX,CX,ICX,SW,IER) IF (IER.EQ.0) GO TO 60 IER = 130 GO TO 9000 !C EXTRACT REAL EIGENSTRUCTURE 60 DO 65 I=1,NV E(I) = -EX(1,I) IS(I) = I 65 CONTINUE !C SORTED INTO DESCENDING ORDER (BY NEGATIVE) CALL VSRTR (E,NV,IS) !C THE NUMBER OF RETAINED EIGENVALUES = NNV NNV = MIN0(NV,NG-1) DO 75 J=1,NNV E(J) = -E(J) K = IS(J) DO 70 I=1,NV C(I,J) = CX(1,I,K) 70 CONTINUE 75 CONTINUE !C COMPUTE FISHER PROJECTION VARIABLES DO 95 J=1,N DO 85 K=1,NNV DSUM = DZERO DO 80 I=1,NV DSUM = DSUM+X(I,J)*C(I,K) 80 CONTINUE SW(K) = DSUM 85 CONTINUE DO 90 K=1,NNV X(K,J) = SW(K) 90 CONTINUE 95 CONTINUE GO TO 9005 9000 CONTINUE CALL UERTST (IER,'ODFISH') 9005 RETURN END !C IMSL ROUTINE NAME - EBALAF !C !C----------------------------------------------------------------------- !C !C COMPUTER - HP1000/SINGLE !C !C LATEST REVISION - JUNE 1, 1982 !C !C PURPOSE - NUCLEUS CALLED ONLY BY IMSL ROUTINE EIGRF !C !C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 !C - SINGLE/H36,H48,H60 !C !C REQD. IMSL ROUTINES - NONE REQUIRED !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C COPYRIGHT - 1982 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- !C SUBROUTINE EBALAF (A,N,IA,D,K,L) !C SPECIFICATIONS FOR ARGUMENTS INTEGER N,IA,K,L REAL A(IA,1),D(1) !C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER L1,K1,K1P1,K11,JJ,J,I,LL,NOCONV REAL R,C,F,G,B,S,B2,ONE,ZERO,P95 DATA B/16.0/,B2/256.0/ DATA ZERO/0.0/,ONE/1.0/,P95/.95/ !C REDUCE NORM A BY DIAGONAL SIMILARITY !C TRANSFORMATION STORED IN D !C FIRST EXECUTABLE STATEMENT L1 = 1 K1 = N !C SEARCH FOR ROWS ISOLATING AN EIGEN- !C VALUE AND PUSH THEM DOWN 5 K1P1 = K1+1 IF (K1.LT.1) GO TO 35 K11=K1 DO 30 JJ=1,K11 J = K1P1-JJ R = ZERO DO 10 I=1,K1 IF (I.EQ.J) GO TO 10 R=R+ABS(A(J,I)) 10 CONTINUE IF (R.NE.ZERO) GO TO 30 D(K1) = J IF (J.EQ.K1) GO TO 25 DO 15 I=1,K1 F = A(I,J) A(I,J) = A(I,K1) A(I,K1) = F 15 CONTINUE DO 20 I=L1,N F = A(J,I) A(J,I) = A(K1,I) A(K1,I) = F 20 CONTINUE 25 K1 = K1-1 GO TO 5 30 CONTINUE !C SEARCH FOR COLUMNS ISOLATING AN !C EIGENVALUE AND PUSH THEM LEFT 35 IF (K1.LT.L1) GO TO 65 LL = L1 DO 60 J=LL,K1 C = ZERO DO 40 I=L1,K1 IF (I.EQ.J) GO TO 40 C = C+ABS(A(I,J)) 40 CONTINUE IF (C.NE.ZERO) GO TO 60 D(L1) = J IF (J.EQ.L1) GO TO 55 DO 45 I=1,K1 F = A(I,J) A(I,J) = A(I,L1) A(I,L1) = F 45 CONTINUE DO 50 I=L1,N F = A(J,I) A(J,I) = A(L1,I) A(L1,I) = F 50 CONTINUE 55 L1 = L1+1 GO TO 35 60 CONTINUE !C NOW BALANCE THE SUBMATRIX IN ROWS !C L1 THROUGH K1 65 K = L1 L = K1 IF (K1.LT.L1) GO TO 75 DO 70 I=L1,K1 D(I) = ONE 70 CONTINUE 75 NOCONV = 0 IF (K1.LT.L1) GO TO 120 DO 115 I=L1,K1 C = ZERO R = ZERO DO 80 J=L1,K1 IF (J.EQ.I) GO TO 80 C = C+ABS(A(J,I)) R = R+ABS(A(I,J)) 80 CONTINUE G = R/B F = ONE S = C+R 85 IF (C.GE.G) GO TO 90 F = F * B C = C*B2 GO TO 85 90 G = R*B 95 IF (C.LT.G) GO TO 100 F = F/B C = C/B2 GO TO 95 !C NOW BALANCE 100 IF ((C+R)/F.GE.P95*S) GO TO 115 G = ONE/F D(I) = D(I)*F NOCONV = 1 DO 105 J=L1,N A(I,J) = A(I,J)*G 105 CONTINUE DO 110 J=1,K1 A(J,I) = A(J,I)*F 110 CONTINUE 115 CONTINUE 120 IF (NOCONV.EQ.1) GO TO 75 RETURN END !C IMSL ROUTINE NAME - EBBCKF !C !C----------------------------------------------------------------------- !C !C COMPUTER - HP1000/SINGLE !C !C LATEST REVISION - JUNE 1, 1982 !C !C PURPOSE - NUCLEUS CALLED ONLY BY IMSL ROUTINE EIGRF !C !C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 !C - DOUBLE/H36,H48,H60 !C !C REQD. IMSL ROUTINES - NONE REQUIRED !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C COPYRIGHT - 1982 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- !C SUBROUTINE EBBCKF (D,Z,K,L,MM,N,IZ) !C SPECIFICATIONS FOR ARGUMENTS INTEGER K,L,MM,N,IZ REAL Z(IZ,1),D(1) !C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER I,J,KM1,II,JJ,LP1 REAL S !C COLUMN SCALE Z BY APPROPRIATE D VALUE !C FIRST EXECUTABLE STATEMENT IF (L.EQ.0) GO TO 15 DO 10 I=K,L S = D(I) DO 5 J=1,MM Z(I,J) = Z(I,J)*S 5 CONTINUE 10 CONTINUE !C INTERCHANGE ROWS IF PERMUTATIONS !C OCCURRED IN EBALAF 15 IF (K.EQ.1) GO TO 30 KM1 = K-1 DO 25 I=1,KM1 II = K-I JJ = D(II) IF (II.EQ.JJ) GO TO 25 DO 20 J=1,MM S = Z(II,J) Z(II,J) = Z(JJ,J) Z(JJ,J) = S 20 CONTINUE 25 CONTINUE 30 IF (L.EQ.N) GO TO 45 LP1 = L+1 DO 40 II=LP1,N JJ = D(II) IF (II.EQ.JJ) GO TO 40 DO 35 J=1,MM S = Z(II,J) Z(II,J) = Z(JJ,J) Z(JJ,J) = S 35 CONTINUE 40 CONTINUE 45 RETURN END !C IMSL ROUTINE NAME - EHBCKF !C !C----------------------------------------------------------------------- !C !C COMPUTER - HP1000/SINGLE !C !C LATEST REVISION - NOVEMBER 1, 1984 !C !C PURPOSE - NUCLEUS CALLED ONLY BY IMSL ROUTINE EIGRF !C !C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 !C - DOUBLE/H36,H48,H60 !C !C REQD. IMSL ROUTINES - NONE REQUIRED !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C COPYRIGHT - 1982 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- !C SUBROUTINE EHBCKF (Z,H,D,N,MM,IZH,K,L) !C SPECIFICATIONS FOR ARGUMENTS INTEGER N,MM,IZH,K,L REAL Z(IZH,1),H(IZH,1),D(1) !C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER LM2,KI,LTEMP,M,MA,MP2,I,J REAL G,ZERO DATA ZERO/0.0/ !C ADAPTED FROM EISPACK ROUTINE ORTBAK !C FIRST EXECUTABLE STATEMENT LM2=L-2 IF(LM2.LT.K) GO TO 9005 LTEMP=LM2+K DO 30 KI=K,LM2 M=LTEMP-KI MA=M+1 IF(H(MA,M).EQ.ZERO) GO TO 30 MP2=M+2 IF(MP2.GT.L) GO TO 10 DO 5 I=MP2,L D(I)=H(I,M) 5 CONTINUE 10 IF(MA.GT.L) GO TO 30 DO 25 J=1,MM G=ZERO DO 15 I=MA,L G=G+D(I)*Z(I,J) 15 CONTINUE !C DOUBLE DIVISION AVOIDS POSSIBLE !C UNDERFLOW G = (G/D(MA))/H(MA,M) DO 20 I=MA,L Z(I,J)=Z(I,J)+G*D(I) 20 CONTINUE 25 CONTINUE 30 CONTINUE 9005 RETURN END !C IMSL ROUTINE NAME - EHESSF !C !C----------------------------------------------------------------------- !C !C COMPUTER - HP1000/SINGLE !C !C LATEST REVISION - JUNE 1, 1982 !C !C PURPOSE - NUCLEUS CALLED ONLY BY IMSL ROUTINE EIGRF !C !C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 !C - SINGLE/H36,H48,H60 !C !C REQD. IMSL ROUTINES - NONE REQUIRED !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C COPYRIGHT - 1982 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- !C SUBROUTINE EHESSF (A,K,L,N,IA,D) !C SPECIFICATIONS FOR ARGUMENTS INTEGER K,L,N,IA REAL A(IA,N),D(N) !C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER LA,KP1,M,I,MP,II,J,JJ REAL F,G,H,SCALE,ZERO DATA ZERO/0.0/ !C FIRST EXECUTABLE STATEMENT LA = L - 1 KP1 = K + 1 IF (LA .LT. KP1) GO TO 50 DO 45 M = KP1, LA H = ZERO D(M) = ZERO SCALE = ZERO !C SCALE COLUMN DO 5 I = M, L SCALE = SCALE + ABS(A(I,M-1)) 5 CONTINUE IF (SCALE .EQ. ZERO ) GO TO 45 MP = M + L !C DO 10 I=L,M,-1 DO 10 II = M, L I = MP - II D(I) = A(I,M-1) / SCALE H = H + D(I) * D(I) 10 CONTINUE G = -SIGN(SQRT(H),D(M)) H = H - D(M) * G D(M) = D(M) - G !C FORM (I-(U*UT)/H) * A DO 25 J = M,N F = ZERO !C DO 15 I=L,M,-1 DO 15 II = M, L I = MP - II F = F + D(I) * A(I,J) 15 CONTINUE F = F / H DO 20 I = M, L A(I,J) = A(I,J) - F * D(I) 20 CONTINUE 25 CONTINUE !C FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) DO 40 I = 1,L F = ZERO !C DO 30 J=L,M,-1 DO 30 JJ = M, L J = MP - JJ F = F + D(J) * A(I,J) 30 CONTINUE F = F / H DO 35 J = M, L A(I,J) = A(I,J) - F * D(J) 35 CONTINUE 40 CONTINUE D(M) = SCALE * D(M) A(M,M-1) = SCALE * G 45 CONTINUE 50 RETURN END !C IMSL ROUTINE NAME - EIGRF !C !C----------------------------------------------------------------------- !C !C COMPUTER - HP1000/SINGLE !C !C LATEST REVISION - JANUARY 1, 1978 !C !C PURPOSE - EIGENVALUES AND (OPTIONALLY) EIGENVECTORS OF !C A REAL GENERAL MATRIX IN FULL STORAGE MODE !C !C USAGE - CALL EIGRF (A,N,IA,IJOB,W,Z,IZ,WK,IER) !C !C ARGUMENTS A - THE INPUT REAL GENERAL MATRIX OF ORDER N !C WHOSE EIGENVALUES AND EIGENVECTORS ARE !C TO BE COMPUTED. INPUT A IS DESTROYED IF !C IJOB IS EQUAL TO 0 OR 1. !C N - THE INPUT ORDER OF THE MATRIX A. !C IA - THE INPUT ROW DIMENSION OF MATRIX A EXACTLY !C AS SPECIFIED IN THE DIMENSION STATEMENT IN !C THE CALLING PROGRAM. !C IJOB - THE INPUT OPTION PARAMETER. WHEN !C IJOB = 0, COMPUTE EIGENVALUES ONLY !C IJOB = 1, COMPUTE EIGENVALUES AND EIGEN- !C VECTORS. !C IJOB = 2, COMPUTE EIGENVALUES, EIGENVECTORS !C AND PERFORMANCE INDEX. !C IJOB = 3, COMPUTE PERFORMANCE INDEX ONLY. !C IF THE PERFORMANCE INDEX IS COMPUTED, IT IS !C RETURNED IN WK(1). THE ROUTINES HAVE !C PERFORMED (WELL, SATISFACTORILY, POORLY) IF !C WK(1) IS (LESS THAN 1, BETWEEN 1 AND 100, !C GREATER THAN 100). !C W - THE OUTPUT COMPLEX VECTOR OF LENGTH N, !C CONTAINING THE EIGENVALUES OF A. !C NOTE - THE ROUTINE TREATS W AS A REAL VECTOR !C OF LENGTH 2*N. AN APPROPRIATE EQUIVALENCE !C STATEMENT MAY BE REQUIRED. SEE DOCUMENT !C EXAMPLE. !C Z - THE OUTPUT N BY N COMPLEX MATRIX CONTAINING !C THE EIGENVECTORS OF A. !C THE EIGENVECTOR IN COLUMN J OF Z CORRES- !C PONDS TO THE EIGENVALUE W(J). !C IF IJOB = 0, Z IS NOT USED. !C NOTE - THE ROUTINE TREATS Z AS A REAL VECTOR !C OF LENGTH 2*N*N. AN APPROPRIATE EQUIVALENCE !C STATEMENT MAY BE REQUIRED. SEE DOCUMENT !C EXAMPLE. !C IZ - THE INPUT ROW DIMENSION OF MATRIX Z EXACTLY !C AS SPECIFIED IN THE DIMENSION STATEMENT IN !C THE CALLING PROGRAM. IZ MUST BE GREATER !C THAN OR EQUAL TO N IF IJOB IS NOT EQUAL TO !C ZERO. !C WK - WORK AREA, THE LENGTH OF WK DEPENDS !C ON THE VALUE OF IJOB, WHEN !C IJOB = 0, THE LENGTH OF WK IS AT LEAST N. !C IJOB = 1, THE LENGTH OF WK IS AT LEAST 2N. !C IJOB = 2, THE LENGTH OF WK IS AT LEAST !C (2+N)N. !C IJOB = 3, THE LENGTH OF WK IS AT LEAST 1. !C IER - ERROR PARAMETER. (OUTPUT) !C TERMINAL ERROR !C IER = 128+J, INDICATES THAT EQRH3F FAILED !C TO CONVERGE ON EIGENVALUE J. EIGENVALUES !C J+1,J+2,...,N HAVE BEEN COMPUTED CORRECTLY !C EIGENVALUES 1,...,J ARE SET TO ZERO. !C IF IJOB = 1 OR 2 EIGENVECTORS ARE SET TO !C ZERO. THE PERFORMANCE INDEX IS SET TO 1000. !C WARNING ERROR (WITH FIX) !C IER = 66, INDICATES IJOB IS LESS THAN 0 OR !C IJOB IS GREATER THAN 3. IJOB SET TO 1. !C IER = 67, INDICATES IJOB IS NOT EQUAL TO !C ZERO, AND IZ IS LESS THAN THE ORDER OF !C MATRIX A. IJOB IS SET TO ZERO. !C !C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 !C - SINGLE/H36,H48,H60 !C !C REQD. IMSL ROUTINES - EBALAF,EBBCKF,EHBCKF,EHESSF,EQRH3F,UERTST, !C UGETIO !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- !C SUBROUTINE EIGRF (A,N,IA,IJOB,W,Z,IZ,WK,IER) !C SPECIFICATIONS FOR ARGUMENTS INTEGER N,IA,IJOB,IZ,IER REAL A(IA,1),WK(N,1),W(1),Z(1) !C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER JER,IZ2,K,L,I,N1,N2,II,JJ,NP1,IIZ,NPI,JW,J,IS,IG,IGZ,LW,LLZ,KKZ,LZ,KZ REAL ANORM,ASUM,PI,SUMZ,SUMR,SUMI,S,TEN,RDELP,ZERO,ONE,THOUS,AN,Z11 DATA RDELP/1.1921E-07/ DATA ZERO,ONE/0.0,1.0/,TEN/10.0/,THOUS/1000.0/ !C INITIALIZE ERROR PARAMETERS !C FIRST EXECUTABLE STATEMENT IER = 0 JER = 0 IZ2 = IZ+IZ IF (IJOB .GE. 0 .AND. IJOB .LE. 3) GO TO 5 !C WARNING ERROR - IJOB IS NOT IN THE !C RANGE IER = 66 IJOB = 1 GO TO 10 5 IF (IJOB .EQ. 0) GO TO 16 10 IF (IZ .GE. N) GO TO 15 !C WARNING ERROR - IZ IS LESS THAN N !C EIGENVECTORS CAN NOT BE COMPUTED, !C IJOB SET TO ZERO IER = 67 IJOB = 0 15 IF (IJOB .EQ. 3) GO TO 95 !C PACK A INTO AN N BY N ARRAY 16 K = 1 L = 1 DO 20 J=1,N DO 20 I=1,N A(K,L) = A(I,J) !C SAVE INPUT A IF IJOB = 2 IF (IJOB .EQ. 2) WK(I,J)=A(I,J) K = K+1 IF (K .GT. IA) K = 1 IF (K .EQ. 1) L = L+1 20 CONTINUE N1 = 1 IF (IJOB .EQ. 2) N1 = N+1 N2 = N1+1 IF (IJOB .EQ. 0) N2 = 1 !C BALANCE THE INPUT A CALL EBALAF (A,N,N,WK(1,N1),K,L) IF (IJOB .EQ. 0 .AND. L .EQ. 0) GO TO 35 !C IF L = 0, A IS ALREADY IN HESSENBERG !C FORM CALL EHESSF (A,K,L,N,N,WK(1,N2)) IF (IJOB .EQ. 0) GO TO 35 !C SET Z IDENTITY MATRIX II = 1 JJ = 1 NP1 = N+1 DO 30 I=1,N DO 25 J=1,N Z(II) = ZERO II = II+1 25 CONTINUE Z(JJ) = ONE JJ = JJ+NP1 30 CONTINUE CALL EHBCKF (Z,A,WK(1,N2),N,N,N,K,L) IIZ = N 35 IF (IJOB .EQ. 0) IIZ = 1 IF(IJOB .EQ. 0 .AND. N .EQ. 1) Z11 = Z(1) CALL EQRH3F (A,N,N,K,L,W(1),W(N+1),Z,IIZ,JER) IF(IJOB .EQ. 0 .AND. N .EQ. 1) Z(1) = Z11 IF (JER .GT. 128 .OR. IJOB .EQ. 0) GO TO 40 CALL EBBCKF (WK(1,N1),Z,K,L,N,N,N) !C CONVERT W (EIGENVALUES) TO COMPLEX !C FORMAT 40 DO 45 I=1,N NPI = N+I WK(I,N1) = W(NPI) 45 CONTINUE JW = N+N J = N DO 50 I=1,N W(JW-1) = W(J) W(JW) = WK(J,N1) JW = JW-2 J = J-1 50 CONTINUE IF (IJOB .EQ. 0) GO TO 9000 !C CONVERT Z (EIGENVECTORS) TO COMPLEX !C FORMAT Z(IZ,N) J = N 60 IF (J .LT. 1) GO TO 85 IF (W(J+J) .EQ. ZERO) GO TO 75 !C MOVE PAIR OF COMPLEX CONJUGATE !C EIGENVECTORS IS = IZ2*(J-1)+1 IG = N*(J-2)+1 IGZ = IG+N !C MOVE COMPLEX CONJUGATE EIGENVECTOR DO 65 I=1,N Z(IS) = Z(IG) Z(IS+1) = -Z(IGZ) IS = IS+2 IG = IG+1 IGZ = IGZ+1 65 CONTINUE !C MOVE COMPLEX EIGENVECTOR IS = IZ2*(J-2)+1 IG = IS+IZ2 DO 70 I=1,N Z(IS) = Z(IG) Z(IS+1) = -Z(IG+1) IS = IS+2 IG = IG+2 70 CONTINUE J = J-2 GO TO 60 !C MOVE REAL EIGENVECTOR 75 IS = IZ2*(J-1)+N+N IG = N*J DO 80 I=1,N Z(IS-1) = Z(IG) Z(IS) = ZERO IS = IS-2 IG = IG-1 80 CONTINUE J = J-1 GO TO 60 !C Z IS NOW IN COMPLEX FORMAT Z(IZ,N). !C NEXT, MOVE ORIGINAL MATRIX BACK !C TO A 85 IF (IJOB .LE. 1) GO TO 9000 DO 90 I=1,N DO 90 J=1,N A(I,J) = WK(I,J) 90 CONTINUE WK(1,1) = THOUS IF (JER .NE. 0) GO TO 9000 !C COMPUTE 1-NORM OF A 95 ANORM = ZERO DO 105 J=1,N ASUM = ZERO DO 100 I=1,N ASUM = ASUM+ABS(A(I,J)) 100 CONTINUE ANORM = AMAX1(ANORM,ASUM) 105 CONTINUE IF (ANORM .EQ. ZERO) ANORM = ONE !C COMPUTE PERFORMANCE INDEX PI = ZERO LW = 1 LLZ = 0 KKZ = 0 DO 120 J=1,N S = ZERO SUMZ = ZERO LZ = LLZ+1 KZ = KKZ+1 LW = J+J-1 DO 115 L=1,N SUMZ = SUMZ+CABS(CMPLX(Z(LZ),Z(LZ+1))) KZ = KKZ+1 SUMR = -W(LW)*Z(LZ)+W(LW+1)*Z(LZ+1) SUMI = -W(LW)*Z(LZ+1)-W(LW+1)*Z(LZ) DO 110 K=1,N SUMR =SUMR+A(L,K)*Z(KZ) SUMI = SUMI+A(L,K)*Z(KZ+1) KZ = KZ+2 110 CONTINUE S = S+CABS(CMPLX(SUMR,SUMI)) LZ = LZ+2 115 CONTINUE PI = AMAX1(PI,S/SUMZ) KKZ = KKZ+IZ2 LLZ = LLZ+IZ2 120 CONTINUE AN = N PI = PI/(ANORM*TEN*AN*RDELP) WK(1,1) = PI 9000 CONTINUE IF (IER .NE. 0) CALL UERTST (IER,'EIGRF ') IF (JER .EQ. 0) GO TO 9005 IER = JER CALL UERTST (IER,'EIGRF ') 9005 RETURN END !C !C----------------------------------------------------------------------- !C !C COMPUTER - HP1000/SINGLE !C !C LATEST REVISION - NOVEMBER 1, 1984 !C !C PURPOSE - NUCLEUS CALLED ONLY BY IMSL ROUTINE EIGRF !C !C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 !C - DOUBLE/H36,H48,H60 !C !C REQD. IMSL ROUTINES - UERTST,UGETIO !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C COPYRIGHT - 1982 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- !C SUBROUTINE EQRH3F (H,N,IH,K,L,WR,WI,Z,IZ,IER) !C SPECIFICATIONS FOR ARGUMENTS INTEGER N,IH,K,L,IZ,IER REAL H(IH,N),WR(N),WI(N),Z(IZ,N) !C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER I,IEN,ITS,IENM2,NPL,LL,LB,NAML,MM,M,MP2,KA,NA,J,JJ REAL T3(2),RDELP,P4,P5,P7,ZERO,ONE,T,X,Y,W,S,ZZ,R,P,Q,RNORM,RA,SA,SCALE,VR,VI COMPLEX Z3 LOGICAL NOTLAS EQUIVALENCE (Z3,T3(1)) DATA RDELP/1.1921E-07/ DATA P4 /0.4375/,P5 /0.5/,P7 /0.75/,ZERO /0.0/,ONE/1.0/ !C FIRST EXECUTABLE STATEMENT IER = 0 !C STORE ROOTS ISOLATED BY EBALAF RNORM = 0.0 KA = 1 DO 10 I=1,N DO 5 J=KA,N 5 RNORM = RNORM+ABS(H(I,J)) KA = I IF (I.GE.K .AND. I.LE.L) GO TO 10 WR(I) = H(I,I) WI(I) = ZERO 10 CONTINUE IEN = L T = ZERO !C SEARCH FOR NEXT EIGENVALUES 15 IF (IEN.LT.K) GO TO 145 ITS = 0 NA = IEN-1 IENM2 = NA-1 !C LOOK FOR SINGLE SMALL SUB-DIAGONAL !C ELEMENT 20 NPL = IEN+K DO 25 LL=K,IEN LB = NPL-LL IF (LB.EQ.K) GO TO 30 S = ABS(H(LB-1,LB-1))+ABS(H(LB,LB)) IF (S.EQ.0.0) S = RNORM IF (ABS(H(LB,LB-1)).LE.RDELP*S) GO TO 30 25 CONTINUE 30 X = H(IEN,IEN) IF (LB.EQ.IEN) GO TO 110 Y = H(NA,NA) W = H(IEN,NA)*H(NA,IEN) IF (LB.EQ.NA) GO TO 115 IF (ITS.EQ.30) GO TO 250 !C FORM SHIFT IF (ITS.NE.10 .AND. ITS.NE.20) GO TO 40 T = T+X DO 35 I=K,IEN H(I,I) = H(I,I)-X 35 CONTINUE S = ABS(H(IEN,NA))+ABS(H(NA,IENM2)) X = P7*S Y = X W = -P4*S*S 40 ITS = ITS+1 !C LOOK FOR TWO CONSECUTIVE SMALL !C SUB-DIAGONAL ELEMENTS NAML = IENM2+LB DO 45 MM=LB,IENM2 M = NAML-MM ZZ = H(M,M) R = X-ZZ S = Y-ZZ P = (R*S-W)/H(M+1,M)+H(M,M+1) Q = H(M+1,M+1)-ZZ-R-S R = H(M+2,M+1) S = ABS(P)+ABS(Q)+ABS(R) P = P/S Q = Q/S R = R/S IF (M.EQ.LB) GO TO 50 IF (ABS(H(M,M-1))*(ABS(Q)+ABS(R)).LE.RDELP*ABS(P)*(ABS(H(M-1,M-1))+ABS(ZZ)+ABS(H(M+1,M+1)))) GO TO 50 45 CONTINUE 50 MP2 = M+2 DO 55 I=MP2,IEN H(I,I-2) = ZERO IF (I.EQ.MP2) GO TO 55 H(I,I-3) = ZERO 55 CONTINUE !C DOUBLE QR STEP INVOLVING ROWS !C L TO EN AND COLUMNS M TO EN DO 105 KA=M,NA NOTLAS = KA.NE.NA IF (KA.EQ.M) GO TO 60 P = H(KA,KA-1) Q = H(KA+1,KA-1) R = ZERO IF (NOTLAS) R = H(KA+2,KA-1) X = ABS(P)+ABS(Q)+ABS(R) IF (X.EQ.ZERO) GO TO 105 P = P/X Q = Q/X R = R/X 60 CONTINUE S = SIGN(SQRT(P*P+Q*Q+R*R),P) IF (KA.EQ.M) GO TO 65 H(KA,KA-1) = -S*X GO TO 70 65 IF (LB.NE.M) H(KA,KA-1) = -H(KA,KA-1) 70 P = P+S X = P/S Y = Q/S ZZ = R/S Q = Q/P R = R/P !C ROW MODIFICATION DO 80 J=KA,N P = H(KA,J)+Q*H(KA+1,J) IF (.NOT.NOTLAS) GO TO 75 P = P+R*H(KA+2,J) H(KA+2,J) = H(KA+2,J)-P*ZZ 75 H(KA+1,J) = H(KA+1,J)-P*Y H(KA,J) = H(KA,J)-P*X 80 CONTINUE J = MIN0(IEN,KA+3) !C COLUMN MODIFICATION DO 90 I=1,J P = X*H(I,KA)+Y*H(I,KA+1) IF (.NOT.NOTLAS) GO TO 85 P = P+ZZ*H(I,KA+2) H(I,KA+2) = H(I,KA+2)-P*R 85 H(I,KA+1) = H(I,KA+1)-P*Q H(I,KA) = H(I,KA)-P 90 CONTINUE IF (IZ.LT.N) GO TO 105 !C ACCUMULATE TRANSFORMATIONS DO 100 I=K,L P = X*Z(I,KA)+Y*Z(I,KA+1) IF (.NOT.NOTLAS) GO TO 95 P = P+ZZ*Z(I,KA+2) Z(I,KA+2) = Z(I,KA+2)-P*R 95 Z(I,KA+1) = Z(I,KA+1)-P*Q Z(I,KA) = Z(I,KA)-P 100 CONTINUE 105 CONTINUE GO TO 20 !C ONE ROOT FOUND 110 H(IEN,IEN) = X+T WR(IEN) = H(IEN,IEN) WI(IEN) = ZERO IEN = NA GO TO 15 !C TWO ROOTS FOUND 115 P = (Y-X)*P5 Q = P*P+W ZZ = SQRT(ABS(Q)) H(IEN,IEN) = X+T X = H(IEN,IEN) H(NA,NA) = Y+T IF (Q.LT.ZERO) GO TO 135 !C REAL PAIR ZZ = P+SIGN(ZZ,P) WR(NA) = X+ZZ WR(IEN) = WR(NA) IF (ZZ.NE.ZERO) WR(IEN) = X-W/ZZ WI(NA) = ZERO WI(IEN) = ZERO X = H(IEN,NA) !C EMPLOY SCALE FACTOR IN CASE X AND !C ZZ ARE VERY SMALL SCALE = ABS(X) + ABS(ZZ) R = SCALE * SQRT( (X/SCALE)**2 + (ZZ/SCALE)**2 ) P = X/R Q = ZZ/R !C ROW MODIFICATION DO 120 J=NA,N ZZ = H(NA,J) H(NA,J) = Q*ZZ+P*H(IEN,J) H(IEN,J) = Q*H(IEN,J)-P*ZZ 120 CONTINUE !C COLUMN MODIFICATION DO 125 I=1,IEN ZZ = H(I,NA) H(I,NA) = Q*ZZ+P*H(I,IEN) H(I,IEN) = Q*H(I,IEN)-P*ZZ 125 CONTINUE IF (IZ.LT.N) GO TO 140 !C ACCUMULATE TRANSFORMATIONS DO 130 I=K,L ZZ = Z(I,NA) Z(I,NA) = Q*ZZ+P*Z(I,IEN) Z(I,IEN) = Q*Z(I,IEN)-P*ZZ 130 CONTINUE GO TO 140 !C COMPLEX PAIR 135 WR(NA) = X+P WR(IEN) = X+P WI(NA) = ZZ WI(IEN) = -ZZ 140 IEN = IENM2 GO TO 15 !C ALL ROOTS FOUND, NOW !C BACKSUBSTITUTE 145 IF (IZ.LT.N) GO TO 9005 IF (RNORM.EQ.ZERO) GO TO 9005 DO 220 NN=1,N IEN = N+1-NN P = WR(IEN) Q = WI(IEN) NA = IEN-1 IF (Q.GT.ZERO) GO TO 220 IF (Q.LT.ZERO) GO TO 180 !C REAL VECTOR M = IEN H(IEN,IEN) = ONE IF (NA.EQ.0) GO TO 220 DO 175 II=1,NA I = IEN-II W = H(I,I)-P R = H(I,IEN) IF (M.GT.NA) GO TO 155 DO 150 J=M,NA R = R+H(I,J)*H(J,IEN) 150 CONTINUE 155 IF (WI(I).GE.ZERO) GO TO 160 ZZ = W S = R GO TO 175 160 M = I IF (WI(I).NE.ZERO) GO TO 165 T = W IF (W.EQ.ZERO) T = RDELP*RNORM H(I,IEN) = -R/T GO TO 175 !C SOLVE REAL EQUATIONS 165 X = H(I,I+1) Y = H(I+1,I) Q = (WR(I)-P)*(WR(I)-P)+WI(I)*WI(I) T = (X*S-ZZ*R)/Q H(I,IEN) = T IF (ABS(X).LE.ABS(ZZ)) GO TO 170 H(I+1,IEN) = (-R-W*T)/X GO TO 175 170 H(I+1,IEN) = (-S-Y*T)/ZZ 175 CONTINUE !C END REAL VECTOR GO TO 220 !C LAST VECTOR COMPONENT CHOSEN !C IMAGINARY SO THAT EIGENVECTOR !C MATRIX IS TRIANGULAR 180 M = NA !C COMPLEX VECTOR IF (ABS(H(IEN,NA)).LE.ABS(H(NA,IEN))) GO TO 185 H(NA,NA) = Q/H(IEN,NA) H(NA,IEN) = -(H(IEN,IEN)-P)/H(IEN,NA) GO TO 190 185 CONTINUE Z3 = CMPLX(ZERO,-H(NA,IEN))/CMPLX(H(NA,NA)-P,Q) H(NA,NA) = T3(1) H(NA,IEN) = T3(2) 190 H(IEN,NA) = ZERO H(IEN,IEN) = ONE IENM2 = NA-1 IF (IENM2.EQ.0) GO TO 220 DO 215 II=1,IENM2 I = NA-II W = H(I,I)-P RA = ZERO SA = H(I,IEN) DO 195 J=M,NA RA = RA+H(I,J)*H(J,NA) SA = SA+H(I,J)*H(J,IEN) 195 CONTINUE IF (WI(I).GE.ZERO) GO TO 200 ZZ = W R = RA S = SA GO TO 215 200 M = I IF (WI(I).NE.ZERO) GO TO 205 Z3 = CMPLX(-RA,-SA)/CMPLX(W,Q) H(I,NA) = T3(1) H(I,IEN) = T3(2) GO TO 215 !C SOLVE COMPLEX EQUATIONS 205 X = H(I,I+1) Y = H(I+1,I) VR = (WR(I)-P)*(WR(I)-P)+WI(I)*WI(I)-Q*Q VI = (WR(I)-P)*Q VI = VI+VI IF (VR.EQ.ZERO .AND. VI.EQ.ZERO) VR = RDELP*RNORM*(ABS(W)+ABS(Q)+ABS(X)+ABS(Y)+ABS(ZZ)) Z3 = CMPLX(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA)/CMPLX(VR,VI) H(I,NA) = T3(1) H(I,IEN) = T3(2) IF (ABS(X).LE.ABS(ZZ)+ABS(Q)) GO TO 210 H(I+1,NA) = (-RA-W*H(I,NA)+Q*H(I,IEN))/X H(I+1,IEN) = (-SA-W*H(I,IEN)-Q*H(I,NA))/X GO TO 215 210 CONTINUE Z3 = CMPLX(-R-Y*H(I,NA),-S-Y*H(I,IEN))/CMPLX(ZZ,Q) H(I+1,NA) = T3(1) H(I+1,IEN) = T3(2) 215 CONTINUE !C END COMPLEX VECTOR 220 CONTINUE !C END BACKSUBSTITUTION !C VECTORS OF ISOLATED ROOTS DO 230 I=1,N IF (I.GE.K .AND. I.LE.L) GO TO 230 DO 225 J=I,N Z(I,J) = H(I,J) 225 CONTINUE 230 CONTINUE IF (L.EQ.0) GO TO 9005 !C MULTIPLY BY TRANSFORMATION MATRIX DO 245 JJ=K,N J = N+K-JJ M = MIN0(J,L) DO 240 I=K,L ZZ = ZERO DO 235 KA=K,M ZZ = ZZ+Z(I,KA)*H(KA,J) 235 CONTINUE Z(I,J) = ZZ 240 CONTINUE 245 CONTINUE GO TO 9005 !C NO CONVERGENCE AFTER 30 ITERATIONS !C SET ERROR INDICATOR TO THE INDEX !C OF THE CURRENT EIGENVALUE 250 IER = 128+IEN DO 255 I=1,IEN WR(I) = ZERO WI(I) = ZERO 255 CONTINUE IF (IZ.LT.N) GO TO 9000 DO 265 I=1,N DO 260 J=1,N Z(I,J) = ZERO 260 CONTINUE 265 CONTINUE 9000 CONTINUE CALL UERTST (IER,'EQRH3F') 9005 RETURN END C IMSL ROUTINE NAME - LEQT1P C C----------------------------------------------------------------------- C C COMPUTER - CDC205/SINGLE C C LATEST REVISION - JANUARY 1, 1978 C C PURPOSE - LINEAR EQUATION SOLUTION - POSITIVE DEFINITE C MATRIX - SYMMETRIC STORAGE MODE - SPACE C ECONOMIZER SOLUTION C C USAGE - CALL LEQT1P (A,M,N,B,IB,IDGT,D1,D2,IER) C C ARGUMENTS A - INPUT VECTOR OF LENGTH N(N+1)/2 CONTAINING THE C N BY N COEFFICIENT MATRIX OF THE EQUATION C AX = B. A IS A POSITIVE DEFINITE SYMMETRIC C MATRIX STORED IN SYMMETRIC STORAGE MODE. C ON OUTPUT, A IS REPLACED BY THE LOWER C TRIANGULAR MATRIX L WHERE A = L*L-TRANSPOSE. C L IS STORED IN SYMMETRIC STORAGE MODE WITH C THE DIAGONAL ELEMENTS OF L IN RECIPROCAL C FORM. !C M - NUMBER OF RIGHT HAND SIDES (COLUMNS IN B). !C (INPUT) !C N - ORDER OF A AND NUMBER OF ROWS IN B. (INPUT) !C B - INPUT MATRIX OF DIMENSION N BY M CONTAINING !C THE RIGHT-HAND SIDES OF THE EQUATION !C AX = B. !C ON OUTPUT, THE N BY M SOLUTION MATRIX X !C REPLACES B. !C IB - ROW DIMENSION OF B EXACTLY AS SPECIFIED IN !C THE DIMENSION STATEMENT IN THE CALLING !C PROGRAM. (INPUT) !C IDGT - THE ELEMENTS OF A ARE ASSUMED TO BE CORRECT !C TO IDGT DECIMAL DIGITS. (CURRENTLY NOT USED) !C D1,D2 - COMPONENTS OF THE DETERMINANT OF A. !C DETERMINANT(A) = D1*2.**D2. (OUTPUT) !C IER - ERROR PARAMETER. (OUTPUT) !C TERMINAL ERROR !C IER = 129 INDICATES THAT THE INPUT MATRIX !C A IS ALGORITHMICALLY NOT POSITIVE !C DEFINITE. (SEE THE CHAPTER L PRELUDE). !C !C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 !C - SINGLE/H36,H48,H60 !C !C REQD. IMSL ROUTINES - LUDECP,LUELMP,UERTST,UGETIO !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- !C SUBROUTINE LEQT1P (A,M,N,B,IB,IDGT,D1,D2,IER) !C DIMENSION A(1),B(IB,1) !C FIRST EXECUTABLE STATEMENT !C INITIALIZE IER IER = 0 !C DECOMPOSE A CALL LUDECP (A,A,N,D1,D2,IER) IF (IER.NE.0) GO TO 9000 !C PERFORM ELIMINATION DO 5 I = 1,M CALL LUELMP (A,B(1,I),N,B(1,I)) 5 CONTINUE GO TO 9005 9000 CONTINUE CALL UERTST(IER,'LEQT1P') 9005 RETURN END !C IMSL ROUTINE NAME - LUDECP !C !C----------------------------------------------------------------------- !C !C COMPUTER - HP1000/SINGLE !C !C LATEST REVISION - JANUARY 1, 1978 !C !C PURPOSE - DECOMPOSITION OF A POSITIVE DEFINITE MATRIX - !C SYMMETRIC STORAGE MODE !C !C USAGE - CALL LUDECP (A,UL,N,D1,D2,IER) !C !C ARGUMENTS A - INPUT VECTOR OF LENGTH N(N+1)/2 CONTAINING !C THE N BY N POSITIVE DEFINITE SYMMETRIC !C MATRIX STORED IN SYMMETRIC STORAGE MODE. !C UL - OUTPUT VECTOR OF LENGTH N(N+1)/2 CONTAINING !C THE DECOMPOSED MATRIX L SUCH THAT A = L* !C L-TRANSPOSE. L IS STORED IN SYMMETRIC !C STORAGE MODE. THE DIAGONAL OF L CONTAINS THE !C RECIPROCALS OF THE ACTUAL DIAGONAL ELEMENTS. !C N - ORDER OF A. (INPUT) !C D1,D2 - COMPONENTS OF THE DETERMINANT OF A. !C DETERMINANT(A) = D1*2.**D2. (OUTPUT) !C IER - ERROR PARAMETER. (OUTPUT) !C TERMINAL ERROR !C IER = 129 INDICATES THAT MATRIX A IS !C ALGORITHMICALLY NOT POSITIVE DEFINITE. !C (SEE THE CHAPTER L PRELUDE). !C !C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 !C - SINGLE/H36,H48,H60 !C !C REQD. IMSL ROUTINES - UERTST,UGETIO !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- !C SUBROUTINE LUDECP (A,UL,N,D1,D2,IER) DIMENSION A(1),UL(1) REAL Q DATA ZERO,ONE,FOUR,SIXTN,SIXTH/0.0,1.,4.,16.,.0625/ !C FIRST EXECUTABLE STATEMENT D1=ONE D2=ZERO RN = ONE/(N*SIXTN) IP = 1 IER=0 DO 45 I = 1,N IQ = IP IR = 1 DO 40 J = 1,I X = A(IP) IF (J .EQ. 1) GO TO 10 DO 5 K=IQ,IP1 X = X - UL(K) * UL(IR) IR = IR+1 5 CONTINUE 10 IF (I.NE.J) GO TO 30 D1 = D1*X Q = A(IP)+X*RN IF (Q .LE. A(IP)) GO TO 50 15 IF (ABS(D1).LE.ONE) GO TO 20 D1 = D1 * SIXTH D2 = D2 + FOUR GO TO 15 20 IF (ABS(D1) .GE. SIXTH) GO TO 25 D1 = D1 * SIXTN D2 = D2 - FOUR GO TO 20 25 UL(IP) = ONE/SQRT(X) GO TO 35 30 UL(IP) = X * UL(IR) 35 IP1 = IP IP = IP+1 IR = IR+1 40 CONTINUE 45 CONTINUE GO TO 9005 50 IER = 129 9000 CONTINUE CALL UERTST(IER,'LUDECP') 9005 RETURN END !C IMSL ROUTINE NAME - LUELMP !C !C----------------------------------------------------------------------- !C !C COMPUTER - HP1000/SINGLE !C !C LATEST REVISION - JANUARY 1, 1978 !C !C PURPOSE - ELIMINATION PART OF THE SOLUTION OF AX=B - !C POSITIVE DEFINITE MATRIX - SYMMETRIC !C STORAGE MODE !C !C USAGE - CALL LUELMP (A,B,N,X) !C !C ARGUMENTS A - INPUT VECTOR OF LENGTH N(N+1)/2 CONTAINING !C THE N BY N MATRIX L WHERE A = L*L-TRANSPOSE. !C L IS A LOWER TRIANGULAR MATRIX STORED IN !C SYMMETRIC STORAGE MODE. THE MAIN DIAGONAL !C ELEMENTS OF L ARE STORED IN RECIPROCAL !C FORM. MATRIX L MAY BE OBTAINED FROM IMSL !C ROUTINE LUDECP. !C B - VECTOR OF LENGTH N CONTAINING THE RIGHT HAND !C SIDE OF THE EQUATION AX = B. (INPUT) !C N - ORDER OF A AND THE LENGTH OF B AND X. (INPUT) !C X - VECTOR OF LENGTH N CONTAINING THE SOLUTION TO !C THE EQUATION AX = B. (OUTPUT) !C !C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32 !C - SINGLE/H36,H48,H60 !C !C REQD. IMSL ROUTINES - NONE REQUIRED !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- !C SUBROUTINE LUELMP (A,B,N,X) !C DIMENSION A(1),B(1),X(1) DATA ZERO/0.0/ !C FIRST EXECUTABLE STATEMENT !C SOLUTION OF LY = B IP=1 IW = 0 DO 15 I=1,N T=B(I) IM1 = I-1 IF (IW .EQ. 0) GO TO 9 IP=IP+IW-1 DO 5 K=IW,IM1 T = T-A(IP)*X(K) IP=IP+1 5 CONTINUE GO TO 10 9 IF (T .NE. ZERO) IW = I IP = IP+IM1 10 X(I)=T*A(IP) IP=IP+1 15 CONTINUE !C SOLUTION OF UX = Y N1 = N+1 DO 30 I = 1,N II = N1-I IP=IP-1 IS=IP IQ=II+1 T=X(II) IF (N.LT.IQ) GO TO 25 KK = N DO 20 K=IQ,N T = T - A(IS) * X(KK) KK = KK-1 IS = IS-KK 20 CONTINUE 25 X(II)=T*A(IS) 30 CONTINUE RETURN END !C IMSL ROUTINE NAME - UERTST !C !C----------------------------------------------------------------------- !C !C COMPUTER - HP1000/SINGLE !C !C LATEST REVISION - JUNE 1, 1982 !C !C PURPOSE - PRINT A MESSAGE REFLECTING AN ERROR CONDITION !C !C USAGE - CALL UERTST (IER,NAME) !C !C ARGUMENTS IER - ERROR PARAMETER. (INPUT) !C IER = I+J WHERE !C I = 128 IMPLIES TERMINAL ERROR MESSAGE, !C I = 64 IMPLIES WARNING WITH FIX MESSAGE, !C I = 32 IMPLIES WARNING MESSAGE. !C J = ERROR CODE RELEVANT TO CALLING !C ROUTINE. !C NAME - A CHARACTER STRING OF LENGTH SIX PROVIDING !C THE NAME OF THE CALLING ROUTINE. (INPUT) !C !C PRECISION/HARDWARE - SINGLE/ALL !C !C REQD. IMSL ROUTINES - UGETIO,USPKD !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C REMARKS THE ERROR MESSAGE PRODUCED BY UERTST IS WRITTEN !C TO THE STANDARD OUTPUT UNIT. THE OUTPUT UNIT !C NUMBER CAN BE DETERMINED BY CALLING UGETIO AS !C FOLLOWS.. CALL UGETIO(1,NIN,NOUT). !C THE OUTPUT UNIT NUMBER CAN BE CHANGED BY CALLING !C UGETIO AS FOLLOWS.. !C NIN = 0 !C NOUT = NEW OUTPUT UNIT NUMBER !C CALL UGETIO(3,NIN,NOUT) !C SEE THE UGETIO DOCUMENT FOR MORE DETAILS. !C !C COPYRIGHT - 1982 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- SUBROUTINE UERTST (IER,NAME) !C SPECIFICATIONS FOR ARGUMENTS INTEGER IER CHARACTER NAME*(*) !C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER I,IEQ,IEQDF,IOUNIT,LEVEL,LEVOLD,NAMEQ(6),NAMSET(6),NAMUPK(6),NIN,NMTB DATA NAMSET/1HU,1HE,1HR,1HS,1HE,1HT/ DATA NAMEQ/6*1H / DATA LEVEL/4/,IEQDF/0/,IEQ/1H=/ !C UNPACK NAME INTO NAMUPK !C FIRST EXECUTABLE STATEMENT CALL USPKD (NAME,6,NAMUPK,NMTB) !C GET OUTPUT UNIT NUMBER CALL UGETIO(1,NIN,IOUNIT) !C CHECK IER IF (IER.GT.999) GO TO 25 IF (IER.LT.-32) GO TO 55 IF (IER.LE.128) GO TO 5 IF (LEVEL.LT.1) GO TO 30 !C PRINT TERMINAL MESSAGE IF (IEQDF.EQ.1) WRITE(IOUNIT,35) IER,NAMEQ,IEQ,NAMUPK IF (IEQDF.EQ.0) WRITE(IOUNIT,35) IER,NAMUPK GO TO 30 5 IF (IER.LE.64) GO TO 10 IF (LEVEL.LT.2) GO TO 30 !C PRINT WARNING WITH FIX MESSAGE IF (IEQDF.EQ.1) WRITE(IOUNIT,40) IER,NAMEQ,IEQ,NAMUPK IF (IEQDF.EQ.0) WRITE(IOUNIT,40) IER,NAMUPK GO TO 30 10 IF (IER.LE.32) GO TO 15 !C PRINT WARNING MESSAGE IF (LEVEL.LT.3) GO TO 30 IF (IEQDF.EQ.1) WRITE(IOUNIT,45) IER,NAMEQ,IEQ,NAMUPK IF (IEQDF.EQ.0) WRITE(IOUNIT,45) IER,NAMUPK GO TO 30 15 CONTINUE !C CHECK FOR UERSET CALL DO 20 I=1,6 IF (NAMUPK(I).NE.NAMSET(I)) GO TO 25 20 CONTINUE LEVOLD = LEVEL LEVEL = IER IER = LEVOLD IF (LEVEL.LT.0) LEVEL = 4 IF (LEVEL.GT.4) LEVEL = 4 GO TO 30 25 CONTINUE IF (LEVEL.LT.4) GO TO 30 !C PRINT NON-DEFINED MESSAGE IF (IEQDF.EQ.1) WRITE(IOUNIT,50) IER,NAMEQ,IEQ,NAMUPK IF (IEQDF.EQ.0) WRITE(IOUNIT,50) IER,NAMUPK 30 IEQDF = 0 RETURN 35 FORMAT(19H *** TERMINAL ERROR,10X,7H(IER = ,I3, 1 20H) FROM IMSL ROUTINE ,6A1,A1,6A1) 40 FORMAT(27H *** WARNING WITH FIX ERROR,2X,7H(IER = ,I3, 1 20H) FROM IMSL ROUTINE ,6A1,A1,6A1) 45 FORMAT(18H *** WARNING ERROR,11X,7H(IER = ,I3, 1 20H) FROM IMSL ROUTINE ,6A1,A1,6A1) 50 FORMAT(20H *** UNDEFINED ERROR,9X,7H(IER = ,I5, 1 20H) FROM IMSL ROUTINE ,6A1,A1,6A1) !C !C SAVE P FOR P = R CASE !C P IS THE PAGE NAMUPK !C R IS THE ROUTINE NAMUPK 55 IEQDF = 1 DO 60 I=1,6 60 NAMEQ(I) = NAMUPK(I) 65 RETURN END !C IMSL ROUTINE NAME - UGETIO !C !C----------------------------------------------------------------------- !C !C COMPUTER - HP1000/SINGLE !C !C LATEST REVISION - JUNE 1, 1981 !C !C PURPOSE - TO RETRIEVE CURRENT VALUES AND TO SET NEW !C VALUES FOR INPUT AND OUTPUT UNIT !C IDENTIFIERS. !C !C USAGE - CALL UGETIO(IOPT,NIN,NOUT) !C !C ARGUMENTS IOPT - OPTION PARAMETER. (INPUT) !C IF IOPT=1, THE CURRENT INPUT AND OUTPUT !C UNIT IDENTIFIER VALUES ARE RETURNED IN NIN !C AND NOUT, RESPECTIVELY. !C IF IOPT=2, THE INTERNAL VALUE OF NIN IS !C RESET FOR SUBSEQUENT USE. !C IF IOPT=3, THE INTERNAL VALUE OF NOUT IS !C RESET FOR SUBSEQUENT USE. !C NIN - INPUT UNIT IDENTIFIER. !C OUTPUT IF IOPT=1, INPUT IF IOPT=2. !C NOUT - OUTPUT UNIT IDENTIFIER. !C OUTPUT IF IOPT=1, INPUT IF IOPT=3. !C !C PRECISION/HARDWARE - SINGLE/ALL !C !C REQD. IMSL ROUTINES - NONE REQUIRED !C !C NOTATION - INFORMATION ON SPECIAL NOTATION AND !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C REMARKS EACH IMSL ROUTINE THAT PERFORMS INPUT AND/OR OUTPUT !C OPERATIONS CALLS UGETIO TO OBTAIN THE CURRENT UNIT !C IDENTIFIER VALUES. IF UGETIO IS CALLED WITH IOPT=2 OR !C IOPT=3, NEW UNIT IDENTIFIER VALUES ARE ESTABLISHED. !C SUBSEQUENT INPUT/OUTPUT IS PERFORMED ON THE NEW UNITS. !C !C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- SUBROUTINE UGETIO(IOPT,NIN,NOUT) !C SPECIFICATIONS FOR ARGUMENTS INTEGER IOPT,NIN,NOUT !C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER NIND,NOUTD DATA NIND/7/,NOUTD/6/ !C FIRST EXECUTABLE STATEMENT IF (IOPT.EQ.3) GO TO 10 IF (IOPT.EQ.2) GO TO 5 IF (IOPT.NE.1) GO TO 9005 NIN = NIND NOUT = NOUTD GO TO 9005 5 NIND = NIN GO TO 9005 10 NOUTD = NOUT 9005 RETURN END SUBROUTINE VSRTR(A,LA,IR) !C IMSL ROUTINE NAME - VSRTR !C !C----------------------------------------------------------------------- !C !C COMPUTER - CDC/SINGLE !C !C LATEST REVISION - JANUARY 1, 1978 !C !C PURPOSE - SORTING OF ARRAYS BY ALGEBRAIC VALUE - !C PERMUTATIONS RETURNED !C !C USAGE - CALL VSRTR (A,LA,IR) !C !C ARGUMENTS A - ON INPUT, A CONTAINS THE ARRAY TO BE SORTED. !C ON OUTPUT, A CONTAINS THE SORTED ARRAY. !C LA - INPUT VARIABLE CONTAINING THE NUMBER OF !C ELEMENTS IN THE ARRAY TO BE SORTED. !C IR - VECTOR OF LENGTH LA. !C ON INPUT, IR CONTAINS THE INTEGER VALUES !C 1,2,...,LA. SEE REMARKS. !C ON OUTPUT, IR CONTAINS A RECORD OF THE !C PERMUTATIONS MADE ON THE VECTOR A. !C !C PRECISION/HARDWARE - SINGLE/ALL !C !C REQD. IMSL ROUTINES - NONE REQUIRED !C !C CONVENTIONS IS AVAILABLE IN THE MANUAL !C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP !C !C REMARKS THE VECTOR IR MUST BE INITIALIZED BEFORE ENTERING !C VSRTR. ORDINARILY, IR(1)=1, IR(2)=2, ..., !C IR(LA)=LA. FOR WIDER APPLICABILITY, ANY INTEGER !C THAT IS TO BE ASSOCIATED WITH A(I) FOR I=1,2,...,LA !C MAY BE ENTERED INTO IR(I). !C COPYRIGHT - 1978 BY IMSL, INC. ALL RIGHTS RESERVED. !C !C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN !C APPLIED TO THIS CODE. NO OTHER WARRANTY, !C EXPRESSED OR IMPLIED, IS APPLICABLE. !C !C----------------------------------------------------------------------- !C DIMENSION A(1),IR(1) !C SPECIFICATIONS FOR ARGUMENTS !C SPECIFICATIONS FOR LOCAL VARIABLES INTEGER IU(21),IL(21),I,M,J,K,IJ,IT,L,ITT REAL T,TT,R !C FIRST EXECUTABLE STATEMENT IF (LA.LE.0) RETURN M = 1 I = 1 J = LA R = .375 5 IF (I.EQ.J) GO TO 45 IF (R.GT..5898437) GO TO 10 R = R+3.90625E-2 GO TO 15 10 R = R-.21875 15 K = I !C SELECT A CENTRAL ELEMENT OF THE !C ARRAY AND SAVE IT IN LOCATION T IJ = I+(J-I)*R T = A(IJ) IT = IR(IJ) !C IF FIRST ELEMENT OF ARRAY IS GREATER !C THAN T, INTERCHANGE WITH T IF (A(I).LE.T) GO TO 20 A(IJ) = A(I) A(I) = T T = A(IJ) IR(IJ) = IR(I) IR(I) = IT IT = IR(IJ) 20 L = J !C IF LAST ELEMENT OF ARRAY IS LESS THAN !C T, INTERCHANGE WITH T IF (A(J).GE.T) GO TO 30 A(IJ) = A(J) A(J) = T T = A(IJ) IR(IJ) = IR(J) IR(J) = IT IT = IR(IJ) !C IF FIRST ELEMENT OF ARRAY IS GREATER !C THAN T, INTERCHANGE WITH T IF (A(I).LE.T) GO TO 30 A(IJ) = A(I) A(I) = T T = A(IJ) IR(IJ) = IR(I) IR(I) = IT IT = IR(IJ) GO TO 30 25 IF (A(L).EQ.A(K)) GO TO 30 TT = A(L) A(L) = A(K) A(K) = TT ITT = IR(L) IR(L) = IR(K) IR(K) = ITT !C FIND AN ELEMENT IN THE SECOND HALF OF !C THE ARRAY WHICH IS SMALLER THAN T 30 L = L-1 IF (A(L).GT.T) GO TO 30 !C FIND AN ELEMENT IN THE FIRST HALF OF !C THE ARRAY WHICH IS GREATER THAN T 35 K = K+1 IF (A(K).LT.T) GO TO 35 !C INTERCHANGE THESE ELEMENTS IF (K.LE.L) GO TO 25 !C SAVE UPPER AND LOWER SUBSCRIPTS OF !C THE ARRAY YET TO BE SORTED IF (L-I.LE.J-K) GO TO 40 IL(M) = I IU(M) = L I = K M = M+1 GO TO 50 40 IL(M) = K IU(M) = J J = L M = M+1 GO TO 50 !C BEGIN AGAIN ON ANOTHER PORTION OF !C THE UNSORTED ARRAY 45 M = M-1 IF (M.EQ.0) RETURN I = IL(M) J = IU(M) 50 IF (J-I.GE.11) GO TO 15 IF (I.EQ.1) GO TO 5 I = I-1 55 I = I+1 IF (I.EQ.J) GO TO 45 T = A(I+1) IT = IR(I+1) IF (A(I).LE.T) GO TO 55 K = I 60 A(K+1) = A(K) IR(K+1) = IR(K) K = K-1 IF (T.LT.A(K)) GO TO 60 A(K+1) = T IR(K+1) = IT GO TO 55 END