femwind/fortran/module_lin_alg.f90

   1 module module_lin_alg
   2
   3 use module_utils
   4
   5 contains
   6
   7 subroutine inv3(M, M_inv)
   8 !Purpose: Calculate the inverse of a 3X3 matrix
   9 !In:
  10 !M  3X3 matrix
  11 !Out:
  12 !M_inv Inverse of M
  13 implicit none
  14 real,intent(in), dimension(3,3):: M
  15 real,intent(out), dimension(3,3):: M_inv
  16
  17 !Local Variables
  18
  19 real :: detM
  20
  21
  22     ! Compute M_inv
  23     detM  = (M(1,1)*M(2,2)*M(3,3) - M(1,1)*M(2,3)*M(3,2)-&
  24              M(1,2)*M(2,1)*M(3,3) + M(1,2)*M(2,3)*M(3,1)+&
  25              M(1,3)*M(2,1)*M(3,2) - M(1,3)*M(2,2)*M(3,1))
  26
  27     if(abs(detM).lt.10.0*tiny(detM))then
  28         print *,'detM=',detM
  29         call crash('The matrix is numerically singular')
  30     endif
  31
  32     detM = 1.0/detM
  33
  34     M_inv(1,1) = +detM * (M(2,2)*M(3,3) - M(2,3)*M(3,2))
  35     M_inv(2,1) = -detM * (M(2,1)*M(3,3) - M(2,3)*M(3,1))
  36     M_inv(3,1) = +detM * (M(2,1)*M(3,2) - M(2,2)*M(3,1))
  37     M_inv(1,2) = -detM * (M(1,2)*M(3,3) - M(1,3)*M(3,2))
  38     M_inv(2,2) = +detM * (M(1,1)*M(3,3) - M(1,3)*M(3,1))
  39     M_inv(3,2) = -detM * (M(1,1)*M(3,2) - M(1,2)*M(3,1))
  40     M_inv(1,3) = +detM * (M(1,2)*M(2,3) - M(1,3)*M(2,2))
  41     M_inv(2,3) = -detM * (M(1,1)*M(2,3) - M(1,3)*M(2,1))
  42     M_inv(3,3) = +detM * (M(1,1)*M(2,2) - M(1,2)*M(2,1))
  43
  44 end subroutine inv3
  45
  46 subroutine print_matrix(name,A)
  47 character(len=*)::name
  48 real :: A(:,:)
  49 integer i,j
  50 print *,'Matrix ',trim(name),lbound(A,1),':',ubound(A,1),' by ', &
  51                           lbound(A,2),':',ubound(A,2)
  52 do i=lbound(A,1),ubound(A,1)
  53     print *,(A(i,j),j=lbound(A,2),ubound(A,2))
  54 enddo
  55 end subroutine print_matrix
  56
  57 end module module_lin_alg