scale_linalgebra.F90

Go to the documentation of this file.

 
#include "scaleFElib.h"
module scale_linalgebra
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_sparsemat, only: sparsemat
 
  !-----------------------------------------------------------------------------
  implicit none
  private
 
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !  
 
  public :: linalgebra_lu
  public :: linalgebra_inv
  interface linalgebra_solvelineq
    module procedure linalgebra_solvelineq_b1d
    module procedure linalgebra_solvelineq_b2d
  end interface
  public :: linalgebra_solvelineq
  public :: linalgebra_solvelineq_gmres
 
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  
  !-----------------------------------------------------------------------------
 
  ! Private procedure
  !
 
  ! Private variable
  !
 
contains
 
!OCL SERIAL
  function linalgebra_inv(A) result(Ainv)
 
    real(rp), intent(in) :: a(:,:)
    real(rp) :: ainv(size(a,1),size(a,2))
 
    real(rp) :: work(size(a,1))
    integer :: ipiv(size(a,1))
    integer :: n, info
 
    !--------------------------------------------------------------------------- 
 
    ainv(:,:) = a
    n = size(a,1)
 
    call dgetrf(n, n, ainv, n, ipiv, info)
    if (info /=0 ) then
      log_error("linalgebra_inv",*) "Matrix is singular"
      call prc_abort
    end if
 
    call dgetri(n, ainv, n, ipiv, work, n, info)
    if (info /=0 ) then
      log_error("linalgebra_inv",*) "Matrix inversion is failed"
      call prc_abort
    end if
 
    return
  end function linalgebra_inv
 
!OCL SERIAL
  subroutine linalgebra_solvelineq_b1d(A, b, x)
 
    real(RP), intent(in) :: A(:,:)
    real(RP), intent(in) :: b(:)
    real(RP), intent(out) :: x(size(b))
 
    real(RP) :: A_lu(size(A,1),size(A,2))
    integer :: ipiv(size(A,1))
    integer :: n, info
 
    !--------------------------------------------------------------------------- 
 
    a_lu = a
    n = size(a,1)
 
    call dgetrf(n, n, a_lu, n, ipiv, info)
    if (info /=0 ) then
      log_error("linalgebra_SolveLinEq",*) "Matrix is singular"
      call prc_abort
    end if
 
    x(:) = b
    call dgetrs('N', n, 1, a_lu, n, ipiv, x, n, info)
    if (info /=0 ) then
      log_error("linalgebra_SolveLinEq",*)  "Matrix inversion is failed"
      call prc_abort
    end if
 
    return
  end subroutine linalgebra_solvelineq_b1d
 
!OCL SERIAL
  subroutine linalgebra_solvelineq_b2d(A, b, x)
 
    real(RP), intent(in) :: A(:,:)
    real(RP), intent(in) :: b(:,:)
    real(RP), intent(out) :: x(size(b,1),size(b,2))
 
    real(RP) :: A_lu(size(A,1),size(A,2))
    integer :: ipiv(size(A,1))
    integer :: n, info
 
    !--------------------------------------------------------------------------- 
 
    a_lu = a
    n = size(a,1)
 
    call dgetrf(n, n, a_lu, n, ipiv, info)
    if (info /=0 ) then
      log_error("linalgebra_SolveLinEq",*) "Matrix is singular"
      call prc_abort
    end if
 
    x(:,:) = b
    call dgetrs('N', n, size(b,2), a_lu, n, ipiv, x, n, info)
    if (info /=0 ) then
      log_error("linalgebra_SolveLinEq",*)  "Matrix inversion is failed"
      call prc_abort
    end if
 
    return
  end subroutine linalgebra_solvelineq_b2d
 
!OCL SERIAL
  subroutine linalgebra_lu(A_lu, ipiv)
 
    real(rp), intent(inout) :: a_lu(:,:)
    integer, intent(out) :: ipiv(size(a_lu,1))
 
    integer :: n
    integer :: info
 
    !--------------------------------------------------------------------------- 
    n = size(a_lu,1)
 
    call dgetrf(n, n, a_lu, n, ipiv, info)
    if (info /=0 ) then
      log_error("linalgebra_LU",*)  "Matrix is singular"
      call prc_abort
    end if
 
    return
  end subroutine linalgebra_lu
 
  subroutine linalgebra_solvelineq_gmres(A, b, x, m, restart_num, CONV_CRIT)
 
    use scale_sparsemat, only: sparsemat_matmul
    implicit none
    
    type(sparsemat), intent(in) :: a
    real(rp), intent(in) :: b(:)
    real(rp), intent(inout) :: x(:)
    integer, intent(in) :: m
    integer, intent(in) ::restart_num
    real(rp), intent(in) :: conv_crit
 
    real(rp) :: x0(size(x))
    real(rp) :: w(size(x))
    real(rp) :: v(size(x),m+1)
    real(rp) :: z(size(x),m+1)
    real(rp) :: av(size(x))
    real(rp) :: g(size(x)+1)
    real(rp) :: r0_l2, r0_l2_new
    real(rp) :: h(m+1,m)
    real(rp) :: r(m,m)
    real(rp) :: c(m), s(m)
    real(rp) :: y(m)
    real(rp) :: tmp, tmp1, tmp2
    integer :: i, j
    integer :: restart_i
 
    type(sparsemat) :: m_pc
 
    !--------------------------------------------------------------------------- 
 
    call precondstep_ilu0_constructmat(a, m_pc)
    x0(:) = x
 
    do restart_i=1, restart_num
 
      if (restart_i == 1) then
       call sparsemat_matmul(a, x, av)
       v(:,1) = b - av
       r0_l2 = sqrt(sum(v(:,1)**2))
      end if
      v(:,1) = v(:,1)/r0_l2
      g(1) = r0_l2
 
      do j=1, m
 
        call precondstep_ilu0_solve(m_pc, v(:,j), z(:,j))
 
        call sparsemat_matmul(a, z(:,j), w)
 
        do i=1, j
          h(i,j) = sum(w*v(:,i))
          w(:) = w - h(i,j)*v(:,i)
        end do
        h(j+1,j) = sqrt(sum(w**2))
        v(:,j+1) = w(:)/h(j+1,j)
 
        r(1,j) = h(1,j)
        do i=1, j-1
          tmp1 =  c(i)*r(i,j) + s(i)*h(i+1,j)
          tmp2 = -s(i)*r(i,j) + c(i)*h(i+1,j)
          r(i,j) = tmp1
          r(i+1,j) = tmp2
        end do
 
        tmp = sqrt(r(j,j)**2 + h(j+1,j)**2)
        c(j) = r(j,j)/tmp
        s(j) = h(j+1,j)/tmp
 
        g(j+1) = -s(j)*g(j)
        g(j) = c(j)*g(j)
        r(j,j) = c(j)*r(j,j) + s(j)*h(j+1,j)
      end do
 
      y(m) = g(m)/r(m,m)
      do i=m-1,1,-1
        y(i) = (g(i) - sum(r(i,i+1:m)*y(i+1:m)))/r(i,i)
      end do
 
      x(:) = x + matmul(z(:,1:m),y)
 
      !write(*,*) "x:", x
      !write(*,*) "x0:", x0
 
      call sparsemat_matmul(a, x, av)
      v(:,1) = b - av
      r0_l2_new = sqrt(sum(v(:,1)**2))
 
      !write(*,*) "r:", r0_l2, "->", r0_l2_new
      if (r0_l2_new < conv_crit) exit
      r0_l2 = r0_l2_new
    end do
 
    return
  end subroutine linalgebra_solvelineq_gmres
 
  subroutine precondstep_ptjacobi(A, b, x)
    implicit none
    type(sparsemat), intent(in) :: a
    real(rp), intent(in) :: b(:)
    real(rp), intent(out) :: x(size(b))
 
    integer :: n
    !--------------------------------------------------------------------------- 
 
    do n=1, size(x)
      x(n) = b(n) / a%GetVal(n, n)
    end do
 
    return
  end subroutine precondstep_ptjacobi
 
  subroutine precondstep_ilu0_constructmat(A, M)
    use scale_const, only: &
      eps => const_eps
    implicit none
 
    type(sparsemat), intent(in) :: a
    type(sparsemat), intent(inout) :: m
 
    integer :: i, j, k, n
    real(rp) :: mij, m_ik, m_kk
    !--------------------------------------------------------------------------- 
 
    n = a%rowPtrSize-1
 
    m = a
    do i=2, n
      do k=1, i-1
        m_ik = m%GetVal(i,k)
        m_kk = m%GetVal(k,k)
        if ( abs(m_ik) < eps .and. abs(m_kk) < eps ) then
          m_ik = m_ik / m_kk
          call m%ReplaceVal( i, k,  m_ik )
 
          do j=k+1, n
            mij = m%GetVal(i,j)
            if ( abs(mij) < eps ) then
              call m%ReplaceVal( i, j, mij - m%GetVal(i,k) * m%GetVal(k,j) )
            end if
          end do
        end if
      end do
    end do
 
    return
  end subroutine precondstep_ilu0_constructmat
 
  subroutine precondstep_ilu0_solve(M, b, x)
    use scale_sparsemat, only: &
      sparsemat_storage_typeid_csr
    implicit none
 
    type(sparsemat), intent(in) :: m
    real(rp), intent(in) :: b(:)
    real(rp), intent(out) :: x(size(b))
 
    integer :: i
    integer :: j
    integer :: n
 
    integer :: j1, j2
    !--------------------------------------------------------------------------- 
 
    if ( m%GetStorageFormatId() /= sparsemat_storage_typeid_csr ) then
      log_error("linalgebra_PreCondStep_ILU0_solve",*)  "The strorge type of specified sparse matrix is not supported. Check!"
      call prc_abort
    end if
 
    n = size(x)
 
    x(1) = b(1)
    do i=2, n
      j1 = m%rowPtr(i)
      j2 = m%rowPtr(i+1)-1
      x(i) = b(i)
      do j=j1, j2
        if (m%colIdx(j) <= i-1) &
          x(i) = x(i) - m%val(j)*x(m%colIdx(j))
      end do
    end do
 
    x(n) = x(n) / m%GetVal(n,n)
    do i=n-1, 1, -1
      j1 = m%rowPtr(i)
      j2 = m%rowPtr(i+1)-1
      do j=j1, j2
        if (m%colIdx(j) >= i+1) &
          x(i) = x(i) - m%val(j)*x(m%colIdx(j))
      end do
      x(i) = x(i)/m%GetVal(i,i)
    end do
 
    return
  end subroutine precondstep_ilu0_solve
 
end module scale_linalgebra