1!------------------------------------------------------------------------------- 2! Copyright (c) 2019 FrontISTR Commons 3! This software is released under the MIT License, see LICENSE.txt 4!------------------------------------------------------------------------------- 5 6module hecmw_solver_las_11 7 8 private 9 10 public :: hecmw_matvec_11 11 public :: hecmw_matresid_11 12 public :: hecmw_rel_resid_L2_11 13 14contains 15 16 !C 17 !C*** 18 !C*** hecmw_matvec_11 19 !C*** 20 !C 21 subroutine hecmw_matvec_11 (hecMESH, hecMAT, X, Y, n, COMMtime) 22 use hecmw_util 23 use m_hecmw_comm_f 24 25 implicit none 26 integer(kind=kint):: n 27 real(kind=kreal), dimension(n) :: X, Y 28 type (hecmwST_matrix) :: hecMAT 29 type (hecmwST_local_mesh) :: hecMESH 30 real(kind=kreal), optional :: COMMtime 31 32 real(kind=kreal) :: START_TIME, END_TIME 33 integer(kind=kint) :: i, j, jS, jE, in 34 real(kind=kreal) :: YV 35 36 START_TIME= HECMW_WTIME() 37 call hecmw_update_1_R (hecMESH, X, hecMAT%NP) 38 END_TIME= HECMW_WTIME() 39 if (present(COMMtime)) COMMtime = COMMtime + END_TIME - START_TIME 40 41 do i= 1, hecMAT%N 42 YV= hecMAT%D(i) * X(i) 43 jS= hecMAT%indexL(i-1) + 1 44 jE= hecMAT%indexL(i ) 45 do j= jS, jE 46 in= hecMAT%itemL(j) 47 YV= YV + hecMAT%AL(j) * X(in) 48 enddo 49 jS= hecMAT%indexU(i-1) + 1 50 jE= hecMAT%indexU(i ) 51 do j= jS, jE 52 in= hecMAT%itemU(j) 53 YV= YV + hecMAT%AU(j) * X(in) 54 enddo 55 Y(i)= YV 56 enddo 57 58 end subroutine hecmw_matvec_11 59 60 !C 61 !C*** 62 !C*** hecmw_matresid_11 63 !C*** 64 !C 65 subroutine hecmw_matresid_11 (hecMESH, hecMAT, X, B, R, COMMtime) 66 use hecmw_util 67 68 implicit none 69 real(kind=kreal) :: X(:), B(:), R(:) 70 type (hecmwST_matrix) :: hecMAT 71 type (hecmwST_local_mesh) :: hecMESH 72 real(kind=kreal), optional :: COMMtime 73 74 integer(kind=kint) :: i 75 real(kind=kreal) :: Tcomm 76 77 Tcomm = 0.d0 78 call hecmw_matvec_11 (hecMESH, hecMAT, X, R, hecMAT%NP, Tcomm) 79 if (present(COMMtime)) COMMtime = COMMtime + Tcomm 80 81 do i = 1, hecMAT%N 82 R(i) = B(i) - R(i) 83 enddo 84 85 end subroutine hecmw_matresid_11 86 87 !C 88 !C*** 89 !C*** hecmw_rel_resid_L2_11 90 !C*** 91 !C 92 function hecmw_rel_resid_L2_11 (hecMESH, hecMAT, COMMtime) 93 use hecmw_util 94 use hecmw_solver_misc 95 96 implicit none 97 real(kind=kreal) :: hecmw_rel_resid_L2_11 98 type ( hecmwST_local_mesh ), intent(in) :: hecMESH 99 type ( hecmwST_matrix ), intent(in) :: hecMAT 100 real(kind=kreal), optional :: COMMtime 101 102 real(kind=kreal) :: r(hecMAT%NDOF*hecMAT%NP) 103 real(kind=kreal) :: bnorm2, rnorm2 104 real(kind=kreal) :: Tcomm 105 106 Tcomm = 0.d0 107 call hecmw_InnerProduct_R(hecMESH, hecMAT%NDOF, hecMAT%B, hecMAT%B, bnorm2, Tcomm) 108 if (bnorm2 == 0.d0) then 109 bnorm2 = 1.d0 110 endif 111 call hecmw_matresid_11(hecMESH, hecMAT, hecMAT%X, hecMAT%B, r, Tcomm) 112 call hecmw_InnerProduct_R(hecMESH, hecMAT%NDOF, r, r, rnorm2, Tcomm) 113 if (present(COMMtime)) COMMtime = COMMtime + Tcomm 114 hecmw_rel_resid_L2_11 = sqrt(rnorm2 / bnorm2) 115 116 end function hecmw_rel_resid_L2_11 117 118end module hecmw_solver_las_11 119