module aln
  implicit none
contains
  subroutine aelup(m,trans)
    !   entree : m de taille nxn contenant A
    !   sortie : m contenant (les coefs utiles de) L et U
    !            trans contenant la suite des (n-1) transpositions de colonnes
    !                  correspondant a la permutation de matrice P
    !   aelup calcul la decomposition A=LUP avec L triangulaire inferieure
    !   avec des 1 sur la diagonale, U triangulaire superieure, P de permutation
    real, dimension(:,:), intent(inout) :: m
    integer, dimension(size(m,2)-1), intent(out) :: trans
    integer, dimension(size(m,2)) :: perm
    integer :: n,j,k,ind_piv
    real :: r,piv
    n=size(m,1)      ! m doit etre carree
    perm=(/(k,k=1,n)/) ! initialise la permutation a l'identite : perm(k)=k
    do k=1,n-1
       ind_piv=k
       piv=abs(m(k,perm(ind_piv)))
       do j=k+1,n
          r=abs(m(k,perm(j)))
          if (r>piv) then
             ind_piv=j
             piv=r
          end if
       end do
       trans(k)=ind_piv        ! enregistre l'indice du j-eme pivot
       if (ind_piv/=k) then    ! echange perm(k) et perm(ind_piv)
          j=perm(k);perm(k)=perm(ind_piv);perm(ind_piv)=j                  
       end if
       m(k+1:n,perm(k))=m(k+1:n,perm(k))/m(k,perm(k))
       do j=k+1,n
          m(k+1:n,perm(j))=m(k+1:n,perm(j))-m(k+1:n,perm(k))*m(k,perm(j))
       end do
    end do
  end subroutine aelup
  subroutine lupxeb(m,trans,b)
    real, dimension(:,:), intent(in) :: m
    integer, dimension(size(m,2)-1), intent(out) :: trans
    real, dimension(size(m,1)), intent(inout) :: b
    integer, dimension(size(m,2)) :: perm
    integer :: n,k,j,ind_piv
    real :: r
    n=size(m,1)      ! m doit etre carree
    perm=(/(k,k=1,n)/)
    do k=1,n-1
       ind_piv=trans(k)       ! enregistre l'indice du k-eme pivot
       if (ind_piv/=k) then   ! echange perm(k) et perm(ind_piv)
          j=perm(k);perm(k)=perm(ind_piv);perm(ind_piv)=j         
       end if
    end do
    ! resolution Lz=b
    do k=2,n 
       b(k:n)=b(k:n)-m(k:n,perm(k-1))*b(k-1)
    end do
    ! resolution Uy=z
    b(n)=b(n)/m(n,perm(n))
    do k=n-1,1,-1
       b(1:k)=b(1:k)-m(1:k,perm(k+1))*b(k+1)
       b(k)=b(k)/m(k,perm(k))
    end do
    ! resolution de Px=y
    do k=n-1,1,-1
       ind_piv=trans(k)          ! enregistre l'indice du k-eme pivot
       if (ind_piv/=k) then      ! echange b(k) et b(ind_piv)
          r=b(k);b(k)=b(ind_piv);b(ind_piv)=r
       end if
    end do
  end subroutine lupxeb
  function pmv(a,x)
    real, dimension(:,:), intent(in) :: a
    real, dimension(size(a,2)), intent(in) :: x
    real, dimension(size(a,1)) :: pmv
    integer :: n,j
    pmv=0.0
    n=size(a,2)
    do j=1,n
       pmv=pmv+a(:,j)*x(j)
    end do
  end function pmv
end module aln
