
      include 'trh.f'
      include 'hax.f'



      subroutine householder(n)
c     ELIMINACJA HOUSEHOLDERA MACIERZY A i TRANSFORMACJA WEKTORA b
c     A i b czytane z dysku: A= file'A', b= file'b' 
c     WEKTORY w (wymiary: n,n-1,...2) SA ZAPAMIETANE JAKO KOLEJNE
c     KOLUMNY MACIERZY WW; H=I-ww^T.
c     ROZWIAZANIE UKLADU Ax=b
c     WYNIK x NA DYSKU: file='x'
c     INFORMACJA O OSOBLIWOSCI
      implicit none
      real*8 A,AA,b,c,x,y,w,alpha,WW,O,E,eps
      integer i,j,k,k1,k2,l,n
      dimension A(n,n),WW(n,n),b(n),c(n),x(n),y(n),w(n)

c     DANE
      O=real(0)
      E=real(1)
      eps=dble(.00000000000001)
      open(unit=10,file='A')
      open(unit=11,file='b')
      rewind(10)
      rewind(11)

c     PO ROZKLADZIE
      open(unit=12,file='HA')
      open(unit=13,file='Hb')
      open(unit=14,file='WW')
      open(unit=15,file='x')
c     CZYTANIE DANYCH
      do 1,k2=1,n,1
      do 2,k1=1,n,1
      read (10,*) A(k1,k2)
    2 continue
    1 continue

      do 3,j=1,n,1
      read (11,*) b(j)
    3 continue

c     START PETLI ELIMINACJI  
c     l - licznik krokow eliminacji
c     ZEROWANIE WW
      do 4,k1=1,n,1
      do 5,k2=1,n,1
      WW(k1,k2)=O
    5 continue
    4 continue
      do 6,l=1,n-1,1
      do 7,k1=l,n,1
      c(k1-l+1)=A(k1,l)
    7 continue

      call trh(n-l+1,c,w,alpha)
      call hax(n-l+1,w,c,y)

c     KOLUMNA NR l
      do 8,k=1,n-l+1,1
      A(l+k-1,l)=y(k)
      WW(l+k-1,l)=w(k)
    8 continue

c     TRANSFORMACJA POZOSTALYCH KOLUMN A I KOLUMNY b
      do 9, k2=l+1,n,1
      do 10,k1=l,n,1
      x(k1-l+1)=A(k1,k2)
   10 continue

c     PRZEKSZTALCENIE I WPISANIE KOLUMNY A
      call hax(n-l+1,w,x,y)
      do 11,k1=l,n,1
      A(k1,k2)=y(k1-l+1)
   11 continue
    9 continue

c     TRANSFORMACJA KOLUMNY b
      do 12 k=l,n,1
      x(k-l+1)=b(k)
   12 continue

c     PRZEKSZTALCENIE I WPISANIE DO KOLUMNY b
      call hax(n-l+1,w,x,y)
      do 13,k=l,n,1
      b(k)=y(k-l+1)
   13 continue
c     ZAMKNIECIE PETLI ELIMINACJI
    6 continue

c     ROZWIAZANIE UKLADU Ax=b
      do 14,k=n,1,-1
      x(k)=b(k)
      do 15,j=k+1,n,1
      x(k)=x(k)-A(k,j)*x(j)
   15 continue
      if(abs(A(k,k)).lt.eps)then
      write(*,*)'ILL CONDITION - WIERSZ:',k
      write(*,*)'WEKTOR x DO KITU!'
      write(*,*)'<ENTER>'
      read (*,*)
      x(k)=x(k)/A(k,k)
      else
      x(k)=x(k)/A(k,k)
      endif
   14 continue

      write(*,*) 'OUTPUT NA DYSK'
      write(*,*) 'MACIERZ A PO ELIMINACJI, KOLUMNAMI: HA'
      write(*,*) 'PRAWA STRONA PO ELIMINACJI: Hb'
      write(*,*) 'WEKTORY W W KOLUMNACH 1 DO n-1 MACIERZY WW'
      write(*,*) 'ROZWIAZANIE UKLADU Ax=b: x'
      do 16,k2=1,n,1
      do 17,k1=1,n,1
      write(12,*) A(k1,k2)
   17 continue
   16 continue
      do 18,k=1,n,1
      write(13,*) b(k)
   18 continue
      do 19,k2=1,n,1
      do 20,k1=1,n,1
      write(14,*) WW(k1,k2)
   20 continue
   19 continue
      do 21,k=1,n,1
      write(15,*) x(k)
   21 continue
      close(10)
      close(11)
      close(12)
      close(13)
      close(14)
      close(15)
      return
      end

