program mat
implicit none
real*8, dimension(:,:), allocatable :: a
real*8, dimension(:), allocatable :: b
integer :: n
  n = 5
  allocate(a(n,n),b(n))
  call matdef(a,b,n)
  call matwrt(a,b,n)
  call gauss_db(a,b,n)
  call matwrt(a,b,n)
  call verif(b,n)
  deallocate(a,b)
end

function xi(a,b,i,n)
! n = numero di punti
! i = 1..n
implicit none
! variabili di scambio
integer :: i,n
real*8 :: xi,a,b
! variabili locali
real*8, parameter :: pi=3.14159265358979d0
  xi=(a+b)/2 + ((b-a)/2)*cos((2*(i-1)+1)*pi/(2*n))
return
end

function runge(x)
implicit none
! variabili di scambio
real*8 :: runge, x
  runge = 1.d0/(1.d0+x*x)
return
end

subroutine matdef(a,b,n)
implicit none
! variabili di scambio
integer :: n
real*8, dimension(1:n,1:n) :: a
real*8, dimension(1:n) :: b
! variabili locali
integer :: i, j
real*8 :: x, runge, xi
!
!  costruisce matrice e vettore noto del sistema da risolvere
!
	do i=1,n
	  x = xi(-5.d0,5.d0,i,n)
	  do j=1,n
	    a(i,j) = x**(j-1)
	  enddo
	  b(i) = runge(x)
	enddo

return
end

subroutine verif(b,n)
implicit none
! variabili di scambio
integer :: n
real*8, dimension(1:n) :: b
! variabili locali
integer :: i, j
real*8 :: pn, x, xi, runge
!
!  stampa scarti calcolati sui nodi
!
	write(*,'(1x,t10,a,t29,a,t48,a,t67,a)')'x','pn','esatta','err rel.'
	do i=1,n
	  x = xi(-5.d0,5.d0,i,n)
	  pn = 0.d0
	  do j=1,n
	    pn = pn + b(j) * x**(j-1)
	  enddo
	  write(*,'(1x,4d19.12)') x,pn,runge(x),abs(runge(x)-pn)/runge(x)
	enddo

return
end


! swap_rows scambia le righe k_1 e k_2 
! della matrice A i corrispondenti elementi del
! vettore B.
! Da utilizzare nel pivoting parziale
subroutine swap_rows(a,b,n,k_1,k_2)
implicit none
! variabili di scambio
integer :: n, k_1, k_2
real*8, dimension(1:n,1:n) :: a
real*8, dimension(1:n) :: b
! variabili locali
real*8, dimension(1:n) :: swap_vector
real*8 :: swap_var
  swap_vector = a(k_1,1:n)
  a(k_1,1:n) = a(k_2,1:n)
  a(k_2,1:n) = swap_vector(1:n)
  swap_var = b(k_1)
  b(k_1) = b(k_2)
  b(k_2) = swap_var
return
end

subroutine gauss(a,b,n)
implicit none
! variabili di scambio
integer :: n
real*8, dimension(1:n,1:n) :: a
real*8, dimension(1:n) :: b
! variabili locali
real*8 :: factor, sum, v_max
integer :: k, i, j, k_1, k_2
!
  do k=1,n-1
!
	do i=k,n
	  if (a(i,k) /= 0.d0) then
		factor = 1.d0 / a(i,k)
		do j = k,n
		  a(i,j) = a(i,j) * factor
		enddo
		b(i) = b(i) * factor
	  endif
	enddo
!
	do i=k+1,n
	  if (a(i,k) /= 0.d0) then
		do j = k,n
		  a(i,j) = a(k,j) - a(i,j)
		enddo
		b(i) = b(k) - b(i)
	  endif
	enddo
!
! visualizza risultati parziali della triangolarizzazione
!
!	call matwrt(a,b,n)
!	pause
!
  enddo ! k=1..n-1
!
! call matwrt(a,b,n)
!
! sostituzione all'indietro
!
	do i = n,1,-1
	  sum = b(i)
	  do j = i+1,n
	    sum = sum - a(i,j)*b(j)
	  enddo
	  b(i) = sum/a(i,i)
	enddo
!
return
end

! eliminazione di gauss con pivoting parziale
! seguito da sostituzione all'indietro
! In ingresso a: matrice di coefficienti, b: vettore dati
! In uscita   a: matrice di coefficienti, b: soluzione 
! NB: algoritmo di triangolarizzazione come in Dahlquist-Bjorck
subroutine gauss_db(a,b,n)
implicit none
! variabili di scambio
integer :: n
real*8, dimension(1:n,1:n) :: a
real*8, dimension(1:n) :: b
! variabili locali
real*8 :: factor, sum, v_max
integer :: k, i, j, k_1, k_2
!
  do k=1,n-1
!
! pivoting parziale: per visualizzare l'evoluzione del pivoting,
! levare i commenti alle successive 5 istruzioni commentate
!
!   write(*,*) 'Visualizza sistema prima del pivoting parziale'
!	call matwrt(a,b,n)
!
    k_1 = k
	k_2 = k
    v_max = abs(a(k,k))
	do i=k+1,n
	  if (abs(a(i,k)) > v_max) then
	    v_max = a(i,k)
		k_2 = i
	  endif
	enddo
	if (k_1 /= k_2) then
	  call swap_rows(a,b,n,k_1,k_2)
	endif
!
!   write(*,*) 'Visualizza sistema dopo il pivoting parziale'
!	call matwrt(a,b,n)
!	pause 
!
	if (a(k,k) /= 0.d0) then
	  do i=k+1,n
		factor = a(i,k) / a(k,k)
		if (factor /= 0.d0) then
		  do j = k,n
			a(i,j) = a(i,j) - a(k,j) * factor
		  enddo
		  b(i) = b(i) - b(k) * factor
		endif
	  enddo
	endif
!
! visualizza risultati parziali della triangolarizzazione
!
 	call matwrt(a,b,n)
 	pause
!
  enddo ! k=1..n
!
! call matwrt(a,b,n)
!
! sostituzione all'indietro
!
	do i = n,1,-1
	  sum = b(i)
	  do j = i+1,n
	    sum = sum - a(i,j)*b(j)
	  enddo
	  b(i) = sum/a(i,i)
	enddo
!
return
end

subroutine matwrt(a,b,n)
implicit none
! variabili di scambio
integer :: n
real*8, dimension(1:n,1:n) :: a
real*8, dimension(1:n) :: b
! variabili locali
integer :: i, j
	
  do i=1,n
	write(*,'(1X,5f12.4)') (a(i,j), j=1,n)
  enddo

  do i=1,n
	write(*,'(1X,f12.4)') b(i)
  enddo

return
end

