25 logical,
save :: async_matvec_flg = .false.
40 real(kind=
kreal),
intent(in) :: x(:)
41 real(kind=
kreal),
intent(out) :: y(:)
42 real(kind=
kreal),
intent(inout),
optional :: time_ax
43 real(kind=
kreal),
intent(inout),
optional :: commtime
45 real(kind=
kreal) :: tcomm
46 real(kind=
kreal),
allocatable :: wk(:)
51 allocate(wk(hecmat%NP * hecmat%NDOF))
55 call hecmw_matvec_44_inner(hecmesh, hecmat, x, y, time_ax, tcomm)
58 if (
present(commtime)) commtime = commtime + tcomm
88 subroutine hecmw_matvec_44_inner (hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
100 real(kind=
kreal),
intent(in) :: x(:)
101 real(kind=
kreal),
intent(out) :: y(:)
102 real(kind=
kreal),
intent(inout) :: time_ax
103 real(kind=
kreal),
intent(inout),
optional :: commtime
105 real(kind=
kreal) :: start_time, end_time, tcomm
106 integer(kind=kint) :: i, j, js, je, in
107 real(kind=
kreal) :: yv1, yv2, yv3, yv4, x1, x2, x3, x4
109 integer(kind=kint) :: n, np
110 integer(kind=kint),
pointer :: indexl(:), iteml(:), indexu(:), itemu(:)
111 real(kind=
kreal),
pointer :: al(:), au(:), d(:)
114 integer,
parameter :: numofblockperthread = 100
115 logical,
save :: isfirst = .true.
116 integer,
save :: numofthread = 1
117 integer,
save,
allocatable :: startpos(:), endpos(:)
118 integer(kind=kint),
save :: sectorcachesize0, sectorcachesize1
119 integer(kind=kint) :: threadnum, blocknum, numofblock
120 integer(kind=kint) :: numofelement, elementcount, blockindex
121 real(kind=
kreal) :: numofelementperblock
129 time_ax = time_ax + end_time - start_time - tcomm
130 if (
present(commtime)) commtime = commtime + tcomm
135 indexl => hecmat%indexL
136 indexu => hecmat%indexU
137 iteml => hecmat%itemL
138 itemu => hecmat%itemU
144 if (.not. isfirst)
then
145 numofblock = numofthread * numofblockperthread
146 if (endpos(numofblock-1) .ne. n-1)
then
147 deallocate(startpos, endpos)
153 numofblock = numofthread * numofblockperthread
154 allocate (startpos(0 : numofblock - 1), endpos(0 : numofblock - 1))
155 numofelement = n + indexl(n) + indexu(n)
156 numofelementperblock = dble(numofelement) / numofblock
159 startpos(blocknum) = 1
161 elementcount = elementcount + 1
162 elementcount = elementcount + (indexl(i) - indexl(i-1))
163 elementcount = elementcount + (indexu(i) - indexu(i-1))
164 if (elementcount > (blocknum + 1) * numofelementperblock)
then
168 blocknum = blocknum + 1
169 startpos(blocknum) = i + 1
170 if (blocknum == (numofblock - 1))
exit
177 do i= blocknum+1, numofblock-1
185 sectorcachesize0, sectorcachesize1)
197 if (
present(commtime)) commtime = commtime + end_time - start_time
212 do blocknum = 0 , numofblockperthread - 1
213 blockindex = blocknum * numofthread + threadnum
214 do i = startpos(blockindex), endpos(blockindex)
219 yv1= d(16*i-15)*x1 + d(16*i-14)*x2 + d(16*i-13)*x3 + d(16*i-12)*x4
220 yv2= d(16*i-11)*x1 + d(16*i-10)*x2 + d(16*i- 9)*x3 + d(16*i- 8)*x4
221 yv3= d(16*i- 7)*x1 + d(16*i- 6)*x2 + d(16*i- 5)*x3 + d(16*i- 4)*x4
222 yv4= d(16*i- 3)*x1 + d(16*i- 2)*x2 + d(16*i- 1)*x3 + d(16*i )*x4
232 yv1= yv1 + al(16*j-15)*x1 + al(16*j-14)*x2 + al(16*j-13)*x3 + al(16*j-12)*x4
233 yv2= yv2 + al(16*j-11)*x1 + al(16*j-10)*x2 + al(16*j- 9)*x3 + al(16*j- 8)*x4
234 yv3= yv3 + al(16*j- 7)*x1 + al(16*j- 6)*x2 + al(16*j- 5)*x3 + al(16*j- 4)*x4
235 yv4= yv4 + al(16*j- 3)*x1 + al(16*j- 2)*x2 + al(16*j- 1)*x3 + al(16*j )*x4
246 yv1= yv1 + au(16*j-15)*x1 + au(16*j-14)*x2 + au(16*j-13)*x3 + au(16*j-12)*x4
247 yv2= yv2 + au(16*j-11)*x1 + au(16*j-10)*x2 + au(16*j- 9)*x3 + au(16*j- 8)*x4
248 yv3= yv3 + au(16*j- 7)*x1 + au(16*j- 6)*x2 + au(16*j- 5)*x3 + au(16*j- 4)*x4
249 yv4= yv4 + au(16*j- 3)*x1 + au(16*j- 2)*x2 + au(16*j- 1)*x3 + au(16*j )*x4
266 time_ax = time_ax + end_time - start_time
270 if (hecmat%cmat%n_val > 0)
then
274 end subroutine hecmw_matvec_44_inner
286 real(kind=
kreal),
intent(in) :: x(:), b(:)
287 real(kind=
kreal),
intent(out) :: r(:)
288 real(kind=
kreal),
intent(inout),
optional :: commtime
290 integer(kind=kint) :: i
291 real(kind=
kreal) :: tcomm
295 if (
present(commtime)) commtime = commtime + tcomm
298 do i = 1, hecmat%N * 4
317 real(kind=
kreal),
intent(inout),
optional :: commtime
319 real(kind=
kreal),
allocatable :: r(:)
320 real(kind=
kreal) :: bnorm2, rnorm2
321 real(kind=
kreal) :: tcomm
323 allocate(r(hecmat%NDOF*hecmat%NP))
327 hecmat%B, hecmat%B, bnorm2, tcomm)
328 if (bnorm2 == 0.d0)
then
335 if (
present(commtime)) commtime = commtime + tcomm
350 real(kind=
kreal),
intent(in) :: x(:)
351 real(kind=
kreal),
intent(out) :: y(:)
352 real(kind=
kreal),
intent(inout) :: commtime
354 real(kind=
kreal) :: start_time, end_time
355 integer(kind=kint) :: i, j, jj, k, kk
360 commtime = commtime + end_time - start_time
364 do i= 1, hecmesh%nn_internal * hecmesh%n_dof
370 outer:
do i= 1, hecmesh%mpc%n_mpc
371 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
372 if (hecmesh%mpc%mpc_dof(j) > 4) cycle outer
374 k = hecmesh%mpc%mpc_index(i-1) + 1
375 kk = 4 * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
377 do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
378 jj = 4 * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
379 y(kk) = y(kk) - hecmesh%mpc%mpc_val(j) * x(jj)
397 real(kind=
kreal),
intent(in) :: x(:)
398 real(kind=
kreal),
intent(out) :: y(:)
399 real(kind=
kreal),
intent(inout) :: commtime
401 real(kind=
kreal) :: start_time, end_time
402 integer(kind=kint) :: i, j, jj, k, kk
407 commtime = commtime + end_time - start_time
411 do i= 1, hecmesh%nn_internal * hecmesh%n_dof
417 outer:
do i= 1, hecmesh%mpc%n_mpc
418 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
419 if (hecmesh%mpc%mpc_dof(j) > 4) cycle outer
421 k = hecmesh%mpc%mpc_index(i-1) + 1
422 kk = 4 * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
424 do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
425 jj = 4 * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
426 y(jj) = y(jj) - hecmesh%mpc%mpc_val(j) * x(kk)
444 real(kind=
kreal),
intent(in) :: x(:)
445 real(kind=
kreal),
intent(out) :: y(:), w(:)
446 real(kind=
kreal),
intent(inout) :: commtime
449 call hecmw_matvec_44_inner(hecmesh, hecmat, y, w, commtime)
466 real(kind=
kreal),
intent(inout),
optional :: commtime
467 real(kind=
kreal),
allocatable :: w(:,:)
468 real(kind=
kreal),
pointer :: d(:)
469 integer(kind=kint) :: ip
470 real(kind=
kreal) :: start_time, end_time
471 allocate(w(4*hecmat%NP,4))
474 w(4*ip-3,1)= d(16*ip-15); w(4*ip-3,2)= d(16*ip-14); w(4*ip-3,3)= d(16*ip-13); w(4*ip-3,4)= d(16*ip-12)
475 w(4*ip-2,1)= d(16*ip-11); w(4*ip-2,2)= d(16*ip-10); w(4*ip-2,3)= d(16*ip- 9); w(4*ip-2,4)= d(16*ip- 8)
476 w(4*ip-1,1)= d(16*ip- 7); w(4*ip-1,2)= d(16*ip- 6); w(4*ip-1,3)= d(16*ip- 5); w(4*ip-1,4)= d(16*ip- 4)
477 w(4*ip ,1)= d(16*ip- 3); w(4*ip ,2)= d(16*ip- 2); w(4*ip ,3)= d(16*ip- 1); w(4*ip ,4)= d(16*ip )
485 if (
present(commtime)) commtime = commtime + end_time - start_time
486 do ip= hecmat%N+1, hecmat%NP
487 d(16*ip-15)= w(4*ip-3,1); d(16*ip-14)= w(4*ip-3,2); d(16*ip-13)= w(4*ip-3,3); d(16*ip-12)= w(4*ip-3,4)
488 d(16*ip-11)= w(4*ip-2,1); d(16*ip-10)= w(4*ip-2,2); d(16*ip- 9)= w(4*ip-2,3); d(16*ip- 8)= w(4*ip-2,4)
489 d(16*ip- 7)= w(4*ip-1,1); d(16*ip- 6)= w(4*ip-1,2); d(16*ip- 5)= w(4*ip-1,3); d(16*ip- 4)= w(4*ip-1,4)
490 d(16*ip- 3)= w(4*ip ,1); d(16*ip- 2)= w(4*ip ,2); d(16*ip- 1)= w(4*ip ,3); d(16*ip )= w(4*ip ,4)
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)
subroutine, public hecmw_tvec_44(hecmesh, x, y, commtime)
real(kind=kreal) function, public hecmw_rel_resid_l2_44(hecmesh, hecmat, commtime)
subroutine, public hecmw_mat_diag_sr_44(hecmesh, hecmat, commtime)
subroutine, public hecmw_matvec_44_unset_async
subroutine, public hecmw_ttmattvec_44(hecmesh, hecmat, x, y, w, commtime)
subroutine, public hecmw_matresid_44(hecmesh, hecmat, x, b, r, commtime)
subroutine, public hecmw_matvec_44_set_async(hecmat)
subroutine, public hecmw_matvec_44(hecmesh, hecmat, x, y, time_ax, commtime)
subroutine, public hecmw_ttvec_44(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_4_r(hecmesh, val, n)