      program ewp
c
c-----------------------------------------------------------------------
c    Hauptprogramm zur
c    Praktikumsaufgabe PA3
c
c    --> (fuer Hauptfachstudenten) <--
c
c    Berechnung der Eigenwerte einer Matrix A durch das QR-Verfahren
c
c    Uebergabeparameter:
c
c    work(ik)........Steifigkeitsmatrix
c    work(im)........Massenmatrix
c    work(ig)........Systemmatrix
c    work(iq)........Transformationsmatrix
c    n...............Dimension der Systemmatrix
c    nd..............Dimenson der bearbeiteten, linken
c                    oberen Blockmatrix
c    delta...........Schranke fuer Deflationskriterium
c    iter............Anzahl an Iterationsschritten
c    sigma...........Shift-Parameter
c 
c    Lokale Variablen:
c
c    work............Arbeitsfeld
c    wmax............Groesse des Arbeitsfelds
c    maxiter.........maximale Anzahl von Iterationsschritten
c    ik, nk ,im ,nm,
c    ig, ng, iq, nq..Adressen im Arbeitsfeld
c    i...............Hilfsvariable
c
c    Benoetigte Unterprogramme:
c
c    SUBROUTINE ausg_ssm(...)
c    SUBROUTINE ausg_mat(...)
c    SUBROUTINE ausg_ew (...)
c    SUBROUTINE defl_tol(...)
c    SUBROUTINE klass_ewp(...)
c    SUBROUTINE ausg_iter(...)
c    SUBROUTINE qr(...)
c    SUBROUTINE rq(...)
c-----------------------------------------------------------------------
c
      IMPLICIT none

      INTEGER          wmax,maxiter, 
     *                 ik,nk,im,nm,ig,ng,iq,nq,
     *                 n,nd, iter,i
      DOUBLE PRECISION delta,sigma

c   Arbeitsfeld, max. Anzahl von Iterationen
      PARAMETER        (wmax=1000, maxiter=25)
      DOUBLE PRECISION work(wmax)

c   Einlesen der Problemdimension
      OPEN(16, FILE='k.dat')
      READ(16,*) n

c   Adressrechnung
c     Steifigkeitsmatrix (in SSM)       
        ik = 1
        nk = (n*(n+1))/2
c     Massenmatrix (nur Diagonalelemente)       
        im = ik+nk
        nm = n
c     Systemmatrix (vollbesetzt)       
        ig = im+nm
        ng = n*n
c     Transformationsmatrix        
        iq = ig+ng
        nq = n*n
      
        IF(iq-1+nq .GT. wmax) 
     *       STOP 'Arbeitsfeld zu klein'

c   Einlesen der Steifigkeitsmatrix im SSM
      DO 100 i=1,nk
 100  READ(16,*) work(ik-1+i)
      CLOSE(16)
c       Kontrollausgabe
          WRITE(*,'(//A/)') 'Steifigkeitsmatrix'
          CALL ausg_ssm(work(ik),n,'(9F10.1)')

c   Einlesen der diagonalen Massenmatrix
      OPEN(16, FILE='m.dat')
      DO 200 i=1,n
 200  READ(16,*) work(im-1+i)
      CLOSE(16)

c   Kontrollausgabe
      WRITE(*,'(//A/)') 'Massenmatrix (Diagonalelemente)'
      CALL ausg_mat(work(im),n,1,'(1F10.3)')

c   Systemmatrix des klassischen EWP
c
c ** Ergaenzen Sie das Programm um die Subroutine **
c    CALL klass_ewp(n,work(ik),work(im),work(ig)) 

c   Kontrollausgabe
      WRITE(*,'(//A/)') 'Systemmatrix'
      CALL ausg_mat(work(ig),n,n,'(9F10.1)')

c   Initialisierung
      iter = 0
      nd   = n
      sigma= 0

c   Deflationstoleranz fuer Matrix A
c
c ** Ergaenzen Sie das Programm um die Subroutine **
c    CALL defl_tol(work(ig),n,delta)

      CALL ausg_iter(work(ig),n,iter,sigma,nd)

c   Iteration
 300  iter = iter + 1
    
c   Eigenwertberechnung durch QR-Zerlegung und
c   RQ-Multiplikation 
c
c ** Ergaenzen Sie das Programm um die Subroutinen **
c    CALL qr(work(ig),work(iq),n,nd)
c    CALL rq(work(ig),work(iq),n,nd)

      CALL ausg_iter(work(ig),n,iter,sigma,nd)
      IF(iter .LT. maxiter) GOTO 300

c   Konvergenzprobleme
      STOP  '*** Abbruch wegen Konvergenzproblemen ***'

c   Eigenwerte
 900  CALL ausg_ew(work(ig),n)

      END




      SUBROUTINE ausg_ew(A,n)
c
c       Ausgabe der Eigenwerte
c       einer Rechtsdreiecksmatrix A
c
c     A(n,n) ... Matrix                (I)
c     n ........ Dimension             (I)

      INTEGER           n, i
      DOUBLE PRECISION  A(n,n)

      WRITE(*   ,'(//A)') 'Eigenwerte:'

      DO i=1,n
        WRITE(*   ,100) A(i,i)
 100    FORMAT(2X,F14.5)
      END DO

      RETURN
      END




      SUBROUTINE ausg_iter(a,n0,iter,s,n)
c
c       Ausgabe einer Matrix A
c
c     A(n0,n0) . Matrix                             (I)
c     n0 ....... Dimension                          (I)
c     iter ..... Iterations-Index                   (I)
c     s ........ Shift-Parameter                    (I)
c     n ........ Ordnung der Restmatrix             (I)

      INTEGER           n0,iter,n, i,j
      DOUBLE PRECISION  A(n0,n0),s

      WRITE(*   ,50) iter
 50   FORMAT(//,'Iterationsschritt: ',I3,/)

      IF(iter .GT. 0) THEN
        WRITE(*   ,51) n,s
 51     FORMAT('Ordnung: ',I3,8X,'Shift: ',G15.6,/)
      END IF 

      DO i=1,n0
        WRITE(*   ,100) (a(i,j), j=1,n0)
 100    FORMAT(9G10.3)
      END DO

      RETURN
      END




      SUBROUTINE ausg_ssm(A,n,fmt)

c       Ausgabe einer symmetrischen
c       n x n - Matrix A gespeichert
c       im SSM
c
c     A.........auszugebende Matrix (I)
c     n.........Dimension           (I)
c     fmt.......Ausgabeformat       (I)

      INTEGER n, izanf,i,j
      DOUBLE PRECISION A(1)
      CHARACTER*20 fmt
      
      izanf=0
      DO i=1,n
        WRITE(*,fmt) (A(izanf+j), j=1,i)
        izanf = izanf+i
      END DO
      
      RETURN
      END




      SUBROUTINE ausg_mat(A,n,m,fmt)

c       Ausgabe einer n x m - Matrix A
c
c     A.........auszugebende Matrix (I)
c     n.........Zeilenzahl          (I)
c     m.........Spaltenzahl         (I)
c     fmt.......Ausgabeformat       (I)

      INTEGER n,m, i,j
      DOUBLE PRECISION A(n,m)
      CHARACTER*20 fmt
      
      DO i=1,n
        WRITE(*,fmt) (A(i,j), j=1,m)
      END DO
      
      RETURN
      END



