        PROGRAM D2R3

C       power method for smallest eigenvalue

        PARAMETER (NP=20,eps=0.5e-6,m=20)

        DIMENSION A(NP,NP),b(np,np),INDX(NP),c(np,np)

        DIMENSION X(NP),y(np)

        WRITE(*,*) 'give # of equations'

        READ(*,*) N

        DO 20 I=1,N

        WRITE(*,'(1X,A,I2,A)') 'give the elements of the row ',I,

     *  '  of the A mat.'

        READ(*,*) (A(I,L), L=1,N)

20      CONTINUE

        DO 21 I=1,N

        WRITE(*,'(1X,A,I2,A)') 'give the elements of the row ',I,

     *  '  of the B mat.'

        READ(*,*) (B(I,L), L=1,N)

21      CONTINUE

        WRITE(*,*) 'give the elements of the initial evc'

        READ(*,*) (X(K), K=1,N)

C       Do LU decomposition

        do 100 i=1,n

        do 100 j=1,n

100     c(i,j)=a(i,j)

        CALL LUDCMP(c,N,NP,INDX,P)

C       start the iterations

        k=0

        xl1=0.

40      do 24 i=1,n

        y(i)=0.

        do 24 j=1,n

        y(i)=y(i)+b(i,j)*x(j)

24      continue

        CALL LUBKSB(c,N,NP,INDX,y)

        xl=0.

        do 22 i=1,n

        xl=xl+y(i)*y(i)

22      continue

        xl=sqrt(xl)

        do 23 i=1,n

        x(i)=y(i)/xl

23      continue

        k=k+1

        if(k.gt.m) then

        write(*,60) m

60      format(5x,'THE ITERATIONS DID NOT CONVERGE IN',i2,'ITEERATIONS')

        go to 50

        endif

        aa=0.

        bb=0.

        do 150 i=1,n

        do 150 j=1,n

        aa=aa+a(i,j)*x(i)*x(j)

150     bb=bb+b(i,j)*x(i)*x(j)

        xl=aa/bb

        if(abs(xl-xl1).le.eps) then

C       xl=1./xl

        WRITE(*,30) k,xl

30      format(5x,'ITER # =',i2,5x,'SMALLEST EVL =',f12.6/)

        DO 25 I=1,N

        WRITE(*,31) i,X(I)

31      format(5x,'X(',I2,')='f12.6)

25      CONTINUE

        else

        xl1=xl

        go to 40

        endif

50      continue

        END

      SUBROUTINE LUDCMP(A,N,NP,INDX,D)

      PARAMETER (NMAX=100,TINY=1.0E-20)

      DIMENSION A(NP,NP),INDX(N),VV(NMAX)

      D=1.

      DO 12 I=1,N

        AAMAX=0.

        DO 11 J=1,N

          IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))

11      CONTINUE

        IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'

        VV(I)=1./AAMAX

12    CONTINUE

      DO 19 J=1,N

        IF (J.GT.1) THEN

          DO 14 I=1,J-1

            SUM=A(I,J)

            IF (I.GT.1)THEN

              DO 13 K=1,I-1

                SUM=SUM-A(I,K)*A(K,J)

13            CONTINUE

              A(I,J)=SUM

            ENDIF

14        CONTINUE

        ENDIF

        AAMAX=0.

        DO 16 I=J,N

          SUM=A(I,J)

          IF (J.GT.1)THEN

            DO 15 K=1,J-1

              SUM=SUM-A(I,K)*A(K,J)

15          CONTINUE

            A(I,J)=SUM

          ENDIF

          DUM=VV(I)*ABS(SUM)

          IF (DUM.GE.AAMAX) THEN

            IMAX=I

            AAMAX=DUM

          ENDIF

16      CONTINUE

        IF (J.NE.IMAX)THEN

          DO 17 K=1,N

            DUM=A(IMAX,K)

            A(IMAX,K)=A(J,K)

            A(J,K)=DUM

17        CONTINUE

          D=-D

          VV(IMAX)=VV(J)

        ENDIF

        INDX(J)=IMAX

        IF(J.NE.N)THEN

          IF(A(J,J).EQ.0.)A(J,J)=TINY

          DUM=1./A(J,J)

          DO 18 I=J+1,N

            A(I,J)=A(I,J)*DUM

18        CONTINUE

        ENDIF

19    CONTINUE

      IF(A(N,N).EQ.0.)A(N,N)=TINY

      RETURN

      END

      SUBROUTINE LUBKSB(A,N,NP,INDX,B)

      DIMENSION A(NP,NP),INDX(N),B(N)

      II=0

      DO 12 I=1,N

        LL=INDX(I)

        SUM=B(LL)

        B(LL)=B(I)

        IF (II.NE.0)THEN

          DO 11 J=II,I-1

            SUM=SUM-A(I,J)*B(J)

11        CONTINUE

        ELSE IF (SUM.NE.0.) THEN

          II=I

        ENDIF

        B(I)=SUM

12    CONTINUE

      DO 14 I=N,1,-1

        SUM=B(I)

        IF(I.LT.N)THEN

          DO 13 J=I+1,N

            SUM=SUM-A(I,J)*B(J)

13        CONTINUE

        ENDIF

        B(I)=SUM/A(I,I)

14    CONTINUE

      RETURN

      END

