SUBROUTINE MINVAL(N,Q,PINIT,R,MMAX,EPS,OP,M,D,X,IECODE) IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,PINIT,R,MMAX REAL*8 EPS EXTERNAL OP DIMENSION D(Q),X(N,Q) INTEGER IECODE DIMENSION E(25),C(25,25) DIMENSION U(1296),V(1296) INTEGER P,S,PS IF (N.LT.2) GO TO 901 IF (N.GT.1296) GO TO 902 IF (R.LT.1) GO TO 903 IF (Q.LE.R) GO TO 904 IF (Q.GT.25) GO TO 905 IF (Q.GT.N) GO TO 906 IMM = 0 P = PINIT IF (P.LT.0) P = -P S = (Q-M)/P IF (S.GT.2) GO TO 100 S = 2 P = Q/2 100 IF (PINIT.LT.0) GO TO 150 DO 120 K = 1,P CALL RANDOM(N,Q,K,X) 120 CONTINUE 150 IF (M.GT.0) GO TO 200 CALL ORTHG(N,Q,M,P,C,X) CALL SECTN(N,Q,M,P,OP,X,C,D,U,V) IMM = IMM+P ERRC = 0.D0 200 ITER = 0 300 IF (M.GE.R) GO TO 900 IF (ITER.GT.MMAX) GO TO 907 ITER = ITER + 1 PS = P*S IMM = IMM + PS CALL BKLANC(N,Q,M,P,S,OP,D,C,X,E,U,V) CALL EIGEN(Q,M,P,PS,C,D) CALL CNVTST(N,Q,M,P,ERRC,EPS,D,E,NCONV) CDBG CALL PCH(N,Q,M,R,NCONV,P,S) CALL ROTATE(N,Q,M,PS,NCONV+P,C,X) CDBG M = M + NCONV PRINT 1001,ITER,IMM,P,PS 1001 FORMAT(' =>ITER,IMM,P,PS =',4I5) GO TO 300 900 IECODE = 0 RETURN 901 IECODE = 1 RETURN 902 IECODE = 2 RETURN 903 IECODE = 3 RETURN 904 IECODE = 4 RETURN 905 IECODE = 5 RETURN 906 IECODE = 6 RETURN 907 IECODE = 7 PINIT = -P RETURN END SUBROUTINE BKLANC(N,Q,M,P,S,OP,D,C,X,E,U,V) IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,M,P,S DIMENSION D(Q),C(Q,Q),X(N,Q) DIMENSION E(Q),U(N),V(N) MP1 = M+1 MPPS = M+PS DO 90 L = 1,S LL = M+(L-1)*P+1 LU = M + L*P DO 70 K = LL,LU DO 10 I = 1,N 10 U(I) = X(I,K) CALL OP(N,U,V) IF (L.GT.1) GO TO 19 DO 12 I=K,LU 12 C(I,K) = 0.D0 C(K,K) = D(K) DO 14 I = 1,N 14 V(I) = V(I) - D(K)*X(I,K) GO TO 61 19 DO 30 I = K,LU T = 0 DO 20 J = 1,N 20 T = T+V(J)*X(J,I) 30 C(I,K) = T IT = K-P DO 60 I = 1,N T = 0 DO 40 J = IT,K 40 T = T + X(I,J)*C(K,J) IF (K.EQ.LU) GO TO 60 KP1 = K+1 DO 50 J = KP1,LU 50 T = T+X(I,J)*C(J,K) 60 V(I) = V(I)-T 61 IF (L.EQ.S) GO TO 70 DO 63 I = 1,N 63 X(I,K+P) = V(I) 70 CONTINUE IF (L.EQ.1) CALL ERR(N,Q,M,P,X,E) IF (L.EQ.S) GO TO 90 CALL ORTHG(N,Q,LU,P,C,X) IL = LU+1 IT = LU DO 80 J = 1,P IT = IT+1 DO 80 I = IL,IT 80 C(I,IT-P) = C(I,IT) 90 CONTINUE RETURN END SUBROUTINE PCH(N,Q,M,R,NCONV,P,S) IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,M,R,NCONV,P,S INTEGER PT,ST MT = M+NCONV PT = R-MT IF (PT.GT.P) PT = P IF (PT.GT.0) GO TO 101 P = 0 RETURN 101 CONTINUE ST = (Q-MT)/PT IF (ST.GT.2) GO TO 110 ST = 2 PT = (Q-MT)/2 110 P = PT S = ST RETURN END SUBROUTINE ERR(N,Q,M,P,X,E) IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,M,P DIMENSION X(N,Q),E(Q) MP1 = M+P+1 MPP = M+P+P DO 200 K = MP1,MPP T = 0.D0 DO 100 I = 1,N 100 T = T+X(I,K)**2 200 E(K-P) = DSQRT(T) RETURN END SUBROUTINE CNVTST(N,Q,M,P,ERRC,EPS,D,E,NCONV) IMPLICIT REAL*8 (A-H,O-Z) INTEGER Q,M,P REAL*8 EPS DIMENSION D(Q),E(Q) REAL*8 MCHEPS/2.22D-16/ K = 0 DO 100 I = 1,P T = DABS(D(M+I)) IF (T.LT.1.D0) T = 1.D0 IF (E(M+I).GT.T*(EPS+10.D0*N*MCHEPS)+ERRC) GO TO 200 100 K = I 200 NCONV = K IF (K.EQ.0) RETURN T = 0.D0 DO 300 I = 1,K T = T+E(M+I)**2 300 CONTINUE ERRC = DSQRT(ERRC**2+T) RETURN END SUBROUTINE EIGEN(Q,M,P,PS,C,D) IMPLICIT REAL*8 (A-H,O-Z) INTEGER M,P,Q,PS,PP1 DIMENSION C(Q,Q),D(Q) DIMENSION DD(25),V(25) DO 150 I=1,PS LIM = I-P IF (I.LE.P) LIM = 1 IF (LIM.LE.1) GO TO 130 LM1 = LIM-1 DO 120 J = 1,LM1 120 C(I,J) = 0.D0 130 DO 140 J = LIM,I 140 C(I,J) = C(I+M,J+M) 150 CONTINUE CALL TRED2(Q,PS,C,DD,V,C) CALL TQL2(Q,PS,DD,V,C,IERR) PRINT 601,PS,(DD(I),I=1,PS) 601 FORMAT(' => ORDER =',I4,/,(' EIGENVALUES =',10F10.5)) PP1 = P+1 DO 155 J = 1,PP1 155 PRINT 602,J,(C(I,J),I=1,PS) 602 FORMAT(' => J =',I4,/,(' EIGENVECTORS =',10F10.5)) DO 160 I = 1,PS 160 D(M+I) = DD(I) RETURN END SUBROUTINE SECTN(N,Q,M,P,OP,X,C,D,U,V) IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,M,P EXTERNAL OP DIMENSION X(N,Q),C(Q,Q),D(Q),U(N),V(N) ICOL1 = M DO 300 J = 1,P ICOL1 = ICOL1 + 1 DO 100 I = 1,N 100 U(I) = X(I,ICOL1) CALL OP(N,U,V) ICOL2 = M DO 300 I = 1,J ICOL2 = ICOL2+1 T = 0.D0 DO 200 K = 1,N 200 T = T + V(K)*X(K,ICOL2) 300 C(ICOL1,ICOL2) = T CALL EIGEN(Q,M,P,P,C,D) CALL ROTATE(N,Q,M,P,P,C,X) RETURN END SUBROUTINE ROTATE(N,Q,M,PS,L,C,X) IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,M,PS,L DIMENSION C(Q,Q),X(N,Q) DIMENSION V(25) DO 300 I = 1,N DO 200 K = 1,L T = 0.D0 DO 100 J = 1,PS 100 T = T+X(I,M+J)*C(J,K) 200 V(K) = T DO 300 K = 1,L 300 X(I,M+K) = V(K) RETURN END SUBROUTINE ORTHG(N,Q,F,P,B,X) IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,F,P DIMENSION B(Q,Q),X(N,Q) INTEGER FP1,FPP LOGICAL ORIG IF (P.EQ.0) RETURN FP1 = F+1 FPP = F+P DO 50 K = FP1,FPP ORIG = .TRUE. KM1 = K-1 10 T = 0.D0 IF (KM1.LT.1) GO TO 25 DO 20 I = 1,KM1 S = 0.D0 DO 15 J = 1,N 15 S = S+X(J,I)*X(J,K) IF (ORIG.AND.I.GT.F) B(I,K) = S T = T + S*S DO 20 J = 1,N 20 X(J,K) = X(J,K) - S*X(J,I) 25 S = 0.D0 DO 30 J = 1,N 30 S = S + X(J,K)*X(J,K) T = T+S IF (S.GT.T/100) GO TO 40 ORIG = .FALSE. GO TO 10 40 S = DSQRT(S) B(K,K) = S IF (S.NE.0) S = 1/S DO 50 J = 1,N 50 X(J,K) = S*X(J,K) RETURN END SUBROUTINE RANDOM(N,Q,L,X) IMPLICIT REAL*8 (A-H,O-Z) INTEGER N,Q,L DIMENSION X(N,Q) DIMENSION T(100) INTEGER F1/71416/,F2/27183/,FT INTEGER A/6821/,C/5327/,X0/5328/,X1 DO 100 I = 1,100 X1 = A *X0+C IF (X1.GE.10000) X1 = X1 - 10000 T(I) = X1/9999.D0 - 0.5D0 100 X0 = X1 DO 200 I = 1,N FT = F1+F2 IF (FT.GE.1000000) FT = FT-1000000 F1 = F2 F2 = FT K = FT/1.D6*100 + 1 X(I,L) = T(K) X1 = A*X0 + C IF (X1.GE.10000) X1 = X1 - 10000 T(K) = X1/9999D0 - 0.5D0 200 X0 = X1 RETURN END SUBROUTINE AX(N,U,V) IMPLICIT REAL*8(A-H,O-Z) DIMENSION U(N),V(N),TMP(040) COMMON/MATRIX/MFLAG,NPOWER,A(040,0030) C DO 3 IROW=1,N 3 V(IROW) = U(IROW) C DO 1 IPOWER=1,NPOWER DO 2 IROW=1,N 2 TMP(IROW) = V(IROW) CALL AXSUB(N,TMP,V) 1 CONTINUE C RETURN END SUBROUTINE AXSUB(N,U,V) IMPLICIT REAL*8(A-H,O-Z) DIMENSION U(N),V(N),TMP(040) COMMON/MATRIX/MFLAG,NPOWER,A(040,0030) C C COMPUTES V = A * U. IF MFLAG .GT. 0 THEN EACH COL OF C ARRAY A IS THE DEFINING VECTOR FOR A HH TRANFORMATION C (FIRST MFLAG COLS ONLY). THE COLS ARE ASSUMED TO BE C ALREADY OF UNIT LENGTH. C IF (MFLAG .NE. 0) GOTO 1 C C COMPUTE A * U NORMALLY DO 2 IROW=1,N SUM = 0D0 DO 3 ISUM=1,N 3 SUM = SUM + A(IROW,ISUM) * U(ISUM) V(IROW) = SUM 2 CONTINUE GOTO 4 C 1 IF (MFLAG .LT. 0) GOTO 9 C C COMPUTE A * U WHERE A IS REPRESENTED BY HH X-FORMATIONS. C I - 2 * A(*,ICOL) * A(*,ICOL) (T), FOR ICOL=1,...,MFLAG C DO 7 IROW=1,N V(IROW) = U(IROW) TMP(IROW) = U(IROW) 7 CONTINUE C DO 5 ICOL=1,MFLAG C SUM = 0D0 DO 6 ISUM=1,N 6 SUM = SUM + A(ISUM,ICOL) * TMP(ISUM) C DO 8 IROW=1,N 8 V(IROW) = V(IROW) - 2 * SUM * A(IROW,ICOL) C DO 14 IROW=1,N 14 TMP(IROW) = V(IROW) C 5 CONTINUE C GOTO 4 C 9 IF (MFLAG .LT. -2) GOTO 21 C C COMPUTE A * U WHERE A IS A SPECIFIC ILL-COND'ED MATRIX C DO 10 ICOL=1,N 10 A(ICOL,1) = 0D0 C DO 11 IROW=1,N C A(IROW,1) = DFLOAT(IROW) IF (IROW .GE. N) GOTO 13 A(IROW+1,1) = DFLOAT(IROW) 13 CONTINUE C SUM = 0D0 DO 12 ISUM=1,N 12 SUM = SUM + A(ISUM,1) * U(ISUM) V(IROW) = SUM C 11 CONTINUE C IF (MFLAG .NE. -2) GOTO 15 C C COMPUTE TWICE THE SYMMETRIC PART OF 'MFLAG=-1' MATRIX... C DO 16 IROW=1,N C IF (IROW .LT. 3) GOTO 17 DO 18 ICOL=3,IROW 18 A(ICOL-2,1) = 0D0 17 CONTINUE IF (IROW .GE. 2) A(IROW-1,1) = DFLOAT(IROW-1) DO 19 ICOL=IROW,N 19 A(ICOL,1) = DFLOAT(IROW) C SUM = 0D0 DO 20 ISUM=1,N 20 SUM = SUM + A(ISUM,1) * U(ISUM) V(IROW) = V(IROW) + SUM 16 CONTINUE C 15 GOTO 4 C 21 CONTINUE C V(1) = 2D0 * U(1) - 1D0 * U(2) C ITMP = N-1 DO 22 IROW=2,ITMP 22 V(IROW) = 2D0*U(IROW) - 1D0*U(IROW-1) - 1D0*U(IROW+1) C V(N) = 2D0*U(N) - 1D0*U(N-1) C 4 CONTINUE C RETURN END IMPLICIT REAL*8 (A-H,O-Z) COMMON/MATRIX/MFLAG,NPOWER,A(040,0030) DIMENSION D(25),X(6000) EXTERNAL AX INTEGER Q,PINIT,R DIMENSION S(32) PI = 3.14159265358979D0 30 READ *,N,PINIT,Q,R,MFLAG,NPOWER PRINT 101,N,PINIT,Q,R,MFLAG,NPOWER 101 FORMAT('1N,PINIT,Q,R,MFLAG,NPOWER:',6I5) IF (N .EQ. 0) STOP DO 10 I = 1,32 10 S(I) = DSIN(I*PI/33)**2 DO 20 I = 1,040 DO 20 J = 1,0030 20 A(I,J) = 0D0 MMAX = 30 EPS = 1.D-14 M = 0 CALL MINVAL(N,Q,PINIT,R,MMAX,EPS,AX,M,D,X,IECODE) PRINT 1001,M,IECODE,(D(I),I=1,M) 1001 FORMAT(/' =>M,IECODE=',I4,1X,I4/(' =>E ',5D23.15)) GOTO 30 END //GO.SYSIN DD * DATA INPUT SET 30 4 12 1 -2 1 30 4 12 1 -2 2 0 0 0 0 0 0 /*