Example of LU factorization using Fortran parameterized derived types

lapack_interface.f90

module lapack

  implicit none

  interface
    subroutine sgetrf(m,n,a,lda,ipiv,info)
      integer, intent(in) :: m, n, lda
      real, intent(inout) :: a(lda,*)
      integer, intent(out) :: ipiv(*)
      integer, intent(out) :: info
    end subroutine
    subroutine dgetrf
      integer, intent(in) :: m, n, lda
      double precision, intent(inout) :: a(lda,*)
      integer, intent(out) :: ipiv(*)
      integer, intent(out) :: info
    end subroutine
  end interface

end module

lu_pdt.f90

module lu_pdt

  implicit none 
  private

  public :: lu_workspace, factorize

  type :: lu_workspace(wp,n)
    integer, kind :: wp
    integer, len :: n
    real(wp) :: a(n,n)
    real(wp) :: b(n)
    integer :: ipiv(n)
    logical :: factorized = .false.
  end type

  interface factorize
    module procedure factorize_sp
    module procedure factorize_dp
  end interface

  integer, parameter :: sp = kind(1.0e0)
  integer, parameter :: dp = kind(1.0d0)

contains

  subroutine factorize_sp(this,info)
    use lapack, only: lapack_factorize => sgetrf
    type(lu_workspace(sp,*)), intent(inout) :: this
    integer, intent(out), optional :: info  
    include "lu_pdt.inc"
  end subroutine

  subroutine factorize_dp(this,info)
    use lapack, only: lapack_factorize => dgetrf
    type(lu_workspace(dp,*)), intent(inout) :: this
    integer, intent(out), optional :: info  
    include "lu_pdt.inc"
  end subroutine

end module

lu_pdt.inc

    integer :: info_

    if (.not. this%factorized) then
      call lapack_factorize(this%n,this%n,this%a,this%n,this%ipiv,info_)
      if (info_ == 0) then
        this%factorized = .true.
      end if
      if (present(info)) info = info_
    else
      return
    end if

main.f90

program main

  use lu_pdt
  implicit none

  integer, parameter :: sp = kind(1.0e0)
  integer, parameter :: dp = kind(1.0d0)

  type(lu_workspace(dp,:)), allocatable :: work
  integer :: info

  allocate(lu_workspace(dp,n=3) :: work)

  work%a = reshape(&
    [real(dp) :: 4, 2, -1, -3, 1, 2, 1, 3, -5], &
    [3,3])

  work%b = [real(dp) :: -10, 0, 17]

  call factorize(work,info)
  print *, "info = ", info

end program

"Examle" should be "Example".

Thanks, noted.