scale_linalgebra Module Reference

module common / Linear algebra More...

Data Types
interface	linalgebra_solvelineq

Functions/Subroutines
real(rp) function, dimension(size(a, 1), size(a, 2)), public	linalgebra_inv (a)
subroutine, public	linalgebra_lu (a_lu, ipiv)
subroutine, public	linalgebra_solvelineq_gmres (a, b, x, m, restart_num, conv_crit)

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

)

Definition at line 66 of file scale_linalgebra.F90.

 
    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

Referenced by scale_element_base::elementbase_construct_massmat(), 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 )

Definition at line 160 of file scale_linalgebra.F90.

 
    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

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::hydrostaic_build_rho_xyz::hydrostaic_build_rho_xyz_moist().

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 )

Definition at line 180 of file scale_linalgebra.F90.

 
    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

References scale_sparsemat::sparsemat_storage_typeid_csr.