InvertMatrix Subroutine

  1. SELF_ColumnOperators_3D_t.f90
  2. SELF_ColumnOperators_3D_t
  3. InvertMatrix

public subroutine InvertMatrix(a, aInv, n)

Inverts the n x n matrix a by Gauss-Jordan elimination with partial pivoting, in double precision. Used once in Init for the (N+1) x (N+1) sweep matrices, which are small and well enough conditioned for a direct solve; this follows the routine of the same name in SELF_RefinementIndicator_3D_t.

Arguments

TypeIntentOptionalAttributesName
real(kind=real64), intent(in) :: a(1:n,1:n)
real(kind=real64), intent(out) :: aInv(1:n,1:n)
integer, intent(in) :: n

Called by

proc~~invertmatrix~3~~CalledByGraph proc~invertmatrix~3 InvertMatrix proc~init_columnoperators3d_t Init_ColumnOperators3D_t proc~init_columnoperators3d_t->proc~invertmatrix~3 proc~init_columnoperators3d Init_ColumnOperators3D proc~init_columnoperators3d->proc~init_columnoperators3d_t

Contents

Source Code

InvertMatrix

Source Code

  subroutine InvertMatrix(a,aInv,n)
    !! Inverts the n x n matrix a by Gauss-Jordan elimination with partial pivoting, in
    !! double precision. Used once in Init for the (N+1) x (N+1) sweep matrices, which
    !! are small and well enough conditioned for a direct solve; this follows the
    !! routine of the same name in SELF_RefinementIndicator_3D_t.
    implicit none
    integer,intent(in) :: n
    real(real64),intent(in) :: a(1:n,1:n)
    real(real64),intent(out) :: aInv(1:n,1:n)
    ! Local
    real(real64) :: b(1:n,1:n)
    integer :: i,j,k,piv
    real(real64) :: pmax,factor,tmp

    b = a
    aInv = 0.0_real64
    do i = 1,n
      aInv(i,i) = 1.0_real64
    enddo

    do k = 1,n
      ! Partial pivot: the largest entry of column k at or below the diagonal
      piv = k
      pmax = abs(b(k,k))
      do i = k+1,n
        if(abs(b(i,k)) > pmax) then
          pmax = abs(b(i,k))
          piv = i
        endif
      enddo
      if(pmax <= 0.0_real64) then
        error stop "ColumnOperators3D : singular sweep matrix"
      endif
      if(piv /= k) then
        do j = 1,n
          tmp = b(k,j); b(k,j) = b(piv,j); b(piv,j) = tmp
          tmp = aInv(k,j); aInv(k,j) = aInv(piv,j); aInv(piv,j) = tmp
        enddo
      endif

      ! Normalize the pivot row
      factor = b(k,k)
      do j = 1,n
        b(k,j) = b(k,j)/factor
        aInv(k,j) = aInv(k,j)/factor
      enddo

      ! Eliminate column k from every other row
      do i = 1,n
        if(i /= k) then
          factor = b(i,k)
          do j = 1,n
            b(i,j) = b(i,j)-factor*b(k,j)
            aInv(i,j) = aInv(i,j)-factor*aInv(k,j)
          enddo
        endif
      enddo
    enddo

  endsubroutine InvertMatrix