27 logical,
save :: async_matvec_flg = .false.
42 real(kind=
kreal),
intent(in) :: x(:)
43 real(kind=
kreal),
intent(out) :: y(:)
44 real(kind=
kreal),
intent(inout) :: time_ax
45 real(kind=
kreal),
intent(inout),
optional :: commtime
47 real(kind=
kreal) :: tcomm
48 real(kind=
kreal),
allocatable :: wk(:)
53 allocate(wk(hecmat%NP * hecmat%NDOF))
57 call hecmw_matvec_nn_inner(hecmesh, hecmat, x, y, tcomm)
60 if (
present(commtime)) commtime = commtime + tcomm
122 subroutine hecmw_matvec_nn_inner (hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
134 real(kind=
kreal),
intent(in) :: x(:)
135 real(kind=
kreal),
intent(out) :: y(:)
136 real(kind=
kreal),
intent(inout) :: time_ax
137 real(kind=
kreal),
intent(inout),
optional :: commtime
139 real(kind=
kreal) :: start_time, end_time, tcomm
140 integer(kind=kint) :: i, j, k, l, js, je, in
141 real(kind=
kreal) :: yv(hecmat%NDOF), xv(hecmat%NDOF)
143 integer(kind=kint) :: n, np, ndof, ndof2
144 integer(kind=kint),
pointer :: indexl(:), iteml(:), indexu(:), itemu(:)
145 real(kind=
kreal),
pointer :: al(:), au(:), d(:)
148 integer,
parameter :: numofblockperthread = 100
149 logical,
save :: isfirst = .true.
150 integer,
save :: numofthread = 1
151 integer,
save,
allocatable :: startpos(:), endpos(:)
152 integer(kind=kint),
save :: sectorcachesize0, sectorcachesize1
153 integer(kind=kint) :: threadnum, blocknum, numofblock
154 integer(kind=kint) :: numofelement, elementcount, blockindex
155 real(kind=
kreal) :: numofelementperblock
163 time_ax = time_ax + end_time - start_time - tcomm
164 if (
present(commtime)) commtime = commtime + tcomm
169 indexl => hecmat%indexL
170 indexu => hecmat%indexU
171 iteml => hecmat%itemL
172 itemu => hecmat%itemU
180 if (.not. isfirst)
then
181 numofblock = numofthread * numofblockperthread
182 if (endpos(numofblock-1) .ne. n-1)
then
183 deallocate(startpos, endpos)
189 numofblock = numofthread * numofblockperthread
190 allocate (startpos(0 : numofblock - 1), endpos(0 : numofblock - 1))
191 numofelement = n + indexl(n) + indexu(n)
192 numofelementperblock = dble(numofelement) / numofblock
195 startpos(blocknum) = 1
197 elementcount = elementcount + 1
198 elementcount = elementcount + (indexl(i) - indexl(i-1))
199 elementcount = elementcount + (indexu(i) - indexu(i-1))
200 if (elementcount > (blocknum + 1) * numofelementperblock)
then
204 blocknum = blocknum + 1
205 startpos(blocknum) = i + 1
206 if (blocknum == (numofblock - 1))
exit
213 do i= blocknum+1, numofblock-1
221 sectorcachesize0, sectorcachesize1)
234 if (
present(commtime)) commtime = commtime + end_time - start_time
249 do blocknum = 0 , numofblockperthread - 1
250 blockindex = blocknum * numofthread + threadnum
251 do i = startpos(blockindex), endpos(blockindex)
253 xv(k) = x(ndof*(i-1)+k)
258 yv(k)=yv(k)+d(ndof2*(i-1)+(k-1)*ndof+l)*xv(l)
266 xv(k) = x(ndof*(in-1)+k)
270 yv(k)=yv(k)+al(ndof2*(j-1)+(k-1)*ndof+l)*xv(l)
280 xv(k) = x(ndof*(in-1)+k)
284 yv(k)=yv(k)+au(ndof2*(j-1)+(k-1)*ndof+l)*xv(l)
289 y(ndof*(i-1)+k) = yv(k)
302 time_ax = time_ax + end_time - start_time
308 if (hecmat%cmat%n_val > 0)
then
311 end subroutine hecmw_matvec_nn_inner
325 real(kind=
kreal),
intent(in) :: x(:), b(:)
326 real(kind=
kreal),
intent(out) :: r(:)
327 real(kind=
kreal),
intent(inout),
optional :: commtime
329 integer(kind=kint) :: i
330 real(kind=
kreal) :: tcomm
334 if (
present(commtime)) commtime = commtime + tcomm
337 do i = 1, hecmat%N * hecmat%NDOF
356 real(kind=
kreal),
intent(inout),
optional :: commtime
358 real(kind=
kreal),
allocatable :: r(:)
359 real(kind=
kreal) :: bnorm2, rnorm2
360 real(kind=
kreal) :: tcomm
362 allocate(r(hecmat%NDOF*hecmat%NP))
366 hecmat%B, hecmat%B, bnorm2, tcomm)
367 if (bnorm2 == 0.d0)
then
374 if (
present(commtime)) commtime = commtime + tcomm
389 integer(kind=kint),
intent(in) :: ndof
390 real(kind=
kreal),
intent(in) :: x(:)
391 real(kind=
kreal),
intent(out) :: y(:)
392 real(kind=
kreal),
intent(inout) :: commtime
394 real(kind=
kreal) :: start_time, end_time
395 integer(kind=kint) :: i, j, jj, k, kk
400 commtime = commtime + end_time - start_time
404 do i= 1, hecmesh%nn_internal * ndof
410 outer:
do i= 1, hecmesh%mpc%n_mpc
411 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
412 if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
414 k = hecmesh%mpc%mpc_index(i-1) + 1
415 kk = ndof * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
417 do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
418 jj = ndof * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
419 y(kk) = y(kk) - hecmesh%mpc%mpc_val(j) * x(jj)
437 integer(kind=kint),
intent(in) :: ndof
438 real(kind=
kreal),
intent(in) :: x(:)
439 real(kind=
kreal),
intent(out) :: y(:)
440 real(kind=
kreal),
intent(inout) :: commtime
442 real(kind=
kreal) :: start_time, end_time
443 integer(kind=kint) :: i, j, jj, k, kk
448 commtime = commtime + end_time - start_time
452 do i= 1, hecmesh%nn_internal * ndof
458 outer:
do i= 1, hecmesh%mpc%n_mpc
459 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
460 if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
462 k = hecmesh%mpc%mpc_index(i-1) + 1
463 kk = ndof * (hecmesh%mpc%mpc_item(k) - 1) + hecmesh%mpc%mpc_dof(k)
465 do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
466 jj = ndof * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
467 y(jj) = y(jj) - hecmesh%mpc%mpc_val(j) * x(kk)
486 real(kind=
kreal),
intent(in) :: x(:)
487 real(kind=
kreal),
intent(out) :: y(:), w(:)
488 real(kind=
kreal),
intent(inout) :: commtime
491 call hecmw_matvec_nn_inner(hecmesh, hecmat, y, w, commtime)
507 real(kind=
kreal),
intent(inout),
optional :: commtime
508 real(kind=
kreal),
allocatable :: w(:,:)
509 real(kind=
kreal),
pointer :: d(:)
510 integer(kind=kint) :: ip, ndof, i, j
511 real(kind=
kreal) :: start_time, end_time
513 allocate(w(ndof*hecmat%NP,ndof))
518 w(ndof*(ip-1)+i,j) = d(ndof*ndof*(ip-1)+(i-1)*ndof+j)
527 if (
present(commtime)) commtime = commtime + end_time - start_time
528 do ip= hecmat%N+1, hecmat%NP
531 d(ndof*ndof*(ip-1)+(i-1)*ndof+j) = w(ndof*(ip-1)+i,j)
542 integer(kind=kint) :: i
544 do i = 1, hecmat1%NP*hecmat1%NDOF*hecmat1%NDOF
545 hecmat3%D(i) = hecmat1%D(i) + hecmat2%D(i)
548 do i = 1, hecmat1%NPU*hecmat1%NDOF*hecmat1%NDOF
549 hecmat3%AU(i) = hecmat1%AU(i) + hecmat2%AU(i)
552 do i = 1, hecmat1%NPL*hecmat1%NDOF*hecmat1%NDOF
553 hecmat3%AL(i) = hecmat1%AL(i) + hecmat2%AL(i)
561 real(kind=
kreal),
intent(in) :: alpha
562 integer(kind=kint) :: i
564 do i = 1, hecmat%NP*hecmat%NDOF*hecmat%NDOF
565 hecmat%D(i) = alpha*hecmat%D(i)
568 do i = 1, hecmat%NPU*hecmat%NDOF*hecmat%NDOF
569 hecmat%AU(i) = alpha*hecmat%AU(i)
572 do i = 1, hecmat%NPL*hecmat%NDOF*hecmat%NDOF
573 hecmat%AL(i) = alpha*hecmat%AL(i)
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_mat_diag_sr_nn(hecMESH, hecMAT, COMMtime)
subroutine, public hecmw_mat_multiple_nn(hecMAT, alpha)
subroutine, public hecmw_ttmattvec_nn(hecMESH, hecMAT, X, Y, W, COMMtime)
subroutine, public hecmw_matresid_nn(hecMESH, hecMAT, X, B, R, COMMtime)
subroutine, public hecmw_ttvec_nn(hecMESH, ndof, X, Y, COMMtime)
subroutine, public hecmw_matvec_nn_unset_async
real(kind=kreal) function, public hecmw_rel_resid_l2_nn(hecMESH, hecMAT, COMMtime)
subroutine, public hecmw_tvec_nn(hecMESH, ndof, X, Y, COMMtime)
subroutine, public hecmw_matvec_nn_set_async(hecMAT)
subroutine, public hecmw_matvec_nn(hecMESH, hecMAT, X, Y, time_Ax, COMMtime)
subroutine, public hecmw_mat_add_nn(hecMAT1, hecMAT2, hecMAT3)
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_m_r(hecMESH, val, n, m)