Source

PetIGA / src / petigainv.f90.in

Full commit
! -*- f90 -*-

pure function Determinant(dim,A) result (detA)
  use PetIGA
  implicit none
  integer(kind=IGA_INTEGER_KIND), intent(in) :: dim
  real   (kind=IGA_REAL_KIND   ), intent(in) :: A(dim,dim)
  real   (kind=IGA_REAL_KIND   )  :: detA
  select case (dim)
  case (1)
     detA = A(1,1)
  case (2)
     detA = + A(1,1)*A(2,2) - A(2,1)*A(1,2)
  case (3)
     detA = + A(1,1) * ( A(2,2)*A(3,3) - A(3,2)*A(2,3) ) &
            - A(2,1) * ( A(1,2)*A(3,3) - A(3,2)*A(1,3) ) &
            + A(3,1) * ( A(1,2)*A(2,3) - A(2,2)*A(1,3) )
  case default
     detA = 0
  end select
end function Determinant

pure function Inverse(dim,detA,A) result (invA)
  use PetIGA
  implicit none
  integer(kind=IGA_INTEGER_KIND), intent(in) :: dim
  real   (kind=IGA_REAL_KIND   ), intent(in) :: detA
  real   (kind=IGA_REAL_KIND   ), intent(in) :: A(dim,dim)
  real   (kind=IGA_REAL_KIND   )  :: invA(dim,dim)
  select case (dim)
  case (1)
     invA = 1/detA
  case (2)
     invA(1,1) = + A(2,2)
     invA(2,1) = - A(2,1)
     invA(1,2) = - A(1,2)
     invA(2,2) = + A(1,1)
     invA = invA/detA
  case (3)
     invA(1,1) = + A(2,2)*A(3,3) - A(2,3)*A(3,2)
     invA(2,1) = - A(2,1)*A(3,3) + A(2,3)*A(3,1)
     invA(3,1) = + A(2,1)*A(3,2) - A(2,2)*A(3,1)
     invA(1,2) = - A(1,2)*A(3,3) + A(1,3)*A(3,2)
     invA(2,2) = + A(1,1)*A(3,3) - A(1,3)*A(3,1)
     invA(3,2) = - A(1,1)*A(3,2) + A(1,2)*A(3,1)
     invA(1,3) = + A(1,2)*A(2,3) - A(1,3)*A(2,2)
     invA(2,3) = - A(1,1)*A(2,3) + A(1,3)*A(2,1)
     invA(3,3) = + A(1,1)*A(2,2) - A(1,2)*A(2,1)
     invA = invA/detA
  case default
     invA = 0
  end select
end function Inverse