Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Example of LU factorization using Fortran parameterized derived types
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
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
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
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
@Beliavsky
Copy link

"Examle" should be "Example".

@ivan-pi
Copy link
Author

ivan-pi commented Mar 22, 2022

Thanks, noted.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment