femwind/fortran/module_lin_alg.f90

   1 module module_lin_alg
   2
   3 contains
   4
   5 subroutine inv3(M, Inv_M)
   6 !Purpose: Calculate the inverse of a 3X3 matrix
   7 !In:
   8 !M  3X3 matrix
   9 !Out:
  10 !M_inv Inverse of M
  11 implicit none
  12 real,intent(in), dimension(3,3):: M
  13 real,intent(out), dimension(3,3):: Inv_M
  14
  15 !Local Variables
  16 real, dimension(3,3)::M_T
  17 real :: det_M, det_M_inv
  18
  19 !!Compute Inverse of M!!
  20 det_M =   M(1,1)*(M(2,2)*M(3,3)-M(3,2)*M(2,3)) -  &
  21           M(1,2)*(M(1,2)*M(3,3) - M(3,1)*M(2,3))+ &
  22           M(1,3)*(M(2,1)*M(2,2) - M(3,1)*M(3,2))
  23
  24 det_M_inv = 1/det_M
  25
  26 M_T   =  transpose(M)
  27
  28 Inv_M(1,1) = det_M_inv*(M_T(2,2)*M_T(3,3) - M_T(2,3)*M_T(3,2))
  29 Inv_M(1,2) = det_M_inv*(M_T(2,1)*M_T(3,3) - M_T(3,1)*M_T(2,3))
  30 Inv_M(1,3) = det_M_inv*(M_T(2,1)*M_T(3,2) - M_T(2,2)*M_T(3,1))
  31 Inv_M(2,1) = det_M_inv*(M_T(1,2)*M_T(3,3) - M_T(3,2)*M_T(1,3))
  32 Inv_M(2,2) = det_M_inv*(M_T(1,1)*M_T(3,3) - M_T(3,1)*M_T(1,3))
  33 Inv_M(2,3) = det_M_inv*(M_T(1,1)*M_T(3,2) - M_T(3,1)*M_T(1,2))
  34 Inv_M(3,1) = det_M_inv*(M_T(1,2)*M_T(2,3) - M_T(1,3)*M_T(2,2))
  35 Inv_M(3,2) = det_M_inv*(M_T(1,1)*M_T(2,3) - M_T(2,1)*M_T(1,3))
  36 Inv_M(3,3) = det_M_inv*(M_T(1,1)*M_T(2,2) - M_T(1,2)*M_T(2,1))
  37
  38 end subroutine inv3
  39 end module module_lin_alg