build_matrixLSF_quadIntp Subroutine

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

This routine builds up the matrix for least square fit used in quadratic interpolation. We extract momentum information completely on the view of the source coordinate system Set the right hand side of the equation system Solve the problem, where b = rhs, x = coefficients Ax = b overdetermined, solve the least Square fit problem (A^T)Ax = (A^T)b x = ((A^T)A)^-1(A^T)b Solve linear system of equation with inverted matrix. Size of matrix: (nCoeffs, QQ) matrix_LSF = ((A^T)A)^-1(A^T)

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_quadintp~~CallsGraph proc~build_matrixlsf_quadintp build_matrixLSF_quadIntp proc~alloc_matrix alloc_matrix proc~build_matrixlsf_quadintp->proc~alloc_matrix proc~invert_matrix invert_matrix proc~build_matrixlsf_quadintp->proc~invert_matrix proc~polyquadratic_1d polyQuadratic_1D proc~build_matrixlsf_quadintp->proc~polyquadratic_1d proc~polyquadratic_2d polyQuadratic_2D proc~build_matrixlsf_quadintp->proc~polyquadratic_2d proc~polyquadratic_3d polyQuadratic_3D proc~build_matrixlsf_quadintp->proc~polyquadratic_3d proc~tem_abort tem_abort proc~build_matrixlsf_quadintp->proc~tem_abort proc~tem_matrix_dump tem_matrix_dump proc~build_matrixlsf_quadintp->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_quadintp~~CalledByGraph proc~build_matrixlsf_quadintp build_matrixLSF_quadIntp proc~append_intpmatrixlsf append_intpMatrixLSF proc~append_intpmatrixlsf->proc~build_matrixlsf_quadintp interface~append~9 append interface~append~9->proc~append_intpmatrixlsf

Source Code

  subroutine build_matrixLSF_quadIntp(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 quadratic'
    ! 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,:) = polyQuadratic_1D( cxDirRK(:,iDir) )
      end do

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

    case (3)
      do iSrc = 1, nSources
        iDir = neighDir(iSrc)
        tmp_matrix%A(iSrc,:) = polyQuadratic_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))
    !write(*,*) 'Matrix_LSF '
    !do iDir = 1, nCoeffs
    !  write(*,*) 'iDir ', iDir, me%A(iDir, :)
    !end do

  end subroutine build_matrixLSF_quadIntp