c      program directsolver
c	programa principal para resolver eq lineares - metodo direto
c	projeto 3.6 do livro Calculo Numerico

      implicit real *8 (a-h,o-z)
	parameter (np=10)
	dimension a(np,np),al(np,np),au(np,np),p(np,np)
	dimension x(np),b(np),z(np),eps(np),temp(np),q(np)
      open(unit=1,file='input.txt',status='old')
      open(unit=2,file='output.txt',status='unknown')        
c
c	Leitura de dados
c
      read(1,*)n,dl
      write(*,*)'n= ',n,' dl= ',dl
	pi=4*atan(1.d0)
	write(*,*)'pi= ',pi
      read(1,*)sigma
      write(*,*)'sigma= ',sigma
      read(1,*)area
      write(*,*)'area= ',area
c
	write(*,*)'o vetor eps(i) eh:'
	write(2,*)'o vetor eps(i) eh:'
	read (1,*) (eps(i),i=1,n)
      write(*,*)'eps(i)=', (eps(i),i=1,n)
	write(2,*)'eps(i)=', (eps(i),i=1,n)

	write(*,*)'o vetor temp(i) eh:'
	write(2,*)'o vetor temp(i) eh:'
	read (1,*) (temp(i),i=1,n)
      write(*,*)'temp(i)=', (temp(i),i=1,n)
	write(2,*)'temp(i)=', (temp(i),i=1,n)

	  ! calculo da matriz associada a(i,j)

	a(1,1)=1.d0
	aux1=-(1.d0-eps(1))/4.d0

	a(1,2)=aux1
	a(1,3)=aux1
	a(1,4)=aux1
	a(1,5)=aux1

	a(2,2)=1.d0
	aux2=-(1.d0-eps(2))/4.d0

	a(2,1)=aux2
	a(2,3)=aux2
	a(2,4)=aux2
	a(2,5)=aux2

	a(3,3)=1.d0
	aux3=-(1.d0-eps(3))/4.d0

	a(3,1)=aux3
	a(3,2)=aux3
	a(3,4)=aux3
	a(3,5)=aux3

	a(4,4)=1.d0
	aux4=-(1.d0-eps(4))/4.d0

	a(4,1)=aux4
	a(4,2)=aux4
	a(4,3)=aux4
	a(4,5)=aux4

	a(5,5)=1.d0
	aux5=-(1.d0-eps(5))/4.d0

	a(5,1)=aux5
	a(5,2)=aux5
	a(5,3)=aux5
	a(5,4)=aux5

c	imprima a matriz a(i,j)
	write(*,*)'a matriz a(i,j) eh'
	write(2,*)'a matriz a(i,j) eh'
	do i=1,n
		write(*,*)'linha ',i
		write(2,*)'linha ',i
		write(2,*)(a(i,j),j=1,n)
		write(*,*)(a(i,j),j=1,n)
	enddo

c	calcule o vetor b(i)
	write(*,*)'o vetor b(i) eh'
	write(2,*)'o vetor b(i) eh'
	do i=1,n
		b(i)=eps(i)*sigma*temp(i)**4
		write(2,*)b(i)
		write(*,*)b(i)
	enddo


	  ! stop
c
c	fatoracao LU
c
	do k=1,n
		al(k,k)=dl
		if (k.gt.1) then
		sum=0.d0
		do is=1,k-1
			sum=sum+al(k,is)*au(is,k)
		enddo
		else
			sum=0.d0
		endif
		au(k,k)=(a(k,k)-sum)/al(k,k)
		do j=k+1,n
			if (k.gt.1) then
				sum=0.d0
				do is=1,k-1
					sum=sum+al(k,is)*au(is,j)
				enddo
			else
				sum=0.d0
			endif
			au(k,j)=(a(k,j)-sum)/al(k,k)
		enddo
		do i=k+1,n
			if (k.gt.1) then
				sum=0.d0
				do is=1,k-1
					sum=sum+al(i,is)*au(is,k)
				enddo
			else
				sum=0.d0
			endif
			al(i,k)=(a(i,k)-sum)/au(k,k)
		enddo
	enddo
c	imprima a matriz l(i,j)
	write(*,*)'a matriz l(i,j) eh'
	write(2,*)'a matriz l(i,j) eh'
	do i=1,n
		write(*,*)'linha ',i
		write(2,*)'linha ',i
		write(2,*)(al(i,j),j=1,n)
		write(*,*)(al(i,j),j=1,n)
	enddo
c	imprima a matriz u(i,j)
	write(*,*)'a matriz u(i,j) eh'
	write(2,*)'a matriz u(i,j) eh'
	do i=1,n
		write(*,*)'linha ',i
		write(2,*)'linha ',i
		write(2,*)(au(i,j),j=1,n)
		write(*,*)(au(i,j),j=1,n)
	enddo
c	verificacao
c	multiplicar l por u e obter a(i,j)
c
	write(*,*)'o produto de l por u eh a matriz p(i,j) abaixo'
	write(2,*)'o produto de l por u eh a matriz p(i,j) abaixo'
	do i=1,n
		do j=1,n
			p(i,j)=0.d0
			do is=1,n
				p(i,j)=p(i,j)+al(i,is)*au(is,j)
			enddo
		enddo
		write(2,*)'linha ',i
		write(*,*)'linha ',i
		write(2,*)(p(i,j),j=1,n)
		write(*,*)(p(i,j),j=1,n)
	enddo
c
c	forward substitution
c
	do i=1,n
		if(i.gt.1) then
			sum=0.d0
			do is=1,i-1
				sum=sum+al(i,is)*z(is)
			enddo
		else
			sum=0.d0
		endif
		z(i)=(b(i)-sum)/al(i,i)
	enddo
c
c	backward substitution
c
	do i=n,1,-1
		if(i.lt.n) then
			sum=0.d0
			do is=i+1,n
				sum=sum+au(i,is)*x(is)
			enddo
		else
			sum=0.d0
		endif
		x(i)=(z(i)-sum)/au(i,i)
	enddo
c
c	solucao
c
	write(*,*)'a solucao eh'
	write(2,*)'a solucao eh'
	do i=1,n
		write(*,*)'x(',i,')=',x(i)
		write(2,*)'x(',i,')=',x(i)
	enddo

	 ! Calculo dos Q's

	write(*,*)'as taxas de trans cal Q_i sao:'
	write(2,*)'as taxas de trans cal Q_i sao:'
	sum=0.d0
	do i=1,n
		q(i)=-area*eps(i)/(1.d0-eps(i))*(sigma*temp(i)**4-x(i))
		write(*,*)'q(',i,')=',q(i)
		write(2,*)'q(',i,')=',q(i)
		sum=sum+q(i)	
	enddo

	  ! verificacao de resultados por conservacao de energia (LEMA)

		write(*,*)'Soma dos Qs =',sum
		write(2,*)'Soma dos Qs =',sum

      stop
      end 
c----------------------------------