C ACM ALGORITHM 570 C C LOPSI: A SIMULTANEOUS ITERATION ALGORITHM FOR REAL MATRICES C C BY W.J. STEWART AND A. JENNINGS C C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, JUNE 1981 C C 00000010 C 00000020 C ** DRIVER PROGRAM FOR LOP-SIDED SIMULTANEOUS ITERATION 00000030 C ALGORITHM, LOPSI. ** 00000040 C 00000050 C 00000060 C NOTE IN THE DATA STATEMENT, IOUNT1 AND IOUNT2 MUST BE INITIAL- 00000070 C **** IZED TO THE INPUT AND OUTPUT STREAM NUMBERS RESPECTIVELY. 00000080 C 00000090 C IOUNT2 MUST ALSO BE INITIALIZED IN THE SUBROUTINE CHOLES. 00000100 C 00000110 C 00000120 C 00000130 DOUBLE PRECISION A,B,ERR,G,H,P,TOL,U,V,W,X,Y 00000140 DIMENSION A(1207),IROW(1207),ICOL(1207),U(163,8),V(163,8),W(163,8)00000150 DIMENSION G(10,10),H(10,10),B(10,10),P(10,10) 00000160 DIMENSION X(10),Y(10),IUSED(10),ERR(10),INT(10) 00000170 DATA IOUNT1/5/,IOUNT2/6/ 00000180 IU = 163 00000190 IG = 10 WRITE(IOUNT2,701) WRITE(IOUNT2,702) C C ** INPUT NUMBER OF MATRICES, NMATX, TO BE ANALYZED. ** C READ(IOUNT1,601)NMATX WRITE(IOUNT2,703)NMATX C C DO 60 LOOP2 = 1,NMATX C C ** FOR EACH MATRIX INPUT ORDER OF MATRIX, N, NUMBER OF C NON-ZERO ELEMENTS, IA, AND NUMBER OF TESTS INVOLVING C THIS MATRIX, NEXS. ** C READ(IOUNT1,601)N,IA,NEXS WRITE(IOUNT2,704)LOOP2,N,IA C C ** INPUT MATRIX IN COMPACT FORM ** C READ(IOUNT1,602)(A(I),I=1,IA) READ(IOUNT1,601)(IROW(I),I=1,IA) READ(IOUNT1,601)(ICOL(I),I=1,IA) C C ** PRINT MATRIX IF ORDER IS LESS THAN 10 ** C IF(N.GT.10)GO TO 4 WRITE(IOUNT2,705) DO 2 I=1,N DO 1 J=1,N G(I,J) = 0.0 1 CONTINUE G(I,I) = 1.0 2 CONTINUE CALL PREMUL(N,A,IA,IROW,ICOL,N,G,IG,H,1,0,0,1) DO 3 I=1,N WRITE(IOUNT2,706)(H(I,J),J=1,N) 3 CONTINUE C C 4 DO 50 LOOP1 = 1,NEXS C C ** FOR EACH TEST ON THIS MATRIX, INPUT THE VALUES OF C PARAMETERS: M, IR, IS, TOL, LMAX, ITMAX AND LR. ** C READ(IOUNT1,603)M,IR,IS,TOL,LMAX,ITMAX,LR WRITE(IOUNT2,707)LOOP1 WRITE(IOUNT2,708) WRITE(IOUNT2,709)M,IR,IS,TOL,LMAX,ITMAX,LR C C ** INPUT INITIAL APPROXIMATIONS TO FIRST IS EIGENVECTORS. ** C IF(IS.GE.1)READ(IOUNT1,604)((U(I,J),I=1,N),J=1,IS) IF(IS.EQ.0.OR.IS.GT.10.OR.N.GT.50)GO TO 25 WRITE(IOUNT2,719) DO 21 I=1,N WRITE(IOUNT2,706)(U(I,J),J=1,IS) 21 CONTINUE C 25 DO 5 I=1,M ERR(I) = 1.23456789 5 CONTINUE C C C C CALL LOPSI(N,A,IA,IROW,ICOL,LR,IR,M,IS,TOL,ITMAX,LMAX,U,IU, 1 V,W,G,H,B,P,IG,INT,X,Y,ERR,ITTOT,INACT,IFAIL) C C C C C ** OUTPUT RESULTS OF CURRENT TEST ** C IF(IFAIL.EQ.1)GO TO 6 WRITE(IOUNT2,710) GO TO 7 6 WRITE(IOUNT2,711) 7 IF(M.EQ.N)WRITE(IOUNT2,712) IF(M.EQ.1)WRITE(IOUNT2,713) WRITE(IOUNT2,714) WRITE(IOUNT2,715) DO 8 I=1,M WRITE(IOUNT2,716)X(I),Y(I),ERR(I) 8 CONTINUE IF(N.GT.50.OR.M.GT.10)GO TO 10 WRITE(IOUNT2,717) DO 9 I=1,N WRITE(IOUNT2,706)(U(I,J),J=1,M) 9 CONTINUE 10 WRITE(IOUNT2,718)ITTOT,INACT C C 50 CONTINUE 60 CONTINUE C C 601 FORMAT(20I4) 602 FORMAT(4F20.12) 603 FORMAT(3I4,F8.5,3I4) 604 FORMAT(10F8.5) 701 FORMAT(1H1/////30X,40H** LOP-SIDED SIMULTANEOUS ITERATION **) 702 FORMAT(//34X,51HFOR PARTIAL EIGENSOLN. OF REAL UNSYMMETRIC MATRICE 1S////) 703 FORMAT(36H NUMBER OF MATRICES TO BE ANALYZED =,I2) 704 FORMAT(//////12H MATRIX NO. ,I2,12H IS OF ORDER,I4,8H AND HAS,I5, 119H NON-ZERO ELEMENTS.//) 705 FORMAT(15H THE MATRIX IS:/) 706 FORMAT(1P10E12.4) 707 FORMAT(//////26H PARAMETER VALUES FOR TEST,I2,4H ARE/) 708 FORMAT(60H M IR IS TOL LMAX ITMAX 1 LR/) 709 FORMAT(3I8,E12.4,3I8) 710 FORMAT(//48H ALL REQUIRED VECTORS SATISFY TOLERANCE CRITERIA) 711 FORMAT(//43H NOT ALL VECTORS SATISFY TOLERANCE CRITERIA) 712 FORMAT(44H RESULTS OBTAINED FROM COMPLETE QR ANALYSIS.) 713 FORMAT(36H RESULTS OBTAINED FROM POWER METHOD.) 714 FORMAT(//61H EIGENVALUE APPROXIMATIONS AND ESTIMATED EIGENVECTOR A 1CCURACY/) 715 FORMAT(15H REAL PART,16X,14HIMAGINARY PART,18X,8HACCURACY/) 716 FORMAT(F20.12,5X,F20.12,12X,F20.12) 717 FORMAT(/27H EIGENVECTOR APPROXIMATIONS/) 718 FORMAT(/34H THESE RESULTS WERE OBTAINED AFTER,I4,23H PREMULTIPLICA 1TIONS AND,I4,22H INTERACTION ANALYSES./) 719 FORMAT(/41H INITIAL APPROXIMATIONS SUPPLIED BY USER:/) C C STOP END SUBROUTINE LOPSI(N,A,IA,IROW,ICOL,LR,IR,M,IS,TOL,ITMAX,LMAX,U,IU, 00000010 1 V,W,G,H,B,P,IG,INT,X,Y,ERR,ITTOT,INACT,IFAIL) C C C C ** LOPSI IS A SIMULTANEOUS ITERATION ALGORITHM WHICH C DETERMINES APPROXIMATIONS TO RIGHT OR LEFT EIGEN- C VECTORS CORRESPONDING TO THE DOMINANT SET OF C EIGENVALUES OF A REAL UNSYMMETRIC MATRIX A. ** C C BY WILLIAM J. STEWART AND ALAN JENNINGS. C C C C FORMAL PARAMETER LIST. C ********************* C C C INPUT PARAMETERS (MUST BE SET BEFORE ENTRY AND ARE UNALTERED C **************** BY THE EXECUTION OF THE ALGORITHM) C C C N AN INTEGER QUANTITY, THE ORDER OF THE MATRIX A FOR WHICH C THE PARTIAL EIGENSOLUTION IS TO BE OBTAINED. C C A A ONE-DIMENSIONAL REAL ARRAY OF LENGTH AT LEAST EQUAL C TO IA. IT CONTAINS IN ARBITRARY ORDER THE NON-ZERO C ELEMENTS OF THE MATRIX A. C C IA AN INTEGER QUANTITY, THE NUMBER OF NON-ZERO ELEMENTS IN C THE MATRIX A. C C IROW A ONE-DIMENSIONAL INTEGER ARRAY OF LENGTH AT LEAST EQUAL TO C IA. THE I-TH COMPONENT DENOTES THE ROW POSITION OF THE C NON-ZERO ELEMENT STORED IN POSITION I OF ARRAY A. C C ICOL A ONE-DIMENSIONAL INTEGER ARRAY OF LENGTH AT LEAST EQUAL TO C IA. THE I-TH COMPONENT DENOTES THE COLUMN POSITION OF THE C NON-ZERO ELEMENT STORED IN POSITION I OF ARRAY A. C C LR AN INTEGER QUANTITY. IF LR=1 ON ENTRY, THEN THE RIGHT EIGEN- C VECTORS CORRESPONDING TO THE DOMINANT EIGENVALUES WILL BE C OBTAINED. OTHERWISE THE LEFT SET WILL BE OBTAINED. C C IR AN INTEGER QUANTITY, THE NO. OF VECTORS REQUIRED ACCURATELY. C C M AN INTEGER QUANTITY, THE NO. OF TRIAL VECTORS EMPLOYED. C M SHOULD NORMALLY BE CHOSEN SO THAT M >= IR+2 C C IS AN INTEGER QUANTITY, THE NUMBER OF TRIAL VECTORS FOR WHICH C INITIAL APPROXIMATIONS ARE SUPPLIED BY THE USER. C C TOL A REAL QUANTITY, THE TOLERANCE DEMANDED OF THE EIGENVECTORS C CORRESPONDING TO THE IR DOMINANT EIGENVALUES. C THE EIGENVALUES ARE NORMALIZED TO HAVE LARGEST COMPONENT C UNITY AND "LOPSI" AIMS FOR EACH COMPONENT TO BE IN ERROR C BY LESS THAN TOL. IN GENERAL IT CAN BE EXPECTED THAT C THE CORRESPONDING EIGENVALUES HAVE A RELATIVE ERROR C LESS THAN TOL. C C ITMAX AN INTEGER QUANTITY, AN UPPER BOUND ON THE NUMBER OF MATRIX C BY VECTOR MULTIPLICATIONS TO BE PERFORMED. C C LMAX AN INTEGER QUANTITY, THE MAXIMUM NUMBER OF PREMULTIPLIC- C ATIONS PERMITTED BETWEEN TWO REORIENTATIONS. C THE VALUE LMAX = 10 SHOULD BE ADEQUATE FOR MOST PROBLEMS. C HOWEVER SET LMAX = 1 WHEN EXPECTING FAST CONVERGENCE. C FOR FURTHER INFORMATION SEE SECTION 3 OF THE PAPER. C C IU AN INTEGER QUANTITY, THE DECLARED FIRST DIMENSION OF C ARRAYS U, V AND W. C C IG AN INTEGER QUANTITY, THE DECLARED FIRST DIMENSION C OF WORK ARRAYS G, H, B, P. C C C INPUT AND OUTPUT PARAMETER C ************************** C C C U A TWO-DIMENSIONAL REAL ARRAY OF AT LEAST (N, M) ELEMENTS. ON C ENTRY IT CONTAINS USER SUPPLIED APPROXIMATIONS TO THE EIGEN- C VECTORS CORRESPONDING TO THE IS DOMINANT EIGENVALUES. C ON EXIT IT CONTAINS THE FINAL ESTIMATES TO THE EIGENVECTORS C CORRESPONDING TO THE M DOMINANT EIGENVALUES. C C C OUTPUT PARAMETERS C ***************** C C C X, Y ONE-DIMENSIONAL REAL ARRAYS OF LENGTH AT LEAST EQUAL TO M. C ON EXIT THE REAL AND IMAGINARY PARTS OF THE J-TH EIGENVALUE C PREDICTION ARE STORED IN X(J) AND Y(J) RESPECTIVELY. C C ERR A ONE-DIMENSIONAL REAL ARRAY OF AT LEAST M ELEMENTS. ON EXIT C ERR(J) DENOTES THE ESTIMATED ACCURACY OF THE J-TH VECTOR. C ERROR ESTIMATES, AS MEASURED BY THE MAXIMUM OF THE ABSOLUTE C DIFFERENCES IN THE COMPONENTS OF SUCCESSIVE APPROXIMATIONS, C ARE GIVEN FOR ALL M VECTORS INCLUDING THOSE THAT HAVE NOT C SATISFIED THE CONVERGENCE CRITERIA SPECIFIED UNDER TOL. C C ITTOT AN INTEGER QUANTITY WHICH ON EXIT DENOTES THE C NUMBER OF PREMULTIPLICATIONS PERFORMED. C C INACT AN INTEGER QUANTITY WHICH ON EXIT DENOTES THE C NUMBER OF INTERACTION ANALYSES PERFORMED. C C IFAIL AN INTEGER QUANTITY. ON EXIT IFAIL=0 IF FIRST IR VECTORS C HAVE SATISFIED THE CONVERGENCE CRITERIA, OTHERWISE IFAIL=1. C C C WORK ARRAYS C *********** C C C V,W TWO-DIMENSIONAL REAL WORK ARRAYS OF AT LEAST (N, M) ELEMENTS C C G,H, FOUR TWO-DIMENSIONAL REAL WORK ARRAYS C B,P OF AT LEAST (M, M) ELEMENTS EACH. C C INT A ONE-DIMENSIONAL INTEGER WORK ARRAY OF AT LEAST M ELEMENTS. C C C DOUBLE PRECISION A,B,DIFF1,DIFF2,DLAM,DLAMM,DLAMR,DLAM1,EMAX DOUBLE PRECISION ERR,FMAX,G,H,P,SIGMA,SUM,TOIT,TOL,TOT,TOT1,TOT2 DOUBLE PRECISION U,V,W,X,Y,EPMACH,EPS DIMENSION A(IA),IROW(IA),ICOL(IA),U(IU,M),V(IU,M),W(IU,M) DIMENSION G(IG,M),H(IG,M),B(IG,M),P(IG,M) DIMENSION X(M),Y(M),INT(M),ERR(M) C C C C ** CALCULATE MACHINE EPSILON, EPMACH. ** C EPMACH = 1.0 1 EPMACH = 0.5*EPMACH EPS = EPMACH+1.0 IF(EPS.GT.1.0)GO TO 1 C C C ** JUMP FOR M = 1, (THE POWER METHOD), C OR M = N, (COMPLETE EIGENSOLUTION) ** C IF(M.EQ.1)GO TO 50 IF(M.EQ.N)GO TO 60 C C C EPS = EPMACH/TOL INACT = 1 ITTOT = 0 ITLOOP = LMAX/2 LOCKED = 0 LOCK = 0 C C C ** INITIALISE TRIAL VECTORS ** C C INITIAL APPROXIMATIONS TO THE FIRST IS VECTORS ARE USER C SUPPLIED. THE COMPONENTS OF THE REMAINING VECTORS ARE INIT- C IALIZED WITH PSEUDO RANDOM NUMBERS IN THE INTERVAL (-1,1), C GENERATED USING THE PROGRAM OF SCHRAGE, C (TOMS, 5, (1979), PP 132-138). C C IZ = 79311503 IF(IS.EQ.M)GO TO 4 IS1 = IS+1 DO 2 J=IS1,M DO 2 I=1,N U(I,J) = RAND(IZ)*2.0-1.0 2 CONTINUE GO TO 4 C C C ** DETERMINE NUMBER OF C PREMULTIPLICATIONS TO PERFORM ** C 100 INACT = INACT+1 C DLAM1 = DSQRT(X(1)*X(1)+Y(1)*Y(1)) DLAMR = DSQRT(X(IR)*X(IR)+Y(IR)*Y(IR)) DLAMM = DSQRT(X(M)*X(M)+Y(M)*Y(M)) L1 = 0 L2 = 0 IF(DLAMM.LE.EPMACH)GO TO 3 L1 = LMAX IF(DLAM1-DLAMM.GT.EPMACH)L1 = 5.0/DLOG10(DLAM1/DLAMM)+0.5 L2 = LMAX DLAM = DLAMR/DLAMM IF(DLAM-1.0.LE.EPMACH)GO TO 3 SIGMA = (DLAM**ITLOOP-1.0)/(DLAM-1.0) IF(ERR(IR).GT.1.0)ERR(IR) = 1.0 L2 = DLOG10(ERR(IR)/(TOL*SIGMA))/DLOG10(DLAM) 3 ITLOOP = MIN0(LMAX,L1,L2) C 4 IF(ITLOOP.LE.0)ITLOOP = 1 ITTOT = ITTOT+ITLOOP C C C ** PREMULTIPLY ITLOOP TIMES ** C LOCK1 = LOCK+1 DO 5 J=LOCK1,M DO 5 I=1,N W(I,J) = U(I,J) 5 CONTINUE CALL PREMUL(N,A,IA,IROW,ICOL,M,U,IU,V,ITLOOP,LOCK,LOCKED,LR) C C C ** MULTIPLY G = U(TRANSPOSE)*U ** C JJ = LOCKED+1 DO 8 I=1,M IF(I.GT.JJ)JJ=I DO 7 J=JJ,M SUM = 0.0 DO 6 K=1,N SUM = SUM+U(K,I)*U(K,J) 6 CONTINUE G(I,J) = SUM G(J,I) = SUM 7 CONTINUE 8 CONTINUE C C C ** MULTIPLY H = U(TRANSPOSE)*V ** C JJ = LOCKED+1 DO 11 I=1,M IF(I.EQ.JJ)JJ=1 DO 10 J=JJ,M SUM = 0.0 DO 9 K=1,N SUM = SUM+U(K,I)*V(K,J) 9 CONTINUE H(I,J) = SUM 10 CONTINUE 11 CONTINUE C C C ** DETERMINATION OF INTERACTION MATRIX ** C CALL CHOLES(M,G,H,B,IG,LOCKED+1,EPS,EPMACH) C C C ** EIGENSOLUTION OF INTERACTION MATRIX ** C USES EISPACK EIGENSYSTEM ROUTINES. C SEE LECTURE NOTES IN COMPUTER SCIENCE, NO. 6. C B.T. SMITH ET AL., SPRINGER-VERLAG, BERLIN (1976). C CALL ELMHES(IG,M,1,M,B,INT) CALL ELTRAN(IG,M,1,M,B,INT,P) CALL HQR2(IG,M,1,M,B,X,Y,P,IERR) C C C ** SORT EIGENVALUES (VECTORS) INTO C DESCENDING ORDER OF MAGNITUDE ** C CALL SORT(M,P,IG,X,Y,LOCK,B) C C C ** MULTIPLY U = V*P ** C DO 14 J=LOCK1,M DO 13 I=1,N SUM = 0.0 DO 12 K=1,M SUM = SUM+V(I,K)*P(K,J) 12 CONTINUE U(I,J) = SUM 13 CONTINUE 14 CONTINUE C DO 15 J=LOCK1,M DO 15 I=1,N V(I,J) = W(I,J) 15 CONTINUE C C C ** NORMALIZATION AND CONVERGENCE TEST ** C C J = LOCK1 LOCKED = LOCK IFAIL = 0 JIR = IR IF(IR.EQ.1)GO TO 20 IF(Y(IR).NE.0.0.AND.Y(IR).EQ.-Y(IR-1))JIR = IR-1 20 IMAG = 0 IF(Y(J).NE.0.0)IMAG = 1 CALL NORMAL(N,M,U,IU,J,IMAG,INDEX,1,IZ) IF(IFAIL.EQ.1.AND.ITTOT.LT.ITMAX.AND.J.NE.JIR)GO TO 32 TOIT = TOL*FLOAT(ITLOOP) IF(IMAG.EQ.1)GO TO 24 C C * REAL VECTOR * C EMAX = 0.0 FMAX = 0.0 DO 21 I=1,M IF(DABS(P(I,J)).LE.EMAX)GO TO 21 EMAX = DABS(P(I,J)) IMAX = I 21 CONTINUE CALL NORMAL(N,M,V,IU,IMAX,IMAG,INDEX,2,IZ) C EMAX = 0.0 DO 23 I=1,N DIFF1 = DABS(U(I,J)-V(I,IMAX)) IF(J.EQ.JIR)GO TO 22 IF(DIFF1.GT.TOIT.AND.J.LE.IR.AND.ITTOT.LT.ITMAX)GO TO 30 22 IF(DIFF1.GT.EMAX)EMAX = DIFF1 23 CONTINUE GO TO 31 C C * COMPLEX PAIR OF VECTORS * C 24 EMAX = 0.0 TOT1 = DABS(P(1,J))+DABS(P(1,J+1)) DO 25 I=2,M TOT2 = DABS(P(I,J))+DABS(P(I,J+1)) TOT = TOT1+TOT2 TOT1 = TOT2 IF(TOT.LE.EMAX)GO TO 25 EMAX = TOT IMAX = I-1 25 CONTINUE CALL NORMAL(N,M,V,IU,IMAX,IMAG,INDEX,2,IZ) C EMAX = 0.0 FMAX = 0.0 J1 = J+1 IMAX1 = IMAX+1 DO 27 I=1,N DIFF1 = DABS(U(I,J)-V(I,IMAX)) DIFF2 = DABS(U(I,J1)-V(I,IMAX1)) IF(J.EQ.JIR)GO TO 26 IF((DIFF1.GT.TOIT.OR.DIFF2.GT.TOIT).AND. 1 J.LE.IR.AND.ITTOT.LT.ITMAX)GO TO 30 26 IF(DIFF1.GT.EMAX)EMAX = DIFF1 IF(DIFF2.GT.FMAX)FMAX = DIFF2 27 CONTINUE GO TO 31 C C 30 IF(J.LE.IR)IFAIL = 1 31 IF(J.LE.JIR.AND.(EMAX.GT.TOIT.OR.FMAX.GT.TOIT))IFAIL = 1 IF(IFAIL.EQ.1.AND.ITTOT.LT.ITMAX)GO TO 32 LOCK = LOCK+1 ERR(J) = EMAX/FLOAT(ITLOOP) IF(IMAG.EQ.0)GO TO 32 LOCK = LOCK+1 ERR(J+1) = FMAX/FLOAT(ITLOOP) 32 IF(IFAIL.EQ.1.AND.J.EQ.JIR)ERR(IR) = EMAX IF(IMAG.EQ.1)J = J+1 J = J+1 IF(J.LE.M)GO TO 20 IF(IFAIL.EQ.1.AND.ITTOT.LT.ITMAX)GO TO 100 LOCK = LOCKED RETURN C C C C ** M = 1 (THE POWER METHOD) ** C C 50 ITTOT = 0 INACT = 0 IF(IS.GE.1)GO TO 52 IZ = 79311503 DO 51 I=1,N U(I,1) = RAND(IZ)*2.0-1.0 51 CONTINUE 52 LASTIN = 0 53 LASTIN = LASTIN + 1 IF(U(LASTIN,1).EQ.0.0.AND.LASTIN.LT.N)GO TO 53 IF(U(LASTIN,1).EQ.0.0)U(N,1) = 1.0 C 54 ITTOT = ITTOT+1 CALL PREMUL(N,A,IA,IROW,ICOL,1,U,IU,V,1,0,0,LR) X(1) =V(LASTIN,1)/U(LASTIN,1) CALL NORMAL(N,1,V,IU,1,0,INDEX,1,IZ) IF(DABS(1.0-U(LASTIN,1)).GT.TOL)CALL NORMAL(N,1,U,IU,1,0,INDEX,2, 1IZ) LASTIN = INDEX C IFAIL = 0 EMAX = 0.0 DO 55 I=1,N DIFF1 = DABS(U(I,1)-V(I,1)) IF(DIFF1.GT.TOL)IFAIL = 1 IF(IFAIL.EQ.1.AND.ITTOT.LT.ITMAX)GO TO 56 IF(EMAX.LT.DIFF1)EMAX = DIFF1 55 CONTINUE ERR(1) = EMAX C 56 DO 57 I=1,N U(I,1) = V(I,1) 57 CONTINUE IF(IFAIL.EQ.1.AND.ITTOT.LT.ITMAX)GO TO 54 Y(1) = 0.0 RETURN C C C ** M = N (COMPLETE EIGENSOLUTION) ** C C 60 ITTOT = 0 INACT = 0 DO 62 I=1,N DO 61 J=1,N U(I,J) = 0.0 61 CONTINUE U(I,I) = 1.0 62 CONTINUE C CALL PREMUL(N,A,IA,IROW,ICOL,N,U,IU,V,1,0,0,LR) CALL ELMHES(IU,N,1,N,V,INT) CALL ELTRAN(IU,N,1,N,V,INT,U) CALL HQR2(IU,N,1,N,V,X,Y,U,IERR) IFAIL = 1 IF(IERR.EQ.0)IFAIL = 0 DO 63 I=1,N ERR(I) = EPMACH 63 CONTINUE IF(IERR.EQ.0)RETURN DO 64 I=IERR,N ERR(I) = 1.0 64 CONTINUE RETURN C C C END SUBROUTINE PREMUL(N,A,IA,IROW,ICOL,M,U,IU,V,ITLOOP,LOCK,LOCKED,LR)00000010 C C C C ** PREMULTIPLICATION V = A*U OR V = A(TRANS)*U ** C C NOT SATISFIED THE TOLERANCE CRITERIA AND CONSEQUENTLY C HAVE NOT BEEN LOCKED. C C C FORMAL PARAMETER LIST C ********************* C C C WITH THE EXCEPTION OF ITLOOP, LOCKED, AND LOCK, ALL C OTHER PARAMETERS ARE EXACTLY AS DEFINED IN THE MAIN C SUBROUTINE, LOPSI. C C C ITLOOP AN INTEGER QUANTITY, THE NUMBER OF TIMES THE C PREMULTIPLICATION IS TO BE CARRIED OUT: C I.E. FORM V = (A**ITLOOP)*U. C C LOCKED AN INTEGER QUANTITY, THE NUMBER OF VECTORS C WHICH SATISFIED THE TOLERANCE CRITERIA PRIOR TO C THE LAST CONVERGENCE TEST. C C LOCK AN INTEGER QUANTITY, THE NUMBER OF VECTORS C WHICH SATISFIED THE TOLERANCE CRITERIA AFTER C THE LAST CONVERGENCE TEST. C C C C SPECIAL NOTE C ************ C C IF LOCKED < LOCK THEN IN ORDER FOR THE MATRIX H = U*V TO BE C CORRECTLY FORMED IN LOPSI, IT IS NECESSARY IN THE FIRST PRE- C MULTIPLICATION TO FORM V(J) = A*U(J) FOR J = LOCKED+1 C THROUGH M. THE REMAINING (ITLOOP-1) PREMULTIPLICATIONS C SHOULD BE CARRIED OUT FOR J = LOCK+1 THRO M. C C C DOUBLE PRECISION A,SCALE,AI,U,V DIMENSION A(IA),IROW(IA),ICOL(IA),U(IU,M),V(IU,M) C C C LOCK1 = LOCKED+1 C DO 100 LOOP = 1,ITLOOP IF(LR.NE.1)GO TO 5 C C C ** PREMULTIPLICATION V = AU ** C DO 1 J=LOCK1,M DO 1 I=1,N V(I,J) = 0.0 1 CONTINUE DO 3 I=1,IA IRO = IROW(I) ICO = ICOL(I) AI = A(I) DO 2 J=LOCK1,M V(IRO,J) = V(IRO,J) + AI*U(ICO,J) 2 CONTINUE 3 CONTINUE GO TO 10 C C C ** PREMULTIPLICATION V = A(TRANS)U ** C C 5 DO 6 J=LOCK1,M DO 6 I=1,N V(I,J) = 0.0 6 CONTINUE DO 8 I=1,IA IRO = IROW(I) ICO = ICOL(I) AI = A(I) DO 7 J=LOCK1,M V(ICO,J) = V(ICO,J) + AI*U(IRO,J) 7 CONTINUE 8 CONTINUE C C 10 IF(LOOP.EQ.ITLOOP)GO TO 100 IF(LOCKED.EQ.LOCK)GO TO 11 LOCK1 = LOCK+1 C C C ** SCALE MATRIX ** C 11 SCALE = 0.0 DO 12 I=1,N IF(DABS(SCALE).LT.DABS(V(I,LOCK1)))SCALE = V(I,LOCK1) 12 CONTINUE SCALE = 1.0/SCALE DO 13 J=LOCK1,M DO 13 I=1,N U(I,J) = V(I,J)*SCALE 13 CONTINUE C C 100 CONTINUE C C C RETURN END SUBROUTINE CHOLES(M,G,H,B,IG,LOCK1,EPS,EPMACH) 00000010 C C C C ** EQUATION SOLVE G*B = H BY METHOD OF CHOLESKI ** C THE THREE MATRICES G, B AND H ARE OF ORDER M. C C C FORMAL PARAMETER LIST C ********************* C C C WITH THE EXCEPTION OF LOCK1 AND EPMACH, ALL OTHER PARAMETERS C ARE EXACTLY AS DEFINED IN THE MAIN SUBROUTINE, LOPSI. C C C LOCK1 ONE MORE THAN THE NUMBER OF VECTORS WHICH SATISFIED C THE TOLERANCE CRITERIA PRIOR TO THE LAST CONVERGENCE C TEST. THE DECOMPOSITION OF G CAN THEREFORE PROCEED C FROM ROW/COLUMN LOCK1. C C EPMACH MACHINE EPSILON, THE SMALLEST NUMBER, EPSILON SUCH THAT C 1+EPSILON > 1. THIS IS CALCULATED IN THE CALLING ROUTINE. C C C SPECIAL NOTE. C ************ C C C THE DECOMPOSITION STAGE INCLUDES A MODIFICATION OF THE PIVOT C ELEMENT(S) IF G IS NEARLY SINGULAR. (SEE PAPER, SECTION 7). C IF THIS HAPPENS THEN A WARNING MESSAGE IS OUTPUT. C C IN THE DATA STATEMENT, IOUNT2 MUST BE INITIALIZED C TO THE OUTPUT STREAM NUMBER. C C DOUBLE PRECISION B,G,H,X,Y,EPS,EPMACH DIMENSION G(IG,M),H(IG,M),B(IG,M) DATA IOUNT2/6/ C C C ** DECOMPOSITION STAGE L*L(TRANS) = G ** C DO 5 I=LOCK1,M DO 4 J=1,I X = G(I,J) IF(J.EQ.1)GO TO 2 J1 = J-1 DO 1 K=1,J1 X = X-G(I,K)*G(J,K) 1 CONTINUE 2 IF(J.EQ.I)GO TO 3 G(I,J) = X/G(J,J) GO TO 4 3 Y = G(I,I)*EPS IF(Y.LT.EPMACH)Y = EPS IF(X.LT.Y)WRITE(IOUNT2,701)Y IF(X.LT.Y)X = Y G(I,I) = DSQRT(X) 4 CONTINUE 5 CONTINUE 701 FORMAT(53H DURING CHOLESKI FACTORIZATION PIVOT WAS INCREASED TO, 1E20.12) C C DO 30 J=1,M C C C ** FORWARD SUBSTITUTION L*X = H ** C DO 12 I=1,M X = H(I,J) IF(I.EQ.1)GO TO 11 I1 = I-1 DO 10 K=1,I1 X = X-G(I,K)*B(K,J) 10 CONTINUE 11 B(I,J) = X/G(I,I) 12 CONTINUE C C C ** BACK SUBSTITUTION L(TRANS)*B = X ** C DO 22 II=1,M I = M-II+1 X = B(I,J) IF(II.EQ.1)GO TO 21 I1 = I+1 DO 20 K=I1,M X = X-G(K,I)*B(K,J) 20 CONTINUE 21 B(I,J) = X/G(I,I) 22 CONTINUE C 30 CONTINUE C C RETURN END SUBROUTINE SORT(M,P,IG,X,Y,LOCK,B) 00000010 C C C C ** SORT EIGENVALUES AND CORRESPONDING EIGENVECTORS C INTO DESCENDING ORDER OF MAGNITUDE USING A STRAIGHT C INSERTION SORT TECHNIQUE. ** C C SEE D.E KNUTH. THE ART OF COMPUTER PROGRAMMING. C VOL 3, SORTING AND SEARCHING 1969, C ADDISON AND WESLEY PUBLISHERS. C C C FORMAL PARAMETER LIST. C ********************* C C C WITH THE EXCEPTION OF LOCK AND B, ALL OTHER PARAMETERS C ARE AS DECLARED IN THE MAIN SUBROUTINE, LOPSI. C NOTE THAT THE REAL AND IMAGINARY PARTS OF THE M C EIGENVALUE APPROXIMATIONS OBTAINED FROM THE INTERACTION C ANALYSIS IN LOPSI ARE STORED IN X AND Y RESPECTIVELY, C AND THAT THE CORRESPONDING EIGENVECTORS ARE STORED IN P. C C C LOCK AN INTEGER QUANTITY, THE NUMBER OF VECTORS WHICH SATISFIED C THE TOLERANCE CRITERIA AFTER THE LAST CONVERGENCE TEST. C C B A ONE-DIMENSIONAL REAL WORK ARRAY OF AT LEAST M ELEMENTS C C C DOUBLE PRECISION B,DIFF,P,RMODI,RMODJ,SX,SY,UNT,X,Y DIMENSION P(IG,M),X(M),Y(M),B(M) C C C ** TEMPORARILY INCREASE MAGNITUDE OF LOCKED EIGENVALUES ** C IF(LOCK.EQ.0)GO TO 20 IMAG = 0 DO 13 I=1,LOCK IF(IMAG.EQ.0)GO TO 10 IMAG = 0 GO TO 13 C C * DETERMINE POSITION OF "UNIT" ELEMENT IN ROW I OF P * C I-TH LOCKED VECTOR IS IN THIS POSITION C 10 UNT = 1.0 DO 11 K=1,M DIFF = DABS(1.0-DABS(P(I,K))) IF(DIFF.GE.UNT)GO TO 11 UNT = DIFF K1 = K 11 CONTINUE X(K1) = X(K1)*2.0 IF(Y(K1).EQ.0.0)GO TO 13 C C * I-TH LOCKED VECTOR IS COMPLEX * C IMAG = 1 IF(K1.EQ.1)GO TO 12 IF(Y(K1).NE.-Y(K1-1))GO TO 12 C C * VECTORS I-1 AND I FORM A COMPLEX PAIR * C X(K1-1) = X(K1) Y(K1) = Y(K1)*2.0 Y(K1-1) = -Y(K1) GO TO 13 C C * VECTORS I AND I+1 FORM A COMPLEX PAIR * C 12 X(K1+1) = X(K1) Y(K1) = Y(K1)*2.0 Y(K1+1) = -Y(K1) 13 CONTINUE C C C ** STRAIGHT INSERTION METHOD. ** C 20 DO 26 J=2,M I=J-1 SX = X(J) SY = Y(J) DO 21 K=1,M B(K) = P(K,J) 21 CONTINUE RMODJ = X(J)*X(J)+Y(J)*Y(J) 22 RMODI = X(I)*X(I)+Y(I)*Y(I) IF(RMODI.GE.RMODJ)GO TO 24 X(I+1)=X(I) Y(I+1)=Y(I) DO 23 K=1,M P(K,I+1) = P(K,I) 23 CONTINUE I=I-1 IF(I.GT.0)GOTO 22 24 X(I+1)=SX Y(I+1)=SY DO 25 K=1,M P(K,I+1) = B(K) 25 CONTINUE 26 CONTINUE C C C ** REDUCE LOCKED EIGENVALUES TO PREVIOUS VALUES ** C IF(LOCK.EQ.0)RETURN DO 30 I=1,LOCK X(I) = X(I)*0.5 Y(I) = Y(I)*0.5 30 CONTINUE C C RETURN END SUBROUTINE NORMAL(N,M,U,IU,J,IMAG,INDEX,ISTAGE,IZ) 00000010 C C C C ** NORMALIZATION OF THE UNLOCKED VECTOR, U(J) ** C C C FORMAL PARAMETER LIST C ********************* C C C WITH THE EXCEPTION OF J, IMAG, INDEX, ISTAGE AND IZ, ALL OTHER C PARAMETERS ARE EXACTLY AS DEFINED IN THE MAIN SUBROUTINE, LOPSI. C C C J AN INTEGER QUANTITY WHICH DENOTES THE COLUMN POSITION C IN ARRAY U OF THE VECTOR TO BE NORMALIZED. C C IMAG AN INTEGER QUANTITY WHICH EQUALS 1 IF THE J-TH C EIGENVALUE IS COMPLEX. C C INDEX, TWO INTEGER QUANTITIES. IF ISTAGE = 2, THE ELEMENT IN THE C ISTAGE INDEX POSITION OF THE VECTOR IS SET EQUAL TO 1 AND THE C REMAINING COMPONENTS ARE NORMALIZED ACCORDINGLY. OTHER- C WISE THE ELEMENT OF MAXIMUM MODULUS IS SET EQUAL TO 1 AND C INDEX IS ASSIGNED THE POSITION OF THIS ELEMENT. C C IZ AN INTEGER QUANTITY WHICH IS USED IN THE GENERATION C OF PSEUDO RANDOM NUMBERS. C C C DOUBLE PRECISION AB,AC,DIV,SS,S1,S2,U DIMENSION U(IU,M) C C IF(IMAG.EQ.1)GO TO 4 C C C ** REAL VECTOR ** C IF(ISTAGE.EQ.2)GO TO 2 C DIV = 0.0 DO 1 I=1,N IF(DABS(DIV).GE.DABS(U(I,J)))GO TO 1 DIV = U(I,J) INDEX = I 1 CONTINUE C 2 IF(U(INDEX,J).EQ.0.0)GO TO 10 DIV = 1.0/U(INDEX,J) DO 3 I=1,N U(I,J) = U(I,J)*DIV 3 CONTINUE RETURN C C C ** COMPLEX PAIR OF VECTORS ** C 4 J1 = J+1 IF(ISTAGE.EQ.2)GO TO 6 C DIV = 0.0 DO 5 I=1,N SS = U(I,J)*U(I,J)+U(I,J1)*U(I,J1) IF(DIV.GE.SS)GO TO 5 DIV = SS INDEX = I 5 CONTINUE C 6 AB = U(INDEX,J) AC = U(INDEX,J1) IF(AB.EQ.0.0.AND.AC.EQ.0.0)GO TO 10 DIV = 1.0/(AB*AB+AC*AC) DO 7 I=1,N S1 = U(I,J) S2 = U(I,J1) U(I,J) = (S1*AB+S2*AC)*DIV U(I,J1) = (-S1*AC+S2*AB)*DIV 7 CONTINUE RETURN C C C ** UNACCEPTABLE ZERO VECTOR. C REPLACE WITH PSEUDO RANDOM VECTOR. ** C 10 IF(ISTAGE.EQ.2)RETURN DO 11 I=1,N U(I,J) = RAND(IZ)*2.0-1.0 11 CONTINUE IF(IMAG.NE.1)RETURN DO 12 I=1,N U(I,J1) = RAND(IZ)*2.0-1.0 12 CONTINUE RETURN C C END FUNCTION RAND(IX) 00000010 C PORTABLE RANDOM NUMBER GENERATOR C USING THE RECURSION C IX = IX*A MOD P C C SOME COMPILERS, E.G., THE HP3000, REQUIRE THE FOLLOWING C DECLARATION TO BE INTEGER*4 INTEGER A,P,IX,B15,B16,XHI,XALO,LEFTLO,FHI,K C C 7**5, 2**15, 2**16, 2**31-1 DATA A/16807/,B15/32768/,B16/65536/,P/2147483647/ C C GET 15 HI ORDER BITS OF IX XHI = IX/B16 C GET 16 LO BITS OF IX AND FORM LO PRODUCT XALO=(IX-XHI*B16)*A C GET 15 HI ORDER BITS OF LO PRODUCT LEFTLO = XALO/B16 C FORM THE 31 HIGHEST BITS OF FULL PRODUCT FHI = XHI*A + LEFTLO C GET OVERFLO PAST 31ST BIT OF FULL PRODUCT K = FHI/B15 C ASSEMBLE ALL THE PARTS AND PRESUBTRACT P C THE PARENTHESES ARE ESSENTIAL IX = (((XALO-LEFTLO*B16) - P) + (FHI-K*B15)*B16) + K C ADD P BACK IN IF NECESSARY IF (IX .LT. 0) IX = IX + P C MULTIPLY BY 1/(2**31-1) RAND = FLOAT(IX)*4.656612875E-10 RETURN END 3 3 6 4 1.0 1.0 1.0 1.0 1.0 1.0 1 1 1 2 3 3 1 2 3 2 1 3 3 3 0 .0001 1 50 1 2 2 0 .0001 1 50 1 2 2 2 .0001 1 50 1 1.0 0.0 0.0 0.0 1.0 0.0 2 2 2 .0001 1 50 1 1.0 0.0 0.0 1.0 0.0 0.0 6 22 13 0.1 0.1 0.2 0.05 0.05 0.25 -0.2236 0.4472 -0.0264 0.2236 -0.2236 0.4472 0.4472 0.4472 0.4472 0.6708 0.3292 0.3292 0.3292 0.3292 0.5528 1.0 1 2 2 3 3 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6 6 1 1 2 1 2 3 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 6 6 6 0 .0001 1 50 1 6 6 0 .0001 1 50 2 1 1 0 .0001 1 50 1 1 1 0 .0001 1 50 2 1 1 1 .0001 1 50 2 1.0 1.0 1.0 1.0 1.0 1.0 4 2 0 .0001 1 50 1 4 3 0 .0001 1 50 1 4 3 3 .0001 1 50 2 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 4 3 3 .0001 1 50 2 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 5 4 0 .0001 1 50 1 5 4 0 .0001 4 50 1 5 4 0 .0001 8 50 1 5 4 0 .0001 8 50 2 1631207 1 .857142862997 .007518796516 .135338340487 .842857148711 .007518796516 .135338340487 .014285714286 .571428577283 .007518796516 .135338340487 .285714285714 .828571434425 .007518796516 .135338340487 .014285714286 .014285714286 .557142862997 .007518796516 .135338340487 .014285714286 .285714285714 .557142862997 .007518796516 .135338340487 .285714285714 .014285714286 .285714291568 .007518796516 .135338340487 .285714285714 .285714285714 .814285714286 .142857142857 .014285714286 .014285714286 .014285714286 .542857142857 .142857142857 .014285714286 .014285714286 .285714285714 .842857148711 .007518796516 .135338340487 .014285714286 .542857142857 .142857142857 .014285714286 .285714285714 .014285714286 .271428571429 .142857142857 .014285714286 .285714285714 .285714285714 .571428577283 .007518796516 .135338340487 .285714285714 .542857142857 .142857142857 .285714285714 .014285714286 .014285714286 .271428571429 .142857142857 .285714285714 .014285714286 .285714285714 .271428571429 .142857142857 .285714285714 .285714285714 .014285714286 .000000000000 .142857142857 .285714285714 .285714285714 .285714285714 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .828571434425 .007518796516 .135338340487 .014285714286 .014285714286 .828571434425 .007518796516 .135338340487 .014285714286 .014285714286 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .557142862997 .007518796516 .135338340487 .014285714286 .285714285714 .557142862997 .007518796516 .135338340487 .014285714286 .285714285714 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .557142862997 .007518796516 .135338340487 .285714285714 .014285714286 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .285714291568 .007518796516 .135338340487 .285714285714 .285714285714 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .557142862997 .007518796516 .135338340487 .285714285714 .014285714286 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .285714291568 .007518796516 .135338340487 .285714285714 .285714285714 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .842857148711 .007518796516 .135338340487 .014285714286 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .571428577283 .007518796516 .135338340487 .285714285714 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .814285716042 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .542857155736 .142857142857 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .542857155736 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .271428595430 .142857142857 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .542857155736 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .271428595430 .142857142857 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .271428595430 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .000000035124 .142857142857 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 .957142858899 .000751879652 .013533834049 .000751879652 .013533834049 .000751879652 .013533834049 .685714298593 .000751879652 .013533834049 .000751879652 .013533834049 .015037593033 .270676680974 .685714298593 .000751879652 .013533834049 .015037593033 .270676680974 .000751879652 .013533834049 .414285738287 .000751879652 .013533834049 .015037593033 .270676680974 .015037593033 .270676680974 .685714298593 .015037593033 .270676680974 .000751879652 .013533834049 .000751879652 .013533834049 .414285738287 .015037593033 .270676680974 .000751879652 .013533834049 .015037593033 .270676680974 .414285738287 .015037593033 .270676680974 .015037593033 .270676680974 .000751879652 .013533834049 .142857177981 .015037593033 .270676680974 .015037593033 .270676680974 .015037593033 .270676680974 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6 7 7 7 7 7 8 8 8 8 8 9 9 9 9 9 10 10 10 10 11 11 11 11 11 12 12 12 12 12 13 13 13 13 14 14 14 14 14 15 15 15 15 15 16 16 16 16 16 17 17 17 17 17 18 18 18 18 18 18 18 18 19 19 19 19 19 20 20 20 20 20 21 21 21 21 21 21 21 21 22 22 22 22 22 23 23 23 23 23 24 24 24 24 24 24 24 24 25 25 25 25 25 26 26 26 26 26 26 26 26 27 27 27 27 27 28 28 28 28 28 28 28 28 29 29 29 29 29 30 30 30 30 30 30 30 30 31 31 31 31 31 32 32 32 32 32 32 32 32 33 33 33 33 33 33 33 33 34 34 34 34 34 34 34 34 35 35 35 35 36 36 36 36 36 36 36 36 37 37 37 37 38 38 38 38 38 38 38 38 39 39 39 39 39 39 39 39 40 40 40 40 40 40 40 40 41 41 41 41 41 41 41 41 42 42 42 42 42 42 42 42 43 43 43 43 43 43 43 43 44 44 44 44 44 44 44 44 45 45 45 45 45 45 45 45 46 46 46 46 46 46 46 46 47 47 47 47 47 47 47 47 48 48 48 48 48 48 48 48 49 49 49 49 49 49 49 49 50 50 50 50 50 50 50 50 51 51 51 51 51 51 51 51 52 52 52 52 52 52 52 52 53 53 53 53 53 53 53 53 54 54 54 54 54 54 54 54 55 55 55 55 55 55 55 55 56 56 56 56 56 56 56 56 57 57 57 57 57 57 57 57 58 58 58 58 58 58 58 58 59 59 59 59 59 59 59 59 60 60 60 60 60 60 60 60 61 61 61 61 61 61 61 61 62 62 62 62 62 62 62 62 63 63 63 63 63 63 63 63 64 64 64 64 64 64 64 64 65 65 65 65 65 65 65 65 66 66 66 66 66 66 66 66 67 67 67 67 67 67 67 67 68 68 68 68 68 68 68 68 69 69 69 69 69 69 69 69 70 70 70 70 70 70 70 70 71 71 71 71 71 71 71 71 72 72 72 72 72 72 72 72 73 73 73 73 73 73 73 73 74 74 74 74 74 74 74 74 75 75 75 75 75 75 75 75 76 76 76 76 76 76 76 76 77 77 77 77 77 77 77 77 78 78 78 78 78 78 78 78 79 79 79 79 79 79 79 79 80 80 80 80 80 80 80 80 81 81 81 81 81 81 81 81 82 82 82 82 82 82 82 82 83 83 83 83 83 83 83 83 84 84 84 84 84 84 84 84 85 85 85 85 85 85 85 85 86 86 86 86 86 86 86 86 87 87 87 87 87 87 87 87 88 88 88 88 88 88 88 88 89 89 89 89 89 89 89 89 90 90 90 90 90 90 90 90 91 91 91 91 91 91 91 91 92 92 92 92 92 92 92 92 93 93 93 93 93 93 93 93 94 94 94 94 94 94 94 94 95 95 95 95 95 95 95 95 96 96 96 96 96 96 96 96 97 97 97 97 97 97 97 97 98 98 98 98 98 98 98 98 99 99 99 99 99 99 99 99 100 100 100 100 100 100 100 100 101 101 101 101 101 101 101 101 102 102 102 102 102 102 102 102 103 103 103 103 103 103 103 103 104 104 104 104 104 104 104 104 105 105 105 105 105 105 105 105 106 106 106 106 106 106 106 106 107 107 107 107 107 107 107 107 108 108 108 108 108 108 108 108 109 109 109 109 109 109 109 109 110 110 110 110 110 110 110 110 111 111 111 111 111 111 111 111 112 112 112 112 112 112 112 112 113 113 113 113 113 113 113 113 114 114 114 114 114 114 114 114 115 115 115 115 115 115 115 115 116 116 116 116 116 116 116 116 117 117 117 117 117 117 117 117 118 118 118 118 118 118 118 118 119 119 119 119 119 119 119 119 120 120 120 120 120 120 120 120 121 121 121 121 121 121 121 121 122 122 122 122 122 122 122 122 123 123 123 123 123 123 123 123 124 124 124 124 124 124 124 124 125 125 125 125 125 125 125 125 126 126 126 126 126 126 126 126 127 127 127 127 127 127 127 127 128 128 128 128 128 128 128 128 129 129 129 129 129 129 129 129 130 130 130 130 130 130 130 130 131 131 131 131 131 131 131 131 132 132 132 132 132 132 132 132 133 133 133 133 133 133 133 133 134 134 134 134 134 134 134 134 135 135 135 135 135 135 135 135 136 136 136 136 136 136 136 136 137 137 137 137 137 137 137 137 138 138 138 138 138 138 138 138 139 139 139 139 139 139 139 139 140 140 140 140 140 140 140 140 141 141 141 141 141 141 141 141 142 142 142 142 142 142 142 142 143 143 143 143 143 143 143 143 144 144 144 144 144 144 144 144 145 145 145 145 145 145 145 145 146 146 146 146 146 146 146 146 147 147 147 147 147 147 147 147 148 148 148 148 148 148 148 148 149 149 149 149 149 149 149 149 150 150 150 150 150 150 150 150 151 151 151 151 151 151 151 151 152 152 152 152 152 152 152 152 153 153 153 153 153 153 153 153 154 154 154 154 154 154 154 154 155 155 155 155 155 155 155 155 156 156 156 156 156 156 156 157 157 157 157 157 157 157 158 158 158 158 158 158 158 159 159 159 159 159 159 159 160 160 160 160 160 160 160 161 161 161 161 161 161 161 162 162 162 162 162 162 162 163 163 163 163 163 163 163 1 2 3 2 4 5 1 3 6 7 1 4 8 9 10 2 5 11 12 13 2 6 14 15 10 3 7 16 17 13 3 8 18 19 20 4 9 21 22 23 4 10 4 6 1 11 24 25 20 5 12 26 27 23 5 13 5 7 1 14 28 19 29 6 15 30 22 31 6 16 32 25 29 7 17 33 27 31 7 18 34 8 14 8 11 8 9 19 8 14 35 10 20 8 11 35 2 21 36 9 15 9 12 8 9 22 9 15 37 10 23 9 12 37 2 24 38 11 16 8 11 11 12 25 11 16 35 13 26 39 12 17 9 12 11 12 27 12 17 37 13 28 40 8 14 14 16 14 15 29 14 16 35 3 30 41 9 15 15 17 14 15 31 15 17 37 3 32 42 11 16 14 16 16 17 33 43 12 17 15 17 16 17 34 44 18 28 18 24 18 21 35 20 29 1 36 45 21 30 21 26 18 21 37 23 31 1 38 46 24 32 18 24 24 26 39 47 26 33 21 26 24 26 40 48 18 28 28 32 28 30 41 49 21 30 30 33 28 30 42 50 24 32 28 32 32 33 43 51 26 33 30 33 32 33 44 52 34 40 34 38 34 36 45 53 36 41 36 39 34 36 46 54 38 42 34 38 38 39 47 55 39 43 36 39 38 39 48 56 34 40 40 42 40 41 49 57 36 41 41 43 40 41 50 58 38 42 40 42 42 43 51 59 39 43 41 43 42 43 52 60 44 48 44 46 44 45 53 61 45 49 45 47 44 45 54 62 46 50 44 46 46 47 55 63 47 51 45 47 46 47 56 64 44 48 48 50 48 49 57 65 45 49 49 51 48 49 58 66 46 50 48 50 50 51 59 67 47 51 49 51 50 51 60 68 52 56 52 54 52 53 61 69 53 57 53 55 52 53 62 70 54 58 52 54 54 55 63 71 55 59 53 55 54 55 64 72 52 56 56 58 56 57 65 73 53 57 57 59 56 57 66 74 54 58 56 58 58 59 67 75 55 59 57 59 58 59 68 76 60 64 60 62 60 61 69 77 61 65 61 63 60 61 70 78 62 66 60 62 62 63 71 79 63 67 61 63 62 63 72 80 60 64 64 66 64 65 73 81 61 65 65 67 64 65 74 82 62 66 64 66 66 67 75 83 63 67 65 67 66 67 76 84 68 72 68 70 68 69 77 85 69 73 69 71 68 69 78 86 70 74 68 70 70 71 79 87 71 75 69 71 70 71 80 88 68 72 72 74 72 73 81 89 69 73 73 75 72 73 82 90 70 74 72 74 74 75 83 91 71 75 73 75 74 75 84 92 76 80 76 78 76 77 85 93 77 81 77 79 76 77 86 94 78 82 76 78 78 79 87 95 79 83 77 79 78 79 88 96 76 80 80 82 80 81 89 97 77 81 81 83 80 81 90 98 78 82 80 82 82 83 91 99 79 83 81 83 82 83 92 100 84 88 84 86 84 85 93 101 85 89 85 87 84 85 94 102 86 90 84 86 86 87 95 103 87 91 85 87 86 87 96 104 84 88 88 90 88 89 97 105 85 89 89 91 88 89 98 106 86 90 88 90 90 91 99 107 87 91 89 91 90 91 100 108 92 96 92 94 92 93 101 109 93 97 93 95 92 93 102 110 94 98 92 94 94 95 103 111 95 99 93 95 94 95 104 112 92 96 96 98 96 97 105 113 93 97 97 99 96 97 106 114 94 98 96 98 98 99 107 115 95 99 97 99 98 99 108 116 100 104 100 102 100 101 109 117 101 105 101 103 100 101 110 118 102 106 100 102 102 103 111 119 103 107 101 103 102 103 112 120 100 104 104 106 104 105 113 121 101 105 105 107 104 105 114 122 102 106 104 106 106 107 115 123 103 107 105 107 106 107 116 124 108 112 108 110 108 109 117 125 109 113 109 111 108 109 118 126 110 114 108 110 110 111 119 127 111 115 109 111 110 111 120 128 108 112 112 114 112 113 121 129 109 113 113 115 112 113 122 130 110 114 112 114 114 115 123 131 111 115 113 115 114 115 124 132 116 120 116 118 116 117 125 133 117 121 117 119 116 117 126 134 118 122 116 118 118 119 127 135 119 123 117 119 118 119 128 136 116 120 120 122 120 121 129 137 117 121 121 123 120 121 130 138 118 122 120 122 122 123 131 139 119 123 121 123 122 123 132 140 124 128 124 126 124 125 133 141 125 129 125 127 124 125 134 142 126 130 124 126 126 127 135 143 127 131 125 127 126 127 136 144 124 128 128 130 128 129 137 145 125 129 129 131 128 129 138 146 126 130 128 130 130 131 139 147 127 131 129 131 130 131 140 148 132 136 132 134 132 133 141 149 133 137 133 135 132 133 142 150 134 138 132 134 134 135 143 151 135 139 133 135 134 135 144 152 132 136 136 138 136 137 145 153 133 137 137 139 136 137 146 154 134 138 136 138 138 139 147 155 135 139 137 139 138 139 148 156 140 144 140 142 140 141 149 157 141 145 141 143 140 141 150 158 142 146 140 142 142 143 151 159 143 147 141 143 142 143 152 160 140 144 144 146 144 145 153 161 141 145 145 147 144 145 154 162 142 146 144 146 146 147 155 163 143 147 145 147 146 147 156 148 152 148 150 148 149 157 149 153 149 151 148 149 158 150 154 148 150 150 151 159 151 155 149 151 150 151 160 148 152 152 154 152 153 161 149 153 153 155 152 153 162 150 154 152 154 154 155 163 151 155 153 155 154 155 5 3 1 .0001 10 500 1 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0