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 the determinant is greater than certain value

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              

          end