28 integer :: nr, nc, nnz, ndof
29 integer(kind=kint),
pointer :: index(:)
30 integer(kind=kint),
pointer :: item(:)
31 real(kind=
kreal),
pointer :: a(:)
35 integer(kind=kint),
parameter :: clid = 1
36 integer(kind=kint),
parameter :: crank = 2
37 integer(kind=kint),
parameter :: cgid = 3
39 integer(kind=kint),
parameter :: debug = 0
40 integer(kind=kint),
parameter :: timer = 0
47 integer(kind=kint),
intent(in) :: iunit
48 integer(kind=kint) :: nr, nc, nnz, ndof, ndof2, i, js, je, j, jj
54 write(iunit,*)
'nr, nc, nnz, ndof', nr, nc, nnz, ndof
55 write(iunit,*)
'i, j, A'
62 write(iunit,*) i, jj, tmat%A(j)
65 write(iunit,*) tmat%A((j-1)*ndof2+1:j*ndof2)
71 subroutine hecmw_localmat_write_ij(Tmat,iunit)
74 integer(kind=kint),
intent(in) :: iunit
75 integer(kind=kint) :: nr, nc, nnz, ndof, i, js, je, j, jj
80 write(iunit,*)
'nr, nc, nnz, ndof', nr, nc, nnz, ndof
90 end subroutine hecmw_localmat_write_ij
95 integer,
intent(in) :: ndof
97 integer,
allocatable :: iw(:)
98 integer :: ndof2, i, icnt, idof, idx, ls, le, l, j, jb, k, lb0, jdof, ks, ke
101 if (mod(tmat%nr, ndof) /= 0 .or. mod(tmat%nc, ndof) /= 0)
then
102 write(0,*) tmat%nr, tmat%nc, ndof
103 stop
'ERROR: blocking_Tmat failed'
105 btmat%nr=tmat%nr/ndof
106 btmat%nc=tmat%nc/ndof
109 allocate(iw(btmat%nc))
110 allocate(btmat%index(0:btmat%nr))
117 ls=tmat%index(idx-1)+1
123 if (iw(k)==jb) cycle lcol
129 btmat%index(i)=btmat%index(i-1)+icnt
132 btmat%nnz=btmat%index(btmat%nr)
133 allocate(btmat%item(btmat%nnz))
134 allocate(btmat%A(btmat%nnz*ndof2))
141 ls=tmat%index(idx-1)+1
147 if (iw(k)==jb) cycle lcol2
157 btmat%item(lb0+k)=iw(k)
161 ls=tmat%index(idx-1)+1
166 jdof=mod((j-1), ndof)+1
167 ks=btmat%index(i-1)+1
170 if (btmat%item(k)==jb)
then
171 btmat%A((k-1)*ndof2+(idof-1)*ndof+jdof)=tmat%A(l)
175 stop
'ERROR: something wrong in blocking Tmat'
184 deallocate(tmat%index)
185 if (
associated(tmat%item))
deallocate(tmat%item)
186 if (
associated(tmat%A))
deallocate(tmat%A)
194 iwS, num_lagrange, hecTKT)
200 integer(kind=kint),
intent(in) :: iws(:)
201 integer(kind=kint),
intent(in) :: num_lagrange
203 if (hecmesh%n_neighbor_pe == 0)
then
204 call hecmw_trimatmul_ttkt_serial(hecmesh, bttmat, hecmat, btmat, &
205 iws, num_lagrange, hectkt)
207 call hecmw_trimatmul_ttkt_parallel(hecmesh, bttmat, hecmat, btmat, &
208 iws, num_lagrange, hectkt)
212 subroutine hecmw_trimatmul_ttkt_serial(hecMESH, BTtmat, hecMAT, BTmat, &
213 iwS, num_lagrange, hecTKT)
219 integer(kind=kint),
intent(in) :: iws(:)
220 integer(kind=kint),
intent(in) :: num_lagrange
223 real(kind=
kreal) :: num
226 call trimatmul_ttkt(bttmat, hecmat, btmat, bttkt)
233 call place_num_on_diag(bttkt, iws, num_lagrange, num)
238 call make_new_hecmat(hecmat, bttkt, hectkt)
240 end subroutine hecmw_trimatmul_ttkt_serial
242 subroutine hecmw_trimatmul_ttkt_parallel(hecMESH, BTtmat, hecMAT, BTmat, &
243 iwS, num_lagrange, hecTKT)
249 integer(kind=kint),
intent(in) :: iws(:)
250 integer(kind=kint),
intent(in) :: num_lagrange
253 real(kind=
kreal) :: num
254 real(kind=
kreal) :: t0, t1
268 if (timer >= 1)
write(0,
'(A,f10.4)')
"#### hecmw_trimatmul_TtKT_parallel (1) : ",t1-t0
272 if (debug >= 2)
write(0,*)
' DEBUG2: multiply Tt and K done'
283 if (timer >= 1)
write(0,
'(A,f10.4)')
"#### hecmw_trimatmul_TtKT_parallel (2) : ",t1-t0
287 if (debug >= 2)
write(0,*)
' DEBUG2: multiply TtK and T done'
298 if (timer >= 1)
write(0,
'(A,f10.4)')
"#### hecmw_trimatmul_TtKT_parallel (3) : ",t1-t0
304 call place_num_on_diag(bttktmat, iws, num_lagrange, num)
317 if (timer >= 1)
write(0,
'(A,f10.4)')
"#### hecmw_trimatmul_TtKT_parallel (4) : ",t1-t0
321 call make_new_hecmat(hecmat, bttktmat, hectkt)
324 if (timer >= 1)
write(0,
'(A,f10.4)')
"#### hecmw_trimatmul_TtKT_parallel (5) : ",t1-t0
325 end subroutine hecmw_trimatmul_ttkt_parallel
327 subroutine trimatmul_ttkt(BTtmat, hecMAT, BTmat, BTtKT)
332 integer :: nr, nc, ndof, ndof2, i, icnt, js, je, j, jj, ks, ke, k, kk
333 integer :: ls, le, l, ll, m, ms, me, mm
334 integer,
allocatable :: iw(:)
335 real(kind=
kreal),
pointer :: ttp(:), kp(:), tp(:), ttktp(:)
346 allocate(bttkt%index(0:nr))
357 js=bttmat%index(i-1)+1
362 ks=hecmat%indexL(jj-1)+1
366 ls=btmat%index(kk-1)+1
371 if (iw(m)==ll) cycle ll1
379 ls=btmat%index(jj-1)+1
384 if (iw(m)==ll) cycle ll2
391 ks=hecmat%indexU(jj-1)+1
395 ls=btmat%index(kk-1)+1
400 if (iw(m)==ll) cycle ll3
408 if (icnt == 0) icnt=1
421 bttkt%index(i)=bttkt%index(i-1)+bttkt%index(i)
425 bttkt%nnz=bttkt%index(nr)
426 allocate(bttkt%item(bttkt%nnz))
427 allocate(bttkt%A(bttkt%nnz*ndof2))
440 ms=bttkt%index(i-1)+1
442 js=bttmat%index(i-1)+1
446 ttp=>bttmat%A((j-1)*ndof2+1:j*ndof2)
448 ks=hecmat%indexL(jj-1)+1
452 kp=>hecmat%AL((k-1)*ndof2+1:k*ndof2)
453 ls=btmat%index(kk-1)+1
457 tp=>btmat%A((l-1)*ndof2+1:l*ndof2)
461 if (bttkt%item(m)==ll) mm=m
469 ttktp=>bttkt%A((mm-1)*ndof2+1:mm*ndof2)
470 call blk_trimatmul_add(ndof, ttp, kp, tp, ttktp)
474 kp=>hecmat%D((jj-1)*ndof2+1:jj*ndof2)
475 ls=btmat%index(jj-1)+1
479 tp=>btmat%A((l-1)*ndof2+1:l*ndof2)
483 if (bttkt%item(m)==ll) mm=m
491 ttktp=>bttkt%A((mm-1)*ndof2+1:mm*ndof2)
492 call blk_trimatmul_add(ndof, ttp, kp, tp, ttktp)
495 ks=hecmat%indexU(jj-1)+1
499 kp=>hecmat%AU((k-1)*ndof2+1:k*ndof2)
500 ls=btmat%index(kk-1)+1
504 tp=>btmat%A((l-1)*ndof2+1:l*ndof2)
508 if (bttkt%item(m)==ll) mm=m
516 ttktp=>bttkt%A((mm-1)*ndof2+1:mm*ndof2)
517 call blk_trimatmul_add(ndof, ttp, kp, tp, ttktp)
528 if (ms-1+icnt /= bttkt%index(i)) stop
'ERROR: trimatmul'
535 end subroutine trimatmul_ttkt
537 subroutine blk_trimatmul_add(ndof, A, B, C, ABC)
539 integer,
intent(in) :: ndof
540 real(kind=
kreal),
intent(in) :: a(:), b(:), c(:)
541 real(kind=
kreal),
intent(inout) :: abc(:)
542 real(kind=
kreal),
allocatable :: ab(:)
543 integer :: ndof2, i, j, k, i0, j0, ij, ik, jk
557 ab(ik)=ab(ik)+a(ij)*b(jk)
570 abc(ik)=abc(ik)+ab(ij)*c(jk)
576 end subroutine blk_trimatmul_add
578 subroutine place_num_on_diag(BTtKT, iwS, num_lagrange, num)
581 integer(kind=kint),
intent(in) :: iws(:)
582 integer(kind=kint),
intent(in) :: num_lagrange
583 real(kind=
kreal),
intent(in) :: num
584 integer(kind=kint) :: ndof, ndof2, ilag, i, idof, js, je, j, jj
585 integer(kind=kint) :: nmissing, k, ks, ke
586 integer(kind=kint),
allocatable :: missing(:), cnt(:)
587 integer(kind=kint),
pointer :: index(:), item(:)
588 real(kind=
kreal),
pointer :: a(:)
594 allocate(missing(num_lagrange))
596 outer1:
do ilag=1,num_lagrange
597 i=(iws(ilag)-1)/ndof+1
598 idof=mod(iws(ilag)-1, ndof)+1
599 js=bttkt%index(i-1)+1
603 if (jj==i) cycle outer1
607 if (missing(k) == i) cycle outer1
609 nmissing = nmissing + 1
610 missing(nmissing) = i
614 if (nmissing > 0)
then
615 allocate(cnt(bttkt%nr))
616 allocate(index(0:bttkt%nr))
618 cnt(i) = bttkt%index(i) - bttkt%index(i-1)
621 cnt(missing(i)) = cnt(missing(i)) + 1
623 call make_index(bttkt%nr, cnt, index)
624 allocate(item(bttkt%nnz + nmissing))
625 allocate(a(ndof2 * (bttkt%nnz + nmissing)))
628 js=bttkt%index(i-1)+1
630 item(ks:ks+(je-js))=bttkt%item(js:je)
631 a(ndof2*(ks-1)+1:ndof2*(ks+(je-js)))=bttkt%A(ndof2*(js-1)+1:ndof2*je)
636 a(ndof2*(ke-1)+1:ndof2*ke)=0.d0
638 deallocate(bttkt%index)
639 deallocate(bttkt%item)
644 bttkt%nnz = index(bttkt%nr)
650 outer:
do ilag=1,num_lagrange
651 i=(iws(ilag)-1)/ndof+1
652 idof=mod(iws(ilag)-1, ndof)+1
653 js=bttkt%index(i-1)+1
659 bttkt%A((j-1)*ndof2+(idof-1)*ndof+idof)=num
664 end subroutine place_num_on_diag
666 subroutine replace_hecmat(hecMAT, BTtKT)
670 integer :: nr, nc, ndof, ndof2, i, nl, nu, js, je, j, jj
671 integer :: ksl, ksu, k
679 if (
associated(hecmat%AL))
deallocate(hecmat%AL)
680 if (
associated(hecmat%AU))
deallocate(hecmat%AU)
681 if (
associated(hecmat%itemL))
deallocate(hecmat%itemL)
682 if (
associated(hecmat%itemU))
deallocate(hecmat%itemU)
693 js=bttkt%index(i-1)+1
714 hecmat%indexL(i)=hecmat%indexL(i-1)+hecmat%indexL(i)
715 hecmat%indexU(i)=hecmat%indexU(i-1)+hecmat%indexU(i)
717 hecmat%NPL=hecmat%indexL(nc)
718 hecmat%NPU=hecmat%indexU(nc)
721 allocate(hecmat%itemL(hecmat%NPL), hecmat%itemU(hecmat%NPU))
722 allocate(hecmat%AL(hecmat%NPL*ndof2), hecmat%AU(hecmat%NPU*ndof2))
736 js=bttkt%index(i-1)+1
738 ksl=hecmat%indexL(i-1)+1
739 ksu=hecmat%indexU(i-1)+1
745 hecmat%AL((k-1)*ndof2+1:k*ndof2)=bttkt%A((j-1)*ndof2+1:j*ndof2)
750 hecmat%AU((k-1)*ndof2+1:k*ndof2)=bttkt%A((j-1)*ndof2+1:j*ndof2)
753 hecmat%D((i-1)*ndof2+1:i*ndof2)=bttkt%A((j-1)*ndof2+1:j*ndof2)
770 end subroutine replace_hecmat
772 subroutine make_new_hecmat(hecMAT, BTtKT, hecTKT)
777 integer(kind=kint) :: nr, nc, ndof, ndof2
793 if (
associated(hectkt%D))
deallocate(hectkt%D)
794 allocate(hectkt%D(nc*ndof2))
796 if (
associated(hectkt%indexL))
deallocate(hectkt%indexL)
797 if (
associated(hectkt%indexU))
deallocate(hectkt%indexU)
798 allocate(hectkt%indexL(0:nc))
799 allocate(hectkt%indexU(0:nc))
801 hectkt%Iarray=hecmat%Iarray
802 hectkt%Rarray=hecmat%Rarray
804 call replace_hecmat(hectkt, bttkt)
805 end subroutine make_new_hecmat
810 real(kind=
kreal),
intent(in),
target :: v(:)
811 real(kind=
kreal),
intent(out),
target :: tv(:)
812 real(kind=
kreal),
pointer :: tvp(:), tp(:), vp(:)
813 integer :: nr, ndof, ndof2, i, js, je, j, jj, k, kl0, l
832 tvp=>tv((i-1)*ndof+1:i*ndof)
833 js=btmat%index(i-1)+1
837 tp=>btmat%A((j-1)*ndof2+1:j*ndof2)
838 vp=>v((jj-1)*ndof+1:jj*ndof)
842 tvp(k)=tvp(k)+tp(kl0+l)*vp(l)
857 integer(kind=kint),
allocatable :: iws(:)
858 integer(kind=kint) :: ndof, n_mpc, i_mpc
859 integer(kind=kint) :: i, j, k, kk, ilag
860 real(kind=
kreal) :: t0, t1
864 outer:
do i=1,hecmesh%mpc%n_mpc
865 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
866 if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
872 outer2:
do i=1,hecmesh%mpc%n_mpc
873 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
874 if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer2
877 k=hecmesh%mpc%mpc_index(i-1)+1
878 kk=ndof*(hecmesh%mpc%mpc_item(k)-1)+hecmesh%mpc%mpc_dof(k)
885 if (timer >= 1)
write(0,
'(A,f10.4)')
"### hecmw_trimatmul_TtKT_mpc (1) : ",t1-t0
887 call make_btmat_mpc(hecmesh, ndof, btmat)
897 if (timer >= 1)
write(0,
'(A,f10.4)')
"### hecmw_trimatmul_TtKT_mpc (2) : ",t1-t0
915 if (timer >= 1)
write(0,
'(A,f10.4)')
"### hecmw_trimatmul_TtKT_mpc (3) : ",t1-t0
919 if (timer >= 1)
write(0,
'(A,f10.4)')
"### hecmw_trimatmul_TtKT_mpc (4) : ",t1-t0
922 if (
associated(hectkt%B))
deallocate(hectkt%B)
923 if (
associated(hectkt%X))
deallocate(hectkt%X)
924 allocate(hectkt%B(ndof*hectkt%NP))
925 allocate(hectkt%X(ndof*hectkt%NP))
926 do i=1,
size(hecmat%B)
927 hectkt%B(i) = hecmat%B(i)
929 do i=1,
size(hecmat%X)
930 hectkt%X(i) = hecmat%X(i)
933 hectkt%X(iws(ilag)) = 0.d0
941 if (timer >= 1)
write(0,
'(A,f10.4)')
"### hecmw_trimatmul_TtKT_mpc (5) : ",t1-t0
944 subroutine make_btmat_mpc(hecMESH, ndof, BTmat)
947 integer(kind=kint),
intent(in) :: ndof
950 integer(kind=kint) :: n_mpc
951 integer(kind=kint) :: i,j,k,js,jj,kk
953 outer:
do i=1,hecmesh%mpc%n_mpc
954 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
955 if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
959 tmat%nr=hecmesh%n_node*ndof
962 allocate(tmat%index(0:tmat%nr))
964 tmat%index(1:tmat%nr)=1
965 outer2:
do i=1,hecmesh%mpc%n_mpc
966 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
967 if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer2
969 k=hecmesh%mpc%mpc_index(i-1)+1
970 kk=ndof*(hecmesh%mpc%mpc_item(k)-1)+hecmesh%mpc%mpc_dof(k)
971 tmat%index(kk)=hecmesh%mpc%mpc_index(i)-hecmesh%mpc%mpc_index(i-1)-1
976 tmat%index(i)=tmat%index(i-1)+tmat%index(i)
978 tmat%nnz=tmat%index(tmat%nr)
979 allocate(tmat%item(tmat%nnz), tmat%A(tmat%nnz))
987 outer3:
do i=1,hecmesh%mpc%n_mpc
988 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
989 if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer3
991 k=hecmesh%mpc%mpc_index(i-1)+1
992 kk=ndof*(hecmesh%mpc%mpc_item(k)-1)+hecmesh%mpc%mpc_dof(k)
993 js=tmat%index(kk-1)+1
994 do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
995 jj = ndof * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
997 tmat%A(js)=-hecmesh%mpc%mpc_val(j)
1005 end subroutine make_btmat_mpc
1007 subroutine make_bttmat_mpc(hecMESH, ndof, BTtmat)
1010 integer(kind=kint),
intent(in) :: ndof
1013 integer(kind=kint) :: n_mpc
1014 integer(kind=kint) :: i,j,k,js,je,jj,kk
1015 integer(kind=kint),
allocatable :: iw(:)
1017 outer:
do i=1,hecmesh%mpc%n_mpc
1018 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
1019 if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer
1023 ttmat%nr=hecmesh%n_node*ndof
1026 allocate(ttmat%index(0:ttmat%nr))
1028 ttmat%index(1:ttmat%nr)=1
1029 outer2:
do i=1,hecmesh%mpc%n_mpc
1030 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
1031 if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer2
1033 k=hecmesh%mpc%mpc_index(i-1)+1
1034 kk=ndof*(hecmesh%mpc%mpc_item(k)-1)+hecmesh%mpc%mpc_dof(k)
1036 do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
1037 jj = ndof * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
1038 ttmat%index(jj)=ttmat%index(jj)+1
1044 ttmat%index(i)=ttmat%index(i-1)+ttmat%index(i)
1046 ttmat%nnz=ttmat%index(ttmat%nr)
1047 allocate(ttmat%item(ttmat%nnz), ttmat%A(ttmat%nnz))
1050 js=ttmat%index(i-1)+1
1058 allocate(iw(ttmat%nr))
1060 outer3:
do i=1,hecmesh%mpc%n_mpc
1061 do j= hecmesh%mpc%mpc_index(i-1) + 1, hecmesh%mpc%mpc_index(i)
1062 if (hecmesh%mpc%mpc_dof(j) > ndof) cycle outer3
1064 k=hecmesh%mpc%mpc_index(i-1)+1
1065 kk=ndof*(hecmesh%mpc%mpc_item(k)-1)+hecmesh%mpc%mpc_dof(k)
1066 do j= hecmesh%mpc%mpc_index(i-1) + 2, hecmesh%mpc%mpc_index(i)
1067 jj = ndof * (hecmesh%mpc%mpc_item(j) - 1) + hecmesh%mpc%mpc_dof(j)
1068 js=ttmat%index(jj-1)+1+iw(jj)
1070 ttmat%A(js)=-hecmesh%mpc%mpc_val(j)
1079 end subroutine make_bttmat_mpc
1085 integer(kind=kint),
allocatable :: iw(:)
1086 integer(kind=kint) :: i, j, jj, ndof, ndof2, k, idof, jdof
1087 allocate(iw(tmat%nc))
1090 do j = tmat%index(i-1)+1, tmat%index(i)
1097 ttmat%nnz = tmat%nnz
1098 ttmat%ndof = tmat%ndof
1101 allocate(ttmat%index(0:ttmat%nr))
1102 allocate(ttmat%item(ttmat%nnz))
1103 allocate(ttmat%A(ttmat%nnz*ndof2))
1106 ttmat%index(i) = ttmat%index(i-1) + iw(i)
1107 iw(i) = ttmat%index(i-1) + 1
1110 do j = tmat%index(i-1)+1, tmat%index(i)
1116 ttmat%A((k-1)*ndof2+(idof-1)*ndof+jdof) = &
1117 tmat%A((j-1)*ndof2+(jdof-1)*ndof+idof)
1125 function hecmw_localmat_equal(Tmat1, Tmat2)
1128 integer(kind=kint) :: hecmw_localmat_equal
1129 integer(kind=kint) :: i, j, k0, k, ndof, ndof2
1130 hecmw_localmat_equal = 0
1131 if (tmat1%nr /= tmat2%nr)
return
1132 if (tmat1%nc /= tmat2%nc)
return
1133 if (tmat1%nnz /= tmat2%nnz)
return
1134 if (tmat1%ndof /= tmat2%ndof)
return
1138 if (tmat1%index(i) /= tmat2%index(i))
return
1139 do j = tmat1%index(i-1)+1, tmat1%index(i)
1140 if (tmat1%item(j) /= tmat2%item(j))
return
1143 if (tmat1%A(k0+k) /= tmat2%A(k0+k))
return
1147 hecmw_localmat_equal = 1
1148 end function hecmw_localmat_equal
1159 integer(kind=kint) :: nn_int, np, ndof, ndof2, nr_ext, nnz_ext
1160 integer(kind=kint),
allocatable :: exp_rows_index(:), exp_cols_index(:)
1161 integer(kind=kint),
allocatable :: exp_rows_item(:,:), exp_cols_item(:,:)
1166 if (debug >= 1)
write(0,*)
'DEBUG: nr,nc,nnz,ndof',btmat%nr,btmat%nc,btmat%nnz,btmat%ndof
1167 if (btmat%nr /= hecmesh%n_node) stop
'ERROR: invalid size in hecmw_localmat_assemble'
1169 nn_int = hecmesh%nn_internal
1174 nr_ext = np - nn_int
1175 nnz_ext = btmat%index(np) - btmat%index(nn_int)
1177 call prepare_bt_ext(btmat, hecmesh, exp_rows_index, exp_rows_item, bt_ext)
1178 if (debug >= 1)
write(0,*)
'DEBUG: prepare_BT_ext done'
1180 call prepare_column_info(hecmesh, bt_ext, exp_cols_index, exp_cols_item)
1181 if (debug >= 1)
write(0,*)
'DEBUG: prepare_column info done'
1183 call send_bt_ext_and_recv_bt_int(hecmesh, exp_rows_index, exp_rows_item, bt_ext, &
1184 exp_cols_index, exp_cols_item, bt_int, hecmeshnew)
1185 if (debug >= 1)
write(0,*)
'DEBUG: send BT_ext and recv BT_int done'
1189 if (debug >= 1)
write(0,*)
'DEBUG: localmat_add done'
1196 btmat%nnz = btnew%nnz
1197 btmat%ndof = btnew%ndof
1198 btmat%index => btnew%index
1199 btmat%item => btnew%item
1219 if (debug >= 1)
write(0,*)
'DEBUG: update BTmat and hecMESH done'
1222 subroutine prepare_bt_ext(BTmat, hecMESH, exp_rows_index, exp_rows_item, BT_ext)
1226 integer(kind=kint),
allocatable,
intent(out) :: exp_rows_index(:)
1227 integer(kind=kint),
allocatable,
intent(out) :: exp_rows_item(:,:)
1229 integer(kind=kint),
allocatable :: incl_nz(:), exp_cols_per_row(:), exp_rows_per_rank(:)
1230 integer(kind=kint) :: nn_int
1231 nn_int = hecmesh%nn_internal
1233 call check_external_nz_blocks(btmat, nn_int, incl_nz)
1235 call count_ext_rows_with_nz(btmat, nn_int, incl_nz, exp_cols_per_row)
1237 call count_exp_rows_per_rank(hecmesh, exp_cols_per_row, exp_rows_per_rank)
1239 allocate(exp_rows_index(0:hecmesh%n_neighbor_pe))
1240 call make_index(hecmesh%n_neighbor_pe, exp_rows_per_rank, exp_rows_index)
1243 deallocate(exp_rows_per_rank)
1245 call make_exp_rows_item(hecmesh, exp_cols_per_row, exp_rows_index, exp_rows_item)
1247 deallocate(exp_cols_per_row)
1249 allocate(bt_ext(hecmesh%n_neighbor_pe))
1250 call extract_bt_ext(hecmesh, btmat, incl_nz, exp_rows_index, exp_rows_item, bt_ext)
1253 end subroutine prepare_bt_ext
1255 subroutine check_external_nz_blocks(BTmat, nn_internal, incl_nz)
1258 integer(kind=kint),
intent(in) :: nn_internal
1259 integer(kind=kint),
allocatable,
intent(out) :: incl_nz(:)
1260 integer(kind=kint) :: ndof2, i0, nnz_ext, i, k, nnz_blk
1261 if (nn_internal > btmat%nr) stop
'ERROR: invalid nn_internal'
1262 ndof2 = btmat%ndof ** 2
1263 i0 = btmat%index(nn_internal)
1264 nnz_ext = btmat%index(btmat%nr) - i0
1265 allocate(incl_nz(nnz_ext))
1270 if (btmat%A(ndof2*(i0+i-1)+k) /= 0.0d0)
then
1272 nnz_blk = nnz_blk + 1
1277 if (debug >= 1)
write(0,*)
'DEBUG: nnz_blk',nnz_blk
1278 end subroutine check_external_nz_blocks
1280 subroutine count_ext_rows_with_nz(BTmat, nn_internal, incl_nz, exp_cols_per_row)
1283 integer(kind=kint),
intent(in) :: nn_internal
1284 integer(kind=kint),
intent(in) :: incl_nz(:)
1285 integer(kind=kint),
allocatable,
intent(out) :: exp_cols_per_row(:)
1286 integer(kind=kint) :: nr_ext, nnz_int, i, irow, js, je, j, jcol
1287 nr_ext = btmat%nr - nn_internal
1288 nnz_int = btmat%index(nn_internal)
1289 allocate(exp_cols_per_row(nr_ext))
1290 exp_cols_per_row(:) = 0
1292 irow = nn_internal+i
1293 js = btmat%index(irow-1)+1
1294 je = btmat%index(irow)
1296 jcol = btmat%item(j)
1297 if (incl_nz(j-nnz_int) == 1) exp_cols_per_row(i) = exp_cols_per_row(i) + 1
1301 end subroutine count_ext_rows_with_nz
1303 subroutine count_exp_rows_per_rank(hecMESH, exp_cols_per_row, exp_rows_per_rank)
1306 integer(kind=kint),
intent(in) :: exp_cols_per_row(:)
1307 integer(kind=kint),
allocatable,
intent(out) :: exp_rows_per_rank(:)
1308 integer(kind=kint) :: nn_int, np, nr_ext, i, irow, exp_rank, idom
1309 allocate(exp_rows_per_rank(hecmesh%n_neighbor_pe))
1310 exp_rows_per_rank(1:hecmesh%n_neighbor_pe) = 0
1311 nn_int = hecmesh%nn_internal
1313 nr_ext = np - nn_int
1315 if (exp_cols_per_row(i) > 0)
then
1317 exp_rank = hecmesh%node_ID(2*irow)
1318 call rank_to_idom(hecmesh, exp_rank, idom)
1319 exp_rows_per_rank(idom) = exp_rows_per_rank(idom) + 1
1323 end subroutine count_exp_rows_per_rank
1325 subroutine rank_to_idom(hecMESH, rank, idom)
1328 integer(kind=kint),
intent(in) :: rank
1329 integer(kind=kint),
intent(out) :: idom
1330 integer(kind=kint) :: i
1331 do i = 1, hecmesh%n_neighbor_pe
1332 if (hecmesh%neighbor_pe(i) == rank)
then
1337 stop
'ERROR: exp_rank not found in neighbor_pe'
1338 end subroutine rank_to_idom
1340 subroutine make_index(len, cnt, index)
1342 integer(kind=kint),
intent(in) :: len
1343 integer(kind=kint),
intent(in) :: cnt(len)
1344 integer(kind=kint),
intent(out) :: index(0:)
1345 integer(kind=kint) :: i
1349 index(i) = index(i-1) + cnt(i)
1351 end subroutine make_index
1353 subroutine make_exp_rows_item(hecMESH, exp_cols_per_row, exp_rows_index, exp_rows_item)
1356 integer(kind=kint),
intent(in) :: exp_cols_per_row(:)
1357 integer(kind=kint),
allocatable,
intent(in) :: exp_rows_index(:)
1358 integer(kind=kint),
allocatable,
intent(out) :: exp_rows_item(:,:)
1359 integer(kind=kint),
allocatable :: cnt(:)
1360 integer(kind=kint) :: nn_int, np, nr_ext, i, irow, exp_rank, idom, idx
1361 allocate(exp_rows_item(2,exp_rows_index(hecmesh%n_neighbor_pe)))
1362 allocate(cnt(hecmesh%n_neighbor_pe))
1364 nn_int = hecmesh%nn_internal
1366 nr_ext = np - nn_int
1368 if (exp_cols_per_row(i) > 0)
then
1370 exp_rank = hecmesh%node_ID(2*irow)
1371 call rank_to_idom(hecmesh, exp_rank, idom)
1372 cnt(idom) = cnt(idom) + 1
1373 idx = exp_rows_index(idom-1) + cnt(idom)
1374 exp_rows_item(1,idx) = irow
1375 exp_rows_item(2,idx) = exp_cols_per_row(i)
1379 do idom = 1, hecmesh%n_neighbor_pe
1380 if (cnt(idom) /= exp_rows_index(idom)-exp_rows_index(idom-1)) stop
'ERROR: make exp_rows_item'
1384 end subroutine make_exp_rows_item
1386 subroutine extract_bt_ext(hecMESH, BTmat, incl_nz, exp_rows_index, exp_rows_item, BT_ext)
1390 integer(kind=kint),
intent(in) :: incl_nz(:)
1391 integer(kind=kint),
allocatable,
intent(in) :: exp_rows_index(:)
1392 integer(kind=kint),
intent(in) :: exp_rows_item(:,:)
1394 integer(kind=kint) :: ndof, ndof2, nn_int, nnz_int, idom, j, idx, ncol, cnt, jrow, ks, ke, k, kcol
1395 allocate(bt_ext(hecmesh%n_neighbor_pe))
1398 nn_int = hecmesh%nn_internal
1399 nnz_int = btmat%index(nn_int)
1400 do idom = 1, hecmesh%n_neighbor_pe
1401 bt_ext(idom)%nr = exp_rows_index(idom) - exp_rows_index(idom-1)
1402 bt_ext(idom)%nc = btmat%nc
1403 bt_ext(idom)%nnz = 0
1404 bt_ext(idom)%ndof = ndof
1405 allocate(bt_ext(idom)%index(0:bt_ext(idom)%nr))
1406 bt_ext(idom)%index(0) = 0
1407 do j = 1, bt_ext(idom)%nr
1408 idx = exp_rows_index(idom-1) + j
1409 ncol = exp_rows_item(2,idx)
1410 bt_ext(idom)%index(j) = bt_ext(idom)%index(j-1) + ncol
1412 bt_ext(idom)%nnz = bt_ext(idom)%index(bt_ext(idom)%nr)
1413 if (debug >= 1)
write(0,*)
'DEBUG: idom,nr,nc,nnz,ndof', &
1414 idom,bt_ext(idom)%nr,bt_ext(idom)%nc,bt_ext(idom)%nnz,bt_ext(idom)%ndof
1415 allocate(bt_ext(idom)%item(bt_ext(idom)%nnz))
1416 allocate(bt_ext(idom)%A(bt_ext(idom)%nnz * ndof2))
1418 do j = 1, bt_ext(idom)%nr
1419 idx = exp_rows_index(idom-1) + j
1420 jrow = exp_rows_item(1,idx)
1421 if (jrow < 1 .or. btmat%nr < jrow) stop
'ERROR: extract BT_ext: jrow'
1422 ks = btmat%index(jrow-1)+1
1423 ke = btmat%index(jrow)
1425 kcol = btmat%item(k)
1426 if (incl_nz(k-nnz_int) == 0) cycle
1428 bt_ext(idom)%item(cnt) = kcol
1429 bt_ext(idom)%A(ndof2*(cnt-1)+1:ndof2*cnt) = btmat%A(ndof2*(k-1)+1:ndof2*k)
1431 if (cnt /= bt_ext(idom)%index(j)) stop
'ERROR: extract BT_ext'
1436 end subroutine extract_bt_ext
1438 subroutine prepare_column_info(hecMESH, BT_ext, exp_cols_index, exp_cols_item)
1442 integer(kind=kint),
allocatable,
intent(out) :: exp_cols_index(:)
1443 integer(kind=kint),
allocatable,
intent(out) :: exp_cols_item(:,:)
1445 call make_exp_cols_index(hecmesh%n_neighbor_pe, bt_ext, exp_cols_index)
1446 if (debug >= 2)
write(0,*)
' DEBUG2: make exp_cols_index done'
1447 if (debug >= 3)
write(0,*)
' DEBUG3: exp_cols_index', exp_cols_index(0:hecmesh%n_neighbor_pe)
1451 call make_exp_cols_item(hecmesh, bt_ext, exp_cols_index, exp_cols_item)
1452 if (debug >= 2)
write(0,*)
' DEBUG2: make exp_cols_item done'
1454 end subroutine prepare_column_info
1456 subroutine make_exp_cols_index(nnb, BT_ext, exp_cols_index)
1458 integer(kind=kint),
intent(in) :: nnb
1460 integer(kind=kint),
allocatable,
intent(out) :: exp_cols_index(:)
1461 integer(kind=kint) :: idom
1462 allocate(exp_cols_index(0:nnb))
1463 exp_cols_index(0) = 0
1465 exp_cols_index(idom) = exp_cols_index(idom-1) + bt_ext(idom)%nnz
1467 end subroutine make_exp_cols_index
1469 subroutine make_exp_cols_item(hecMESH, BT_ext, exp_cols_index, exp_cols_item)
1473 integer(kind=kint),
allocatable,
intent(in) :: exp_cols_index(:)
1474 integer(kind=kint),
allocatable,
intent(out) :: exp_cols_item(:,:)
1475 integer(kind=kint) :: cnt, idom, j, jcol
1476 allocate(exp_cols_item(
cncol_item,exp_cols_index(hecmesh%n_neighbor_pe)))
1478 do idom = 1, hecmesh%n_neighbor_pe
1479 do j = 1, bt_ext(idom)%nnz
1481 jcol = bt_ext(idom)%item(j)
1487 exp_cols_item(clid,cnt) = hecmesh%node_ID(2*jcol-1)
1488 exp_cols_item(crank,cnt) = hecmesh%node_ID(2*jcol)
1489 if (
cncol_item >= 3) exp_cols_item(cgid,cnt) = hecmesh%global_node_ID(jcol)
1492 if (cnt /= exp_cols_index(idom)) stop
'ERROR: make exp_cols_item'
1494 end subroutine make_exp_cols_item
1496 subroutine send_bt_ext_and_recv_bt_int(hecMESH, exp_rows_index, exp_rows_item, BT_ext, &
1497 exp_cols_index, exp_cols_item, BT_int, hecMESHnew)
1500 integer(kind=kint),
allocatable,
intent(inout) :: exp_rows_index(:), exp_cols_index(:)
1501 integer(kind=kint),
allocatable,
intent(inout) :: exp_rows_item(:,:), exp_cols_item(:,:)
1505 integer(kind=kint),
allocatable :: imp_rows_index(:), imp_cols_index(:)
1506 integer(kind=kint),
allocatable :: imp_rows_item(:,:), imp_cols_item(:,:)
1507 real(kind=
kreal),
allocatable :: imp_vals_item(:)
1508 integer(kind=kint),
allocatable :: map(:), add_nodes(:,:)
1509 integer(kind=kint) :: ndof, ndof2, idom, n_add_node, i0
1510 if (hecmesh%n_neighbor_pe == 0)
return
1511 ndof = bt_ext(1)%ndof
1514 call convert_rowid_to_remote_localid(hecmesh, exp_rows_index(hecmesh%n_neighbor_pe), exp_rows_item)
1515 if (debug >= 2)
write(0,*)
' DEBUG2: convert rowID to remote localID done'
1517 call send_recv_bt_ext_nr_nnz(hecmesh, bt_ext, imp_rows_index, imp_cols_index)
1518 if (debug >= 2)
write(0,*)
' DEBUG2: send recv BT_ext nr and nnz done'
1520 call send_recv_bt_ext_contents(hecmesh, bt_ext, &
1521 exp_rows_index, exp_cols_index, exp_rows_item, exp_cols_item, &
1522 imp_rows_index, imp_cols_index, &
1523 imp_rows_item, imp_cols_item, imp_vals_item)
1524 if (debug >= 2)
write(0,*)
' DEBUG2: send recv BT_ext contents done'
1526 do idom = 1, hecmesh%n_neighbor_pe
1531 call allocate_bt_int(hecmesh, ndof, imp_rows_index, imp_rows_item, bt_int)
1532 if (debug >= 2)
write(0,*)
' DEBUG2: allocate BT_int done'
1537 call map_imported_cols(hecmeshnew, imp_cols_index(hecmesh%n_neighbor_pe), &
1538 imp_cols_item, n_add_node, add_nodes, map, i0)
1539 if (debug >= 2)
write(0,*)
' DEBUG2: map imported cols done'
1541 call update_comm_table(hecmeshnew, n_add_node, add_nodes, i0)
1542 if (debug >= 2)
write(0,*)
' DEBUG2: update comm_table done'
1544 bt_int%nc = hecmeshnew%n_node
1546 call copy_vals_to_bt_int(hecmesh%n_neighbor_pe, imp_rows_index, imp_cols_index, &
1547 imp_rows_item, map, ndof2, imp_vals_item, bt_int)
1548 if (debug >= 2)
write(0,*)
' DEBUG2: copy vals to BT_int done'
1550 deallocate(imp_rows_index)
1551 deallocate(imp_cols_index)
1552 deallocate(imp_rows_item)
1553 deallocate(imp_cols_item)
1554 deallocate(imp_vals_item)
1557 call sort_and_uniq_rows(bt_int)
1558 if (debug >= 2)
write(0,*)
' DEBUG2: sort and uniq rows of BT_int done'
1559 end subroutine send_bt_ext_and_recv_bt_int
1561 subroutine convert_rowid_to_remote_localid(hecMESH, len, exp_rows_item)
1564 integer(kind=kint),
intent(in) :: len
1565 integer(kind=kint),
intent(out) :: exp_rows_item(:,:)
1566 integer(kind=kint) :: i
1568 exp_rows_item(1,i) = hecmesh%node_ID(2 * exp_rows_item(1,i) - 1)
1570 end subroutine convert_rowid_to_remote_localid
1572 subroutine send_recv_bt_ext_nr_nnz(hecMESH, BT_ext, imp_rows_index, imp_cols_index)
1577 integer(kind=kint),
allocatable,
intent(out) :: imp_rows_index(:), imp_cols_index(:)
1578 integer(kind=kint) :: nnb, idom, irank, tag, recvbuf(2)
1579 integer(kind=kint),
allocatable :: sendbuf(:,:)
1580 integer(kind=kint),
allocatable :: requests(:)
1581 integer(kind=kint),
allocatable :: statuses(:,:)
1582 nnb = hecmesh%n_neighbor_pe
1583 allocate(imp_rows_index(0:nnb))
1584 allocate(imp_cols_index(0:nnb))
1585 allocate(requests(nnb))
1587 allocate(sendbuf(2,nnb))
1589 irank = hecmesh%neighbor_pe(idom)
1591 sendbuf(1,idom) = bt_ext(idom)%nr
1592 sendbuf(2,idom) = bt_ext(idom)%nnz
1594 call hecmw_isend_int(sendbuf(1,idom), 2, irank, tag, hecmesh%MPI_COMM, &
1597 imp_rows_index(0) = 0
1598 imp_cols_index(0) = 0
1600 irank = hecmesh%neighbor_pe(idom)
1603 hecmesh%MPI_COMM, statuses(:,1))
1604 imp_rows_index(idom) = imp_rows_index(idom-1) + recvbuf(1)
1605 imp_cols_index(idom) = imp_cols_index(idom-1) + recvbuf(2)
1608 deallocate(requests)
1609 deallocate(statuses)
1611 end subroutine send_recv_bt_ext_nr_nnz
1613 subroutine send_recv_bt_ext_contents(hecMESH, BT_ext, &
1614 exp_rows_index, exp_cols_index, exp_rows_item, exp_cols_item, &
1615 imp_rows_index, imp_cols_index, &
1616 imp_rows_item, imp_cols_item, imp_vals_item)
1621 integer(kind=kint),
allocatable,
intent(inout) :: exp_rows_index(:), exp_cols_index(:)
1622 integer(kind=kint),
allocatable,
intent(inout) :: exp_rows_item(:,:), exp_cols_item(:,:)
1623 integer(kind=kint),
allocatable,
intent(in) :: imp_rows_index(:), imp_cols_index(:)
1624 integer(kind=kint),
allocatable,
intent(out) :: imp_rows_item(:,:), imp_cols_item(:,:)
1625 real(kind=
kreal),
allocatable,
intent(out) :: imp_vals_item(:)
1626 integer(kind=kint) :: nnb, ndof2, n_send, idom, irank, tag, nr, nnz
1627 integer(kind=kint),
allocatable :: requests(:)
1628 integer(kind=kint),
allocatable :: statuses(:,:)
1629 nnb = hecmesh%n_neighbor_pe
1631 ndof2 = bt_ext(1)%ndof ** 2
1632 allocate(imp_rows_item(2,imp_rows_index(nnb)))
1633 allocate(imp_cols_item(
cncol_item,imp_cols_index(nnb)))
1634 allocate(imp_vals_item(ndof2*imp_cols_index(nnb)))
1635 allocate(requests(3*nnb))
1639 irank = hecmesh%neighbor_pe(idom)
1640 if (bt_ext(idom)%nr > 0)
then
1644 2*bt_ext(idom)%nr, irank, tag, hecmesh%MPI_COMM, &
1649 cncol_item*bt_ext(idom)%nnz, irank, tag, hecmesh%MPI_COMM, &
1653 call hecmw_isend_r(bt_ext(idom)%A, ndof2*bt_ext(idom)%nnz, irank, &
1654 tag, hecmesh%MPI_COMM, requests(n_send))
1658 irank = hecmesh%neighbor_pe(idom)
1659 nr = imp_rows_index(idom) - imp_rows_index(idom-1)
1660 nnz = imp_cols_index(idom) - imp_cols_index(idom-1)
1664 2*nr, irank, tag, hecmesh%MPI_COMM, statuses(:,1))
1667 cncol_item*nnz, irank, tag, hecmesh%MPI_COMM, statuses(:,1))
1669 call hecmw_recv_r(imp_vals_item(ndof2*imp_cols_index(idom-1)+1), &
1670 ndof2*nnz, irank, tag, hecmesh%MPI_COMM, statuses(:,1))
1674 deallocate(exp_rows_index)
1675 deallocate(exp_rows_item)
1676 deallocate(exp_cols_index)
1677 deallocate(exp_cols_item)
1678 end subroutine send_recv_bt_ext_contents
1680 subroutine allocate_bt_int(hecMESH, ndof, imp_rows_index, imp_rows_item, BT_int)
1683 integer(kind=kint),
intent(in) :: ndof
1684 integer(kind=kint),
allocatable,
intent(in) :: imp_rows_index(:), imp_rows_item(:,:)
1686 integer(kind=kint),
allocatable :: cnt(:)
1687 integer(kind=kint) :: idom, is, ie, i, irow, ncol, ndof2
1689 bt_int%nr = hecmesh%nn_internal
1690 bt_int%nc = hecmesh%n_node
1693 allocate(cnt(bt_int%nr))
1695 do idom = 1, hecmesh%n_neighbor_pe
1696 is = imp_rows_index(idom-1)+1
1697 ie = imp_rows_index(idom)
1699 irow = imp_rows_item(1,i)
1700 ncol = imp_rows_item(2,i)
1701 if (irow < 1 .or. bt_int%nr < irow) stop
'ERROR: allocate BT_int'
1702 cnt(irow) = cnt(irow) + ncol
1706 allocate(bt_int%index(0:bt_int%nr))
1707 call make_index(bt_int%nr, cnt, bt_int%index)
1709 bt_int%nnz = bt_int%index(bt_int%nr)
1710 allocate(bt_int%item(bt_int%nnz))
1711 allocate(bt_int%A(bt_int%nnz * ndof2))
1713 end subroutine allocate_bt_int
1715 subroutine copy_mesh(src, dst)
1720 dst%MPI_COMM = src%MPI_COMM
1721 dst%PETOT = src%PETOT
1722 dst%PEsmpTOT = src%PEsmpTOT
1723 dst%my_rank = src%my_rank
1724 dst%n_subdomain = src%n_subdomain
1725 dst%n_node = src%n_node
1726 dst%nn_internal = src%nn_internal
1727 dst%n_dof = src%n_dof
1728 dst%n_neighbor_pe = src%n_neighbor_pe
1729 allocate(dst%neighbor_pe(dst%n_neighbor_pe))
1730 dst%neighbor_pe(:) = src%neighbor_pe(:)
1731 allocate(dst%import_index(0:dst%n_neighbor_pe))
1732 allocate(dst%export_index(0:dst%n_neighbor_pe))
1733 dst%import_index(:)= src%import_index(:)
1734 dst%export_index(:)= src%export_index(:)
1735 allocate(dst%import_item(dst%import_index(dst%n_neighbor_pe)))
1736 dst%import_item(:) = src%import_item(:)
1737 allocate(dst%export_item(dst%export_index(dst%n_neighbor_pe)))
1738 dst%export_item(:) = src%export_item(:)
1739 allocate(dst%node_ID(2*dst%n_node))
1740 dst%node_ID(1:2*dst%n_node) = src%node_ID(1:2*src%n_node)
1741 allocate(dst%global_node_ID(dst%n_node))
1742 dst%global_node_ID(1:dst%n_node) = src%global_node_ID(1:src%n_node)
1744 dst%node => src%node
1745 end subroutine copy_mesh
1747 subroutine map_imported_cols(hecMESHnew, ncols, cols, n_add_node, add_nodes, map, i0)
1750 integer(kind=kint),
intent(in) :: ncols
1751 integer(kind=kint),
intent(in) :: cols(
cncol_item,ncols)
1752 integer(kind=kint),
allocatable,
intent(out) :: map(:)
1753 integer(kind=kint),
intent(out) :: n_add_node
1754 integer(kind=kint),
allocatable,
intent(out) :: add_nodes(:,:)
1755 integer(kind=kint),
intent(out) :: i0
1756 allocate(map(ncols))
1758 call map_present_nodes(hecmeshnew, ncols, cols, map, n_add_node)
1762 call extract_add_nodes(ncols, cols, map, n_add_node, add_nodes)
1764 call append_nodes(hecmeshnew, n_add_node, add_nodes, i0)
1766 call map_additional_nodes(ncols, cols, n_add_node, add_nodes, i0, map)
1767 end subroutine map_imported_cols
1769 subroutine map_present_nodes(hecMESH, ncols, cols, map, n_add_node)
1772 integer(kind=kint),
intent(in) :: ncols
1773 integer(kind=kint),
intent(in) :: cols(
cncol_item,ncols)
1774 integer(kind=kint),
intent(out) :: map(ncols)
1775 integer(kind=kint),
intent(out) :: n_add_node
1776 integer(kind=kint) :: i, j, lid, rank, llid, n_ext_node, idx
1777 integer(kind=kint),
allocatable :: ext_node(:)
1781 do i = hecmesh%nn_internal + 1, hecmesh%n_node
1788 allocate(ext_node(ncols))
1796 rank = cols(crank,i)
1798 if (rank == hecmesh%my_rank)
then
1802 n_ext_node = n_ext_node + 1
1810 do j = 1, n_ext_node
1813 rank = cols(crank,i)
1820 n_add_node = n_add_node + 1
1825 deallocate(ext_node)
1828 end subroutine map_present_nodes
1830 subroutine extract_add_nodes(ncols, cols, map, n_add_node, add_nodes)
1832 integer(kind=kint),
intent(in) :: ncols
1833 integer(kind=kint),
intent(in) :: cols(
cncol_item,ncols), map(ncols)
1834 integer(kind=kint),
intent(inout) :: n_add_node
1835 integer(kind=kint),
allocatable,
intent(out) :: add_nodes(:,:)
1836 integer(kind=kint) :: cnt, i
1840 if (map(i) == -1)
then
1845 if (cnt /= n_add_node) stop
'ERROR: extract add_nodes'
1846 call sort_and_uniq_add_nodes(n_add_node, add_nodes)
1847 end subroutine extract_add_nodes
1849 subroutine sort_and_uniq_add_nodes(n_add_node, add_nodes)
1851 integer(kind=kint),
intent(inout) :: n_add_node
1852 integer(kind=kint),
intent(inout) :: add_nodes(
cncol_item,n_add_node)
1853 integer(kind=kint) :: ndup
1854 call sort_add_nodes(add_nodes, 1, n_add_node)
1855 call uniq_add_nodes(add_nodes, n_add_node, ndup)
1856 n_add_node = n_add_node - ndup
1857 end subroutine sort_and_uniq_add_nodes
1859 recursive subroutine sort_add_nodes(add_nodes, id1, id2)
1861 integer(kind=kint),
intent(inout) :: add_nodes(:,:)
1862 integer(kind=kint),
intent(in) :: id1, id2
1863 integer(kind=kint) :: center, left, right
1865 if (id1 >= id2)
return
1866 center = (id1 + id2) / 2
1871 do while ((add_nodes(crank,left) < pivot(crank)) .or. &
1872 (add_nodes(crank,left) == pivot(crank) .and. add_nodes(clid,left) < pivot(clid)))
1875 do while ((pivot(crank) < add_nodes(crank,right)) .or. &
1876 (pivot(crank) == add_nodes(crank,right) .and. pivot(clid) < add_nodes(clid,right)))
1879 if (left >= right)
exit
1886 if (id1 < left-1)
call sort_add_nodes(add_nodes, id1, left-1)
1887 if (right+1 < id2)
call sort_add_nodes(add_nodes, right+1, id2)
1889 end subroutine sort_add_nodes
1891 subroutine uniq_add_nodes(add_nodes, len, ndup)
1893 integer(kind=kint),
intent(inout) :: add_nodes(:,:)
1894 integer(kind=kint),
intent(in) :: len
1895 integer(kind=kint),
intent(out) :: ndup
1896 integer(kind=kint) :: i
1899 if (add_nodes(clid,i) == add_nodes(clid,i-1-ndup) .and. &
1900 add_nodes(crank,i) == add_nodes(crank,i-1-ndup))
then
1902 else if (ndup > 0)
then
1906 end subroutine uniq_add_nodes
1908 subroutine search_add_nodes(n_add_node, add_nodes, rank, lid, idx)
1910 integer(kind=kint),
intent(in) :: n_add_node
1911 integer(kind=kint),
intent(in) :: add_nodes(
cncol_item,n_add_node)
1912 integer(kind=kint),
intent(in) :: rank
1913 integer(kind=kint),
intent(in) :: lid
1914 integer(kind=kint),
intent(out) :: idx
1915 integer(kind=kint) :: left, right, center
1918 do while (left <= right)
1919 center = (left + right) / 2
1920 if ((rank == add_nodes(crank,center)) .and. (lid == add_nodes(clid,center)))
then
1923 else if ((rank < add_nodes(crank,center)) .or. &
1924 (rank == add_nodes(crank,center) .and. lid < add_nodes(clid,center)))
then
1926 else if ((add_nodes(crank,center) < rank) .or. &
1927 (add_nodes(crank,center) == rank .and. add_nodes(clid,center) < lid))
then
1932 end subroutine search_add_nodes
1934 subroutine append_nodes(hecMESHnew, n_add_node, add_nodes, i0)
1937 integer(kind=kint),
intent(in) :: n_add_node
1938 integer(kind=kint),
intent(in) :: add_nodes(:,:)
1939 integer(kind=kint),
intent(out) :: i0
1940 integer(kind=kint) :: n_node, i, ii
1941 integer(kind=kint),
pointer :: node_id(:), global_node_id(:)
1942 i0 = hecmeshnew%n_node
1943 n_node = hecmeshnew%n_node + n_add_node
1944 allocate(node_id(2*n_node))
1945 allocate(global_node_id(n_node))
1946 do i = 1, hecmeshnew%n_node
1947 node_id(2*i-1) = hecmeshnew%node_ID(2*i-1)
1948 node_id(2*i ) = hecmeshnew%node_ID(2*i )
1949 global_node_id(i) = hecmeshnew%global_node_ID(i)
1951 do i = 1, n_add_node
1952 ii = hecmeshnew%n_node + i
1953 node_id(2*ii-1) = add_nodes(clid,i)
1954 node_id(2*ii ) = add_nodes(crank,i)
1956 global_node_id(i) = add_nodes(cgid,i)
1958 global_node_id(i) = -1
1961 deallocate(hecmeshnew%node_ID)
1962 deallocate(hecmeshnew%global_node_ID)
1963 hecmeshnew%n_node = n_node
1964 hecmeshnew%node_ID => node_id
1965 hecmeshnew%global_node_ID => global_node_id
1966 end subroutine append_nodes
1968 subroutine map_additional_nodes(ncols, cols, n_add_node, add_nodes, i0, map)
1970 integer(kind=kint),
intent(in) :: ncols
1971 integer(kind=kint),
intent(in) :: cols(
cncol_item,ncols)
1972 integer(kind=kint),
intent(in) :: n_add_node
1973 integer(kind=kint),
intent(in) :: add_nodes(
cncol_item,n_add_node)
1974 integer(kind=kint),
intent(in) :: i0
1975 integer(kind=kint),
intent(inout) :: map(ncols)
1976 integer(kind=kint) :: i, j
1978 if (map(i) > 0) cycle
1979 call search_add_nodes(n_add_node, add_nodes, cols(crank,i), cols(clid,i), j)
1980 if (j == -1) stop
'ERROR: map_additional_nodes'
1983 end subroutine map_additional_nodes
1985 subroutine update_comm_table(hecMESHnew, n_add_node, add_nodes, i0)
1989 integer(kind=kint),
intent(in) :: n_add_node
1990 integer(kind=kint),
allocatable,
intent(inout) :: add_nodes(:,:)
1991 integer(kind=kint),
intent(in) :: i0
1992 integer(kind=kint),
allocatable :: n_add_imp(:), add_imp_index(:)
1993 integer(kind=kint),
allocatable :: add_imp_item_remote(:), add_imp_item_local(:)
1994 integer(kind=kint),
allocatable :: n_add_exp(:), add_exp_index(:), add_exp_item(:)
1995 integer(kind=kint),
allocatable :: n_new_imp(:), n_new_exp(:)
1996 integer(kind=kint) :: npe, nnb, comm, new_nnb
1997 integer(kind=kint),
pointer :: nbpe(:), new_nbpe(:)
1998 integer(kind=kint),
pointer :: import_index(:), export_index(:), import_item(:), export_item(:)
1999 integer(kind=kint),
pointer :: new_import_index(:), new_export_index(:)
2000 integer(kind=kint),
pointer :: new_import_item(:), new_export_item(:)
2001 npe = hecmeshnew%PETOT
2002 nnb = hecmeshnew%n_neighbor_pe
2003 comm = hecmeshnew%MPI_COMM
2004 nbpe => hecmeshnew%neighbor_pe
2005 import_index => hecmeshnew%import_index
2006 export_index => hecmeshnew%export_index
2007 import_item => hecmeshnew%import_item
2008 export_item => hecmeshnew%export_item
2010 call count_add_imp_per_rank(n_add_node, add_nodes, npe, n_add_imp)
2011 if (debug >= 3)
write(0,*)
' DEBUG3: count add_imp per rank done'
2013 allocate(add_imp_index(0:npe))
2014 call make_index(npe, n_add_imp, add_imp_index)
2015 if (debug >= 3)
write(0,*)
' DEBUG3: make add_imp_index done'
2017 call make_add_imp_item(n_add_node, add_nodes, npe, i0, add_imp_index, &
2018 add_imp_item_remote, add_imp_item_local)
2019 if (debug >= 3)
write(0,*)
' DEBUG3: make add_imp_item done'
2021 deallocate(add_nodes)
2025 allocate(n_add_exp(npe))
2027 if (debug >= 3)
write(0,*)
' DEBUG3: alltoall n_add_imp to n_add_exp done'
2029 allocate(add_exp_index(0:npe))
2030 call make_index(npe, n_add_exp, add_exp_index)
2031 if (debug >= 3)
write(0,*)
' DEBUG3: make add_exp_index done'
2033 call send_recv_add_imp_exp_item(npe, add_imp_index, add_imp_item_remote, &
2034 add_exp_index, add_exp_item, comm)
2035 if (debug >= 3)
write(0,*)
' DEBUG3: send recv add_imp/exp_item done'
2039 call count_new_comm_nodes(npe, nnb, nbpe, import_index, n_add_imp, n_new_imp)
2040 if (debug >= 3)
write(0,*)
' DEBUG3: count new comm_nodes (import) done'
2044 call count_new_comm_nodes(npe, nnb, nbpe, export_index, n_add_exp, n_new_exp)
2045 if (debug >= 3)
write(0,*)
' DEBUG3: count new comm_nodes (export) done'
2047 call update_neighbor_pe(npe, n_new_imp, n_new_exp, new_nnb, new_nbpe)
2048 if (debug >= 3)
write(0,*)
' DEBUG3: update neighbor_pe done'
2052 call merge_comm_table(npe, nnb, nbpe, import_index, import_item, &
2053 new_nnb, new_nbpe, add_imp_index, add_imp_item_local, n_add_imp, n_new_imp, &
2054 new_import_index, new_import_item)
2055 if (debug >= 3)
write(0,*)
' DEBUG3: merge comm_table (import) done'
2057 deallocate(n_add_imp)
2058 deallocate(add_imp_index)
2059 deallocate(add_imp_item_remote, add_imp_item_local)
2060 deallocate(n_new_imp)
2064 call merge_comm_table(npe, nnb, nbpe, export_index, export_item, &
2065 new_nnb, new_nbpe, add_exp_index, add_exp_item, n_add_exp, n_new_exp, &
2066 new_export_index, new_export_item)
2067 if (debug >= 3)
write(0,*)
' DEBUG3: merge comm_table (export) done'
2069 deallocate(n_add_exp)
2070 deallocate(add_exp_index)
2071 deallocate(add_exp_item)
2072 deallocate(n_new_exp)
2075 deallocate(import_index,import_item)
2076 deallocate(export_index,export_item)
2077 hecmeshnew%n_neighbor_pe = new_nnb
2078 hecmeshnew%neighbor_pe => new_nbpe
2079 hecmeshnew%import_index => new_import_index
2080 hecmeshnew%export_index => new_export_index
2081 hecmeshnew%import_item => new_import_item
2082 hecmeshnew%export_item => new_export_item
2083 end subroutine update_comm_table
2085 subroutine count_add_imp_per_rank(n_add_node, add_nodes, npe, n_add_imp)
2087 integer(kind=kint),
intent(in) :: n_add_node
2088 integer(kind=kint),
intent(in) :: add_nodes(
cncol_item,n_add_node)
2089 integer(kind=kint),
intent(in) :: npe
2090 integer(kind=kint),
allocatable,
intent(out) :: n_add_imp(:)
2091 integer(kind=kint) :: i, rank
2092 allocate(n_add_imp(npe))
2094 do i = 1, n_add_node
2095 rank = add_nodes(crank,i)
2096 n_add_imp(rank+1) = n_add_imp(rank+1) + 1
2098 end subroutine count_add_imp_per_rank
2100 subroutine make_add_imp_item(n_add_node, add_nodes, npe, i0, add_imp_index, &
2101 add_imp_item_remote, add_imp_item_local)
2103 integer(kind=kint),
intent(in) :: n_add_node
2104 integer(kind=kint),
intent(in) :: add_nodes(
cncol_item,n_add_node)
2105 integer(kind=kint),
intent(in) :: npe, i0
2106 integer(kind=kint),
allocatable,
intent(in) :: add_imp_index(:)
2107 integer(kind=kint),
allocatable,
intent(out) :: add_imp_item_remote(:), add_imp_item_local(:)
2108 integer(kind=kint),
allocatable :: cnt(:)
2109 integer(kind=kint) :: i, lid, rank, ipe
2110 allocate(add_imp_item_remote(add_imp_index(npe)))
2111 allocate(add_imp_item_local(add_imp_index(npe)))
2114 do i = 1, n_add_node
2115 lid = add_nodes(clid,i)
2116 rank = add_nodes(crank,i)
2118 cnt(ipe) = cnt(ipe) + 1
2119 add_imp_item_remote(add_imp_index(ipe-1) + cnt(ipe)) = lid
2120 add_imp_item_local(add_imp_index(ipe-1) + cnt(ipe)) = i0 + i
2123 end subroutine make_add_imp_item
2125 subroutine send_recv_add_imp_exp_item(npe, add_imp_index, add_imp_item_remote, &
2126 add_exp_index, add_exp_item, mpi_comm)
2129 integer(kind=kint),
intent(in) :: npe
2130 integer(kind=kint),
allocatable,
intent(in) :: add_imp_index(:), add_imp_item_remote(:)
2131 integer(kind=kint),
allocatable,
intent(in) :: add_exp_index(:)
2132 integer(kind=kint),
allocatable,
intent(out) :: add_exp_item(:)
2133 integer(kind=kint),
intent(in) :: mpi_comm
2134 integer(kind=kint) :: n_send, i, irank, is, ie, len, tag
2135 integer(kind=kint),
allocatable :: requests(:)
2136 integer(kind=kint),
allocatable :: statuses(:,:)
2137 allocate(add_exp_item(add_exp_index(npe)))
2138 allocate(requests(npe))
2143 is = add_imp_index(i-1)+1
2144 ie = add_imp_index(i)
2150 mpi_comm, requests(n_send))
2155 is = add_exp_index(i-1)+1
2156 ie = add_exp_index(i)
2161 mpi_comm, statuses(:,1))
2164 end subroutine send_recv_add_imp_exp_item
2166 subroutine count_new_comm_nodes(npe, org_nnb, org_nbpe, org_index, n_add, n_new)
2168 integer(kind=kint),
intent(in) :: npe, org_nnb
2170 integer(kind=kint),
pointer,
intent(in) :: org_nbpe(:), org_index(:)
2171 integer(kind=kint),
intent(in) :: n_add(:)
2172 integer(kind=kint),
allocatable,
intent(out) :: n_new(:)
2173 integer(kind=kint) :: i, irank, n_org
2174 allocate(n_new(npe))
2178 n_org = org_index(i) - org_index(i-1)
2179 n_new(irank+1) = n_new(irank+1) + n_org
2181 end subroutine count_new_comm_nodes
2183 subroutine update_neighbor_pe(npe, n_new_imp, n_new_exp, &
2186 integer(kind=kint),
intent(in) :: npe
2187 integer(kind=kint),
intent(in) :: n_new_imp(npe), n_new_exp(npe)
2188 integer(kind=kint),
intent(out) :: new_nnb
2189 integer(kind=kint),
pointer,
intent(out) :: new_nbpe(:)
2190 integer(kind=kint) :: i
2193 if (n_new_imp(i) > 0 .or. n_new_exp(i) > 0) new_nnb = new_nnb+1
2195 allocate(new_nbpe(new_nnb))
2198 if (n_new_imp(i) > 0 .or. n_new_exp(i) > 0)
then
2200 new_nbpe(new_nnb) = i-1
2203 end subroutine update_neighbor_pe
2205 subroutine merge_comm_table(npe, org_nnb, org_nbpe, org_index, org_item, &
2206 new_nnb, new_nbpe, add_index, add_item, n_add, n_new, new_index, new_item)
2208 integer(kind=kint),
intent(in) :: npe, org_nnb
2210 integer(kind=kint),
pointer,
intent(in) :: org_nbpe(:), org_index(:), org_item(:)
2211 integer(kind=kint),
intent(in) :: new_nnb
2213 integer(kind=kint),
pointer,
intent(in) :: new_nbpe(:)
2214 integer(kind=kint),
allocatable,
intent(in) :: add_index(:), add_item(:)
2215 integer(kind=kint),
intent(in) :: n_add(npe), n_new(npe)
2216 integer(kind=kint),
pointer,
intent(out) :: new_index(:), new_item(:)
2217 integer(kind=kint),
allocatable :: cnt(:)
2218 integer(kind=kint) :: i, irank, j, jrank, i0, j0, len
2221 allocate(new_index(0:new_nnb))
2225 new_index(i) = new_index(i-1) + n_new(irank+1)
2227 allocate(new_item(new_index(new_nnb)))
2233 if (org_index(i) - org_index(i-1) == 0) cycle
2235 do while (jrank < irank)
2237 if (j > new_nnb)
exit
2240 if (jrank /= irank) stop
'ERROR: merging comm table: org into new'
2242 len = org_index(i) - i0
2244 new_item(j0+1:j0+len) = org_item(i0+1:i0+len)
2250 if (n_add(i) == 0) cycle
2252 do while (jrank < irank)
2256 if (jrank /= irank) stop
'ERROR: merging comm table: add into new'
2258 len = add_index(i) - i0
2259 j0 = new_index(j-1) + cnt(jrank+1)
2260 new_item(j0+1:j0+len) = add_item(i0+1:i0+len)
2261 cnt(jrank+1) = cnt(jrank+1) + len
2262 if (cnt(jrank+1) /= new_index(j)-new_index(j-1)) stop
'ERROR: merging comm table'
2265 end subroutine merge_comm_table
2267 subroutine copy_vals_to_bt_int(nnb, imp_rows_index, imp_cols_index, &
2268 imp_rows_item, map, ndof2, imp_vals_item, BT_int)
2270 integer(kind=kint),
intent(in) :: nnb
2271 integer(kind=kint),
allocatable,
intent(in) :: imp_rows_index(:), imp_cols_index(:)
2272 integer(kind=kint),
intent(in) :: imp_rows_item(:,:), map(:)
2273 integer(kind=kint),
intent(in) :: ndof2
2274 real(kind=
kreal),
intent(in) :: imp_vals_item(:)
2276 integer(kind=kint),
allocatable :: cnt(:)
2277 integer(kind=kint) :: idom, is, ie, ic0, i, irow, ncol, j0, j
2278 allocate(cnt(bt_int%nr))
2281 is = imp_rows_index(idom-1)+1
2282 ie = imp_rows_index(idom)
2283 ic0 = imp_cols_index(idom-1)
2285 irow = imp_rows_item(1,i)
2286 ncol = imp_rows_item(2,i)
2287 if (irow < 1 .or. bt_int%nr < irow) stop
'ERROR: copy vals to BT_int: irow'
2288 j0 = bt_int%index(irow-1) + cnt(irow)
2290 bt_int%item(j0+j) = map(ic0+j)
2291 bt_int%A(ndof2*(j0+j-1)+1:ndof2*(j0+j)) = imp_vals_item(ndof2*(ic0+j-1)+1:ndof2*(ic0+j))
2293 cnt(irow) = cnt(irow) + ncol
2296 if (ic0 /= imp_cols_index(idom)) stop
'ERROR: copy vals to BT_int: ic0'
2299 end subroutine copy_vals_to_bt_int
2301 subroutine sort_and_uniq_rows(BTmat)
2305 integer(kind=kint) :: nr, ndof, ndof2
2306 integer(kind=kint) :: irow, is, ie, is_new, ie_new, i, i_new
2307 integer(kind=kint) :: ndup, ndup_tot
2308 integer(kind=kint) :: js, je, js_new, je_new
2309 integer(kind=kint) :: new_nnz
2310 integer(kind=kint),
allocatable :: cnt(:)
2311 integer(kind=kint),
pointer :: sort_item(:), new_index(:), new_item(:)
2312 real(kind=
kreal),
pointer :: new_a(:)
2314 real(kind=
kreal) :: t0, t1
2319 outer:
do irow = 1, nr
2320 is = btmat%index(irow-1)+1
2321 ie = btmat%index(irow)
2323 if (btmat%item(i) >= btmat%item(i+1))
then
2330 if (timer >= 4)
write(0,
'(A,f10.4,L2)')
"####### sort_and_uniq_rows (1) : ",t1-t0,sorted
2337 allocate(sort_item(btmat%nnz))
2339 sort_item(i) = btmat%item(i)
2350 is = btmat%index(irow-1)+1
2351 ie = btmat%index(irow)
2354 cnt(irow) = (ie-is+1) - ndup
2355 ndup_tot = ndup_tot + ndup
2359 if (timer >= 4)
write(0,
'(A,f10.4,I5)')
"####### sort_and_uniq_rows (2) : ",t1-t0,ndup_tot
2362 if (ndup_tot == 0)
then
2363 new_index => btmat%index
2365 new_item => sort_item
2367 allocate(new_index(0:nr))
2368 call make_index(nr, cnt, new_index)
2369 new_nnz = new_index(nr)
2370 allocate(new_item(new_nnz))
2372 is = btmat%index(irow-1)+1
2374 is_new = new_index(irow-1)+1
2375 ie_new = is_new+cnt(irow)-1
2376 new_item(is_new:ie_new) = sort_item(is:ie)
2378 deallocate(sort_item)
2382 if (timer >= 4)
write(0,
'(A,f10.4)')
"####### sort_and_uniq_rows (3) : ",t1-t0
2385 allocate(new_a(ndof2*new_nnz))
2393 is = btmat%index(irow-1)+1
2394 ie = btmat%index(irow)
2395 is_new = new_index(irow-1)+1
2396 ie_new = new_index(irow)
2401 if (i_new == -1) stop
'ERROR: sort_and_uniq_rows'
2404 js_new = ndof2*(i_new-1)+1
2405 je_new = ndof2*i_new
2406 new_a(js_new:je_new) = new_a(js_new:je_new) + btmat%A(js:je)
2411 if (timer >= 4)
write(0,
'(A,f10.4)')
"####### sort_and_uniq_rows (4) : ",t1-t0
2414 if (ndup_tot == 0)
then
2415 deallocate(btmat%item)
2416 btmat%item => new_item
2421 deallocate(btmat%index)
2422 btmat%index => new_index
2423 deallocate(btmat%item)
2424 btmat%item => new_item
2428 end subroutine sort_and_uniq_rows
2435 integer(kind=kint) :: ndof, ndof2, nr, nc, i, icnt, js, je, j, jcol, idx, i0, k
2436 integer(kind=kint),
allocatable :: iw(:)
2437 if (amat%ndof /= bmat%ndof) stop
'ERROR: hecmw_localmat_add: non-matching ndof'
2440 nr = min(amat%nr, bmat%nr)
2441 nc = max(amat%nc, bmat%nc)
2446 allocate(cmat%index(0:nr))
2452 js = amat%index(i-1)+1
2460 js = bmat%index(i-1)+1
2465 if (iw(k) == jcol) cycle lj1
2470 cmat%index(i) = cmat%index(i-1) + icnt
2472 cmat%nnz = cmat%index(nr)
2473 allocate(cmat%item(cmat%nnz))
2474 allocate(cmat%A(ndof2*cmat%nnz))
2476 i0 = cmat%index(i-1)
2479 js = amat%index(i-1)+1
2485 cmat%item(idx) = jcol
2486 cmat%A(ndof2*(idx-1)+1:ndof2*idx) = amat%A(ndof2*(j-1)+1:ndof2*j)
2489 js = bmat%index(i-1)+1
2495 if (cmat%item(idx) == jcol)
then
2496 cmat%A(ndof2*(idx-1)+1:ndof2*idx) = &
2497 cmat%A(ndof2*(idx-1)+1:ndof2*idx) + bmat%A(ndof2*(j-1)+1:ndof2*j)
2503 cmat%item(idx) = jcol
2504 cmat%A(ndof2*(idx-1)+1:ndof2*idx) = bmat%A(ndof2*(j-1)+1:ndof2*j)
2506 if (i0 + icnt /= cmat%index(i)) stop
'ERROR: merge localmat'
2508 call sort_and_uniq_rows(cmat)
2558 integer(kind=kint),
optional,
intent(in) :: num_lagrange
2559 integer(kind=kint) :: ndof, ndof2, i, idx, idx2, js, je, j, k
2560 integer(kind=kint),
allocatable :: incl_nz(:), cnt(:)
2561 logical :: check_nonzero
2563 check_nonzero = .true.
2568 bkmat%nr = hecmat%NP
2569 bkmat%nc = hecmat%NP
2572 if (
present(num_lagrange))
then
2573 if (num_lagrange == 0)
then
2575 allocate(bkmat%index(0:bkmat%nr))
2577 bkmat%item => null()
2581 check_nonzero = .true.
2584 if (check_nonzero)
then
2585 allocate(incl_nz(hecmat%NPL + hecmat%NPU + hecmat%NP))
2586 allocate(cnt(bkmat%nr))
2593 idx = hecmat%indexL(i-1) + (i-1) + hecmat%indexU(i-1)
2596 js = hecmat%indexL(i-1)+1
2597 je = hecmat%indexL(i)
2601 if (hecmat%AL(ndof2*(j-1)+k) /= 0.0d0)
then
2611 if (hecmat%D(ndof2*(i-1)+k) /= 0.0d0)
then
2618 js = hecmat%indexU(i-1)+1
2619 je = hecmat%indexU(i)
2623 if (hecmat%AU(ndof2*(j-1)+k) /= 0.0d0)
then
2630 if (idx /= hecmat%indexL(i) + i + hecmat%indexU(i)) stop
'ERROR: hecmw_localmat_init_with_hecmat: count'
2635 allocate(bkmat%index(0:bkmat%nr))
2636 call make_index(bkmat%nr, cnt, bkmat%index)
2638 bkmat%nnz = bkmat%index(bkmat%nr)
2640 allocate(bkmat%item(bkmat%nnz))
2641 allocate(bkmat%A(ndof2 * bkmat%nnz))
2647 idx = hecmat%indexL(i-1) + (i-1) + hecmat%indexU(i-1)
2648 idx2 = bkmat%index(i-1)
2650 js = hecmat%indexL(i-1)+1
2651 je = hecmat%indexL(i)
2654 if (incl_nz(idx) == 1)
then
2656 bkmat%item(idx2) = hecmat%itemL(j)
2657 bkmat%A(ndof2*(idx2-1)+1:ndof2*idx2) = hecmat%AL(ndof2*(j-1)+1:ndof2*j)
2662 if (incl_nz(idx) == 1)
then
2664 bkmat%item(idx2) = i
2665 bkmat%A(ndof2*(idx2-1)+1:ndof2*idx2) = hecmat%D(ndof2*(i-1)+1:ndof2*i)
2668 js = hecmat%indexU(i-1)+1
2669 je = hecmat%indexU(i)
2672 if (incl_nz(idx) == 1)
then
2674 bkmat%item(idx2) = hecmat%itemU(j)
2675 bkmat%A(ndof2*(idx2-1)+1:ndof2*idx2) = hecmat%AU(ndof2*(j-1)+1:ndof2*j)
2678 if (idx /= hecmat%indexL(i) + i + hecmat%indexU(i)) stop
'ERROR: hecmw_localmat_init_with_hecmat: copy'
2679 if (idx2 /= bkmat%index(i)) stop
'ERROR: hecmw_localmat_init_with_hecmat: index'
2685 bkmat%nnz = hecmat%NPL + hecmat%NP + hecmat%NPU
2686 allocate(bkmat%index(0:bkmat%nr))
2687 allocate(bkmat%item(bkmat%nnz))
2688 allocate(bkmat%A(ndof2 * bkmat%nnz))
2694 idx = hecmat%indexL(i-1) + (i-1) + hecmat%indexU(i-1)
2696 js = hecmat%indexL(i-1)+1
2697 je = hecmat%indexL(i)
2700 bkmat%item(idx) = hecmat%itemL(j)
2701 bkmat%A(ndof2*(idx-1)+1:ndof2*idx) = hecmat%AL(ndof2*(j-1)+1:ndof2*j)
2706 bkmat%A(ndof2*(idx-1)+1:ndof2*idx) = hecmat%D(ndof2*(i-1)+1:ndof2*i)
2708 js = hecmat%indexU(i-1)+1
2709 je = hecmat%indexU(i)
2712 bkmat%item(idx) = hecmat%itemU(j)
2713 bkmat%A(ndof2*(idx-1)+1:ndof2*idx) = hecmat%AU(ndof2*(j-1)+1:ndof2*j)
2715 bkmat%index(i) = idx
2716 if (idx /= hecmat%indexL(i) + i + hecmat%indexU(i)) stop
'ERROR: hecmw_localmat_init_with_hecmat: copy'
2729 if (debug >= 3)
then
2731 if (debug == 3)
then
2742 bkmat%nnz = w2mat%nnz
2743 bkmat%ndof = w2mat%ndof
2744 bkmat%index => w2mat%index
2745 bkmat%item => w2mat%item
2759 integer(kind=kint),
allocatable :: exp_cols_index(:)
2760 integer(kind=kint),
allocatable :: exp_cols_item(:,:)
2761 real(kind=
kreal) :: t0, t1
2764 call make_comm_table(bkmat, hecmesh, heccomm)
2765 if (debug >= 1)
write(0,*)
'DEBUG: hecmw_localmat_multmat: make_comm_table done'
2767 if (timer >= 2)
write(0,
'(A,f10.4)')
'##### hecmw_localmat_multmat (1) : ',t1-t0
2770 if (btmat%nr > hecmesh%nn_internal)
then
2772 if (debug >= 1)
write(0,
'(A)')
'DEBUG: hecmw_localmat_multmat: ignore external part of BTmat'
2773 btmat%nr = hecmesh%nn_internal
2774 btmat%nnz = btmat%index(btmat%nr)
2777 call extract_bt_exp(btmat, heccomm, bt_exp)
2778 if (debug >= 1)
write(0,*)
'DEBUG: hecmw_localmat_multmat: extract_BT_exp done'
2780 if (timer >= 2)
write(0,
'(A,f10.4)')
'##### hecmw_localmat_multmat (2) : ',t1-t0
2783 call prepare_column_info(hecmesh, bt_exp, exp_cols_index, exp_cols_item)
2784 if (debug >= 1)
write(0,*)
'DEBUG: hecmw_localmat_multmat: prepare column info done'
2786 if (timer >= 2)
write(0,
'(A,f10.4)')
'##### hecmw_localmat_multmat (3) : ',t1-t0
2789 call send_bt_exp_and_recv_bt_imp(hecmesh, heccomm, bt_exp, exp_cols_index, exp_cols_item, bt_imp, hecmeshnew)
2790 if (debug >= 1)
write(0,*)
'DEBUG: hecmw_localmat_multmat: send BT_exp and recv BT_imp done'
2792 if (timer >= 2)
write(0,
'(A,f10.4)')
'##### hecmw_localmat_multmat (4) : ',t1-t0
2794 call free_comm_table(heccomm)
2796 call concat_btmat_and_bt_imp(btmat, bt_imp, bt_all)
2797 if (debug >= 1)
write(0,*)
'DEBUG: hecmw_localmat_multmat: concat BTmat and BT_imp into BT_all done'
2799 if (timer >= 2)
write(0,
'(A,f10.4)')
'##### hecmw_localmat_multmat (5) : ',t1-t0
2803 call multiply_mat_mat(bkmat, bt_all, bktmat)
2804 if (debug >= 1)
write(0,*)
'DEBUG: hecmw_localmat_multmat: multiply BKmat and BT_all into BKTmat done'
2806 if (timer >= 2)
write(0,
'(A,f10.4)')
'##### hecmw_localmat_multmat (6) : ',t1-t0
2810 if (hecmesh%n_neighbor_pe > 0)
then
2811 hecmesh%n_node = hecmeshnew%n_node
2812 hecmesh%n_neighbor_pe = hecmeshnew%n_neighbor_pe
2813 deallocate(hecmesh%neighbor_pe)
2814 deallocate(hecmesh%import_index)
2815 deallocate(hecmesh%export_index)
2816 deallocate(hecmesh%import_item)
2817 deallocate(hecmesh%export_item)
2818 deallocate(hecmesh%node_ID)
2819 deallocate(hecmesh%global_node_ID)
2820 hecmesh%neighbor_pe => hecmeshnew%neighbor_pe
2821 hecmesh%import_index => hecmeshnew%import_index
2822 hecmesh%export_index => hecmeshnew%export_index
2823 hecmesh%import_item => hecmeshnew%import_item
2824 hecmesh%export_item => hecmeshnew%export_item
2825 hecmesh%node_ID => hecmeshnew%node_ID
2826 hecmesh%global_node_ID => hecmeshnew%global_node_ID
2827 if (debug >= 1)
write(0,*)
'DEBUG: hecmw_localmat_multmat: update hecMESH done'
2829 if (timer >= 2)
write(0,
'(A,f10.4)')
'##### hecmw_localmat_multmat (7) : ',t1-t0
2833 subroutine make_comm_table(BKmat, hecMESH, hecCOMM)
2839 integer(kind=kint) :: nn_int, nn_ext, nnb, i, icol, irank, idom, idx, n_send, tag, js, je, len
2840 integer(kind=kint),
allocatable :: is_nz_col(:), imp_cnt(:), exp_cnt(:), import_item_remote(:)
2841 integer(kind=kint),
allocatable :: requests(:), statuses(:,:)
2842 heccomm%zero = hecmesh%zero
2843 heccomm%HECMW_COMM = hecmesh%MPI_COMM
2844 heccomm%PETOT = hecmesh%PETOT
2845 heccomm%PEsmpTOT = hecmesh%PEsmpTOT
2846 heccomm%my_rank = hecmesh%my_rank
2847 heccomm%errnof = hecmesh%errnof
2848 heccomm%n_subdomain = hecmesh%n_subdomain
2849 heccomm%n_neighbor_pe = hecmesh%n_neighbor_pe
2850 allocate(heccomm%neighbor_pe(heccomm%n_neighbor_pe))
2851 heccomm%neighbor_pe(:) = hecmesh%neighbor_pe(:)
2853 nn_int = hecmesh%nn_internal
2854 nn_ext = hecmesh%n_node - hecmesh%nn_internal
2855 nnb = heccomm%n_neighbor_pe
2858 allocate(is_nz_col(nn_ext))
2860 do i = 1, bkmat%index(nn_int)
2861 icol = bkmat%item(i)
2862 if (icol > nn_int) is_nz_col(icol - nn_int) = 1
2866 allocate(imp_cnt(nnb))
2869 if (is_nz_col(i) == 1)
then
2870 irank = hecmesh%node_ID(2*(nn_int+i))
2871 call rank_to_idom(hecmesh, irank, idom)
2872 imp_cnt(idom) = imp_cnt(idom) + 1
2875 if (debug >= 3)
write(0,*)
' DEBUG3: imp_cnt',imp_cnt(:)
2878 allocate(heccomm%import_index(0:nnb))
2879 call make_index(nnb, imp_cnt, heccomm%import_index)
2880 if (debug >= 3)
write(0,*)
' DEBUG3: import_index',heccomm%import_index(:)
2883 allocate(heccomm%import_item(heccomm%import_index(nnb)))
2886 if (is_nz_col(i) == 1)
then
2887 irank = hecmesh%node_ID(2*(nn_int+i))
2888 call rank_to_idom(hecmesh, irank, idom)
2889 imp_cnt(idom) = imp_cnt(idom) + 1
2890 idx = heccomm%import_index(idom-1)+imp_cnt(idom)
2891 heccomm%import_item(idx) = nn_int+i
2894 if (debug >= 3)
write(0,*)
' DEBUG3: import_item',heccomm%import_item(:)
2896 allocate(import_item_remote(heccomm%import_index(nnb)))
2897 do i = 1, heccomm%import_index(nnb)
2898 import_item_remote(i) = hecmesh%node_ID(2*heccomm%import_item(i)-1)
2900 if (debug >= 3)
write(0,*)
' DEBUG3: import_item_remote',import_item_remote(:)
2902 allocate(requests(2*nnb))
2908 irank = heccomm%neighbor_pe(idom)
2911 call hecmw_isend_int(imp_cnt(idom), 1, irank, tag, heccomm%HECMW_COMM, requests(n_send))
2912 if (imp_cnt(idom) > 0)
then
2913 js = heccomm%import_index(idom-1)+1
2914 je = heccomm%import_index(idom)
2919 heccomm%HECMW_COMM, requests(n_send))
2924 allocate(exp_cnt(nnb))
2926 irank = heccomm%neighbor_pe(idom)
2928 call hecmw_recv_int(exp_cnt(idom), 1, irank, tag, heccomm%HECMW_COMM, statuses(:,1))
2930 allocate(heccomm%export_index(0:nnb))
2931 call make_index(nnb, exp_cnt, heccomm%export_index)
2932 if (debug >= 3)
write(0,*)
' DEBUG3: export_index',heccomm%export_index(:)
2935 allocate(heccomm%export_item(heccomm%export_index(nnb)))
2937 if (exp_cnt(idom) <= 0) cycle
2938 irank = heccomm%neighbor_pe(idom)
2939 js = heccomm%export_index(idom-1)+1
2940 je = heccomm%export_index(idom)
2943 call hecmw_recv_int(heccomm%export_item(js:je), len, irank, tag, &
2944 heccomm%HECMW_COMM, statuses(:,1))
2946 if (debug >= 3)
write(0,*)
' DEBUG3: export_item',heccomm%export_item(:)
2951 deallocate(import_item_remote)
2952 end subroutine make_comm_table
2954 subroutine free_comm_table(hecCOMM)
2957 deallocate(heccomm%neighbor_pe)
2958 deallocate(heccomm%import_index)
2959 deallocate(heccomm%import_item)
2960 deallocate(heccomm%export_index)
2961 deallocate(heccomm%export_item)
2962 end subroutine free_comm_table
2964 subroutine extract_bt_exp(BTmat, hecCOMM, BT_exp)
2969 integer(kind=kint) :: ndof, ndof2, idom, idx_0, idx_n, j, jrow, nnz_row, idx, ks, ke, k
2970 if (heccomm%n_neighbor_pe == 0)
return
2971 allocate(bt_exp(heccomm%n_neighbor_pe))
2974 do idom = 1, heccomm%n_neighbor_pe
2975 idx_0 = heccomm%export_index(idom-1)
2976 idx_n = heccomm%export_index(idom)
2977 bt_exp(idom)%nr = idx_n - idx_0
2978 bt_exp(idom)%nc = btmat%nc
2979 bt_exp(idom)%nnz = 0
2980 bt_exp(idom)%ndof = ndof
2981 allocate(bt_exp(idom)%index(0:bt_exp(idom)%nr))
2982 bt_exp(idom)%index(0) = 0
2983 do j = 1, bt_exp(idom)%nr
2984 jrow = heccomm%export_item(idx_0 + j)
2985 nnz_row = btmat%index(jrow) - btmat%index(jrow-1)
2986 bt_exp(idom)%index(j) = bt_exp(idom)%index(j-1) + nnz_row
2988 bt_exp(idom)%nnz = bt_exp(idom)%index(bt_exp(idom)%nr)
2989 allocate(bt_exp(idom)%item(bt_exp(idom)%nnz))
2990 allocate(bt_exp(idom)%A(ndof2 * bt_exp(idom)%nnz))
2992 do j = 1, bt_exp(idom)%nr
2993 jrow = heccomm%export_item(idx_0 + j)
2994 ks = btmat%index(jrow-1) + 1
2995 ke = btmat%index(jrow)
2998 bt_exp(idom)%item(idx) = btmat%item(k)
2999 bt_exp(idom)%A(ndof2*(idx-1)+1:ndof2*idx) = btmat%A(ndof2*(k-1)+1:ndof2*k)
3001 if (idx /= bt_exp(idom)%index(j)) stop
'ERROR: extract BT_exp'
3004 end subroutine extract_bt_exp
3006 subroutine send_bt_exp_and_recv_bt_imp(hecMESH, hecCOMM, BT_exp, exp_cols_index, exp_cols_item, BT_imp, hecMESHnew)
3012 integer(kind=kint),
allocatable,
intent(inout) :: exp_cols_index(:)
3013 integer(kind=kint),
allocatable,
intent(inout) :: exp_cols_item(:,:)
3016 integer(kind=kint),
allocatable :: nnz_imp(:), cnt(:), index_imp(:)
3017 integer(kind=kint),
allocatable :: imp_cols_index(:)
3018 integer(kind=kint),
allocatable :: imp_cols_item(:,:)
3019 real(kind=
kreal),
allocatable :: imp_vals_item(:)
3020 integer(kind=kint) :: nnb, ndof, ndof2, idom, irank, nr, n_send, tag, idx_0, idx_n, j, jj, nnz
3021 integer(kind=kint),
allocatable :: requests(:)
3022 integer(kind=kint),
allocatable :: statuses(:,:)
3023 integer(kind=kint),
allocatable :: map(:), add_nodes(:,:)
3024 integer(kind=kint) :: n_add_node, i0
3025 nnb = heccomm%n_neighbor_pe
3031 allocate(bt_imp%index(0:0))
3035 ndof = bt_exp(1)%ndof
3037 allocate(requests(nnb*3))
3041 irank = heccomm%neighbor_pe(idom)
3042 nr = bt_exp(idom)%nr
3046 call hecmw_isend_int(bt_exp(idom)%index(0:bt_exp(idom)%nr), bt_exp(idom)%nr + 1, &
3047 irank, tag, heccomm%HECMW_COMM, requests(n_send))
3048 if (bt_exp(idom)%nnz == 0) cycle
3052 cncol_item * bt_exp(idom)%nnz, irank, tag, heccomm%HECMW_COMM, requests(n_send))
3055 call hecmw_isend_r(bt_exp(idom)%A, ndof2 * bt_exp(idom)%nnz, &
3056 irank, tag, heccomm%HECMW_COMM, requests(n_send))
3060 bt_imp%nr = hecmesh%n_node - hecmesh%nn_internal
3065 allocate(nnz_imp(nnb))
3066 allocate(cnt(bt_imp%nr))
3070 irank = heccomm%neighbor_pe(idom)
3071 idx_0 = heccomm%import_index(idom-1)
3072 idx_n = heccomm%import_index(idom)
3078 allocate(index_imp(0:nr))
3081 heccomm%HECMW_COMM, statuses(:,1))
3082 nnz_imp(idom) = index_imp(nr)
3084 jj = heccomm%import_item(idx_0 + j) - hecmesh%nn_internal
3085 if (jj < 1 .or. bt_imp%nr < jj) stop
'ERROR: jj out of range'
3086 if (cnt(jj) /= 0) stop
'ERROR: duplicate import rows?'
3087 cnt(jj) = index_imp(j) - index_imp(j-1)
3089 deallocate(index_imp)
3092 allocate(imp_cols_index(0:nnb))
3093 call make_index(nnb, nnz_imp, imp_cols_index)
3096 allocate(bt_imp%index(0:bt_imp%nr))
3097 call make_index(bt_imp%nr, cnt, bt_imp%index)
3100 bt_imp%nnz = bt_imp%index(bt_imp%nr)
3101 if (bt_imp%nnz /= imp_cols_index(nnb)) &
3102 stop
'ERROR: total num of nonzero of BT_imp'
3104 allocate(imp_cols_item(
cncol_item, bt_imp%nnz))
3105 allocate(imp_vals_item(ndof2 * bt_imp%nnz))
3108 irank = heccomm%neighbor_pe(idom)
3109 idx_0 = imp_cols_index(idom-1)
3110 idx_n = imp_cols_index(idom)
3115 irank, tag, heccomm%HECMW_COMM, statuses(:,1))
3117 call hecmw_recv_r(imp_vals_item(ndof2*idx_0 + 1), ndof2 * nnz, &
3118 irank, tag, heccomm%HECMW_COMM, statuses(:,1))
3121 if (debug >= 2)
write(0,*)
' DEBUG2: send BT_imp and recv into temporary data done'
3123 deallocate(requests)
3124 deallocate(statuses)
3130 deallocate(exp_cols_index)
3131 deallocate(exp_cols_item)
3133 call copy_mesh(hecmesh, hecmeshnew)
3135 call map_imported_cols(hecmeshnew, imp_cols_index(nnb), imp_cols_item, n_add_node, add_nodes, map, i0)
3136 if (debug >= 2)
write(0,*)
' DEBUG2: map imported cols done'
3138 call update_comm_table(hecmeshnew, n_add_node, add_nodes, i0)
3139 if (debug >= 2)
write(0,*)
' DEBUG2: update comm_table done'
3141 bt_imp%nc = hecmeshnew%n_node
3143 allocate(bt_imp%item(bt_imp%nnz))
3144 allocate(bt_imp%A(ndof2 * bt_imp%nnz))
3145 call copy_vals_to_bt_imp(heccomm, hecmesh%nn_internal, imp_cols_index, map, imp_vals_item, bt_imp)
3146 if (debug >= 2)
write(0,*)
' DEBUG2: copy vals to BT_imp done'
3148 deallocate(imp_cols_index)
3149 deallocate(imp_cols_item)
3150 deallocate(imp_vals_item)
3152 end subroutine send_bt_exp_and_recv_bt_imp
3154 subroutine copy_vals_to_bt_imp(hecCOMM, nn_internal, imp_cols_index, map, imp_vals_item, BT_imp)
3157 integer(kind=kint),
intent(in) :: nn_internal
3158 integer(kind=kint),
allocatable,
intent(in) :: imp_cols_index(:)
3159 integer(kind=kint),
intent(in) :: map(:)
3160 real(kind=
kreal),
intent(in) :: imp_vals_item(:)
3162 integer(kind=kint) :: nnb, ndof2, idx, idom, idx_0, idx_n, nr, j, jrow, ks, ke, k
3163 nnb = heccomm%n_neighbor_pe
3164 ndof2 = bt_imp%ndof ** 2
3167 idx_0 = heccomm%import_index(idom-1)
3168 idx_n = heccomm%import_index(idom)
3172 jrow = heccomm%import_item(idx_0 + j) - nn_internal
3173 ks = bt_imp%index(jrow-1)+1
3174 ke = bt_imp%index(jrow)
3177 bt_imp%item(k) = map(idx)
3178 bt_imp%A(ndof2*(k-1)+1:ndof2*k) = imp_vals_item(ndof2*(idx-1)+1:ndof2*idx)
3181 if (idx /= imp_cols_index(idom)) stop
'ERROR: copy vals to BT_imp'
3183 end subroutine copy_vals_to_bt_imp
3185 subroutine concat_btmat_and_bt_imp(BTmat, BT_imp, BT_all)
3190 integer(kind=kint) :: ndof, ndof2, i, ii
3192 if (bt_imp%nr > 0 .and. bt_imp%ndof /= ndof) stop
'ERROR: concat BTmat and BT_imp: ndof'
3194 bt_all%nr = btmat%nr + bt_imp%nr
3195 bt_all%nc = max(btmat%nc, bt_imp%nc)
3196 bt_all%nnz = btmat%nnz + bt_imp%nnz
3198 allocate(bt_all%index(0:bt_all%nr))
3199 allocate(bt_all%item(bt_all%nnz))
3200 allocate(bt_all%A(ndof2 * bt_all%nnz))
3203 bt_all%index(i) = btmat%index(i)
3206 bt_all%index(btmat%nr+i) = bt_all%index(btmat%nr+i-1) + &
3207 bt_imp%index(i) - bt_imp%index(i-1)
3210 bt_all%item(i) = btmat%item(i)
3211 bt_all%A(ndof2*(i-1)+1:ndof2*i) = btmat%A(ndof2*(i-1)+1:ndof2*i)
3213 do i = 1, bt_imp%nnz
3215 bt_all%item(ii) = bt_imp%item(i)
3216 bt_all%A(ndof2*(ii-1)+1:ndof2*ii) = bt_imp%A(ndof2*(i-1)+1:ndof2*i)
3218 end subroutine concat_btmat_and_bt_imp
3220 subroutine multiply_mat_mat(Amat, Bmat, Cmat)
3225 integer(kind=kint) :: ndof, ndof2, nr, nc, nnz, i, icnt
3226 integer(kind=kint) :: js, je, j, jj, ks, ke, k, kk, l, ll, l0
3227 integer(kind=kint),
allocatable :: iw(:)
3228 real(kind=
kreal),
pointer :: ap(:), bp(:), cp(:)
3229 real(kind=
kreal) :: t0, t1
3231 if (amat%ndof /= bmat%ndof) stop
'ERROR: multiply_mat_mat: unmatching ndof'
3236 if (amat%nc /= bmat%nr)
then
3237 write(0,*)
'Amat: nr, nc = ', amat%nr, amat%nc
3238 write(0,*)
'Bmat: nr, nc = ', bmat%nr, bmat%nc
3239 stop
'ERROR: multiply_mat_mat: unmatching size'
3244 allocate(cmat%index(0:nr))
3253 js = amat%index(i-1)+1
3257 ks = bmat%index(jj-1)+1
3262 if (iw(l) == kk) cycle kl1
3268 cmat%index(i) = icnt
3274 cmat%index(i) = cmat%index(i-1) + cmat%index(i)
3276 nnz = cmat%index(nr)
3280 if (timer >= 3)
write(0,
'(A,f10.4)')
"###### multiply_mat_mat (1) : ",t1-t0
3282 allocate(cmat%item(nnz))
3283 allocate(cmat%A(ndof2 * nnz))
3291 l0 = cmat%index(i-1)
3293 js = amat%index(i-1)+1
3297 ap => amat%A(ndof2*(j-1)+1:ndof2*j)
3298 ks = bmat%index(jj-1)+1
3302 bp => bmat%A(ndof2*(k-1)+1:ndof2*k)
3305 if (cmat%item(l0+l) == kk)
then
3315 cp => cmat%A(ndof2*(ll-1)+1:ndof2*ll)
3316 call blk_matmul_add(ndof, ap, bp, cp)
3320 if (l0+icnt /= cmat%index(i)) stop
'ERROR: multiply_mat_mat: unknown error'
3325 if (timer >= 3)
write(0,
'(A,f10.4)')
"###### multiply_mat_mat (2) : ",t1-t0
3327 call sort_and_uniq_rows(cmat)
3329 if (timer >= 3)
write(0,
'(A,f10.4)')
"###### multiply_mat_mat (3) : ",t1-t0
3330 end subroutine multiply_mat_mat
3332 subroutine blk_matmul_add(ndof, A, B, AB)
3334 integer,
intent(in) :: ndof
3335 real(kind=
kreal),
intent(in) :: a(:), b(:)
3336 real(kind=
kreal),
intent(inout) :: ab(:)
3337 integer :: ndof2, i, j, k, i0, j0, ij, ik, jk
3348 ab(ik)=ab(ik)+a(ij)*b(jk)
3352 end subroutine blk_matmul_add
3359 call make_new_hecmat(hecmat, bttktmat, hectkt)
3367 call make_comm_table(bkmat, hecmesh, heccomm)
3368 deallocate(hecmesh%import_index)
3369 deallocate(hecmesh%import_item)
3370 deallocate(hecmesh%export_index)
3371 deallocate(hecmesh%export_item)
3372 hecmesh%import_index => heccomm%import_index
3373 hecmesh%import_item => heccomm%import_item
3374 hecmesh%export_index => heccomm%export_index
3375 hecmesh%export_item => heccomm%export_item
3376 deallocate(heccomm%neighbor_pe)
subroutine, public hecmw_bsearch_int_array(array, istart, iend, val, idx)
recursive subroutine, public hecmw_qsort_int_array(array, istart, iend)
subroutine, public hecmw_uniq_int_array(array, istart, iend, ndup)
subroutine, public hecmw_localmat_write(tmat, iunit)
subroutine, public hecmw_localmat_add(amat, bmat, cmat)
subroutine, public hecmw_localmat_transpose(tmat, ttmat)
subroutine, public hecmw_trimatmul_ttkt(hecmesh, bttmat, hecmat, btmat, iws, num_lagrange, hectkt)
subroutine, public hecmw_localmat_init_with_hecmat(bkmat, hecmat, num_lagrange)
subroutine, public hecmw_localmat_assemble(btmat, hecmesh, hecmeshnew)
subroutine, public hecmw_localmat_free(tmat)
subroutine, public hecmw_localmat_blocking(tmat, ndof, btmat)
subroutine, public hecmw_localmat_make_hecmat(hecmat, bttktmat, hectkt)
subroutine, public hecmw_localmat_mulvec(btmat, v, tv)
subroutine, public hecmw_localmat_multmat(bkmat, btmat, hecmesh, bktmat)
subroutine, public hecmw_localmat_add_hecmat(bkmat, hecmat)
subroutine, public hecmw_localmat_shrink_comm_table(bkmat, hecmesh)
integer(kind=kint), parameter cncol_item
num of column items to be migrated (2 or 3)
subroutine, public hecmw_trimatmul_ttkt_mpc(hecmesh, hecmat, hectkt)
subroutine, public hecmw_pair_array_append(parray, id, i1, i2)
subroutine, public hecmw_pair_array_finalize(parray)
integer(kind=kint) function, public hecmw_pair_array_find_id(parray, i1, i2)
subroutine, public hecmw_pair_array_init(parray, max_num)
subroutine, public hecmw_pair_array_sort(parray)
integer(kind=4), parameter kreal
integer(kind=kint), parameter hecmw_status_size
integer(kind=kint) function hecmw_comm_get_rank()
real(kind=kreal) function hecmw_wtime()
subroutine hecmw_isend_int(sbuf, sc, dest, tag, comm, req)
subroutine hecmw_recv_int(rbuf, rc, source, tag, comm, stat)
subroutine hecmw_isend_r(sbuf, sc, dest, tag, comm, req)
subroutine hecmw_waitall(cnt, reqs, stats)
subroutine hecmw_recv_r(rbuf, rc, source, tag, comm, stat)
subroutine hecmw_alltoall_int(sbuf, sc, rbuf, rc, comm)