test_solvers Program

program test_solvers
use omp_lib,   only :   omp_get_wtime
use data_arch, only :   I4, R4, R8, get_unit

use num_param, only :   OPU,          & ! *default output unit*
                        IPU,          & ! *default input unit*
                        SOLV_MESS,    & ! *Solver message: yes=```PRINT_MESS```, no=```NO_MESS```*
                        NO_MESS,      & ! *no output message*
                        PRINT_MESS      ! *solver infos output*

use solver,    only :   MAT_SOLV,     & ! *system data type*
                        solve_syst,   & ! *subroutine for the system resolution*
                        MUMP,         & ! *integer*
                        UMFP,         & ! *integer*
                        SULU,         & ! *integer*
                        convert_matrice_format ! *various matrix format conversion*
implicit none

integer(kind=I4) :: i, ii, uu, size_a_elt, state
real(kind=R8)    :: time1, time2, error
type(MAT_SOLV)   :: slv_struct

real(kind=R8),    dimension(:), allocatable :: b            !! *system right hand side*
real(kind=R8),    dimension(:), allocatable :: a_elt        !! *matrix non-zero entries in CC*
real(kind=R8),    dimension(:), allocatable :: a_elt_ref    !! *initial ```a_elt```*
integer(kind=I4), dimension(:), allocatable :: eltptr       !! *elemental matrices position in ```eltvar```*
integer(kind=I4), dimension(:), allocatable :: eltvar       !! *elemental matrix global lines*

   SOLV_MESS = NO_MESS!PRINT_MESS!
   
   write(OPU, *) 'small matrix? 0:n, 1:y'
   read( IPU, *) i
   
   ! ======================== SOLVER TYPE ===========================================
   do
      write(OPU, *) 'data read, choose solver n° : 0-MA48, 1-SULU, 2-MUMP (ref), 3-UMFP'
      read( IPU, *) ii
      slv_struct%slv_t = ii
      if (ii>=0 .and. ii<=3) exit
   enddo

   ! ======================== SYST DATA READ ===========================================
   call get_unit(uu)
   
   if (i==0) then
      open(uu, file="mat/big_syst.sys")
   else
      open(uu, file="mat/small_syst.sys")
   endif

   ! ======================== SYST INFO ===========================================
   read(uu, *) slv_struct%nn, slv_struct%ne, slv_struct%nvar, slv_struct%nt
   write(OPU, *) '*************** INFO ********************'
   write(OPU, *) 'system size:                  ', slv_struct%nn
   write(OPU, *) 'number of elemental matrices: ', slv_struct%ne
   write(OPU, *) 'number of nnz entries:        ', slv_struct%nt
   write(OPU, *) '*****************************************'
   slv_struct%first = .true.

   allocate( eltvar(slv_struct%nvar  ) )
   allocate( a_elt(     slv_struct%nt    ), stat = state ) ; if (state/=0) stop 'Memory allocation problem in PRG'
   allocate( a_elt_ref( slv_struct%nt    ), stat = state ) ; if (state/=0) stop 'Memory allocation problem in PRG'
   allocate( eltptr(slv_struct%ne +1 ) )
   allocate(      b(slv_struct%nn    ) )

   allocate( slv_struct%eltvar(slv_struct%nvar  ) )
   allocate( slv_struct%a_elt( slv_struct%nt    ), stat = state ) ; if (state/=0) stop 'Memory allocation problem in PRG'
   allocate( slv_struct%eltptr(slv_struct%ne +1 ) )

   call solve_syst(mat = slv_struct, step = 'ini')

   do ii = 1, slv_struct%nvar
      read(uu, *) slv_struct%eltvar(ii)
   enddo
   do ii = 1, slv_struct%ne +1
      read(uu, *) slv_struct%eltptr(ii)
   enddo
   do ii = 1, slv_struct%nt
      read(uu, *) slv_struct%a_elt(ii)
   enddo
   do ii = 1, slv_struct%nn
      read(uu, *) slv_struct%b(ii)
   enddo
   
   close(uu)

   ! The matrices are in CC HB format. Only MUMPS accepts elemental entries, so the
   ! following subroutine converts elemental marices to assembled vectors.
   ! If MUMPS is chosen, nothing is done.
   call convert_matrice_format(mat = slv_struct)
   ! ======================== backup ============================================
   eltvar     = slv_struct%eltvar
   a_elt_ref  = slv_struct%a_elt
   eltptr     = slv_struct%eltptr
   b          = slv_struct%b

   ! ======================== PROCESS ===========================================
   time1 = omp_get_wtime()
      call solve_syst(mat = slv_struct, step = 'ana') ; write(OPU, *) 'system analyzed' ; slv_struct%first = .false.
      call solve_syst(mat = slv_struct, step = 'fac') ; write(OPU, *) 'system factorized'
      call solve_syst(mat = slv_struct, step = 'sol') ; write(OPU, *) 'system solved'
      call solve_syst(mat = slv_struct, step = 'fre') ; write(OPU, *) 'system freed'
   time2 = omp_get_wtime()

   call verif_solution(slv_struct = slv_struct, a_elt = a_elt_ref, b = b, error = error)

   write(OPU,*) 'max error    = ', error
   write(OPU,*) 'elapsed time = ', real(abs(time1 -time2), kind=R4)

   ! Here, the matrix coefficients are modified, but the sparsity is conserved. It gives a look to the ability
   ! of the solver to exploit the symbolic calculations performed before.
   do

      if (slv_struct%slv_t==MUMP) then
         size_a_elt   = size(slv_struct%a_elt)
      else
         size_a_elt   = slv_struct%nt
      endif

      ! the original entries are retrieved, then modified
      a_elt = a_elt_ref
      call modify_a_elt(tab = a_elt, nz = size_a_elt)
      slv_struct%a_elt  = a_elt
      slv_struct%b      = b
      slv_struct%x      = 0._R8

      time1 = omp_get_wtime()
         call solve_syst(mat = slv_struct, step = 'fac')
         call solve_syst(mat = slv_struct, step = 'sol')
         call solve_syst(mat = slv_struct, step = 'fre')
      time2 = omp_get_wtime()

      call verif_solution(slv_struct = slv_struct, a_elt = a_elt, b = b, error = error)

      write(OPU,*) 'max error    = ', error
      write(OPU,*) 'elapsed time = ', real(abs(time1 -time2), kind=R4)

      write(OPU, *) 'STOP? (y=1, n=0)'
      read( IPU, *) ii
      if (ii==1) exit
   enddo

   deallocate(eltvar, a_elt, a_elt_ref, eltptr, b)

   call solve_syst(mat = slv_struct, step = 'end')
   
stop

contains

   !=========================================================================================
   !< @note The product \(\{y\}\) of the system matrix \([A]\) by the solution \(\{x\}\), is
   !        calculated, and compared to the right hand side \(\{b\}\).
   !        The calculated error is the absolute error in %.
   !-----------------------------------------------------------------------------------------
   subroutine verif_solution(slv_struct, a_elt, b, error)
   implicit none
   type(MAT_SOLV), intent(in) :: slv_struct
   real(kind=R8), dimension(:), intent(in) :: a_elt
   real(kind=R8), dimension(:), intent(in) :: b
   real(kind=R8), intent(out) :: error
      real(kind=R8), dimension(slv_struct%nn) :: y
      ! to assess the accuracy, the solution is applied to the
      ! system matrix and compared to the rhs.
      if (slv_struct%slv_t==MUMP) then
         call prod_elemental_x(n      = slv_struct%nn,      &
                               nz     = slv_struct%nt,      &
                               nelt   = slv_struct%ne,      &
                               nvar   = slv_struct%nvar,    &
                               x      = slv_struct%x,       &
                               y      = y,                  &
                               a_elt  = a_elt,              &
                               eltptr = slv_struct%eltptr,  &
                               eltvar = slv_struct%eltvar)
      else
         call prod_a_x(n  = slv_struct%nn,         &
                       nz = slv_struct%nz,         &
                       x  = slv_struct%x,          &
                       y  = y,                     &
                       a_elt = a_elt,              &
                       irow  = slv_struct%irow,    &
                       jptr  = slv_struct%jptr,    &
                       slvt  = slv_struct%slv_t)
      endif

      error = 100*maxval( abs(y(1:slv_struct%nn) -b(1:slv_struct%nn)) )/ &
                  maxval( abs(y(1:slv_struct%nn) +b(1:slv_struct%nn)) )
   return
   endsubroutine verif_solution
   

   !=========================================================================================
   !> @note multiplication of the system coefficient by a random factor
   !-----------------------------------------------------------------------------------------
   subroutine modify_a_elt(tab, nz)
   implicit none
   integer(kind=I4), intent(in) :: nz
   real(kind=R8), intent(inout) :: tab(1:nz)
      real(kind=R8), allocatable :: tmp(:)
      allocate(tmp(1:nz))
      call random_number(harvest=tmp(1:nz))
      tmp(1:nz) = 2*tmp(1:nz) -1.0_R8
      tab(1:nz) = tab(1:nz)*tmp(1:nz)
      deallocate(tmp)
   return
   endsubroutine modify_a_elt

   
   !=========================================================================================
   !> @note \([A] \{x\}\), assembled CC format
   !-----------------------------------------------------------------------------------------
   subroutine prod_a_x(n, nz, x, y, a_elt, irow, jptr, slvt)
   implicit none
   integer(kind=I4), intent(in) :: n, nz, slvt
   real(kind=R8), dimension(nz), intent(in)     :: a_elt
   integer(kind=I4), dimension(nz ), intent(in) :: irow
   integer(kind=I4), dimension(n+1), intent(in) :: jptr
   real(kind=R8), dimension(n), intent(in)      :: x
   real(kind=R8), dimension(n), intent(out)     :: y
      integer(kind=I4) :: i, k, dec
      y(1:n) = 0._R8
      dec = 0
      if (slvt==UMFP .or. slvt==SULU) dec = 1
      do i = 1, n
         do k = jptr(i), jptr(i+1)-1
            y(irow(k +dec) +dec) = y(irow(k +dec) +dec) + x(i)*a_elt(k +dec)
         enddo
      enddo
   return
   endsubroutine prod_a_x

   
   !=========================================================================================
   !> @note \([A] \{x\}\), elemental CC format
   !-----------------------------------------------------------------------------------------
   subroutine prod_elemental_x(n, nz, nelt, nvar, x, y, a_elt, eltptr, eltvar)
   implicit none
   integer(kind=I4), intent(in) :: n, nz, nelt, nvar
   real(kind=R8),    dimension(nz     ), intent(in) :: a_elt
   integer(kind=I4), dimension(nelt +1), intent(in) :: eltptr
   integer(kind=I4), dimension(nvar   ), intent(in) :: eltvar
   real(kind=R8), dimension(n), intent(in)      :: x
   real(kind=R8), dimension(n), intent(out)     :: y
      integer(kind=I4) :: i, j, k, kk, inc_nz, inc_nn, n_elem, i_elem, max_n_elem
      real(kind=R8), dimension(:), allocatable :: a_elt_tmp, x_tmp, y_tmp
      inc_nz = 0
      inc_nn = 0
      n_elem = 0
      y      = 0._R8
      
      max_n_elem = 0
      do i_elem = 1, nelt
         max_n_elem = max(max_n_elem, eltptr(i_elem +1) -eltptr(i_elem))
      enddo

      allocate( a_elt_tmp(1:max_n_elem**2) )
      allocate(     x_tmp(1:max_n_elem   ) )
      allocate(     y_tmp(1:max_n_elem   ) )

      do i_elem = 1, nelt                                               ! browse all elemental matrices
         inc_nn = inc_nn +n_elem                                        ! step in eltvar for matrix number i_elem
         inc_nz = inc_nz +n_elem**2                                     ! step in a_elt for matrix number i_elem
         n_elem = eltptr(i_elem +1) -eltptr(i_elem)                     ! elemental matrix size
         a_elt_tmp(1:n_elem**2) = a_elt(inc_nz +1 : inc_nz +n_elem**2)  ! elemental matrix coefficients
         ! --- elemental rhs
         kk = 0
         do k = inc_nn +1, inc_nn +n_elem
            kk = kk +1
            x_tmp(kk) = x(eltvar(k))
         enddo
         ! --- elemental product
         y_tmp(1:n_elem) = 0._R8
         do i = 1, n_elem
            do j = 1, n_elem
               y_tmp(i) = y_tmp(i) +a_elt_tmp(i +n_elem*(j-1))*x_tmp(j)
            enddo
         enddo
         ! --- elemental product in global vector y
         kk = 0
         do k = inc_nn +1, inc_nn +n_elem
            kk = kk +1
            y(eltvar(k)) = y(eltvar(k)) +y_tmp(kk)
         enddo
      enddo
      deallocate(a_elt_tmp, x_tmp, y_tmp)
   return
   endsubroutine prod_elemental_x

endprogram test_solvers