SUBROUTINE lu(a)
  ! Calcola la fattorizzazione LU di A
  IMPLICIT NONE
  INTEGER                          :: n, i, k
  REAL(KIND(0.d0)), DIMENSION(:,:) :: a
  n=SIZE(a,1)
  elimina: DO k=1,n-1
     ! Aggiorna la matrice L
     a(k+1:n,k)=a(k+1:n,k)/a(k,k)
     ! Aggiorna la matrice U
     DO i=k+1,n
        a(i,k+1:n)=a(i,k+1:n)-a(i,k)*a(k,k+1:n)
     END DO
  END DO elimina
END SUBROUTINE lu


SUBROUTINE lusolve(a,b)
  ! A contiene le matrici L e U
  ! b contiene il termine noto (input) e la soluzione (output)
  IMPLICIT NONE
  INTEGER :: n, i, j
  REAL(KIND(0.d0)), DIMENSION(:,:):: a
  REAL(KIND(0.d0)), DIMENSION(:)  :: b
  REAL(KIND(0.d0)), DIMENSION(:),allocatable  :: y
  n=SIZE(a,1)
  ALLOCATE(y(n))
  y(1)=b(1)
  lower: DO i=2,n
     y(i)=b(i)
     DO j=1,i-1
        y(i)=y(i)-a(i,j)*y(j)
     END DO
  END DO lower
  b(n)=y(n)/a(n,n)
  upper: DO i=n-1,1,-1
     b(i)=y(i)
     DO j=i+1,n
        b(i)=b(i)-a(i,j)*b(j)
     END DO
     b(i)=b(i)/a(i,i)
  END DO upper
END SUBROUTINE lusolve

SUBROUTINE lupp(a,r)
  ! Calcola la fattorizzazione LU di A con la strategia
  ! del pivot parziale
  IMPLICIT NONE
  INTEGER                          :: n, i, k, piv
  REAL(KIND(0.d0)), DIMENSION(:,:) :: a
  INTEGER, DIMENSION(:)            :: r ! Vettore per lo scambio
  ! delle righe
  n=SIZE(a,1)
  DO i=1,n
     r(i)=i
  END DO

  elimina: DO k=1,n-1

     ! Cerca l'indice di pivot
     piv=MINVAL(MAXLOC(ABS(a(r(k:n),k))))+k-1

     ! Scambia le righe piv e k
     IF(piv/=k)THEN
        i=r(piv)
        r(piv)=r(k)
        r(k)=i
     END IF

     ! Aggiorna la matrice L
     a(r(k+1:n),k)=a(r(k+1:n),k)/a(r(k),k)

     ! Aggiorna la matrice U
     DO i=k+1,n
        a(r(i),k+1:n)=a(r(i),k+1:n)-a(r(i),k)*a(r(k),k+1:n)
     END DO

  END DO elimina

END SUBROUTINE lupp

SUBROUTINE luppsolve(a,r,b)
  IMPLICIT NONE
  INTEGER :: n, i, j
  REAL(KIND(0.d0)), DIMENSION(:,:):: a
  REAL(KIND(0.d0)), DIMENSION(:)  :: b
  REAL(KIND(0.d0)), DIMENSION(:), allocatable  :: y
  INTEGER, DIMENSION(:)           :: r
  n=SIZE(a,1)
  ALLOCATE(y(n))
  y(1)=b(r(1))
  lower: DO i=2,n
     y(i)=b(r(i))
     DO j=1,i-1
        y(i)=y(i)-a(r(i),j)*y(j)
     END DO
  END DO lower
  b(n)=y(n)/a(r(n),n)
  upper: DO i=n-1,1,-1
     b(i)=y(i)
     DO j=i+1,n
        b(i)=b(i)-a(r(i),j)*b(j)
     END DO
     b(i)=b(i)/a(r(i),i)
  END DO upper
END SUBROUTINE luppsolve


SUBROUTINE lupt(a,r,c)
  ! Calcola la fattorizzazione LU di A con la strategia
  ! del pivot totale
  IMPLICIT NONE
  INTEGER                          :: n, i, k, j
  INTEGER, DIMENSION(2)            :: piv
  REAL(KIND(0.d0)), DIMENSION(:,:) :: a
  REAL(KIND(0.d0))                 :: xpiv
  INTEGER, DIMENSION(:)            :: r, c ! Vettori per lo scambio
  ! delle righe e colonne
  n=SIZE(a,1)
  DO i=1,n
     r(i)=i
     c(i)=i
  END DO

  elimina: DO k=1,n-1

     ! Cerca l'indice di pivot
     !     piv=minval(MAXLOC(ABS(a(r(k:n),c(k:n)))))+k-1
     piv(1)=k; piv(2)=k
     xpiv=ABS(a(r(k),c(k)))
     DO i=k,n
        DO j=k,n
           IF(ABS(a(r(i),c(j)))>xpiv)THEN
              piv(1)=i
              piv(2)=j
              xpiv=ABS(a(r(i),c(j)))
           END IF
        END DO
     END DO

     ! Scambia righe e colonne
     IF(piv(1)/=k)THEN
        i=r(piv(1))
        r(piv(1))=r(k)
        r(k)=i
     END IF
     IF(piv(2)/=k)THEN
        i=c(piv(2))
        c(piv(2))=c(k)
        c(k)=i
     END IF

     ! Aggiorna la matrice L
     a(r(k+1:n),c(k))=a(r(k+1:n),c(k))/a(r(k),c(k))

     ! Aggiorna la matrice U
     DO i=k+1,n
        a(r(i),c(k+1:n))=a(r(i),c(k+1:n))-a(r(i),c(k))*a(r(k),c(k+1:n))
     END DO

  END DO elimina

END SUBROUTINE lupt

SUBROUTINE luptsolve(a,r,c,b)
  IMPLICIT NONE
  INTEGER :: n, i, j
  REAL(KIND(0.d0)), DIMENSION(:,:):: a
  REAL(KIND(0.d0)), DIMENSION(:)  :: b
  REAL(KIND(0.d0)), DIMENSION(:), allocatable  :: y
  INTEGER, DIMENSION(:)           :: r,c
  n=SIZE(a,1)
  ALLOCATE(y(n))
  y(1)=b(r(1))
  lower: DO i=2,n
     y(i)=b(r(i))
     DO j=1,i-1
        y(i)=y(i)-a(r(i),c(j))*y(j)
     END DO
  END DO lower
  b(c(n))=y(n)/a(r(n),c(n))
  upper: DO i=n-1,1,-1
     b(c(i))=y(i)
     DO j=i+1,n
        b(c(i))=b(c(i))-a(r(i),c(j))*b(c(j))
     END DO
     b(c(i))=b(c(i))/a(r(i),c(i))
  END DO upper
END SUBROUTINE luptsolve









