      PROGRAM INVPOW
C     
C.. THISPROGRAM WILL EXTRACT NUMEIG LOWEST EIGENVALUES AND EIGENVECTORS
C     OF SYMMETRIC POSITIVE DEFINITE MATRIX A USING INVERSE POWER METHOD
C     
C     Original version by Sami Saarinen, CSC for Finnish Supercomputer
C     Evaluation Project (1988)
C     
C     Version 2 by Sami Saarinen, CSC for OPMVAX Replacement Project
C     (1991). Major changes in release:
C     
C     - ABS -> DABS, SQRT -> DSQRT
C     - source code indenting for better understanding
C     - matrix input file read in packed format (subr. READI)
C     - removal of I/O -operations and introdution of tmp-tables
c
c     Version 3
c     - I/O added back (in PUTSCR) and atmp(,) removed
C     
CDEP  
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (N=605)
c--      PARAMETER (N=2015)
      DIMENSION A(N,N),ZK(N),YK(N),WORK2(N),WORK1(1),EIGV(N)
c--      dimension atmp(n,n)
c--      common /base/ a,zk,atmp
      common /base/ a,zk
      CHARACTER*40 FILNAM                           
      EQUIVALENCE (WORK1(1),YK(1))
C     
      CALL MESAGE('** JOB STARTED **   ',1,1)
C     
      WRITE(6,*)'ENTER NUMBER OF EIGENVALUES TO BE FOUND:'
      READ(5,*) NUMEIG
      WRITE(6,*)'NUMEIG=',NUMEIG
      WRITE(6,*)'ENTER CONVERGENCE TOLERANCE:'
      READ(5,*) CTOL
      CTOL=DABS(CTOL)
      WRITE(6,*)'CTOL=',CTOL
      WRITE(6,*)'ENTER NUMBER OF ITERATIONS ALLOWED:'
      READ(5,*) MAXIT  
      WRITE(6,*)'MAXIT=',MAXIT
C     WRITE(6,*)'ENTER INPUT MATRIX FILE NAME:'
C     READ(5,FMT='(A)') FILNAM
C     WRITE(6,*)'FILENAME='//FILNAM
CDEP  
      OPEN(1,file='invpow.mat',STATUS='OLD')
C     
      CALL MESAGE('** READ MATRIX A    ',6,1)
C     
      READ(1,*) NN      
      WRITE(6,*)'MATRIX SIZE=',NN
      CALL READI(1,A,ZK,NN)
C     
      CALL MESAGE('** READ MATRIX DONE ',6,0)
C     
      CLOSE(1)     
C     
C.. FINDNUMEIG EIGENVALUES AND -VECTORS
C     
CDEP  
      OPEN(20,STATUS='UNKNOWN',FORM='UNFORMATTED')
C     
      NDOF=NN
      DO 100 IEIG=1,NUMEIG
         WRITE(6,*)'----- EIGENVALUE NUMBER ',IEIG,' -----'
C     
         CALL MESAGE('** MAIN LOOP 100 ** ',2,1)
C     
C.. WRITEMATRIX A TO A SCRATCH FILE
         CALL PUTSCR(20,A,NN,NN,1)
c--         CALL PUTSCR(atmp,A,NN,NN,1)
C.. COPYDIAGONALS TO INITIAL GUESS EIGENVECTOR
         CALL COPDIA(A,NN,ZK)                   
C.. FINDTHE MOST SUITABLE SHIFT VALUE
         CALL PUTSHF(A,YK,NN,ESHIFT)
C.. DECOMPOSEMATRIX  A-ESHIFT*I  ONCE PER EIGENVALUE
C     
         CALL MESAGE('** LU-DECOMP STARTED',4,1)
C     
         CALL DECOM3(1,A,YK,NN)            
C     
         CALL MESAGE('** LU-DECOMP DONE   ',4,0)
C     
C.. EXTRACTEIGENVALUE BY USING INVERSE POWER METHOD
C     
         CALL MESAGE('** EIGV. EXTRACTION ',5,1)
C     
         CALL EIGENV(A,ZK,YK,NN,CTOL,MAXIT,IEIG,ESHIFT,EIGV(IEIG))
C     
         CALL MESAGE('** EXTRACTION DONE  ',5,0)
C     
C.. PRINTOUT THE SOLUTION VECTOR
C     CALL WRITEO(ZK,NN)
C.. GETMATRIX A FROM AUXILIARY STORAGE
         CALL GETSCR(20,A,NN,NN)      
c--         CALL GETSCR(atmp,A,NN,NN)      
C.. PERFORMSOLUTION ERROR CHECK
         CALL ECHECK(A,ZK,YK,NN,EIGV(IEIG))
         IF (IEIG.EQ.NUMEIG) GOTO 105
         IF (IEIG.EQ.1)                               
     +        OPEN(21,STATUS='UNKNOWN',FORM='UNFORMATTED') 
C.. DO DEFLATION
C     
         CALL MESAGE('** DEFLATION STARTED',3,1)
C     
         CALL DEFLAT(A,ZK,NN,WORK1,WORK2,21)
c--         CALL DEFLAT(A,ZK,NN,WORK1,WORK2)
C.. REDUCETHE ORDER OF PROBLEM
         NN=NN-1
         WRITE(6,*)
     +        '** THE ORDER OF PROBLEM HAS BEEN REDUCED TO ',NN
C     
         CALL MESAGE('** DEFLATION DONE   ',3,0)
C     
 105     CONTINUE
C     
         CALL MESAGE('** MAIN LOOP 100 ** ',2,0)
C     
 100  CONTINUE                             
C     
      WRITE(6,2010)'NUMBER','EIGENVALUE'
 2010 FORMAT(5X,A12,5X,A15)
      DO 200 IEIG=1,NUMEIG
         WRITE(6,2000) IEIG,EIGV(IEIG)
 200  CONTINUE
 2000 FORMAT(5X,I12,5X,1P,G20.10)            
C     
      CALL MESAGE('** END OF JOB **    ',1,0)
C     
      STOP 'OK EXIT'
      END
