sll_m_pastix.F90

Go to the documentation of this file.

module sll_m_pastix
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#define MPI_MASTER 0
#include "pastix_fortran.h"
#include "sll_working_precision.h"
#include "sll_memory.h"
 
   use sll_m_collective, only: &
      sll_s_boot_collective, &
      sll_f_get_collective_rank, &
      sll_f_get_collective_size, &
      sll_v_world_collective
 
   implicit none
 
   public :: &
      pastix_solver, &
      initialize, &
      factorize, &
      solve, &
      delete
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   type :: pastix_solver
      pastix_data_ptr_t        :: pastix_data 
      pastix_int_t             :: ncols       
      pastix_int_t, pointer :: colptr(:)   
      pastix_int_t, pointer :: row(:)      
      pastix_float_t, pointer :: avals(:)    
      pastix_int_t, pointer :: perm(:)     
      pastix_int_t, pointer :: invp(:)     
      pastix_float_t, pointer :: rhs(:)      
      pastix_int_t             :: nrhs        
      pastix_int_t             :: iparm(iparm_size) ! sll_int32 parameters
      double precision         :: dparm(dparm_size) ! Floating poin parameters
      sll_int32                :: nbthread    
      sll_int32                :: verbose     
   end type pastix_solver
 
   interface initialize
      module procedure init_pastix
   end interface initialize
 
   interface solve
      module procedure solve_pastix_with_rhs
      module procedure solve_pastix_without_rhs
   end interface solve
 
   interface factorize
      module procedure factorize_pastix
   end interface factorize
 
   interface delete
      module procedure free_pastix
   end interface delete
 
contains
 
   subroutine init_pastix(self, n, nnzeros, row_ptr, col_ind, val)
 
      type(pastix_solver)              :: self
      sll_int32, intent(in)         :: n
      pastix_int_t, intent(in)         :: nnzeros
      sll_int32, intent(in), target :: row_ptr(:)
      sll_int32, intent(in), target :: col_ind(:)
      sll_real64, intent(in), target :: val(:)
      sll_int32                        :: error
      sll_int32                        :: comm
      sll_int32                        :: prank
      sll_int32                        :: psize
 
      if (.not. associated(sll_v_world_collective)) then
         call sll_s_boot_collective()
      end if
      self%ncols = n
      prank = sll_f_get_collective_rank(sll_v_world_collective)
      psize = sll_f_get_collective_size(sll_v_world_collective)
      comm = sll_v_world_collective%comm
 
      ! Get options ftom command line
      self%nbthread = 1
      self%verbose = api_verbose_no
 
      ! Allocating
      self%colptr => row_ptr
      self%row => col_ind
      self%avals => val
      sll_clear_allocate(self%rhs(nnzeros), error)
 
      ! First PaStiX call to initiate parameters
 
      self%pastix_data = 0
      self%nrhs = 1
      self%iparm(iparm_modify_parameter) = api_no
      self%iparm(iparm_start_task) = api_task_init
      self%iparm(iparm_end_task) = api_task_init
 
      sll_allocate(self%perm(n), error)
      sll_allocate(self%invp(n), error)
 
      call pastix_fortran(self%pastix_data, comm, n, self%colptr, self%row, &
                          self%avals, self%perm, self%invp, self%rhs, &
                          self%nrhs, self%iparm, self%dparm)
 
      ! Customize some parameters
 
      self%iparm(iparm_thread_nbr) = self%nbthread
      self%iparm(iparm_verbose) = self%verbose
 
      self%iparm(iparm_sym) = api_sym_no !API_SYM_YES !API_SYM_NO
      self%iparm(iparm_factorization) = api_fact_lu !API_FACT_LDLT
 
      self%iparm(iparm_matrix_verification) = api_yes
 
      !The matrix is in CSR format, we transpose to get CSC
      self%iparm(iparm_transpose_solve) = api_yes
 
   end subroutine init_pastix
 
   subroutine factorize_pastix(self)
 
      type(pastix_solver) :: self
      sll_int32           :: comm
 
      comm = sll_v_world_collective%comm
      ! Call PaStiX first steps (Scotch - Fax - Blend)
      self%iparm(iparm_start_task) = api_task_ordering
      self%iparm(iparm_end_task) = api_task_analyse
 
      call pastix_fortran(self%pastix_data, &
                          comm, &
                          self%ncols, &
                          self%colptr, &
                          self%row, &
                          self%avals, &
                          self%perm, &
                          self%invp, &
                          self%rhs, &
                          self%nrhs, &
                          self%iparm, &
                          self%dparm)
 
      ! Call PaStiX factorization
      self%iparm(iparm_start_task) = api_task_numfact
      self%iparm(iparm_end_task) = api_task_numfact
 
      call pastix_fortran(self%pastix_data, &
                          comm, &
                          self%ncols, &
                          self%colptr, &
                          self%row, &
                          self%avals, &
                          self%perm, &
                          self%invp, &
                          self%rhs, &
                          self%nrhs, &
                          self%iparm, &
                          self%dparm)
 
   end subroutine factorize_pastix
 
   subroutine solve_pastix_without_rhs(self, sol)
 
      type(pastix_solver)              :: self
      sll_real64, dimension(:), target :: sol
      sll_int32                        :: comm
 
      comm = sll_v_world_collective%comm
 
      self%rhs => sol
 
      ! Call PaStiX updown and refinement
      self%iparm(iparm_start_task) = api_task_solve
      self%iparm(iparm_end_task) = api_task_refine
      ! rhs will be changed to solution
      call pastix_fortran(self%pastix_data, comm, &
                          self%ncols, self%colptr, self%row, &
                          self%avals, self%perm, self%invp, &
                          self%rhs, self%nrhs, self%iparm, self%dparm)
 
      sol = self%rhs
 
   end subroutine solve_pastix_without_rhs
 
   subroutine solve_pastix_with_rhs(self, rhs, sol)
 
      type(pastix_solver)      :: self
      sll_real64, dimension(:) :: rhs
      sll_real64, dimension(:) :: sol
      sll_int32                :: comm
 
      comm = sll_v_world_collective%comm
 
      self%rhs = rhs
 
      ! Call PaStiX updown and refinement
      self%iparm(iparm_start_task) = api_task_solve
      self%iparm(iparm_end_task) = api_task_refine
      ! rhs will be changed to solution
      call pastix_fortran(self%pastix_data, comm, &
                          self%ncols, self%colptr, self%row, &
                          self%avals, self%perm, self%invp, &
                          self%rhs, self%nrhs, self%iparm, self%dparm)
 
      sol = self%rhs
 
   end subroutine solve_pastix_with_rhs
 
   subroutine free_pastix(self)
 
      type(pastix_solver) :: self
      sll_int32           :: comm
 
      comm = sll_v_world_collective%comm
 
      ! Call PaStiX clean
      self%iparm(iparm_start_task) = api_task_clean
      self%iparm(iparm_end_task) = api_task_clean
 
      call pastix_fortran(self%pastix_data, comm, &
                          self%ncols, self%colptr, self%row, &
                          self%avals, self%perm, self%invp, &
                          self%rhs, self%nrhs, self%iparm, self%dparm)
 
      deallocate (self%colptr)
      deallocate (self%row)
      deallocate (self%avals)
      deallocate (self%perm)
      deallocate (self%invp)
      deallocate (self%rhs)
 
   end subroutine free_pastix
 
end module sll_m_pastix