PROGRAM EIG_L
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      INTEGER K,LK
      DIMENSION U(33000),U1(33000)
     *,S1(33000),S2(67000),IT1(33000),IT2(67000)
      COMMON /AB/N,H,PI,RO,H2,HPI,N1,N2,RR1
      COMMON /C/D1Y(33000),D2Y(67000)
      OPEN(1,FILE='RES_L')
      WRITE(*,*)'EIGENVALUES OF THE PROBLEM -Lu=$Mu IN THE L-SHAPED '
      WRITE(*,*)' REGION WHERE  M  IS PIECE-WISE CONSTANT FUNCTION'
      WRITE(*,*)'Skorokhodov&Knyazev knyazev@na-net.ornl.gov 1990'
      WRITE(1,*)'EIGENVALUES OF THE PROBLEM -Lu=$Mu IN THE L-SHAPED '
      WRITE(1,*)' REGION WHERE "M" IS PIECE-WISE CONSTANT FUNCTION'
      WRITE(1,*)'Skorokhodov&Knyazev knyazev@na-net.ornl.gov 1990'
      WRITE(*,*)'PLEASE ENTER "n",WHERE 2**n-1 IS THE '
     *,'NUMBER OF POINTS ON THE INTERFACE'
      WRITE(*,*)'!!!!! REMEMBER THAT 2<n<16 !!!!!'
      READ(*,*)N
      WRITE(*,*)'n=',N
      WRITE(*,*)'PLEASE ENTER THE COEFFICIENT "M1" WHERE  M1 IS '
      WRITE(*,*)'THE VALUE OF THE FUNCTION  "M" IN THE QUADRAT'
      WRITE(*,*)' (SMALLER PART OF THE REGION)'
     *,'NUMBER OF POINTS ON THE INTERFACE'
      WRITE(1,*)'n=',N
      WRITE(1,*)' '
      WRITE(1,*)'PLEASE ENTER THE COEFFICIENT "M1" WHERE  M1 IS '
      WRITE(1,*)'THE VALUE OF THE FUNCTION  "M" IN THE QUADRAT'
      WRITE(1,*)' (SMALLER PART OF THE REGION)'
      READ(*,*)RR1
      WRITE(1,*)' '
      WRITE(1,*)'M1=',RR1
      WRITE(*,*)'M1=',RR1
      WRITE(1,*)' '
      WRITE(1,*)'THE LANCOSZ METHOD IS IMPLEMENTED'
      PI=4.D+00*DATAN(1.D+00)
      N1=2**N-1
      HH1=1.D0/N1
      N=N1+1
      H=1.D+00/N
      HPI=PI*H
      H2=H*H
      RO=H2/6.
      N2=2*N-1
      RN3=3.D00*N1*N1+2.D00*N1
      WRITE(1,*)' '
      WRITE(1,*)'DIMENSION OF THE SUBSPACE IS',N1
      WRITE(1,*)' '
      WRITE(1,*)'DIMENSION OF THE WHOLE PROBLEM IS ',RN3
      WRITE(1,*)' '
      CALL INGS(U)
      CALL CRCOM(R)
      CALL TRANS(N,S1,IT1)
      CALL TRANS(2*N,S2,IT2)
      CALL CRRAY(RAY0,U,S1,IT1,S2,IT2)
      WRITE(1,*)'LAMBDA DESH=',R
      WRITE(1,*)' '
      WRITE(1,*)'INITIAL RAY=',RAY0
      WRITE(1,*)'   '
      LK=2
      DO K=1,6
      WRITE(*,*)'             ON ITERATION NUMBER ',K
*      WRITE(1,*)'             ON ITERATION NUMBER ',K
      CALL LANC(U,RAY0,U1,LK,S1,IT1,S2,IT2)
      S=SCH(U1,RAY0,S1,IT1,S2,IT2)
      R=SM(U1,RAY0,S1,IT1,S2,IT2)
      RAY1=RAY0+S/R
      DO I=1,N1
      U(I)=U1(I)
      END DO
      RAY0=RAY1
      WRITE(1,*)'        '
      WRITE(*,*)'             THE RAYLEIGH QUOTIENT IS'
      WRITE(*,*)'                ',RAY0
      WRITE(*,*)'        '
      END DO
      RY=4.D+00*DSIN(PI*.25*HH1)**2/HH1**2
      RY=(D1Y(1)+RY)/RR1
      WRITE(1,*)'             THE RAYLEIGH QUOTIENT IS'
      WRITE(1,*)'                ',RAY0
*      WRITE(1,*)'EIGV WITH NEUMAN CONDITIONS=',RY
      WRITE(1,*)'   '
      CLOSE(1)
      STOP
      END
*
      SUBROUTINE CRCOM(RAY)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      COMMON /B/DS(33000)
      COMMON /C/D1Y(33000),D2Y(67000)
      COMMON /AB/N,H,PI,RO,H2,HPI,N1,N2,RR1
      DO K=1,N1
      R=4.D+00*DSIN(K*HPI*.5)**2
      RR=R/H2
      D1Y(K)=RR
      D2Y(2*K)=RR
      DS(K)=DSQRT(R+R*R*.25)/H2
      END DO
      DO K=1,N
      R=4.D+00/H2*DSIN((2*K-1)*HPI*.25)**2
      D2Y(2*K-1)=R
      END DO
      RAY=D1Y(1)+D2Y(1)
      RAY1=(D1Y(1)+D1Y(1))/(RR1+1.D-111)
      IF(DABS(RAY1).LT.RAY)RAY=RAY1
      RETURN
      END
*
      SUBROUTINE INGS(U)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      DIMENSION U(33000)
      COMMON /AB/N,H,PI,RO,H2,HPI,N1,N2,RR1
      B=.5D00
      A=DSQRT(2.D+00/N)
      DO I=1,N1
      U(I)=A*DSIN(HPI*I)/(N-I)**(1./3.)
*     U(I)=A*DSIN(2.*HPI*I)
*      WRITE(1,*)U(I)
      END DO
      RETURN
      END
*
      SUBROUTINE CRRAY(RAY,U,S1,IT1,S2,IT2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      INTEGER IERR
      DIMENSION U(33000),U1(33000),EV(3),A(3,3),RAB(6)
     *,S1(33000),S2(67000),IT1(33000),IT2(67000)
      COMMON /AB/N,H,PI,RO,H2,HPI,N1,N2,RR1
      COMMON /C/D1Y(33000),D2Y(67000)
      RL1=D1Y(1)+D2Y(1)
      RL2=(D1Y(1)+D1Y(1))
      C=0.
      RL3=SCH(U,C,S1,IT1,S2,IT2)
      R1=0.D00
      R2=0.D00
      AA=DSQRT(1.D00/N)
      B=DSQRT(2.D00/N)
      SS2=B*DSIN(HPI)
      DO I=1,N1
      R1=R1+U(I)*SS2*AA*DSIN(HPI*.5*I)
      R2=R2+U(I)*SS2*B*DSIN(HPI*I)
      END DO
      RM13=R1/RL1/H2
      RM23=R2*RR1/RL2/H2
      RM3=SM(U,C,S1,IT1,S2,IT2)
      WRITE(*,*)'RL3=',RL3,'   RM3=',RM3
      A(1,1)=1.D00/RL1
      A(1,2)=0.D00
      A(2,1)=0.D00
      A(2,2)=RR1/RL2
      A(1,3)=RM13/DSQRT(RL3*RL1)
      A(2,3)=RM23/DSQRT(RL3*RL2)
      A(3,1)=A(1,3)
      A(3,2)=A(2,3)
      A(3,3)=RM3/RL3
      CALL AEH2D(3,A,EV,RAB,IERR)
      RAY=1.D00/EV(3)
*     WRITE(1,*)EV
      RETURN
      END
*
      SUBROUTINE SCHUR(U,RAY,SU,S1,IT1,S2,IT2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      DIMENSION U(33000),SU(33000)
     *,S1(33000),S2(67000),IT1(33000),IT2(67000),U2(67000),U1(33000)
      COMMON /C/D1Y(33000),D2Y(67000)
      COMMON /AB/N,H,PI,RO,H2,HPI,N1,N2,RR1
      CALL TRIDIAG(U,SU,-1.D+00,4.D+00)
      DO I=1,N1
      SU(I)=SU(I)/H2-RAY*U(I)*(1.D0+RR1)/2
      END DO
      DO I=1,N1
      U2(I)=U(I)
      U1(I)=U(I)
      U2(I+N)=0.D+00
      END DO
      U2(N)=0.D+00
      CALL FFT(2*N,U2,S2,IT2)
      CALL FFT(N,U1,S1,IT1)
      DO I=1,N1
      A=RAY*RR1-D1Y(I)
                  IF(A.GT.0.D00)THEN
      ALF=2.*DASIN(.5*H*DSQRT(A))
      G=DCOS(ALF)-DSIN(ALF)/DTAN(ALF*N)
      U1(I)=U1(I)*G
                   ELSE
      ALF=2.*DASINH(.5*H*DSQRT(-A))
      G=DCOSH(ALF)-DSINH(ALF)/DTANH(ALF*N)
      U1(I)=U1(I)*G
                   END IF
      END DO
      DO I=1,N2
      A=RAY-D2Y(I)
                  IF(A.GT.0.D00)THEN
      ALF=2.*DASIN(.5*H*DSQRT(A))
      G=DCOS(ALF)-DSIN(ALF)/DTAN(ALF*N)
      U2(I)=U2(I)*G
                   ELSE
      ALF=2.*DASINH(.5*H*DSQRT(-A))
      G=DCOSH(ALF)-DSINH(ALF)/DTANH(ALF*N)
      U2(I)=U2(I)*G
                   END IF
      END DO
      CALL FFT(2*N,U2,S2,IT2)
      CALL FFT(N,U1,S1,IT1)
      DO I=1,N1
      SU(I)=SU(I)-(U2(I)+U1(I))/H2
      END DO
      RETURN
      END
*
      SUBROUTINE SB(U,V,ML,S1,IT1,S2,IT2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      INTEGER ML
      COMMON /AB/N,H,PI,SO,H2,HPI,N1,N2,RR1
      DIMENSION U(33000),V(33000)
     *,S1(33000),S2(67000),IT1(33000),IT2(67000)
      COMMON /B/DS(33000)
      DO I=1,N1
      V(I)=U(I)
      END DO
      CALL FFT(N,V,S1,IT1)
      DO I=1,N1
      V(I)=V(I)*DS(I)**ML
      END DO
      CALL FFT(N,V,S1,IT1)
      RETURN
      END
*
      SUBROUTINE LANC(U,RAY,W,LK,S1,IT1,S2,IT2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      INTEGER IRAB,LK,K
      DIMENSION U(33000),R(33000,0:2),V(33000),W(33000),A(2,2),VV(2,2),
     *EV(2),IRAB(2),RAB(2)
     *,S1(33000),S2(67000),IT1(33000),IT2(67000)
      COMMON /AB/N,H,PI,RO,H2,HPI,N1,N2,RR1
      RLB=0.D+00
      A(1,1)=0.D+00
      RUB=0.D+00
      BETA=0.D+00
      CALL SB(U,V,1,S1,IT1,S2,IT2)
      UN=DSQRT(DSCAL(U,V,N1))
      DO I=1,N1
      U(I)=U(I)/UN
      R(I,1)=U(I)
      END DO
      DO K=1,LK
      CALL SCHUR(U,RAY,V,S1,IT1,S2,IT2)
      ALFA=DSCAL(U,V,N1)
      CALL SB(V,W,-1,S1,IT1,S2,IT2)
      DO I=1,N1
      V(I)=W(I)-ALFA*R(I,K)-BETA*R(I,K-1)
      END DO
      GAMMA=DSQRT(DSCAL(V,V,N1))
      CALL SB(V,U,1,S1,IT1,S2,IT2)
      BETA=DSQRT(DSCAL(U,V,N1))
      A(K,2)=ALFA
      IF(K.LT.LK)THEN
      A(K+1,1)=BETA
      DO I=1,N1
      U(I)=V(I)/BETA
      R(I,K+1)=U(I)
      END DO
      END IF
      END DO
      CALL AEE1D(LK,A,EV,VV,RAB,IERR)
*     WRITE(1,*) 'LANC, VV(1,1)=',VV(1,1) ,'  VV(2,1)=',VV(2,1)
*      VPAR=VV(2,1)/VV(1,1)/A(2,1)
      VPAR=VV(2,1)/VV(1,1)/GAMMA
*      VPAR=.5
      WRITE(1,*) 'IH FORMULA v=u+@*su PARAM @=',VPAR
*     WRITE(1,*) 'LANC, EV(1)=',EV(1) ,' EV(2)=',EV(2)
      WRITE(*,*) 'LANC, EV(1)=',EV(1) ,' EV(2)=',EV(2)
      IF(VV(1,1).EQ.0.D+00.AND.VV(2,1).EQ.0.D+00)THEN
      WRITE(1,*)'   '
      WRITE(1,*)'WARNING!!!!!!!!!'
      END IF
      DO I=1,N1
      W(I)=0.D+00
      DO J=1,LK
      W(I)=W(I)+VV(J,1)*R(I,J)
*      W(I)=W(I)+R(I,1)-.5D-00*A(2,1)*R(I,2)
      END DO
      END DO
      RETURN
      END
*
      DOUBLE PRECISION FUNCTION SCH(U,RAY,S1,IT1,S2,IT2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      DIMENSION U(33000),U2(67000),U1(33000)
     *,S1(33000),S2(67000),IT1(33000),IT2(67000)
      COMMON /AB/N,H,PI,RO,H2,HPI,N1,N2,RR1
      COMMON /C/D1Y(33000),D2Y(67000)
      CALL TRIDIAG(U,U1,-1.D+00,4.D+00)
      RL=DSCAL(U,U1,N1)/H2-RAY*(1.D0+RR1)/2*DSCAL(U,U,N1)
*     WRITE(*,*)'RL=',RL
      DO I=1,N1
      U2(I)=U(I)
      U1(I)=U(I)
      U2(I+N)=0.D+00
      END DO
      U2(N)=0.D+00
      CALL FFT(2*N,U2,S2,IT2)
      CALL FFT(N,U1,S1,IT1)
      A1=0.D+00
      A2=0.D+00
      DO I=N1,1,-1
      A=D1Y(I)-RAY*RR1
                  IF(A.GT.0.D00)THEN
      ALF=2.*DASINH(.5*H*DSQRT(A))
      G=DCOSH(ALF)-DSINH(ALF)/DTANH(ALF*N)
      A1=A1+U1(I)*U1(I)*G
                  ELSE
      ALF=2.*DASIN(.5*H*DSQRT(-A))
      G=DCOS(ALF)-DSIN(ALF)/DTAN(ALF*N)
      A1=A1+U1(I)*U1(I)*G
                  END IF
      END DO
      DO I=N2,1,-1
      A=D2Y(I)-RAY
                  IF(A.GT.0.D00)THEN
      ALF=2.*DASINH(.5*H*DSQRT(A))
      G=DCOSH(ALF)-DSINH(ALF)/DTANH(ALF*N)
      A2=A2+U2(I)*U2(I)*G
                  ELSE
      ALF=2.*DASIN(.5*H*DSQRT(-A))
      G=DCOS(ALF)-DSIN(ALF)/DTAN(ALF*N)
      A2=A2+U2(I)*U2(I)*G
                  END IF
      END DO
*     WRITE(*,*)'A1=',A1,'        A2=',A2
      SCH=RL-(A1+A2)/H2
*     WRITE(*,*)'SCH=',SCH
      RETURN
      END
*
      DOUBLE PRECISION FUNCTION SM(U,RAY,S1,IT1,S2,IT2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      DIMENSION U(33000),U1(33000),U2(67000)
     *,S1(33000),S2(67000),IT1(33000),IT2(67000)
      COMMON /AB/N,H,PI,RO,H2,HPI,N1,N2,RR1
      COMMON /C/D1Y(33000),D2Y(67000)
      RM=(1.D0+RR1)/2*DSCAL(U,U,N1)
      DO I=1,N1
      U1(I)=U(I)
      U2(I)=U(I)
      U2(I+N)=0.D+00
      END DO
      U2(N)=0.D+00
      CALL FFT(2*N,U2,S2,IT2)
      CALL FFT(N,U1,S1,IT1)
      B2=0.D+00
      B1=0.D+00
      DO I=N1,1,-1
      A=RAY*RR1-D1Y(I)
                  IF(A.GT.0.D00)THEN
      ALF=2.*DASIN(.5*H*DSQRT(A))
      C=(DSIN(2.*ALF*(DBLE(N)-.5D00))+DSIN(ALF))/DSIN(ALF)
      G=(.5*DBLE(N)-.25*C)/DSIN(ALF*N)**2
      B1=B1+U1(I)*U1(I)*G
                  ELSE
      ALF=2.D00*DASINH(.5*H*DSQRT(-A))
      C=DEXP(-2.D00*ALF*N)
      C1=DEXP(-2.D00*ALF*N1)
      C2=EXP(2.D00*ALF)
      G=(1.D00-C1*C)/(C2-1.D00)+C-2.D00*N*C
      G=G/(1.D00-2.D00*C+C*C)
      B1=B1+U1(I)*U1(I)*G
                  END IF
      END DO
      DO I=N2,1,-1
      A=RAY-D2Y(I)
                  IF(A.GT.0.D00)THEN
      ALF=2.*DASIN(.5*H*DSQRT(A))
      C=(DSIN(2.*ALF*(DBLE(N)-.5D00))+DSIN(ALF))/DSIN(ALF)
      G=(.5*DBLE(N)-.25*C)/DSIN(ALF*N)**2
      B2=B2+U2(I)*U2(I)*G
                  ELSE
      ALF=2.D00*DASINH(.5*H*DSQRT(-A))
      C=DEXP(-2.D00*ALF*N)
      C1=DEXP(-2.D00*ALF*N1)
      C2=EXP(2.D00*ALF)
      G=(1.D00-C1*C)/(C2-1.D00)+C-2.D00*N*C
      G=G/(1.D00-2.D00*C+C*C)
      B2=B2+U2(I)*U2(I)*G
                  END IF
      END DO
      SM=B1*RR1+B2+RM
*     WRITE(*,*)'SM=',SM
*     WRITE(*,*)'B1=',B1,'     B2=',B2
      RETURN
      END
*
*
      SUBROUTINE TRIDIAG(U,V,A,B)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      DIMENSION U(33000),V(33000)
      COMMON /AB/N,H,PI,RO,H2,HPI,N1,N2,RR1
      NN=N1-1
      V(1)=B*U(1)+A*U(2)
      DO I=2,NN
      V(I)=A*(U(I-1)+U(I+1))+B*U(I)
      END DO
      V(N1)=A*U(N1-1)+B*U(N1)
      RETURN
      END
*
      SUBROUTINE FFT(M,U,S,IT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      DIMENSION S(67000),IT(67000),X(67000),B(67000),U(67000)
      IB=1
      M1=M-1
      DO I=1,M1
      X(IB+I+1)=U(I)
      END DO
      CALL SIN816(M,IB,X,B,S,IT)
      DO I=1,M1
      U(I)=X(IB+I+1)
      END DO
      RETURN
      END

      SUBROUTINE TRANS(M,S,ITRAN)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      DIMENSION S(67000),ITRAN(67000)
      DATA RCON5/.5D00/
      RCON3=4.*DATAN(1.D00)
      MH=M/2
      MQ=M/4
      MC=M/8
      ITRAN(MQ)=MQ
      ITRAN(MC)=MC
      I1=MH-MC
      ITRAN(I1)=I1
      JH=MH
      JQ=MQ
      JC=MC
   50 IF(JC.EQ.1)GOTO 70
      JS=JH
      JH=JQ
      JQ=JC
      JC=JC/2
      JT=JH-JC
      JE=MH-JH
      DO 60 J=JH,JE,JS
      ITRAN(J-JT)=ITRAN(J+JQ)/2
      ITRAN(J-JC)=MH-ITRAN(J-JT)
      ITRAN(J+JC)=ITRAN(J-JQ)/2
   60 ITRAN(J+JT)=MH-ITRAN(J+JC)
      GOTO 50
   70 DO 80 J=2,MH
      S(MH+J)=RCON5/DSIN(ITRAN(J-1)*RCON3/M)
   80 CONTINUE
      JE=MQ
   90 DO 100 J=2,JE
      I1=J+JE
  100 S(I1)=S(2*I1-1)
      IF(JE.EQ.2)GOTO 110
      JE=JE/2
      GOTO 90
  110 CONTINUE
      ITRAN(M)=2
      JH=MQ
      JS=MH
      JQ=6
  120 JB=JS+JH
      DO 130 J=JB,M,JS
  130 ITRAN(J)=JQ-ITRAN(J+JH)
      ITRAN(JS)=JQ/2
      IF(JH.EQ.1)GOTO 170
      JS=JH
      JH=JH/2
      JQ=JQ+JQ-2
      GOTO 120
 170  CONTINUE
      RETURN
      END
*
      SUBROUTINE SIN816(M,IB,X,B,S,IT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      DIMENSION X(67000),B(67000),S(67000),IT(67000)
      JS=M
      JH=M/2
   10 JL=JS+2
      LB=JL+IB
      DO 20 JJ=2,JH,1
      J1=IB+JJ
      J2=JL-JJ
      J3=LB-JJ
      B(J2)=X(J1)+X(J3)
      X(J1)=X(J1)-X(J3)
   20 CONTINUE
      B(JH+1)=X(IB+JH+1)
      JS=JH
      JH=JH/2
      JB=JS+JH+1
      JE=M-JH+1
      IF(JH.EQ.2)GOTO 85
      JQ=JH/2
      IE=JQ-1
      DO 80 JJ=JB,JE,JS
      IP=JJ+JQ
      IM=JJ-JQ
      I5=JJ-JH
      RC=S(JJ)
      R1=B(JJ)*RC
      B(JJ)=R1+B(I5)
      B(I5)=R1-B(I5)
      R1=(B(IP)+B(IM))*RC
      B(IP)=R1+B(IM)
      B(IM)=R1-B(IM)
      DO 70 I=1,IE,1
      I1=IP+I
      I2=IP-I
      I3=IM+I
      I4=IM-I
      R1=(B(I1)+B(I4))*RC
      R2=(B(I2)+B(I3))*RC
      B(I1)=R1+B(I3)
      B(I3)=R1-B(I3)
      B(I2)=R2+B(I4)
      B(I4)=R2-B(I4)
   70 CONTINUE
   80 CONTINUE
      GO TO 10
   85 R6=X(IB+2)
      R7=X(IB+3)
      R8=X(IB+4)
      DO 95 JJ=JB,JE,JS
      J1=JJ-2
      J2=JJ-1
      J3=JJ+1
      I1=IB+IT(J1)
      I2=IB+IT(J2)
      I3=IB+IT(JJ)
      I4=IB+IT(J3)
      R1=B(JJ)*S(JJ)
      R2=R1+B(J1)
      R1=R1-B(J1)
      R3=(B(J3)+B(J2))*S(JJ)
      R4=(R3+B(J2))*S(J3)
      R5=(R3-B(J2))*S(J2)
      X(I4)=R4+R2
      X(I3)=R4-R2
      X(I2)=R5+R1
      X(I1)=R5-R1
   95 CONTINUE
      MQ=M/4
      JB=IB+1
      X(2*MQ+JB)=R6-R8
      R5=(R6+R8)*S(4)
      X(MQ+JB)=R5+R7
      X(3*MQ+JB)=R5-R7
   90 CONTINUE
      A=DSQRT(2.D00/M)
      M2=M+IB
      DO I=1,M2
      X(I)=X(I)*A
      END DO
      RETURN
      END
*
      DOUBLE PRECISION FUNCTION DSCAL(U,V,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z),INTEGER*4 (I-N)
      DIMENSION U(67000),V(67000)
      A=0.D00
      DO I=1,N
      A=A+U(I)*V(I)
      END DO
      DSCAL=A
      RETURN
      END
*
      DOUBLE PRECISION FUNCTION DASINH(Y)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DASINH=DLOG(Y+DSQRT(Y*Y+1.D0))
      RETURN
      END
*
      SUBROUTINE AEE1D(N,A,EV,V,RAB1,IERR)
      INTEGER N,IERR,I,J,L,M,O,P,Q,R,S
      DOUBLE PRECISION A,EV,V,RAB1,B,C,D,E,F,G,H,K,DSQRT,DABS,DSIGN,
     *SYS051
      DIMENSION A(N,2),EV(N),V(N,N),RAB1(N)
      DATA SYS051/2.2204460492503D-16/
      K=SYS051
      DO 2 I=1,N
      DO 1 J=1,N
      V(I,J)=0.
    1 CONTINUE
      V(I,I)=1.
      EV(I)=A(I,2)
      RAB1(I)=A(I,1)
    2 CONTINUE
      IERR=0
      IF(N.EQ.1) GO TO 16
      DO 3 L=2,N
    3 RAB1(L-1)=RAB1(L)
      RAB1(N)=0.0D0
      DO 11 P=1,N
      M=0
    4 DO 5 Q=P,N
      IF(Q.EQ.N) GO TO 6
      IF(DABS(RAB1(Q)).LE.K*(DABS(EV(Q))+DABS(EV(Q+1)))) GO TO 6
    5 CONTINUE
    6 F=EV(P)
      IF(Q.EQ.P) GO TO 11
      IF(M.EQ.30) GO TO 15
      M=M+1
      E=(EV(P+1)-F)/(2.0*RAB1(P))
      IF(DABS(E).LE.1.0D0) G=DSQRT(1.0D0+E*E)
      IF(DABS(E).GT.1.0D0) G=DABS(E)*DSQRT(1.0D0+(1.0D0/E)**2)
      E=EV(Q)-F+RAB1(P)/(E+DSIGN(G,E))
      H=1.0D0
      C=1.0D0
      F=0.0D0
      S=Q-P
      DO 10 R=1,S
      L=Q-R
      D=H*RAB1(L)
      B=C*RAB1(L)
      IF(DABS(D).LT.DABS(E)) GO TO 7
      C=E/D
      G=DSQRT(C*C+1.0D0)
      RAB1(L+1)=D*G
      H=1.0D0/G
      C=C*H
      GO TO 8
    7 H=D/E
      G=DSQRT(H*H+1.0D0)
      RAB1(L+1)=E*G
      C=1.0D0/G
      H=H*C
    8 E=EV(L+1)-F
      G=(EV(L)-E)*H+2.0*C*B
      F=H*G
      EV(L+1)=E+F
      E=C*G-B
      DO 9 O=1,N
      D=V(O,L+1)
      V(O,L+1)=H*V(O,L)+C*D
      V(O,L)=C*V(O,L)-H*D
    9 CONTINUE
   10 CONTINUE
      EV(P)=EV(P)-F
      RAB1(P)=E
      RAB1(Q)=0.0D0
      GO TO 4
   11 CONTINUE
      DO 14 R=2,N
      L=R-1
      O=L
      F=EV(L)
      DO 12 M=R,N
      IF(EV(M).GE.F) GO TO 12
      O=M
      F=EV(M)
   12 CONTINUE
      IF(O.EQ.L) GO TO 14
      EV(O)=EV(L)
      EV(L)=F
      DO 13 M=1,N
      F=V(M,L)
      V(M,L)=V(M,O)
      V(M,O)=F
   13 CONTINUE
   14 CONTINUE
      GO TO 16
   15 IERR=P
   16 CONTINUE
      RETURN
      END
      SUBROUTINE AEH2D(N,A,EV,RAB,IERR)
      DIMENSION A(N,N),EV(N),RAB(*)
      INTEGER N,IERR,B,C,D,E,F,G,L,M,O,P,Q,R,S
      INTEGER NB,NL,NO,NP
      DOUBLE PRECISION A,EV,RAB,H,I,J,K,DSQRT,DABS,DSIGN,T,U,V,W,
     *X,Y,Z,BA,BB,SYS051
      DATA SYS051 /8.2718061255D-16/
       BB=SYS051
      DO 1 B=1,N
   1  EV(B)=A(B,B)
      DO 14 F=1,N
      B=N+1-F
      NB=N+B
      E=B-1
      J=0.0D0
      K=0.0D0
      IF(E.LT.1) GO TO 3
      DO 2 D=1,E
   2  K=K+DABS(A(B,D))
      IF(K.NE.0.0D0) GO TO 4
   3  RAB(B)=0.0
      RAB(NB)=0.0
      GO TO 13
   4  DO 5 D=1,E
      A(B,D)=A(B,D)/K
      J=J+A(B,D)*A(B,D)
   5  CONTINUE
      RAB(NB)=K*K*J
      H=A(B,E)
      I=-DSIGN(DSQRT(J),H)
      RAB(B)=K*I
      J=J-H*I
      A(B,E)=H-I
      IF(E.EQ.1) GO TO 11
      H=0.0D0
      DO 9 C= 1,E
      I=0.0D0
      DO 6 D=1,C
   6  I=I+A(C,D)*A(B,D)
      G=C+1
      IF(E.LT.G) GO TO 8
      DO 7 D=G,E
   7  I=I+A(D,C)*A(B,D)
   8  RAB(C)=I/J
      H=H+RAB(C)*A(B,C)
   9  CONTINUE
      J=H/(J+J)
      DO 10 C=1,E
      H=A(B,C)
      I=RAB(C)-J*H
      RAB(C)=I
      DO 10 D=1,C
      A(C,D)=A(C,D)-H*RAB(D)-I*A(B,D)
  10  CONTINUE
  11  DO 12 D=1,E
  12  A(B,D)=K*A(B,D)
  13  J=EV(B)
      EV(B)=A(B,B)
      A(B,B)=J
  14  CONTINUE
      IERR=0
      IF(N.EQ.1) GO TO 28
      DO 15 L=2,N
      NL=N+L
  15  RAB(NL-1)=RAB(NL)
      V=0.0D0
      T=0.0D0
      RAB(2*N)=0.0
      DO 26 O=1,N
      NO=N+O
      M=0
      X=BB*(DABS(EV(O))+DSQRT(RAB(NO)))
      IF(T.GT.X) GO TO 16
      T=X
      U=T*T
  16  DO 17 P=O,N
      NP=N+P
      IF(RAB(NP).LE.U) GO TO 18
  17  CONTINUE
  18  IF(P.EQ.O) GO TO 22
  19  IF(M.EQ.30) GO TO 27
      M=M+1
      R=O+1
      BA=DSQRT(RAB(NO))
      W=EV(O)
      Y=(EV(R)-W)/(2.0D0*BA)
      IF(DABS(Y).GT.1.0D0) Z=DABS(Y)*DSQRT(1.0+(1.0/Y)**2)
      IF(DABS(Y).LE.1.0D0) Z=DSQRT(Y*Y+1.0)
      EV(O)=BA/(Y+DSIGN(Z,Y))
      X=W-EV(O)
      DO 20 L=R,N
  20  EV(L)=EV(L)-X
      V=V+X
      W=EV(P)
      IF(W.EQ.0.0D0) W=T
      X=W
      BA=0.0D0
      S=P-O
      DO 21 Q=1,S
      L=P-Q
      Y=W*X
      NL=N+L
      Z=Y+RAB(NL)
      RAB(NL+1)=BA*Z
      BA=RAB(NL)/Z
      EV(L+1)=X+BA*(X+EV(L))
      W=EV(L)-RAB(NL)/W
      IF(W.EQ.0.0D0) W=T
      X=W*Y/Z
  21  CONTINUE
      RAB(NO)=BA*W
      EV(O)=X
      IF(X.EQ.0.0D0) GO TO 22
      IF(DABS(RAB(NO)).LE.DABS(U/X)) GO TO 22
      RAB(NO)=X*RAB(NO)
      IF(RAB(NO).NE.0.0) GO TO 19
  22  Y=EV(O)+V
      IF(O.EQ.1) GO TO 24
      DO 23 Q=2,O
      L=O+2-Q
      IF(Y.GE.EV(L-1)) GO TO 25
      EV(L)=EV(L-1)
  23  CONTINUE
  24  L=1
  25  EV(L)=Y
  26  CONTINUE
      GO TO 28
  27  IERR=O
  28  CONTINUE
      RETURN
      END