scale_linalgebra.F90

Go to the documentation of this file.

!> Module common / Linear algebra
!!
!! @par Description
!!      A module to provide utilities for linear algebra
!!
!! @par Reference
!!
!! @author Yuta Kawai, Team SCALE
!!
#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
  !> Solve a linear equation Ax=b using a direct solver 
  interface linalgebra_solvelineq
    module procedure linalgebra_solvelineq_b1d
    module procedure linalgebra_solvelineq_b2d
  end interface
  public :: linalgebra_solvelineq
  public :: linalgebra_solvelineq_bndmat
  public :: linalgebra_solvelineq_gmres
 
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  
  !-----------------------------------------------------------------------------
 
  ! Private procedure
  !
  private :: my_dger
  private :: my_idamax
  private :: my_swap
 
  private :: precondstep_ilu0_constructmat
  private :: precondstep_ilu0_solve
  private :: precondstep_ptjacobi
  
  ! Private variable
  !
 
contains
 
!> Calculate a inversion of matrix A
!OCL SERIAL
  function linalgebra_inv(A) result(Ainv)
    implicit none
    real(rp), intent(in) :: a(:,:)         !< Matrix A
    real(rp) :: ainv(size(a,1),size(a,2))  !< Result of matrix inversion
 
    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
 
!> Solve a linear equation Ax=b using a direct solver
!OCL SERIAL
  subroutine linalgebra_solvelineq_b1d(A, b, x)
    implicit none
    real(RP), intent(in) :: A(:,:)        !< Coefficient matrix
    real(RP), intent(in) :: b(:)          !< Vector in the right-hand side
    real(RP), intent(out) :: x(size(b))   !< Solution vector of linear equation Ax=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
 
!> Solve a linear equation Ax=b using a direct solver where b is a matrix
!OCL SERIAL
  subroutine linalgebra_solvelineq_b2d(A, b, x)
    implicit none
    real(RP), intent(in) :: A(:,:)                   !< Coefficient matrix
    real(RP), intent(in) :: b(:,:)                   !< Matrix in the right-hand side
    real(RP), intent(out) :: x(size(b,1),size(b,2))  !< Matrix storing solution of linear equation Ax=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, 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
 
!> Perform LU factorization
!OCL SERIAL
  subroutine linalgebra_lu(A_lu, ipiv)
    implicit none
    real(rp), intent(inout) :: a_lu(:,:)       !> Matrix performed LU factorization. Note that original matrix is overwritten.
    integer, intent(out) :: ipiv(size(a_lu,1)) !> Array storing the pivoting information
 
    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
 
 
!> Calculate the solution of linear equations with a band matrix 
!!
!! In this module, we treat a linear equations system, A * X = B
!! where A is a band matrix of order N with KL subdiagonals and KU superdiagonals, 
!! and X and B are matrices whose size is N * NRHS.  
!! 
!! @param A Matrix A in band storage same as in dgbsv of LAPACK
!! @param b Right hand side matrix B with N-by-RHS 
!! @param ipiv Pivot indices that define the permutation matrix
!! @param KL Number of subdiagonals within the band of A
!! @param KU Number of superdiagonals within the band of A
!! @param NRHS Number of right hand sides
!! @param use_lapack Flag whether LAPACK is used
!OCL SERIAL
  subroutine linalgebra_solvelineq_bndmat(A, b, ipiv, &
    N, KL, KU, NRHS, use_lapack )
    implicit none
    integer, intent(in) :: n
    integer, intent(in) :: kl
    integer, intent(in) :: ku
    integer, intent(in) :: nrhs
    real(rp), intent(inout) :: a(2*kl+ku+1,n)
    real(rp), intent(inout) :: b(n,nrhs)
    integer, intent(out) :: ipiv(n)
    logical, intent(in), optional :: use_lapack
 
    integer :: i, j, l, lm, p
    integer :: km
    integer :: ju, jp
    integer :: kv
    real(rp) :: tmp
    integer :: lda, ldb
    logical :: use_lapack_
 
    integer :: info
    !--------------------------------------------------------------------------- 
 
    if (present(use_lapack)) then
      use_lapack_ = use_lapack
    else
      use_lapack_ = .false.
    end if
 
    if ( use_lapack_ ) then
      call dgbsv( n, kl, ku, nrhs, a, 2*kl+ku+1, ipiv, b, n, info)
      return
    end if
 
    !--
    kv = kl + ku
    lda = size(a,1)
    ldb = size(b,1)
 
    do j=ku+2, min(kv, n)
    do i=kv-j+2, kl
      a(i,j) = 0.0_rp
    end do
    end do
    
    ju = 1
    do j=1, n
      if ( j+kv <= n ) then
        do i=1, kl
          a(i,j+kv) = 0.0_rp
        end do
      end if
 
      km = min(kl, n-j)
      call my_idamax( km+1, a(kv+1,j), &
        jp )
      ipiv(j) = jp + j - 1
 
      if ( a(kv+jp,j) /= 0.0_rp ) then
        ju = max( ju, min(j+ku+jp-1,n) )
 
        if ( jp /= 1 ) then
          call my_swap( ju-j+1, a(kv+jp,j), lda-1, &
                                a(kv +1,j), lda-1  )
        end if
 
        if ( km > 0 ) then
          tmp = 1.0_rp / a(kv+1,j)
          do l=kv+2, kv+1+km
            a(l,j) = tmp * a(l,j)
          end do
          call my_dger( km, ju-j, &
            a(kv+2,j  ),          &
            a(kv  ,j+1), lda-1,   &
            a(kv+1,j+1), lda-1    )
        end if
      end if
    end do
    
    !--
    ! Solve L*X = B, overwriting B with X
 
    if ( kl > 0 ) then
      do j=1, n-1
        lm = min( kl, n-j )
        l = ipiv(j)
        if ( l /= j ) call my_swap( nrhs, b(l,1), ldb, b(j,1), ldb )
        call my_dger( lm, nrhs, &
          a(kv+2,j),     &
          b(j  ,1), ldb, &
          b(j+1,1), ldb  )
      end do
    end if
 
    !--
    ! Solve U*X = B, overwriting B with X
    do p=1, nrhs
      do j=n, 1, -1
        l = kv + 1 - j
        b(j,p) = b(j,p) / a(kv+1,j)
        tmp = b(j,p)
        do i=j-1, max(1,j-kv),-1
          b(i,p) = b(i,p) - tmp * a(l+i,j)
        end do
      end do 
    end do
 
    return
  end subroutine linalgebra_solvelineq_bndmat
 
  !> Solve a linear equation Ax=b using the Generalized Minimal Residual solver (GMRES)
  !!
  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   !< Object managing sparse matrix A
    real(rp), intent(in) :: b(:)       !< Right hand side vector
    real(rp), intent(inout) :: x(:)    !< Solution vector
    integer, intent(in) :: m           !<
    integer, intent(in) ::restart_num  !< Number of times by which restarting is performed if convergence condition is not satisfied
    real(rp), intent(in) :: conv_crit  !< Threshold for finishing the iteration solver
 
    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
 
!- private ----------------------------------------------------------------
 
!OCL SERIAL
  subroutine my_idamax( N, DX, ind )
    implicit none
    integer, intent(in) :: n
    real(rp), intent(in) :: dx(*)
    integer, intent(out) :: ind
 
    integer :: i
    real(rp) :: dmax
    !----------------------------
    dmax = abs(dx(1))
    ind = 1
    do i=2, n
      if ( abs(dx(i)) > dmax ) then
        ind = i
        dmax = abs(dx(i))
      end if
    end do
    return
  end subroutine my_idamax
 
!> Calculate A_ij = A_ij - x_i y_j 
!OCL SERIAL
  subroutine my_dger( m, n, x, y, incy, A, LDA)
    implicit none
    integer, intent(in) :: m
    integer, intent(in) :: n
    integer, intent(in) :: lda    
    real(rp), intent(inout) :: a(lda,*)
    real(rp), intent(in) :: x(*)
    real(rp), intent(in) :: y(*)
    integer, intent(in) :: incy
 
    integer :: i, j
    integer :: jy
    real(rp) :: tmp
    !---------------------------------------------------
    jy = 1
    do j=1, n
      tmp = y(jy)
      if ( tmp /= 0.0_rp ) then
        do i=1, m
          a(i,j) = a(i,j) - tmp * x(i)
        end do
      end if
      jy = jy + incy
    end do
    return
  end subroutine my_dger
 
!> Interchange two vectors
!!
!! Note that incx, incy should be > 0
!OCL SERIAL
  subroutine my_swap( n, dx, incx, dy, incy )
    implicit none
    integer, intent(in) :: n
    real(rp), intent(inout) :: dx(*)
    integer, intent(in) :: incx
    real(rp), intent(inout) :: dy(*)
    integer, intent(in) :: incy
 
    integer :: i, ix, iy
    integer :: m, mp1
    real(rp) :: dtemp
    !---------------------------------------------------
    if ( incx==1 .and. incy==1 ) then
      m = mod(n,3)
      if (m /= 0) then
        do i=1, m
          dtemp = dx(i)
          dx(i) = dy(i)
          dy(i) = dtemp
        end do
        mp1 = m + 1
        do i=mp1, n, 3
          dtemp = dx(i)
          dx(i) = dy(i)
          dy(i) = dtemp
          dtemp = dx(i+1)
          dx(i+1) = dy(i+1)
          dy(i+1) = dtemp
          dtemp = dx(i+2)
          dx(i+2) = dy(i+2)
          dy(i+2) = dtemp
        end do
      end if
    else
      ix = 1; iy = 1
      do i=1, n
        dtemp = dx(ix)
        dx(ix) = dy(iy)
        dy(iy) = dtemp
        ix = ix + incx
        iy = iy + incy
      end do
    end if
    return
  end subroutine my_swap
 
!OCL SERIAL
  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
 
!OCL SERIAL
  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
 
!OCL SERIAL
  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