1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5
6module hecmw_solver
7contains
8
9  subroutine hecmw_solve (hecMESH, hecMAT)
10
11    use hecmw_util
12    use hecmw_solver_las_nn
13    use hecmw_solver_iterative
14    use hecmw_solver_direct
15    use hecmw_solver_direct_parallel
16    use hecmw_solver_direct_MUMPS
17    use hecmw_solver_direct_MKL
18    use hecmw_solver_direct_clusterMKL
19    use hecmw_matrix_misc
20    use m_hecmw_comm_f
21    implicit none
22
23    type (hecmwST_matrix), target :: hecMAT
24    type (hecmwST_local_mesh) :: hecMESH
25
26    real(kind=kreal) :: resid
27    integer(kind=kint) :: i, myrank, NDOF
28    integer(kind=kint) :: imsg = 51
29    NDOF=hecMAT%NDOF
30
31    !C ERROR CHECK
32    if(hecmw_solve_check_zerorhs(hecMESH, hecMAT))then
33      hecMAT%X = 0.0d0
34      return
35    endif
36
37    select case(hecMAT%Iarray(99))
38        !* Call Iterative Solver
39      case (1)
40        call hecmw_solve_iterative(hecMESH,hecMAT)
41        !* Call Direct Solver
42      case(2:)
43        !C
44        !* Please note the following:
45        !* Flag to activate symbolic factorization: 1(yes) 0(no)  hecMESH%Iarray(98)
46        !* Flag to activate numeric  factorization: 1(yes) 0(no)  hecMESH%Iarray(97)
47
48        if (hecMAT%Iarray(97) .gt. 1) hecMAT%Iarray(97)=1
49
50        call hecmw_mat_set_flag_converged(hecMAT, 0)
51        call hecmw_mat_set_flag_diverged(hecMAT, 0)
52
53        if (hecMAT%Iarray(2) .eq. 102) then
54          if(hecMESH%PETOT.GT.1) then
55            call hecmw_solve_direct_ClusterMKL(hecMESH, hecMAT)
56          else
57            call hecmw_solve_direct_MKL(hecMESH,hecMAT)
58          endif
59        elseif (hecMAT%Iarray(2) .eq. 104) then
60          call hecmw_solve_direct_MUMPS(hecMESH, hecMAT)
61        else
62          if(hecMESH%PETOT.GT.1) then
63            call hecmw_solve_direct_parallel(hecMESH,hecMAT,imsg)
64          else
65            call hecmw_solve_direct(hecMESH,hecMAT,imsg)
66          endif
67        endif
68
69        resid=hecmw_rel_resid_L2_nn(hecMESH,hecMAT)
70        myrank=hecmw_comm_get_rank()
71        if (myrank==0) then
72          if (hecMAT%Iarray(21) > 0 .or. hecMAT%Iarray(22) > 0) then
73            write(*,"(a,1pe12.5)")'### Relative residual =', resid
74          endif
75          if( resid >= 1.0d-8) then
76            write(*,"(a)")'### Relative residual exceeded 1.0d-8---Direct Solver### '
77            !            stop
78          endif
79        endif
80        if (resid < hecmw_mat_get_resid(hecMAT)) then
81          call hecmw_mat_set_flag_converged(hecMAT, 1)
82        endif
83        !C
84    end select
85
86  end subroutine hecmw_solve
87
88  subroutine hecmw_substitute_solver(hecMESH, hecMATorig, NDOF)
89
90    use hecmw_util
91    use hecmw_solver_iterative
92    use hecmw_solver_direct
93    use hecmw_solver_direct_parallel
94    use hecmw_solver_direct_MUMPS
95    use hecmw_solver_direct_clusterMKL
96    use hecmw_matrix_contact
97    type (hecmwST_local_mesh)     :: hecMESH
98    type (hecmwST_matrix)         :: hecMATorig
99    type (hecmwST_matrix),pointer :: hecMAT => null()
100    integer(kind=kint) NDOF
101    if (NDOF == hecMATorig%NDOF) then
102      call hecmw_clone_matrix(hecMATorig,hecMAT)
103    else if (NDOF < hecMATorig%NDOF) then
104      call hecmw_abort( hecmw_comm_get_comm() )
105    else
106      call hecmw_blockmatrix_expand(hecMATorig,hecMAT,NDOF)
107      call hecmw_cmat_init(hecMAT%cmat)
108    end if
109    !  select case(NDOF)
110    !    case(1)
111    !      call hecmw_solve_11(hecMESH,hecMAT)
112    !    case(2)
113    !      call hecmw_solve_22(hecMESH,hecMAT)
114    !    case(3)
115    !      call hecmw_solve_33(hecMESH,hecMAT)
116    !    case(4)
117    !      call hecmw_solve_iterative(hecMESH,hecMAT)
118    !    case(5)
119    !      call hecmw_solve_iterative(hecMESH,hecMAT)
120    !    case(6)
121    !      call hecmw_solve_66(hecMESH,hecMAT)
122    !    case(7:)
123    !      call hecmw_solve_iterative(hecMESH,hecMAT)
124    !  end select
125    !call hecmw_solve_direct_MUMPS(hecMESH, hecMAT)
126    call hecmw_solve_iterative(hecMESH,hecMAT)
127    if (NDOF /= hecMATorig%NDOF) then
128      call hecmw_vector_contract(hecMATorig,hecMAT,NDOF)
129    end if
130  end subroutine hecmw_substitute_solver
131
132end module hecmw_solver
133