!===================================================================================== module det ! �s�v�Z���W���[���u���l�v�Z�̂��߂�Fortran90/95�v���O���~���O����v�����p implicit none contains recursive function det_mat(a, n) result(det) integer, intent(in) :: n real(8), intent(in) :: a(n, n) real(8) det, b(n - 1, n - 1) integer i if (n > 1) then det = 0.0d0 do i = 1, n b(1:i - 1, 1:n - 1) = a(1:i - 1, 2:n) b(i:n - 1, 1:n - 1) = a(i + 1:n, 2:n) det = det + (-1.0d0)**(i + 1) & *a(i, 1)*det_mat(b, n - 1) end do else det = a(1, 1) end if end function det_mat !===================================================================================== subroutine trans_rank_2(A, A_T) real(8), intent(in)::A(:, :) real(8), allocatable, intent(out)::A_T(:, :) integer n, m, i, j n = size(A, 1) m = size(A, 2) if (.not. allocated(A_T)) allocate (A_T(m, n)) do i = 1, n do j = 1, m A_T(j, i) = A(i, j) end do end do end subroutine trans_rank_2 !================================================================================== function trans1(A) result(A_T) real(8), intent(in)::A(:) real(8), allocatable::A_T(:, :) integer n, m, i, j n = size(A) if (.not. allocated(A_T)) allocate (A_T(1, n)) do i = 1, n A_T(1, i) = A(i) end do end function trans1 !================================================================================== function trans2(A) result(A_T) real(8), intent(in)::A(:, :) real(8), allocatable::A_T(:, :) integer n, m, i, j n = size(A, 1) m = size(A, 2) if (.not. allocated(A_T)) allocate (A_T(m, n)) do i = 1, n do j = 1, m A_T(j, i) = A(i, j) end do end do end function trans2 !================================================================================== subroutine inverse_rank_2(A, A_inv) real(8), intent(in)::A(:, :) real(8), allocatable::A_inv(:, :) real(8) detA, detA_1 integer m, n m = size(A, 1) n = size(A, 2) if (.not. allocated(A_inv)) allocate (A_inv(m, n)) detA = det_mat(A, n) if (detA == 0.0d0) stop "ERROR: inverse, detA=0" detA_1 = 1.0d0/detA if (n == 2) then A_inv(1, 1) = detA_1*A(2, 2) A_inv(1, 2) = -detA_1*A(1, 2) A_inv(2, 1) = -detA_1*A(2, 1) A_inv(2, 2) = detA_1*A(1, 1) elseif (n == 3) then A_inv(1, 1) = detA_1*(A(2, 2)*A(3, 3) - A(2, 3)*A(3, 2)) A_inv(1, 2) = detA_1*(A(1, 3)*A(3, 2) - A(1, 2)*A(3, 3)) A_inv(1, 3) = detA_1*(A(1, 2)*A(2, 3) - A(1, 3)*A(2, 2)) A_inv(2, 1) = detA_1*(A(2, 3)*A(3, 1) - A(2, 1)*A(3, 3)) A_inv(2, 2) = detA_1*(A(1, 1)*A(3, 3) - A(1, 3)*A(3, 1)) A_inv(2, 3) = detA_1*(A(1, 3)*A(2, 1) - A(1, 1)*A(2, 3)) A_inv(3, 1) = detA_1*(A(2, 1)*A(3, 2) - A(2, 2)*A(3, 1)) A_inv(3, 2) = detA_1*(A(1, 2)*A(3, 1) - A(1, 1)*A(3, 2)) A_inv(3, 3) = detA_1*(A(1, 1)*A(2, 2) - A(1, 2)*A(2, 1)) else print *, "ERROR: Aij with i=j=", n, "/=2or3" end if end subroutine inverse_rank_2 !================================================================================== subroutine tensor_exponential(A, expA, TOL, itr_tol) real(8), intent(in)::A(:, :), TOL real(8), allocatable, intent(inout)::expA(:, :) real(8), allocatable::increA(:, :) integer, intent(in)::itr_tol real(8) increment, NN integer i, j, n if (.not. allocated(expA)) allocate (expA(size(A, 1), size(A, 2))) allocate (increA(size(A, 1), size(A, 2))) if (size(A, 1) /= size(A, 2)) stop "ERROR:tensor exp is not a square matrix" expA(:, :) = 0.0d0 do n = 1, size(expA, 1) expA(n, n) = 1.0d0 end do NN = 1.0d0 increA(:, :) = expA(:, :) do n = 1, itr_tol if (n > 1) then NN = NN*(NN + 1.0d0) end if increA(:, :) = matmul(increA, A) expA(:, :) = expA(:, :) + 1.0d0/NN*increA(:, :) increment = 0.0d0 do i = 1, size(A, 1) do j = 1, size(A, 2) increment = increment + 1.0d0/NN*increA(i, j)*increA(i, j) end do end do if (increment <= TOL) then exit else if (n >= itr_tol) then stop "tensor exponential is not converged" end if cycle end if end do deallocate (increA) end subroutine tensor_exponential !================================================================================== subroutine tensor_expo_der(A, expA_A, TOL, itr_tol) real(8), intent(in)::A(:, :), TOL real(8), allocatable, intent(inout)::expA_A(:, :, :, :) real(8), allocatable::increA_1(:, :), increA_2(:, :), increA_3(:, :, :, :), I_ij(:, :), A_inv(:, :) integer, intent(in)::itr_tol real(8) increment, NN integer i, j, k, l, n, m, o if (.not. allocated(expA_A)) allocate (expA_A(size(A, 1), size(A, 1), size(A, 1), size(A, 1))) allocate (I_ij(size(A, 1), size(A, 1))) allocate (increA_1(size(A, 1), size(A, 1))) allocate (increA_2(size(A, 1), size(A, 1))) allocate (increA_3(size(A, 1), size(A, 1), size(A, 1), size(A, 1))) if (size(A, 1) /= size(A, 2)) stop "ERROR:tensor exp is not a square matrix" call inverse_rank_2(A, A_inv) I_ij(:, :) = 0.0d0 do n = 1, size(expA_A, 1) I_ij(n, n) = 1.0d0 end do NN = 1.0d0 do i = 1, size(A, 1) do j = 1, size(A, 1) do k = 1, size(A, 1) do l = 1, size(A, 1) expA_A(i, j, k, l) = I_ij(i, k)*I_ij(l, j) end do end do end do end do increA_1(:, :) = I_ij(:, :) increA_2(:, :) = I_ij(:, :) do n = 1, itr_tol if (n > 2) then NN = NN*(NN + 1.0d0) end if increA_1(:, :) = A_inv(:, :) increA_2(:, :) = matmul(increA_2, A) increA_3(:, :, :, :) = 0.0d0 do m = 1, n increA_1(:, :) = matmul(increA_1, A) increA_2(:, :) = matmul(increA_2, A_inv) do i = 1, size(A, 1) do j = 1, size(A, 1) do k = 1, size(A, 1) do l = 1, size(A, 1) increA_3(i, j, k, l) = increA_3(i, j, k, l) + increA_1(i, k)*increA_2(l, j) expA_A(i, j, k, l) = expA_A(i, j, k, l) + 1.0d0/NN*increA_3(i, j, k, l) end do end do end do end do end do do i = 1, size(A, 1) do j = 1, size(A, 1) do k = 1, size(A, 1) do l = 1, size(A, 1) increment = increment + 1.0d0/NN*increA_3(i, j, k, l) & *increA_3(i, j, k, l) & *increA_3(i, j, k, l) & *increA_3(i, j, k, l) end do end do end do end do if (increment <= TOL) then exit else if (n >= itr_tol) then stop "tensor exponential is not converged" end if cycle end if end do deallocate (increA_1, increA_2, increA_3, I_ij, A_inv) end subroutine tensor_expo_der !================================================================================== end module det