scale_linalgebra Module Reference

Module common / Linear algebra. More...

Data Types
interface	linalgebra_solvelineq
	Solve a linear equation Ax=b using a direct solver. More...

Functions/Subroutines
real(rp) function, dimension(size(a, 1), size(a, 2)), public	linalgebra_inv (a)
	Calculate a inversion of matrix A.
subroutine, public	linalgebra_lu (a_lu, ipiv)
	Perform LU factorization.
subroutine, public	linalgebra_solvelineq_bndmat (a, b, ipiv, n, kl, ku, nrhs, use_lapack)
	Calculate the solution of linear equations with a band matrix.
subroutine, public	linalgebra_solvelineq_gmres (a, b, x, m, restart_num, conv_crit)
	Solve a linear equation Ax=b using the Generalized Minimal Residual solver (GMRES)

Detailed Description

Module common / Linear algebra.

Description

A module to provide utilities for linear algebra

Reference

Author

Yuta Kawai, Team SCALE

Function/Subroutine Documentation

linalgebra_inv()

real(rp) function, dimension(size(a,1),size(a,2)), public scale_linalgebra::linalgebra_inv

(

real(rp), dimension(:,:), intent(in)

a

)

Calculate a inversion of matrix A.

Parameters

[in]

a

Matrix A

Returns

Result of matrix inversion

Definition at line 76 of file scale_linalgebra.F90.

    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

Referenced by scale_element_base::elementbase_construct_massmat(), scale_element_operation_tensorprod3d_interp::elementoptrtensorprod3dinterp_init(), get_pmatd_lu(), scale_element_hexahedral::hexhedralelement_init(), scale_element_line::lineelement_init(), and scale_element_quadrilateral::quadrilateralelement_init().

linalgebra_lu()

subroutine, public scale_linalgebra::linalgebra_lu	(	real(rp), dimension(:,:), intent(inout)	a_lu,
		integer, dimension(size(a_lu,1)), intent(out)	ipiv )

Perform LU factorization.

Definition at line 173 of file scale_linalgebra.F90.

    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

Referenced by scale_atm_dyn_dgm_nonhydro3d_hevi_gmres::atm_dyn_dgm_nonhydro3d_hevi_cal_vi(), get_pmatd_lu(), and scale_atm_dyn_dgm_hydrostatic::hydrostatic_build_rho_xyz::hydrostatic_build_rho_xyz_moist().

linalgebra_solvelineq_bndmat()

subroutine, public scale_linalgebra::linalgebra_solvelineq_bndmat	(	real(rp), dimension(2*kl+ku+1,n), intent(inout)	a,
		real(rp), dimension(n,nrhs), intent(inout)	b,
		integer, dimension(n), intent(out)	ipiv,
		integer, intent(in)	n,
		integer, intent(in)	kl,
		integer, intent(in)	ku,
		integer, intent(in)	nrhs,
		logical, intent(in), optional	use_lapack )

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.

Parameters

A	Matrix A in band storage same as in dgbsv of LAPACK
b	Right hand side matrix B with N-by-RHS
ipiv	Pivot indices that define the permutation matrix
KL	Number of subdiagonals within the band of A
KU	Number of superdiagonals within the band of A
NRHS	Number of right hand sides
use_lapack	Flag whether LAPACK is used

Definition at line 208 of file scale_linalgebra.F90.

    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

Referenced by scale_atm_dyn_dgm_globalnonhydro3d_rhot_hevi::atm_dyn_dgm_globalnonhydro3d_rhot_hevi_cal_vi_asis(), scale_atm_dyn_dgm_nonhydro3d_rhot_hevi::atm_dyn_dgm_nonhydro3d_rhot_hevi_cal_vi_asis(), and scale_atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform::atm_dyn_dgm_nonhydro3d_rhot_hevi_splitform_cal_vi().

linalgebra_solvelineq_gmres()

subroutine, public scale_linalgebra::linalgebra_solvelineq_gmres	(	type(sparsemat), intent(in)	a,
		real(rp), dimension(:), intent(in)	b,
		real(rp), dimension(:), intent(inout)	x,
		integer, intent(in)	m,
		integer, intent(in)	restart_num,
		real(rp), intent(in)	conv_crit )

Solve a linear equation Ax=b using the Generalized Minimal Residual solver (GMRES)

Parameters

[in]	a	Object managing sparse matrix A
[in]	b	Right hand side vector
[in,out]	x	Solution vector
[in]	restart_num	Number of times by which restarting is performed if convergence condition is not satisfied
[in]	conv_crit	Threshold for finishing the iteration solver

Definition at line 320 of file scale_linalgebra.F90.

    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

References scale_sparsemat::sparsemat_storage_typeid_csr.