1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5!> This module provides linear equation solver interface for Pardiso
6
7module hecmw_solver_direct_MKL
8  use hecmw_util
9  use m_sparse_matrix
10  use m_sparse_matrix_hec
11  use m_hecmw_MKL_wrapper
12  use hecmw_matrix_ass
13  use hecmw_matrix_dump
14
15  private
16  public :: hecmw_solve_direct_MKL
17
18  logical, save :: INITIALIZED = .false.
19  type (sparse_matrix), save :: spMAT
20
21contains
22
23  subroutine hecmw_solve_direct_MKL(hecMESH,hecMAT)
24    implicit none
25    type (hecmwST_local_mesh), intent(in) :: hecMESH
26    type (hecmwST_matrix    ), intent(inout) :: hecMAT
27    integer(kind=kint) :: spmat_type
28    integer(kind=kint) :: spmat_symtype
29    integer(kind=kint) :: phase_start
30    integer(kind=kint) :: istat,myrank
31    real(kind=kreal) :: t1,t2
32
33    t1=hecmw_wtime()
34
35    call hecmw_mat_dump(hecMAT, hecMESH)
36
37    myrank=hecmw_comm_get_rank()
38
39    if (INITIALIZED .and. hecMAT%Iarray(98) .eq. 1) then
40      call sparse_matrix_finalize(spMAT)
41      INITIALIZED = .false.
42    endif
43
44    if (.not. INITIALIZED) then
45      spmat_type = SPARSE_MATRIX_TYPE_CSR
46      if (hecMAT%symmetric) then
47        spmat_symtype = SPARSE_MATRIX_SYMTYPE_SYM
48      else
49        spmat_symtype = SPARSE_MATRIX_SYMTYPE_ASYM
50      end if
51      call sparse_matrix_set_type(spMAT, spmat_type, spmat_symtype)
52
53      INITIALIZED = .true.
54      hecMAT%Iarray(98) = 1
55    endif
56
57    !* Flag to activate symbolic factorization: 1(yes) 0(no)  hecMESH%Iarray(98)
58    !* Flag to activate numeric  factorization: 1(yes) 0(no)  hecMESH%Iarray(97)
59
60    if (hecMAT%Iarray(98) .eq. 1) then
61      ! ANALYSIS and FACTORIZATION
62      call sparse_matrix_hec_init_prof(spMAT, hecMAT, hecMESH)
63      call sparse_matrix_hec_set_vals(spMAT, hecMAT)
64      !call sparse_matrix_dump(spMAT)
65      phase_start = 1
66      hecMAT%Iarray(98) = 0
67      hecMAT%Iarray(97) = 0
68    endif
69    if (hecMAT%Iarray(97) .eq. 1) then
70      ! FACTORIZATION
71      call sparse_matrix_hec_set_vals(spMAT, hecMAT)
72      !call sparse_matrix_dump(spMAT)
73      phase_start = 2
74      hecMAT%Iarray(97) = 0
75    endif
76    call sparse_matrix_hec_set_rhs(spMAT, hecMAT)
77
78    t2=hecmw_wtime()
79    if (myrank==0 .and. (spMAT%iterlog > 0 .or. spMAT%timelog > 0)) &
80       write(*,'(A,f10.3)') ' [Pardiso]: Setup completed.          time(sec)=',t2-t1
81
82    ! SOLVE
83    call hecmw_mkl_wrapper(spMAT, phase_start, hecMAT%X, istat)
84
85    call hecmw_mat_dump_solution(hecMAT)
86
87  end subroutine hecmw_solve_direct_MKL
88
89end module hecmw_solver_direct_MKL
90