PROGRAM sistema
  ! Programma principale per il calcolo del determinante di una matrice
  ! e per la risoluzione di un sistema con la regola di Cramer
  ! versione che usa l' interface
  IMPLICIT NONE

  ! Interfaccia, che permette di passare implicitamente la dimensione 
  ! delle variabili multidimensionali

  INTERFACE

     RECURSIVE FUNCTION det(a) RESULT(risultato) 
       REAL(KIND(0.d0)) :: risultato
       REAL(KIND(0.d0)), DIMENSION(:,:) :: a
     END FUNCTION det

     SUBROUTINE cramer(a,b,x)
       REAL(KIND(0.d0)), DIMENSION(:,:) :: a
       REAL(KIND(0.d0)), DIMENSION(:)   :: b, x
     END SUBROUTINE cramer

  END INTERFACE

  INTEGER :: n, i, j
  REAL(KIND(0.d0)), ALLOCATABLE, DIMENSION(:,:) :: a
  REAL(KIND(0.d0)), ALLOCATABLE, DIMENSION(:) :: b, x
  REAL(KIND(0.d0)) :: d, err

  WRITE(*,*)"Dimensione della matrice"
  READ(*,*)n
  ALLOCATE(a(n,n))
  ALLOCATE(b(n),x(n))

  ! Assegna i valori alla matrice
  a=0.d0
  DO i=1,n
     DO j=1,n
        IF(i==j)a(i,j)=2.d0
        IF(ABS(i-j)==1)a(i,j)=1.d0
     END DO
  END DO

  ! Assegna i valori al termine noto in modo che la soluzione
  ! del sitema Ax=b sia il vettore x=(1,1,...1)'
  b = SUM(a, dim=2)

  ! Calcola il determinante
  d=det(a)
  WRITE(*,*)"Determinante=", d

  ! Risolve il sistema Ax=b
  CALL cramer(a,b,x)
  WRITE(*,*)"Soluzione=",x

  ! Calcola l'errore relativo della soluzione in norma 2
  err = SQRT( SUM( (x-1.d0)**2  )/(1.d0*n))
  WRITE(*,*)"Errore relativo=", err
  
END PROGRAM sistema


RECURSIVE FUNCTION det(a) RESULT(risultato) 
! Calcola il determinante della matrice a mediante la regola di Laplace,
! sviluppando il determinante lungo la prima riga
  IMPLICIT NONE
  INTEGER :: n
  REAL(KIND(0.d0)), DIMENSION(:,:) :: a
  REAL(KIND(0.d0)), DIMENSION(:,:), ALLOCATABLE :: b
  INTEGER :: i
  REAL(KIND(0.d0)) :: risultato,s

  n=SIZE(a,1)
  ALLOCATE(b(n-1,n-1))
  IF(n==1)THEN
     risultato=a(1,1)
  ELSE
     s=0.d0
     DO i=1,n
        b(1:n-1,1:i-1)=a(2:n,1:i-1)
        b(1:n-1,i:n-1)=a(2:n,i+1:n)
        s=s+(-1)**(i+1)*a(1,i)*det(b)
     END DO
     risultato=s
  END IF
END FUNCTION det

SUBROUTINE cramer(a,b,x)
  ! Risolve il sistema Ax=b con la regola di Cramer
  IMPLICIT NONE

  INTERFACE
     RECURSIVE FUNCTION det(a) RESULT(risultato) 
       REAL(KIND(0.d0)) :: risultato
       REAL(KIND(0.d0)), DIMENSION(:,:) :: a
     END FUNCTION det
  END INTERFACE

  REAL(KIND(0.d0)), DIMENSION(:,:) :: a
  REAL(KIND(0.d0)), DIMENSION(:)   :: b, x
  REAL(KIND(0.d0)), DIMENSION(:,:), ALLOCATABLE :: w
  REAL(KIND(0.d0)) :: d
  INTEGER :: n, i

  n=SIZE(a,1)
  ALLOCATE(w(n,n))
  d=det(a)
  DO i=1,n
     w=a
     w(:,i)=b
     x(i)=det(w)/d
  END DO

END SUBROUTINE cramer

