submodule(direct_solver_m) direct_solver_MKL_s !! Direct solver based on MKL PARDISO library use error_handling_m, only : handle_error, error_info_t use status_codes_m, only : PARALLAX_ERR_SOLVER2D implicit none contains module subroutine create_MKL(self, a) class(direct_solver_MKL_t), intent(inout) :: self type(csrmat_t), intent(inout) :: a self%msglvl = 0 ! set by default to no output self%mtype = 11 ! type of matrix: real non-symmetric #ifdef ENABLE_MKL call pardisoinit( self%pt, self%mtype, self%iparm) #endif #ifndef DOUBLE_PREC ! set pardiso to single precision self%iparm(28) = 1 #endif ! solver created self%phase = 1 self%ndim = a%ndim allocate(self%perm(a%ndim)) end subroutine module subroutine factorise_MKL(self, a, task) class(direct_solver_MKL_t), intent(inout) :: self type(csrmat_t), intent(in) :: a character(len=*), intent(in) :: task integer :: err, i real(FP), dimension(a%ndim) :: rdmy1, rdmy2 select case(task) case('ANALYSE') self%phase = 11 case('NUMFAC') self%phase = 22 case('ANALYSE_AND_NUMFAC') self%phase = 12 case default call handle_error('invalid task', & PARALLAX_ERR_SOLVER2D, __LINE__, __FILE__) end select #ifdef ENABLE_MKL call pardiso( self%pt, & maxfct = 1, & mnum = 1, & mtype = self%mtype, & phase = self%phase, & n = a%ndim, & a = a%val, & ia = a%i, & ja = a%j, & perm = self%perm, & nrhs = 1, & iparm = self%iparm, & msglvl = self%msglvl,& b = rdmy1, & ! dummy array not accessed x = rdmy2, & ! dummy array not accessed error = err ) if (err /= 0) then call handle_error('factorisation failed', & PARALLAX_ERR_SOLVER2D, __LINE__, __FILE__, & additional_info=& error_info_t('pardiso_error=',[err])) endif #endif end subroutine module subroutine backsub_MKL(self, a, b, x) class(direct_solver_MKL_t), intent(inout) :: self type(csrmat_t), intent(in) :: a real(FP), dimension(a%ndim), intent(inout) :: b real(FP), dimension(a%ndim), intent(out) :: x integer :: err self%phase = 33 #ifdef ENABLE_MKL call pardiso( self%pt, & maxfct = 1, & mnum = 1, & mtype = self%mtype, & phase = self%phase, & n = a%ndim, & a = a%val, & ia = a%i, & ja = a%j, & perm = self%perm, & nrhs = 1, & iparm = self%iparm, & msglvl = self%msglvl,& b = b, & x = x, & error = err ) if (err /= 0) then call handle_error('backsubstitution failed', & PARALLAX_ERR_SOLVER2D, __LINE__, __FILE__, & additional_info=& error_info_t('pardiso_error=',[err])) endif #endif end subroutine module subroutine destructor_MKL(self) type(direct_solver_MKL_t), intent(inout) :: self integer :: err real(FP), dimension(1) :: rdmy1, rdmy2, rdmy3 integer, dimension(1) :: idmy1, idmy2 if ((self%phase /= 1) .and. (self%phase /= 12) .and. (self%phase /= 33)) then ! Pardiso not initialised self%phase = -1 else self%phase = -1 rdmy1 = 0.0_FP rdmy2 = 0.0_FP rdmy3 = 0.0-FP #ifdef ENABLE_MKL call pardiso( self%pt, & maxfct = 1, & mnum = 1, & mtype = self%mtype, & phase = self%phase, & n = self%ndim, & a = rdmy1, & ! dummy array not accessed ia = idmy1, & ! dummy array not accessed ja = idmy2, & ! dummy array not accessed perm = self%perm, & nrhs = 1, & iparm = self%iparm, & msglvl = self%msglvl,& b = rdmy2, & ! dummy array not accessed x = rdmy3, & ! dummy array not accessed error = err ) #endif endif if (allocated(self%perm)) then deallocate(self%perm) endif self%ndim = 0 end subroutine end submodule