program triangolare
  implicit none
  ! Programma che calcola le soluzioni di un sistema triangolare

  real(8),allocatable,dimension(:,:) :: a
  real(8),allocatable,dimension(:)   :: x,b
  real(8) :: r
  integer :: n,i,j

  write(*,'(a)',advance='no') "Inserisci la dimensione del sistema "
  read (*,*) n
  allocate(a(n,n),b(n),x(n))

  ! Inizializzo la matrice 
  a=0.d0
  do i=1,n
     do j=1,i  ! ALT: j=i,n
        a(i,j)=1.d0/(i+j)
     end do
  end do
  ! Inizializzo un termine noto che dia come soluzione (1,1,...,1)^T
  b=sum(a,dim=2)

  ! Calcolo la soluzione
  call lsolve(n,a,x,b) ! ALT: call usolve(n,a,x,b)

  ! Scrivo il risultato
  write(*,*) "L'errore e` ",maxval(abs(x-1.d0))

end program triangolare

subroutine lsolve (n,a,x,b)
  implicit none
  
  integer :: n,i,j
  real(8),dimension(n,n) :: a
  real(8),dimension(n)   :: x,b
  real(8) :: s

  do i=1,n
!     s=0                     ! Versione lenta
!     do j=1,i-1              ! lenta
!        s=s+x(j)*a(i,j)      ! lenta
!     end do                  ! lenta
     s=sum(x(1:i-1)*a(i,1:i-1))  ! Versione veloce
     x(i)=1/a(i,i)*(b(i)-s)
  end do
  
end subroutine lsolve

subroutine usolve (n,a,x,b)
  implicit none
  
  integer :: n,i,j
  real(8),dimension(n,n) :: a
  real(8),dimension(n)   :: x,b
  real(8) :: s

  do i=n,1,-1
     s=0
     do j=i+1,n
        s=s+x(j)*a(i,j)
     end do
     x(i)=1/a(i,i)*(b(i)-s)
  end do

end subroutine usolve
