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