pqm_solve_diag_dominant Subroutine

private pure subroutine pqm_solve_diag_dominant(n, al, ac, au, r, x)

Diagonally-dominant tridiagonal solve; central diagonal supplied as the OFFSET ac from al + au (full pivot = ac + al + au). Never divides by zero for positive-definite ac, al, au (White & Adcroft 2008).

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n

Number of unknowns (= number of edges = nz+1)
real(kind=wp), intent(in) :: al(n)

Lower diagonal (al(1) unused)
real(kind=wp), intent(in) :: ac(n)

Central-diagonal OFFSET from al+au (full diagonal = ac+al+au)
real(kind=wp), intent(in) :: au(n)

Upper diagonal (au(n) unused)
real(kind=wp), intent(in) :: r(n)

Right-hand side
real(kind=wp), intent(out) :: x(n)

Solution vector

Called by

proc~~pqm_solve_diag_dominant~~CalledByGraph proc~pqm_solve_diag_dominant pqm_solve_diag_dominant proc~remap_column_pqm remap_column_pqm proc~remap_column_pqm->proc~pqm_solve_diag_dominant proc~remap_column remap_column proc~remap_column->proc~remap_column_pqm proc~ocean_remap_tracer_column ocean_remap_tracer_column proc~ocean_remap_tracer_column->proc~remap_column proc~ocean_remap_tracer_field ocean_remap_tracer_field proc~ocean_remap_tracer_field->proc~remap_column proc~remap_layer_to_density_impl remap_layer_to_density_impl proc~remap_layer_to_density_impl->proc~remap_column proc~remap_layer_to_vcoord_impl remap_layer_to_vcoord_impl proc~remap_layer_to_vcoord_impl->proc~remap_column proc~remap_tracer_grounded remap_tracer_grounded proc~remap_tracer_grounded->proc~remap_column proc~remap_x_face_grounded remap_x_face_grounded proc~remap_x_face_grounded->proc~remap_column proc~remap_x_face_velocity remap_x_face_velocity proc~remap_x_face_velocity->proc~remap_column proc~remap_y_face_grounded remap_y_face_grounded proc~remap_y_face_grounded->proc~remap_column proc~remap_y_face_velocity remap_y_face_velocity proc~remap_y_face_velocity->proc~remap_column proc~ocean_apply_ale_remap_centres ocean_apply_ale_remap_centres proc~ocean_apply_ale_remap_centres->proc~ocean_remap_tracer_field proc~ocean_apply_ale_remap_faces ocean_apply_ale_remap_faces proc~ocean_apply_ale_remap_faces->proc~remap_x_face_velocity proc~ocean_apply_ale_remap_faces->proc~remap_y_face_velocity proc~ocean_apply_ale_remap_step ocean_apply_ale_remap_step proc~ocean_apply_ale_remap_step->proc~ocean_remap_tracer_field proc~ocean_apply_ale_remap_step->proc~remap_x_face_velocity proc~ocean_apply_ale_remap_step->proc~remap_y_face_velocity proc~ocean_apply_conservative_min_thickness ocean_apply_conservative_min_thickness proc~ocean_apply_conservative_min_thickness->proc~remap_tracer_grounded proc~ocean_apply_conservative_min_thickness->proc~remap_x_face_grounded proc~ocean_apply_conservative_min_thickness->proc~remap_y_face_grounded proc~remap_layer_to_density remap_layer_to_density proc~remap_layer_to_density->proc~remap_layer_to_density_impl proc~remap_layer_to_sigma remap_layer_to_sigma proc~remap_layer_to_sigma->proc~remap_layer_to_vcoord_impl proc~remap_layer_to_z remap_layer_to_z proc~remap_layer_to_z->proc~remap_layer_to_vcoord_impl proc~remap_layer_to_zstar remap_layer_to_zstar proc~remap_layer_to_zstar->proc~remap_layer_to_vcoord_impl proc~ocean_dyn_step_split ocean_dyn_step_split proc~ocean_dyn_step_split->proc~ocean_apply_ale_remap_step proc~run_continuity_chain run_continuity_chain proc~run_continuity_chain->proc~ocean_apply_conservative_min_thickness

Variables

Type Visibility Attributes Name Initial
real(kind=wp), private :: c1(NZ_STACK_MAX+1)
real(kind=wp), private :: d1
real(kind=wp), private :: denom_t1
real(kind=wp), private :: i_pivot
integer, private :: k

Source Code

   pure subroutine pqm_solve_diag_dominant(n, al, ac, au, r, x)
      !$acc routine seq
      !! Diagonally-dominant tridiagonal solve; central diagonal supplied as the
      !! OFFSET `ac` from `al + au` (full pivot = ac + al + au). Never divides by
      !! zero for positive-definite ac, al, au (White & Adcroft 2008).
      integer, intent(in) :: n
         !! Number of unknowns (= number of edges = nz+1)
      real(wp), intent(in) :: al(n)
         !! Lower diagonal (al(1) unused)
      real(wp), intent(in) :: ac(n)
         !! Central-diagonal OFFSET from al+au (full diagonal = ac+al+au)
      real(wp), intent(in) :: au(n)
         !! Upper diagonal (au(n) unused)
      real(wp), intent(in) :: r(n)
         !! Right-hand side
      real(wp), intent(out) :: x(n)
         !! Solution vector

      real(wp) :: c1(NZ_STACK_MAX + 1)
      real(wp) :: d1, i_pivot, denom_t1
      integer :: k

      i_pivot = 1.0_wp/(ac(1) + au(1))
      d1 = ac(1)*i_pivot
      c1(1) = au(1)*i_pivot
      x(1) = r(1)*i_pivot
      do k = 2, n - 1
         denom_t1 = ac(k) + d1*al(k)
         i_pivot = 1.0_wp/(denom_t1 + au(k))
         d1 = denom_t1*i_pivot
         c1(k) = au(k)*i_pivot
         x(k) = (r(k) - al(k)*x(k - 1))*i_pivot
      end do
      i_pivot = 1.0_wp/(ac(n) + d1*al(n))
      x(n) = (r(n) - al(n)*x(n - 1))*i_pivot
      do k = n - 1, 1, -1
         x(k) = x(k) - c1(k)*x(k + 1)
      end do
   end subroutine pqm_solve_diag_dominant