linear_algebra.f90

! linear_algebra.f90 --- basic tools from linear algebra
!
! Copyright  (C)  2010  Florian Wilhelm <Florian.Wilhelm@kit.edu>
!
! Version: 1.0
! Keywords: scales ppm solver linear algebra residual
! Author: Florian Wilhelm <Florian.Wilhelm@kit.edu>
! Maintainer: Florian Wilhelm <Florian.Wilhelm@kit.edu>
! URL: https://www.dkrz.de/redmine/projects/show/scales-ppm
!
! Redistribution and use in source and binary forms, with or without
! modification, are  permitted provided that the following conditions are
! met:
!
! Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! Redistributions in binary form must reproduce the above copyright
! notice, this list of conditions and the following disclaimer in the
! documentation and/or other materials provided with the distribution.
!
! Neither the name of the DKRZ GmbH nor the names of its contributors
! may be used to endorse or promote products derived from this software
! without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
! Commentary:
!
! Code:
!

MODULE linear_algebra
    USE solver_internal, ONLY: config
    USE mo_kind, ONLY: wp, dp
    USE ppm_extents, ONLY: extent_type => extent, extent_start, &
         extent_end, extent_size
    USE ppm_base, ONLY: abort_ppm

    IMPLICIT NONE

    PRIVATE

    !-------------------------------------------------------------------
    ! Module Function & Procedure Interface
    !-------------------------------------------------------------------

    INTERFACE arr_dotproduct
        MODULE PROCEDURE arr_dotproduct1
        MODULE PROCEDURE arr_dotproduct2
    END INTERFACE

    PUBLIC:: calc_abs_res, calc_rel_res, arr_dotproduct, global_sum, arr_norm_2

    CONTAINS

     ! Array-wise dotproduct for one matrix with itself
    FUNCTION arr_dotproduct1(x, global_opt) RESULT(ans)
        REAL(wp), INTENT(IN) :: x(:,:)
        LOGICAL, INTENT(IN), OPTIONAL :: global_opt
        REAL(wp) :: ans
        LOGICAL :: global

        IF ( PRESENT(global_opt) ) THEN
            global = global_opt
        ELSE
            global = config%do_exchange
        ENDIF

        ans = SUM(x**2)

        IF (global) ans = global_sum(ans)

    END FUNCTION arr_dotproduct1

    ! Array-wise dotproduct for two different matrices
    FUNCTION arr_dotproduct2(x, y, global_opt) RESULT(ans)
#ifdef BLAS
        REAL(wp) :: DDOT
#endif
        REAL(wp), INTENT(IN) :: x(:,:), y(:,:)
        LOGICAL, INTENT(IN), OPTIONAL :: global_opt
        REAL(wp) :: ans
        REAL(wp), ALLOCATABLE :: xy(:)
        LOGICAL :: global
        INTEGER :: j, ie, je

        IF ( PRESENT(global_opt) ) THEN
            global = global_opt
        ELSE
            global = config%do_exchange
        ENDIF

        ie = SIZE(x,1)
        je = SIZE(x,2)

#ifndef BLAS
        ALLOCATE(xy(je))

        DO j=1,je
            xy(j) = DOT_PRODUCT(x(:,j), y(:,j))
        ENDDO

        ans = SUM(xy)
        DEALLOCATE(xy)
#else
        ans = DDOT(ie*je, x, 1, y, 1)
#endif
        IF (global) ans = global_sum(ans)

    END FUNCTION arr_dotproduct2

    ! Array-wise 2-norm
    FUNCTION arr_norm_2(x, global_opt) RESULT(ans)
        REAL(wp), INTENT(IN) :: x(:,:)
        LOGICAL, INTENT(IN), OPTIONAL :: global_opt
        REAL(wp) :: ans
        LOGICAL :: global

        IF ( PRESENT(global_opt) ) THEN
            global = global_opt
        ELSE
            global = config%do_exchange
        ENDIF

        ans = SQRT(arr_dotproduct(x, global))

    END FUNCTION arr_norm_2

    ! Sums a variable globally up
    FUNCTION global_sum(summand, comm_opt) RESULT(all_sum)
#ifdef USE_MPI
        USE solver_internal, ONLY: mpi_real_dp, mpi_comm_world, MPI_SUM
#endif

        REAL(wp), INTENT(in) :: summand
        INTEGER, OPTIONAL, INTENT(in) :: comm_opt
        REAL(wp) :: all_sum
#ifdef USE_MPI
        INTEGER :: comm, ierror

        IF (PRESENT(comm_opt)) THEN
            comm = comm_opt
        ELSE
            comm = mpi_comm_world
        ENDIF

        CALL MPI_ALLREDUCE (summand, all_sum, 1, mpi_real_dp, MPI_SUM, comm, ierror)
#else
        all_sum = summand
#endif

    END FUNCTION global_sum

    ! Calculate absolute residual of system
    ! If global_opt is .TRUE. then calculate the global absolute residual
    FUNCTION calc_abs_res(A, b, x, ext_x, global_opt) RESULT(abs_res)
        REAL(wp), INTENT(IN) :: x(:,:)
        REAL(wp), INTENT(IN) :: b(:,:)
        TYPE(extent_type), INTENT(IN) :: ext_x(:)      ! extent of x
        LOGICAL, INTENT(IN), OPTIONAL :: global_opt
        REAL(wp), ALLOCATABLE :: r(:,:), Ax(:,:)
        INTEGER :: ie, je, ib, jb, il, jl
        REAL(wp) :: abs_res
        LOGICAL :: global

        INTERFACE
            SUBROUTINE A(field, res_field)
                USE mo_kind, ONLY: wp, dp
                REAL(wp), INTENT(IN) :: field(:,:)
                REAL(wp), INTENT(OUT) :: res_field(:,:)
            END SUBROUTINE
        END INTERFACE

        IF ( PRESENT(global_opt) ) THEN
            global = global_opt
        ELSE
            global = config%do_exchange
        ENDIF

        ie = SIZE(x,1)
        je = SIZE(x,2)
        ib = extent_start(ext_x(1))
        jb = extent_start(ext_x(2))
        il = extent_end(ext_x(1))
        jl = extent_end(ext_x(2))

        ALLOCATE(r(ie,je), Ax(ie,je))

        CALL A(x, Ax)
        r(ib:il,jb:jl) = b(ib:il,jb:jl) - Ax(ib:il,jb:jl)

        abs_res = arr_norm_2(r(ib:il,jb:jl), global)

        DEALLOCATE(Ax, r)

    END FUNCTION calc_abs_res

    ! Calculate relative residual of system
    ! If global_opt is .TRUE. then calculate the global relative residual
    FUNCTION calc_rel_res(A, b, x, ext_x, global_opt, x0_opt) RESULT(rel_res)
        REAL(wp), INTENT(IN) :: x(:,:)
        REAL(wp), INTENT(IN) :: b(:,:)
        TYPE(extent_type), INTENT(IN) :: ext_x(:)      ! extent of x
        LOGICAL, INTENT(IN), OPTIONAL :: global_opt
        REAL(wp), INTENT(IN), OPTIONAL :: x0_opt(:,:)
        REAL(wp), ALLOCATABLE :: r(:,:), Ax(:,:), x0(:,:)
        INTEGER :: ie, je
        REAL(wp) :: rel_res, numer, denom

        INTERFACE
            SUBROUTINE A(field, res_field)
                USE mo_kind, ONLY: wp, dp
                REAL(wp), INTENT(IN) :: field(:,:)
                REAL(wp), INTENT(OUT) :: res_field(:,:)
            END SUBROUTINE
        END INTERFACE

        IF ( PRESENT(x0_opt) ) THEN
            x0 = x0_opt
        ELSE
            ie = SIZE(x,1)
            je = SIZE(x,2)
            ALLOCATE(x0(ie,je))
            x0 = 0._wp
        ENDIF

        IF ( PRESENT(global_opt) ) THEN
            numer = calc_abs_res(A, b, x, ext_x, global_opt)
            denom = calc_abs_res(A, b, x0, ext_x, global_opt)
        ELSE
            numer = calc_abs_res(A, b, x, ext_x)
            denom = calc_abs_res(A, b, x0, ext_x)
        ENDIF

        rel_res = numer/denom

        IF (.NOT. PRESENT(x0_opt)) THEN
            DEALLOCATE(x0)
        ENDIF

    END FUNCTION calc_rel_res

END MODULE linear_algebra