build_matrixLSF_linearIntp Subroutine

private subroutine build_matrixLSF_linearIntp(me, QQ, nDims, nSources, cxDirRK, neighDir, nCoeffs, success)

This routine builds up the matrix for least square fit used in linear interpolation.

Arguments

Type IntentOptional Attributes Name
type(tem_matrix_type), intent(out) :: me

Matrix to fill
integer, intent(in) :: QQ

Number of stencil directions
integer, intent(in) :: nDims

Number of dimensions
integer, intent(in) :: nSources

Number of sources from coarser found
real(kind=rk), intent(in) :: cxDirRK(3,QQ)

Stencil directions
integer, intent(in) :: neighDir(nSources)

direction in which sources are found
integer, intent(in) :: nCoeffs

nUnknown coeffs
logical, intent(out) :: success

success if false if matrix is singular reduce interpolation order

Calls

proc~~build_matrixlsf_linearintp~~CallsGraph proc~build_matrixlsf_linearintp build_matrixLSF_linearIntp proc~alloc_matrix alloc_matrix proc~build_matrixlsf_linearintp->proc~alloc_matrix proc~invert_matrix invert_matrix proc~build_matrixlsf_linearintp->proc~invert_matrix proc~polylinear_1d polyLinear_1D proc~build_matrixlsf_linearintp->proc~polylinear_1d proc~polylinear_2d polyLinear_2D proc~build_matrixlsf_linearintp->proc~polylinear_2d proc~polylinear_3d polyLinear_3D proc~build_matrixlsf_linearintp->proc~polylinear_3d proc~tem_abort tem_abort proc~build_matrixlsf_linearintp->proc~tem_abort proc~tem_matrix_dump tem_matrix_dump proc~build_matrixlsf_linearintp->proc~tem_matrix_dump proc~alloc_matrix->proc~tem_abort proc~invert_matrix->proc~tem_abort dgetrf dgetrf proc~invert_matrix->dgetrf dgetri dgetri proc~invert_matrix->dgetri mpi_abort mpi_abort proc~tem_abort->mpi_abort

Called by

proc~~build_matrixlsf_linearintp~~CalledByGraph proc~build_matrixlsf_linearintp build_matrixLSF_linearIntp proc~append_intpmatrixlsf append_intpMatrixLSF proc~append_intpmatrixlsf->proc~build_matrixlsf_linearintp interface~append~9 append interface~append~9->proc~append_intpmatrixlsf

Source Code

  subroutine build_matrixLSF_linearIntp(me, QQ, nDims, nSources, cxDirRK, &
    &                                 neighDir, nCoeffs, success)
    ! --------------------------------------------------------------------------
    !> Matrix to fill
    type(tem_matrix_type), intent(out) :: me
    !> Number of stencil directions
    integer, intent(in) :: QQ
    !> Number of dimensions
    integer, intent(in) :: nDims
    !> Number of sources from coarser found
    integer, intent(in) :: nSources
    !> Stencil directions
    real(kind=rk), intent(in) :: cxDirRK(3,QQ)
    !> direction in which sources are found
    integer, intent(in) :: neighDir(nSources)
    !> nUnknown coeffs
    integer, intent(in) :: nCoeffs
    !> success if false if matrix is singular reduce interpolation order
    logical, intent(out) :: success
    ! --------------------------------------------------------------------------
    integer :: iDir, iSrc
    !> Each row represents a polynomial evaluated at coord of elements in
    ! stencil directions
    type(tem_matrix_type) :: tmp_matrix
    real(kind=rk) :: inv_AtA(nCoeffs,nCoeffs), AtA(nCoeffs,nCoeffs)
    integer :: errCode
    ! --------------------------------------------------------------------------
write(dbgUnit(4),"(A)") 'Inside build least square fit matrix for linear'
    ! me%A = ((A^T)A)^-1*(A^T)
    ! inv_AtA = ((A^T)A)^-1
    call alloc_matrix(tmp_matrix, nSources, nCoeffs)
    select case(nDims)
    case (1)
      do iSrc = 1, nSources
        iDir = neighDir(iSrc)
        tmp_matrix%A(iSrc,:) = polyLinear_1D( cxDirRK(:,iDir) )
      end do

    case (2)
      do iSrc = 1, nSources
        iDir = neighDir(iSrc)
        tmp_matrix%A(iSrc,:) = polyLinear_2D( cxDirRK(:,iDir) )
      end do

    case (3)
      do iSrc = 1, nSources
        iDir = neighDir(iSrc)
        tmp_matrix%A(iSrc,:) = polyLinear_3D( cxDirRK(:,iDir) )
      end do

    case default
      write(logUnit(1),*) 'Unknown nDims for quadratic interpolation'
      call tem_abort()
    end select

write(dbgUnit(4),"(A)") ' tmp matrix:'
call tem_matrix_dump(tmp_matrix, dbgUnit(4))
flush(dbgUnit(4))

    AtA =  matmul( transpose(tmp_matrix%A), tmp_matrix%A)
    inv_AtA = invert_matrix(AtA, errCode)
    if (errCode == 0) then
      ! matrix_LSF size is transpose of tmp_matrix
      call alloc_matrix(me, nCoeffs, nSources)
      me%A = matmul( inv_AtA, transpose(tmp_matrix%A) )
      success = .true.
    else
      ! singular matrix, reduce interpolation order
      call alloc_matrix(me, 1, 1)
      success = .false.
      return
    end if

write(dbgUnit(4),"(A)") ' Matrix LSF:'
call tem_matrix_dump(me, dbgUnit(4))
flush(dbgUnit(4))

  end subroutine build_matrixLSF_linearIntp