solve_small_linear_system Subroutine

public subroutine solve_small_linear_system(n, a_in, b_in, x, ok)

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n
real(kind=dp), intent(in) :: a_in(n,n)
real(kind=dp), intent(in) :: b_in(n)
real(kind=dp), intent(out) :: x(n)
logical, intent(out) :: ok

Called by

proc~~solve_small_linear_system~~CalledByGraph proc~solve_small_linear_system solve_small_linear_system proc~try_newton_solve try_newton_solve proc~try_newton_solve->proc~solve_small_linear_system proc~solve_nonlinear_system solve_nonlinear_system proc~solve_nonlinear_system->proc~try_newton_solve proc~try_solve_zhao_branch_a try_solve_zhao_branch_a proc~try_solve_zhao_branch_a->proc~solve_nonlinear_system proc~try_solve_zhao_branch_b try_solve_zhao_branch_b proc~try_solve_zhao_branch_b->proc~solve_nonlinear_system proc~try_solve_zhao_branch_c try_solve_zhao_branch_c proc~try_solve_zhao_branch_c->proc~solve_nonlinear_system proc~solve_zhao_branch_a solve_zhao_branch_a proc~solve_zhao_branch_a->proc~try_solve_zhao_branch_a proc~solve_zhao_branch_b solve_zhao_branch_b proc~solve_zhao_branch_b->proc~try_solve_zhao_branch_b proc~solve_zhao_branch_c solve_zhao_branch_c proc~solve_zhao_branch_c->proc~try_solve_zhao_branch_c proc~solve_zhao_unknowns solve_zhao_unknowns proc~solve_zhao_unknowns->proc~try_solve_zhao_branch_a proc~solve_zhao_unknowns->proc~try_solve_zhao_branch_b proc~solve_zhao_unknowns->proc~try_solve_zhao_branch_c proc~solve_zhao_unknowns->proc~solve_zhao_branch_a proc~solve_zhao_unknowns->proc~solve_zhao_branch_b proc~solve_zhao_unknowns->proc~solve_zhao_branch_c proc~resolve_sheath_injection_context resolve_sheath_injection_context proc~resolve_sheath_injection_context->proc~solve_zhao_unknowns proc~init_particle_batch_from_config init_particle_batch_from_config proc~init_particle_batch_from_config->proc~resolve_sheath_injection_context

Source Code

  subroutine solve_small_linear_system(n, a_in, b_in, x, ok)
    integer, intent(in) :: n
    real(dp), intent(in) :: a_in(n, n), b_in(n)
    real(dp), intent(out) :: x(n)
    logical, intent(out) :: ok

    integer :: i, j, k, pivot_row
    real(dp) :: a(n, n), b(n), factor, pivot_abs, tmp_row(n), tmp_val

    a = a_in
    b = b_in
    ok = .true.

    do k = 1, n
      pivot_row = k
      pivot_abs = abs(a(k, k))
      do i = k + 1, n
        if (abs(a(i, k)) > pivot_abs) then
          pivot_abs = abs(a(i, k))
          pivot_row = i
        end if
      end do
      if (pivot_abs <= 1.0d-18) then
        ok = .false.
        x = 0.0d0
        return
      end if
      if (pivot_row /= k) then
        tmp_row = a(k, :)
        a(k, :) = a(pivot_row, :)
        a(pivot_row, :) = tmp_row
        tmp_val = b(k)
        b(k) = b(pivot_row)
        b(pivot_row) = tmp_val
      end if
      do i = k + 1, n
        factor = a(i, k)/a(k, k)
        a(i, k:n) = a(i, k:n) - factor*a(k, k:n)
        b(i) = b(i) - factor*b(k)
      end do
    end do

    x = 0.0d0
    do i = n, 1, -1
      x(i) = b(i)
      do j = i + 1, n
        x(i) = x(i) - a(i, j)*x(j)
      end do
      x(i) = x(i)/a(i, i)
    end do
  end subroutine solve_small_linear_system