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) :: 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_33_inner(hecmesh, hecmat, x, y, time_ax, tcomm)
58 if (
present(commtime)) commtime = commtime + tcomm
120 subroutine hecmw_matvec_33_inner (hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
132 real(kind=
kreal),
intent(in) :: x(:)
133 real(kind=
kreal),
intent(out) :: y(:)
134 real(kind=
kreal),
intent(inout) :: time_ax
135 real(kind=
kreal),
intent(inout),
optional :: commtime
137 real(kind=
kreal) :: start_time, end_time, tcomm
138 integer(kind=kint) :: i, j, js, je, in
139 real(kind=
kreal) :: yv1, yv2, yv3, x1, x2, x3
141 integer(kind=kint) :: n, np
142 integer(kind=kint),
pointer :: indexl(:), iteml(:), indexu(:), itemu(:)
143 real(kind=
kreal),
pointer :: al(:), au(:), d(:)
146 integer,
parameter :: numofblockperthread = 100
147 logical,
save :: isfirst = .true.
148 integer,
save :: numofthread = 1
149 integer,
save,
allocatable :: startpos(:), endpos(:)
150 integer(kind=kint),
save :: sectorcachesize0, sectorcachesize1
151 integer(kind=kint) :: threadnum, blocknum, numofblock
152 integer(kind=kint) :: numofelement, elementcount, blockindex
153 real(kind=
kreal) :: numofelementperblock
161 time_ax = time_ax + end_time - start_time - tcomm
162 if (
present(commtime)) commtime = commtime + tcomm
167 indexl => hecmat%indexL
168 indexu => hecmat%indexU
169 iteml => hecmat%itemL
170 itemu => hecmat%itemU
176 if (.not. isfirst)
then
177 numofblock = numofthread * numofblockperthread
178 if (endpos(numofblock-1) .ne. n-1)
then
179 deallocate(startpos, endpos)
185 numofblock = numofthread * numofblockperthread
186 allocate (startpos(0 : numofblock - 1), endpos(0 : numofblock - 1))
187 numofelement = n + indexl(n) + indexu(n)
188 numofelementperblock = dble(numofelement) / numofblock
191 startpos(blocknum) = 1
193 elementcount = elementcount + 1
194 elementcount = elementcount + (indexl(i) - indexl(i-1))
195 elementcount = elementcount + (indexu(i) - indexu(i-1))
196 if (elementcount > (blocknum + 1) * numofelementperblock)
then
200 blocknum = blocknum + 1
201 startpos(blocknum) = i + 1
202 if (blocknum == (numofblock - 1))
exit
209 do i= blocknum+1, numofblock-1
217 sectorcachesize0, sectorcachesize1)
230 if (
present(commtime)) commtime = commtime + end_time - start_time
245 do blocknum = 0 , numofblockperthread - 1
246 blockindex = blocknum * numofthread + threadnum
247 do i = startpos(blockindex), endpos(blockindex)
251 yv1= d(9*i-8)*x1 + d(9*i-7)*x2 + d(9*i-6)*x3
252 yv2= d(9*i-5)*x1 + d(9*i-4)*x2 + d(9*i-3)*x3
253 yv3= d(9*i-2)*x1 + d(9*i-1)*x2 + d(9*i )*x3
262 yv1= yv1 + al(9*j-8)*x1 + al(9*j-7)*x2 + al(9*j-6)*x3
263 yv2= yv2 + al(9*j-5)*x1 + al(9*j-4)*x2 + al(9*j-3)*x3
264 yv3= yv3 + al(9*j-2)*x1 + al(9*j-1)*x2 + al(9*j )*x3
274 yv1= yv1 + au(9*j-8)*x1 + au(9*j-7)*x2 + au(9*j-6)*x3
275 yv2= yv2 + au(9*j-5)*x1 + au(9*j-4)*x2 + au(9*j-3)*x3
276 yv3= yv3 + au(9*j-2)*x1 + au(9*j-1)*x2 + au(9*j )*x3
292 time_ax = time_ax + end_time - start_time
329 if (hecmat%cmat%n_val > 0)
then
333 end subroutine hecmw_matvec_33_inner
345 real(kind=
kreal),
intent(in) :: x(:), b(:)
346 real(kind=
kreal),
intent(out) :: r(:)
347 real(kind=
kreal),
intent(inout),
optional :: commtime
349 integer(kind=kint) :: i
350 real(kind=
kreal) :: tcomm
354 if (
present(commtime)) commtime = commtime + tcomm
357 do i = 1, hecmat%N * 3
376 real(kind=
kreal),
intent(inout),
optional :: commtime
378 real(kind=
kreal),
allocatable :: r(:)
379 real(kind=
kreal) :: bnorm2, rnorm2
380 real(kind=
kreal) :: tcomm
382 allocate(r(hecmat%NDOF*hecmat%NP))
386 hecmat%B, hecmat%B, bnorm2, tcomm)
387 if (bnorm2 == 0.d0)
then
394 if (
present(commtime)) commtime = commtime + tcomm
409 real(kind=
kreal),
intent(in) :: x(:)
410 real(kind=
kreal),
intent(out) :: y(:)
411 real(kind=
kreal),
intent(inout) :: commtime
413 real(kind=
kreal) :: start_time, end_time
414 integer(kind=kint) :: i, j, jj, k, kk
419 commtime = commtime + end_time - start_time
423 do i= 1, hecmesh%nn_internal * hecmesh%n_dof
429 outer:
do i= 1, hecmesh%mpc%n_mpc
430 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
431 if (hecmesh%mpc%mpc_dof(j) > 3) cycle outer
433 k = hecmesh%mpc%mpc_index(i-1) + 1
434 kk = 3 * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
436 do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
437 jj = 3 * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
438 y(kk) = y(kk) - hecmesh%mpc%mpc_val(j) * x(jj)
456 real(kind=
kreal),
intent(in) :: x(:)
457 real(kind=
kreal),
intent(out) :: y(:)
458 real(kind=
kreal),
intent(inout) :: commtime
460 real(kind=
kreal) :: start_time, end_time
461 integer(kind=kint) :: i, j, jj, k, kk
466 commtime = commtime + end_time - start_time
470 do i= 1, hecmesh%nn_internal * hecmesh%n_dof
476 outer:
do i= 1, hecmesh%mpc%n_mpc
477 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
478 if (hecmesh%mpc%mpc_dof(j) > 3) cycle outer
480 k = hecmesh%mpc%mpc_index(i-1) + 1
481 kk = 3 * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
483 do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
484 jj = 3 * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
485 y(jj) = y(jj) - hecmesh%mpc%mpc_val(j) * x(kk)
503 real(kind=
kreal),
intent(in) :: x(:)
504 real(kind=
kreal),
intent(out) :: y(:), w(:)
505 real(kind=
kreal),
intent(inout) :: commtime
508 call hecmw_matvec_33_inner(hecmesh, hecmat, y, w, commtime)
525 real(kind=
kreal),
intent(inout),
optional :: commtime
526 real(kind=
kreal),
allocatable :: w(:,:)
527 real(kind=
kreal),
pointer :: d(:)
528 integer(kind=kint) :: ip
529 real(kind=
kreal) :: start_time, end_time
530 allocate(w(3*hecmat%NP,3))
533 w(3*ip-2,1)= d(9*ip-8); w(3*ip-2,2)= d(9*ip-7); w(3*ip-2,3)= d(9*ip-6)
534 w(3*ip-1,1)= d(9*ip-5); w(3*ip-1,2)= d(9*ip-4); w(3*ip-1,3)= d(9*ip-3)
535 w(3*ip ,1)= d(9*ip-2); w(3*ip ,2)= d(9*ip-1); w(3*ip ,3)= d(9*ip )
542 if (
present(commtime)) commtime = commtime + end_time - start_time
543 do ip= hecmat%N+1, hecmat%NP
544 d(9*ip-8)= w(3*ip-2,1); d(9*ip-7)= w(3*ip-2,2); d(9*ip-6)= w(3*ip-2,3)
545 d(9*ip-5)= w(3*ip-1,1); d(9*ip-4)= w(3*ip-1,2); d(9*ip-3)= w(3*ip-1,3)
546 d(9*ip-2)= w(3*ip ,1); d(9*ip-1)= w(3*ip ,2); d(9*ip )= w(3*ip ,3)
Jagged Diagonal Matrix storage for vector processors. Original code was provided by JAMSTEC.
subroutine, public hecmw_jad_matvec(hecMESH, hecMAT, X, Y, COMMtime)
integer(kind=kint) function, public hecmw_jad_is_initialized()
integer(kind=kint) function, public hecmw_mat_get_flag_mpcmatvec(hecMAT)
subroutine, public hecmw_matvec_33(hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
subroutine, public hecmw_mat_diag_sr_33(hecMESH, hecMAT, COMMtime)
subroutine, public hecmw_matvec_33_unset_async
subroutine, public hecmw_matresid_33(hecMESH, hecMAT, X, B, R, COMMtime)
real(kind=kreal) function, public hecmw_rel_resid_l2_33(hecMESH, hecMAT, COMMtime)
subroutine, public hecmw_ttmattvec_33(hecMESH, hecMAT, X, Y, W, COMMtime)
subroutine, public hecmw_ttvec_33(hecMESH, X, Y, COMMtime)
subroutine, public hecmw_tvec_33(hecMESH, X, Y, COMMtime)
subroutine, public hecmw_matvec_33_set_async(hecMAT)
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_3_r(hecMESH, val, n)