1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5!> Jagged Diagonal Matrix storage for vector processors.
6!> Original code was provided by JAMSTEC.
7
8module hecmw_JAD_TYPE_44
9  use hecmw_util
10  use m_hecmw_comm_f
11  implicit none
12
13  private
14
15  public :: hecmw_JAD_INIT_44
16  public :: hecmw_JAD_FINALIZE_44
17  public :: hecmw_JAD_IS_INITIALIZED_44
18  public :: hecmw_JAD_MATVEC_44
19
20  !C---------------------- AU&AL
21  real(kind=kreal), allocatable      :: AJAD(:)
22  integer(kind=kint), allocatable    :: JAJAD(:)
23  integer(kind=kint), allocatable    :: JADORD(:)
24  integer(kind=kint), allocatable    :: IAJAD(:)
25  integer(kind=kint) :: MJAD
26  real(kind=kreal), allocatable  :: WP1(:), WP2(:), WP3(:), WP4(:)
27  integer(kind=kint) :: INITIALIZED = 0
28
29contains
30
31  subroutine hecmw_JAD_INIT_44(hecMAT)
32    type(hecmwST_matrix) :: hecMAT
33    allocate(WP1(hecMAT%NP), WP2(hecMAT%NP), WP3(hecMAT%NP), WP4(hecMAT%NP))
34    allocate(AJAD((hecMAT%NPL+hecMAT%NPU)*16))
35    allocate(JAJAD(hecMAT%NPL+hecMAT%NPU))
36    allocate(JADORD(hecMAT%NP))
37    allocate(IAJAD(hecMAT%NP+1))
38    call REPACK(hecMAT%N, hecMAT, MJAD, AJAD, JAJAD, IAJAD, JADORD)
39    INITIALIZED = 1
40  end subroutine hecmw_JAD_INIT_44
41
42  subroutine hecmw_JAD_FINALIZE_44()
43    deallocate(AJAD)
44    deallocate(JAJAD)
45    deallocate(JADORD)
46    deallocate(IAJAD)
47    deallocate(WP1,WP2,WP3,WP4)
48    INITIALIZED = 0
49  end subroutine hecmw_JAD_FINALIZE_44
50
51  function hecmw_JAD_IS_INITIALIZED_44()
52    integer(kind=kint) :: hecmw_JAD_IS_INITIALIZED_44
53    hecmw_JAD_IS_INITIALIZED_44 = INITIALIZED
54  end function hecmw_JAD_IS_INITIALIZED_44
55
56  subroutine hecmw_JAD_MATVEC_44(hecMESH, hecMAT, X, Y, COMMtime)
57    type(hecmwST_local_mesh), intent(in) :: hecMESH
58    type(hecmwST_matrix), intent(in), target :: hecMAT
59    real(kind=kreal), intent(in) :: X(:)
60    real(kind=kreal), intent(out) :: Y(:)
61    real(kind=kreal), intent(inout) :: COMMtime
62    real(kind=kreal) :: START_TIME, END_TIME
63    real(kind=kreal), pointer :: D(:)
64    integer(kind=kint) :: i
65    real(kind=kreal) :: X1, X2, X3, X4
66
67    START_TIME= HECMW_WTIME()
68    call hecmw_update_4_R (hecMESH, X, hecMAT%NP)
69    END_TIME= HECMW_WTIME()
70    COMMtime = COMMtime + END_TIME - START_TIME
71
72    D => hecMAT%D
73
74    !$OMP PARALLEL PRIVATE(i)
75    !$OMP DO
76    do i= 1, hecMAT%N
77      X1= X(4*i-3)
78      X2= X(4*i-2)
79      X3= X(4*i-1)
80      X4= X(4*i  )
81      Y(4*i-3)= D(16*i-15)*X1 + D(16*i-14)*X2 + D(16*i-13)*X3 + D(16*i-12)*X4
82      Y(4*i-2)= D(16*i-11)*X1 + D(16*i-10)*X2 + D(16*i- 9)*X3 + D(16*i- 8)*X4
83      Y(4*i-1)= D(16*i- 7)*X1 + D(16*i- 6)*X2 + D(16*i- 5)*X3 + D(16*i- 4)*X4
84      Y(4*i  )= D(16*i- 3)*X1 + D(16*i- 2)*X2 + D(16*i- 1)*X3 + D(16*i- 0)*X4
85    enddo
86    !$OMP END DO
87    !$OMP END PARALLEL
88    call MATJAD(hecMAT%N, MJAD, IAJAD, JAJAD, AJAD, JADORD, X, Y, WP1, WP2, WP3, WP4)
89  end subroutine hecmw_JAD_MATVEC_44
90
91  subroutine REPACK(N, hecMAT, MJAD,  AJAD, JAJAD, IAJAD, JADORD)
92    use hecmw_util
93    !C---------------------------------
94    type (hecmwST_matrix)     :: hecMAT
95    !C----------------------
96    integer(kind = kint) :: N, MJAD
97    real(kind = kreal), dimension(*) :: AJAD
98    integer(kind = kint), dimension(*) :: JAJAD
99    integer(kind = kint), dimension(*) :: IAJAD
100    integer(kind = kint), dimension(*) :: JADORD
101
102    integer(kind = kint) :: IJAD, MAXNZ, MINNZ
103    integer(kind = kint) :: I, J, JS, JE, in, JC
104    integer(kind = kint), allocatable :: len(:), LENZ(:), JADREORD(:)
105
106    allocate(len(N))
107    allocate(JADREORD(N))
108    do I=1,N
109      len(I)= hecMAT%indexL(I) - hecMAT%indexL(I-1) &
110        &        + hecMAT%indexU(I) - hecMAT%indexU(I-1)
111    end do
112    MAXNZ=maxval(len(1:N))
113    MINNZ=minval(len(1:N))
114    MJAD =MAXNZ
115    allocate(LENZ(0:MJAD))
116    LENZ = 0
117    do I=1,N
118      LENZ(len(I))=LENZ(len(I))+1
119    enddo
120    do I=MAXNZ-1,MINNZ,-1
121      LENZ(I)=LENZ(I)+LENZ(I+1)
122    enddo
123    do I=1,N
124      JADORD(I)=LENZ(len(I))
125      LENZ(len(I))=LENZ(len(I))-1
126    enddo
127    do I=1,N
128      JADREORD(JADORD(I))=I
129    enddo
130    do I=1,N
131      LENZ(len(JADREORD(I)))=I
132    enddo
133    do I=MAXNZ-1,1,-1
134      LENZ(I)=max(LENZ(I+1),LENZ(I))
135    enddo
136    IAJAD(1)=1
137    do I=1,MAXNZ
138      IAJAD(I+1)=IAJAD(I)+LENZ(I)
139    enddo
140    len=0
141    do I= 1, N
142      IJAD=JADORD(I)
143      JS= hecMAT%indexL(I-1) + 1
144      JE= hecMAT%indexL(I  )
145      do J=JS,JE
146        in  = hecMAT%itemL(J)
147        len(IJAD)=len(IJAD)+1
148        JC=IAJAD(len(IJAD))+IJAD-1
149        AJAD(JC*16-15:JC*16) = hecMAT%AL(16*J-15:16*J)
150        JAJAD(JC) = in
151      end do
152    end do
153    do I= 1, N
154      IJAD=JADORD(I)
155      JS= hecMAT%indexU(I-1) + 1
156      JE= hecMAT%indexU(I  )
157      do J=JS,JE
158        in  = hecMAT%itemU(J)
159        len(IJAD)=len(IJAD)+1
160        JC=IAJAD(len(IJAD))+IJAD-1
161        AJAD(JC*16-15:JC*16) = hecMAT%AU(16*J-15:16*J)
162        JAJAD(JC) = in
163      end do
164    end do
165    deallocate(len)
166    deallocate(JADREORD)
167    deallocate(LENZ)
168  end subroutine REPACK
169
170  subroutine MATJAD(N, MJAD, IAJAD, JAJAD, AJAD, JADORD, X, Y, W1, W2, W3, W4)
171    use hecmw_util
172    integer(kind=kint) :: N, MJAD
173    integer(kind=kint) :: IAJAD(*), JAJAD(*), JADORD(*)
174    real(kind=kreal)   :: AJAD(*), X(*), Y(*), W1(*), W2(*), W3(*), W4(*)
175
176    integer(kind=kint) :: I, K, NZ, IXX
177    real(kind=kreal)   :: X1, X2, X3, X4
178
179    !$OMP PARALLEL PRIVATE(I,K,X1,X2,X3,X4,IXX)
180    !$OMP DO
181    do I=1,N
182      W1(I)=0.D0
183      W2(I)=0.D0
184      W3(I)=0.D0
185      W4(I)=0.D0
186    enddo
187    !$OMP END DO
188
189    do NZ=1,MJAD
190      !$OMP DO
191      do K=IAJAD(NZ),IAJAD(NZ+1)-1
192        X1=X(JAJAD(K)*4-3)
193        X2=X(JAJAD(K)*4-2)
194        X3=X(JAJAD(K)*4-1)
195        X4=X(JAJAD(K)*4  )
196        IXX = K-IAJAD(NZ)+1
197        W1(IXX)=W1(IXX) + AJAD(K*16-15)*X1 + AJAD(K*16-14)*X2 + AJAD(K*16-13)*X3 + AJAD(K*16-12)*X4
198        W2(IXX)=W2(IXX) + AJAD(K*16-11)*X1 + AJAD(K*16-10)*X2 + AJAD(K*16- 9)*X3 + AJAD(K*16- 8)*X4
199        W3(IXX)=W3(IXX) + AJAD(K*16- 7)*X1 + AJAD(K*16- 6)*X2 + AJAD(K*16- 5)*X3 + AJAD(K*16- 4)*X4
200        W4(IXX)=W4(IXX) + AJAD(K*16- 3)*X1 + AJAD(K*16- 2)*X2 + AJAD(K*16- 1)*X3 + AJAD(K*16- 0)*X4
201      enddo
202      !$OMP END DO
203    enddo
204
205    !$OMP DO
206    do I=1,N
207      Y(4*I-3)=Y(4*I-3)+W1(JADORD(I))
208      Y(4*I-2)=Y(4*I-2)+W2(JADORD(I))
209      Y(4*I-1)=Y(4*I-1)+W3(JADORD(I))
210      Y(4*I  )=Y(4*I  )+W4(JADORD(I))
211    enddo
212    !$OMP END DO
213    !$OMP END PARALLEL
214  end subroutine MATJAD
215
216end module hecmw_JAD_TYPE_44
217