23 logical,
save :: async_matvec_flg = .false.
37 real(kind=
kreal),
intent(in) :: x(:)
38 real(kind=
kreal),
intent(out) :: y(:)
39 real(kind=
kreal),
intent(inout) :: time_ax
40 real(kind=
kreal),
intent(inout),
optional :: commtime
42 real(kind=
kreal) :: tcomm
43 real(kind=
kreal),
allocatable :: wk(:)
48 allocate(wk(hecmat%NP * hecmat%NDOF))
52 call hecmw_matvec_22_inner(hecmesh, hecmat, x, y, time_ax, tcomm)
55 if (
present(commtime)) commtime = commtime + tcomm
72 subroutine hecmw_matvec_22_inner (hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
84 real(kind=
kreal),
intent(in) :: x(:)
85 real(kind=
kreal),
intent(out) :: y(:)
86 real(kind=
kreal),
intent(inout) :: time_ax
87 real(kind=
kreal),
intent(inout),
optional :: commtime
89 real(kind=
kreal) :: start_time, end_time, tcomm
90 integer(kind=kint) :: i, j, js, je, in
91 real(kind=
kreal) :: yv1, yv2, x1, x2
93 integer(kind=kint) :: n, np
94 integer(kind=kint),
pointer :: indexl(:), iteml(:), indexu(:), itemu(:)
95 real(kind=
kreal),
pointer :: al(:), au(:), d(:)
98 integer,
parameter :: numofblockperthread = 100
99 logical,
save :: isfirst = .true.
100 integer,
save :: numofthread = 1
101 integer,
save,
allocatable :: startpos(:), endpos(:)
102 integer(kind=kint),
save :: sectorcachesize0, sectorcachesize1
103 integer(kind=kint) :: threadnum, blocknum, numofblock
104 integer(kind=kint) :: numofelement, elementcount, blockindex
105 real(kind=
kreal) :: numofelementperblock
113 time_ax = time_ax + end_time - start_time - tcomm
114 if (
present(commtime)) commtime = commtime + tcomm
119 indexl => hecmat%indexL
120 indexu => hecmat%indexU
121 iteml => hecmat%itemL
122 itemu => hecmat%itemU
128 if (.not. isfirst)
then
129 numofblock = numofthread * numofblockperthread
130 if (endpos(numofblock-1) .ne. n-1)
then
131 deallocate(startpos, endpos)
137 numofblock = numofthread * numofblockperthread
138 allocate (startpos(0 : numofblock - 1), endpos(0 : numofblock - 1))
139 numofelement = n + indexl(n) + indexu(n)
140 numofelementperblock = dble(numofelement) / numofblock
143 startpos(blocknum) = 1
145 elementcount = elementcount + 1
146 elementcount = elementcount + (indexl(i) - indexl(i-1))
147 elementcount = elementcount + (indexu(i) - indexu(i-1))
148 if (elementcount > (blocknum + 1) * numofelementperblock)
then
150 blocknum = blocknum + 1
151 startpos(blocknum) = i + 1
152 if (blocknum == (numofblock - 1))
exit
157 do i= blocknum+1, numofblock-1
163 sectorcachesize0, sectorcachesize1)
172 if (
present(commtime)) commtime = commtime + end_time - start_time
187 do blocknum = 0 , numofblockperthread - 1
188 blockindex = blocknum * numofthread + threadnum
189 do i = startpos(blockindex), endpos(blockindex)
192 yv1= d(4*i-3)*x1 + d(4*i-2)*x2
193 yv2= d(4*i-1)*x1 + d(4*i )*x2
201 yv1= yv1 + al(4*j-3)*x1 + al(4*j-2)*x2
202 yv2= yv2 + al(4*j-1)*x1 + al(4*j )*x2
210 yv1= yv1 + au(4*j-3)*x1 + au(4*j-2)*x2
211 yv2= yv2 + au(4*j-1)*x1 + au(4*j )*x2
225 time_ax = time_ax + end_time - start_time
230 if (hecmat%cmat%n_val > 0)
then
234 end subroutine hecmw_matvec_22_inner
250 real(kind=
kreal) :: x(:), b(:), r(:)
253 real(kind=
kreal),
optional :: commtime
255 integer(kind=kint) :: i
256 real(kind=
kreal) :: tcomm
260 if (
present(commtime)) commtime = commtime + tcomm
263 do i = 1, hecmat%N * 2
283 real(kind=
kreal),
optional :: commtime
285 real(kind=
kreal) :: r(hecmat%NDOF*hecmat%NP)
286 real(kind=
kreal) :: bnorm2, rnorm2
287 real(kind=
kreal) :: tcomm
291 if (bnorm2 == 0.d0)
then
296 if (
present(commtime)) commtime = commtime + tcomm
311 real(kind=
kreal),
intent(in) :: x(:)
312 real(kind=
kreal),
intent(out) :: y(:)
313 real(kind=
kreal),
intent(inout) :: commtime
315 real(kind=
kreal) :: start_time, end_time
316 integer(kind=kint) :: i, j, jj, k, kk
321 commtime = commtime + end_time - start_time
325 do i= 1, hecmesh%nn_internal * hecmesh%n_dof
331 outer:
do i= 1, hecmesh%mpc%n_mpc
332 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
333 if (hecmesh%mpc%mpc_dof(j) > 2) cycle outer
335 k = hecmesh%mpc%mpc_index(i-1) + 1
336 kk = 2 * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
338 do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
339 jj = 2 * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
340 y(kk) = y(kk) - hecmesh%mpc%mpc_val(j) * x(jj)
358 real(kind=
kreal),
intent(in) :: x(:)
359 real(kind=
kreal),
intent(out) :: y(:)
360 real(kind=
kreal),
intent(inout) :: commtime
362 real(kind=
kreal) :: start_time, end_time
363 integer(kind=kint) :: i, j, jj, k, kk
368 commtime = commtime + end_time - start_time
372 do i= 1, hecmesh%nn_internal * hecmesh%n_dof
378 outer:
do i= 1, hecmesh%mpc%n_mpc
379 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
380 if (hecmesh%mpc%mpc_dof(j) > 2) cycle outer
382 k = hecmesh%mpc%mpc_index(i-1) + 1
383 kk = 2 * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
385 do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
386 jj = 2 * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
387 y(jj) = y(jj) - hecmesh%mpc%mpc_val(j) * x(kk)
405 real(kind=
kreal),
intent(in) :: x(:)
406 real(kind=
kreal),
intent(out) :: y(:), w(:)
407 real(kind=
kreal),
intent(inout) :: commtime
410 call hecmw_matvec_22_inner(hecmesh, hecmat, y, w, commtime)
427 real(kind=
kreal),
intent(inout),
optional :: commtime
428 real(kind=
kreal),
allocatable :: w(:,:)
429 real(kind=
kreal),
pointer :: d(:)
430 integer(kind=kint) :: ip
431 real(kind=
kreal) :: start_time, end_time
432 allocate(w(2*hecmat%NP,2))
435 w(2*ip-1,1)= d(4*ip-3); w(2*ip-1,2)= d(4*ip-2);
436 w(2*ip-0,1)= d(4*ip-1); w(2*ip-0,2)= d(4*ip-0);
442 if (
present(commtime)) commtime = commtime + end_time - start_time
443 do ip= hecmat%N+1, hecmat%NP
444 d(4*ip-3)= w(2*ip-1,1); d(4*ip-2)= w(2*ip-1,2);
445 d(4*ip-1)= w(2*ip-0,1); d(4*ip-0)= w(2*ip-0,2);
Jagged Diagonal Matrix storage for vector processors. Original code was provided by JAMSTEC.
integer(kind=kint) function, public hecmw_jad_is_initialized()
subroutine, public hecmw_jad_matvec(hecmesh, hecmat, x, y, commtime)
integer(kind=kint) function, public hecmw_mat_get_flag_mpcmatvec(hecmat)
real(kind=kreal) function, public hecmw_rel_resid_l2_22(hecmesh, hecmat, commtime)
subroutine, public hecmw_matvec_22(hecmesh, hecmat, x, y, time_ax, commtime)
subroutine, public hecmw_ttvec_22(hecmesh, x, y, commtime)
subroutine, public hecmw_ttmattvec_22(hecmesh, hecmat, x, y, w, commtime)
subroutine, public hecmw_matresid_22(hecmesh, hecmat, x, b, r, commtime)
subroutine, public hecmw_matvec_22_set_async(hecmat)
subroutine, public hecmw_mat_diag_sr_22(hecmesh, hecmat, commtime)
subroutine, public hecmw_matvec_22_unset_async
subroutine, public hecmw_tvec_22(hecmesh, x, y, commtime)
subroutine hecmw_innerproduct_r(hecmesh, ndof, x, y, sum, commtime)
subroutine, public hecmw_tuning_fx_calc_sector_cache(n, ndof, sectorcachesize0, sectorcachesize1)
integer(kind=4), parameter kreal
real(kind=kreal) function hecmw_wtime()
subroutine hecmw_update_2_r(hecmesh, val, n)
subroutine hecmw_update_3_r(hecmesh, val, n)