FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
hecmw_jadm_44.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
7 
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
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 
29 contains
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 
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_r (hecmesh, x, hecmat%NP, 4)
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)
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 
216 end module hecmw_jad_type_44
hecmw_util::hecmw_wtime
real(kind=kreal) function hecmw_wtime()
Definition: hecmw_util_f.F90:549
hecmw_jad_type_44::hecmw_jad_is_initialized_44
integer(kind=kint) function, public hecmw_jad_is_initialized_44()
Definition: hecmw_jadm_44.f90:52
hecmw_jad_type_44::hecmw_jad_init_44
subroutine, public hecmw_jad_init_44(hecMAT)
Definition: hecmw_jadm_44.f90:32
hecmw_jad_type_44::matjad
subroutine matjad(N, MJAD, IAJAD, JAJAD, AJAD, JADORD, X, Y, W1, W2, W3, W4)
Definition: hecmw_jadm_44.f90:171
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
hecmw_util::hecmwst_local_mesh
Definition: hecmw_util_f.F90:234
m_hecmw_comm_f
Definition: hecmw_comm_f.F90:6
hecmw_jad_type_44
Jagged Diagonal Matrix storage for vector processors. Original code was provided by JAMSTEC.
Definition: hecmw_jadm_44.f90:8
hecmw_jad_type_44::hecmw_jad_matvec_44
subroutine, public hecmw_jad_matvec_44(hecMESH, hecMAT, X, Y, COMMtime)
Definition: hecmw_jadm_44.f90:57
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
hecmw_jad_type_44::hecmw_jad_finalize_44
subroutine, public hecmw_jad_finalize_44()
Definition: hecmw_jadm_44.f90:43
m_hecmw_comm_f::hecmw_update_r
subroutine hecmw_update_r(hecMESH, val, n, m)
Definition: hecmw_comm_f.F90:683
hecmw_util::hecmwst_matrix
Definition: hecmw_util_f.F90:444