PROGRAM sistema
  ! Risolve un sistema mediante fattorizzazione LU
  use lu_interface
  IMPLICIT NONE
  REAL(KIND(0.d0)), DIMENSION(:,:), ALLOCATABLE :: a,aux
  REAL(KIND(0.d0)), DIMENSION(:), ALLOCATABLE   :: b, baux
  REAL(KIND(0.d0))                              :: err
  INTEGER                                       :: n,i,j, test
  INTEGER,DIMENSION(:), ALLOCATABLE             :: r ,c

  write(*,*) "numero del test (1-4)"
  write(*,*)" 1-  Hankel"
  write(*,*)" 2-  Crescita esponenziale"
  write(*,*)" 3-  Hilbert"
  write(*,*)" 4-  Albero"
  read(*,*) test

  WRITE(*,*)"Dimensione matrice="
  READ(*,*)n

  ALLOCATE(a(n,n))
  ALLOCATE(b(n))


  if(test==1)then
     DO i=1,n
        DO j=1,n-i+1
           a(i,j)=2.d0**(1.d0/(n-j-i+2))
        END DO
        DO j=n-i+2,n
           a(i,j)=2.d0**(i+j-n)
        END DO
     END DO
  elseif(test==2) then
     do i=1,n
        do j=1,n
           if(i>j)then
              a(i,j)=-1.d0
           elseif(j==n)then
              a(i,j)=1.d0
           elseif(i==j)then
              a(i,j)=1.d0
           else
              a(i,j)=0.d0
           end if
        end do
     end do
  elseif(test==3) then
     do i=1,n
        do j=1,n
           a(i,j)=1.d0/(i+j-1)
        end do
     end do
  elseif (test==4) then
     a=0
     do i=1,n
        a(i,i)=2
     end do
     do i=2,n
        a(1,i)=i*2
        a(i,1)=n-i+1
     end do
     a(1,1)=1.d-4
  end if


  b=SUM(a,dim=2)
  ALLOCATE(aux(n,n))
  ALLOCATE(baux(n))
  aux=a
  baux=b

  WRITE(*,*)"Fattorizzazione LU"
  CALL lu(aux)
  CALL lusolve(aux,baux)
  err=0.d0
  DO i=1,n
     err=err+(baux(i)-1.d0)**2
  END DO
  err=SQRT(err/n)
  WRITE(*,*)"Errore=",err
  WRITE(*,*)

  WRITE(*,*)"Strategia del pivot parziale"
  aux=a
  baux=b
  ALLOCATE(r(n))
  CALL lupp(aux,r)
  CALL luppsolve(aux,r,baux)
  err=0.d0
  DO i=1,n
     err=err+(baux(i)-1.d0)**2
  END DO
  err=SQRT(err/n)
  WRITE(*,*)"Errore=",err
  WRITE(*,*)


  WRITE(*,*)"Strategia del pivot totale"
  ALLOCATE(c(n))
  aux=a
  baux=b
  CALL lupt(aux,r,c)
  CALL luptsolve(aux,r,c,baux)
  err=0.d0
  DO i=1,n
     err=err+(baux(i)-1.d0)**2
  END DO
  err=SQRT(err/n)
  WRITE(*,*)"Errore=",err

END PROGRAM sistema




