umat子程序编写常用的fortran函数分享（三）

计算3*3矩阵的逆矩阵

subroutine inv3x3(A,invA,det)
implicit none
real(8), intent(in)  :: A(3,3)
real(8), intent(out) :: invA(3,3), det
integer :: i,j
call deter3x3(A,det)
if (abs(det) < 1e-6) then
    invA=0.0d+0
else
    invA(1,1)=((A(2,2)*A(3,3))-(A(2,3)*A(3,2)))/det
    invA(2,1)=-((A(2,1)*A(3,3))-(A(2,3)*A(3,1)))/det
    invA(3,1)=((A(2,1)*A(3,2))-(A(2,2)*A(3,1)))/det
    invA(1,2)=-((A(1,2)*A(3,3))-(A(1,3)*A(3,2)))/det
    invA(2,2)=((A(1,1)*A(3,3))-(A(1,3)*A(3,1)))/det
    invA(3,2)=-((A(1,1)*A(3,2))-(A(1,2)*A(3,1)))/det
    invA(1,3)=((A(1,2)*A(2,3))-(A(1,3)*A(2,2)))/det
    invA(2,3)=-((A(1,1)*A(2,3))-(A(2,1)*A(1,3)))/det
    invA(3,3)=((A(1,1)*A(2,2))-(A(1,2)*A(2,1)))/det
endif
return
end subroutine inv3x3

计算2*2矩阵的逆

subroutine inv2x2(A,invA,det)
implicit none
real(8), intent(in)  :: A(2,2)
real(8), intent(out) :: invA(2,2), det
integer :: i,j
call deter2x2(A,det)
if (abs(det) < 1e-6) then
  invA=0.0d+0
else
  invA(1,1) = A(2,2)/det
  invA(1,2) =-A(1,2)/det
  invA(2,1) =-A(2,1)/det
  invA(2,2) = A(1,1)/det
endif
return
end subroutine inv2x2

向量单位化

subroutine normalize_vector(x,y,z)
include 'ABA_PARAM.INC'           
xlength = sqrt(x*x+y*y+z*z)      
x = x / xlength     
y = y / xlength     
z = z / xlength     
return              
end

矩阵转置

subroutine transpose(n,a,b)      
include 'ABA_PARAM.INC'           
dimension a(n,n), b(n,n)         
do i = 1,n          
  do j = 1,n       
    b(i,j) = a(j,i)            
  end do           
end do              
return              
en