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