c      program directsolver
c	programa principal para resolver eq lineares - metodo direto
c
      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)
      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
	write(*,*)'a matriz a(i,j) eh:'
	do i=1,n
		read(1,*) (a(i,j),j=1,n)
		write(*,*)(a(i,j),j=1,n)
	enddo
c
	write(*,*)'o vetor b(i) eh:'
	write(2,*)'o vetor b(i) eh:'
	read (1,*) (b(i),i=1,n)
      write(*,*)'b(i)=', (b(i),i=1,n)
	write(2,*)'b(i)=', (b(i),i=1,n)
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)
		au(k,k)=sqrt((a(k,k)-sum)) 
		al(k,k)=au(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
      stop
      end 
c----------------------------------