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