40 public :: hecmwst_mf_factor
53 real(kind=
kreal),
parameter :: mf_pivot_u = 0.01d0
55 real(kind=
kreal),
parameter :: mf_pivot_alpha = (1.0d0 + sqrt(17.0d0)) / 8.0d0
57 real(kind=
kreal),
parameter :: mf_pivot_zero = 1.0d-14
62 real(kind=
kreal),
parameter :: mf_blr_eps = 1.0d-8
64 integer(kind=kint),
parameter :: mf_par_rows = 512
73 integer(kind=kint) :: ncol = 0
74 integer(kind=kint) :: npiv = 0
75 integer(kind=kint) :: nt = 0
76 integer(kind=kint) :: ntc = 0
77 integer(kind=8) :: cbsize = 0
78 integer(kind=8) :: pwords = 0
79 integer(kind=kint),
allocatable :: tbnd(:)
80 integer(kind=kint),
allocatable :: fsdof(:)
81 integer(kind=kint),
allocatable :: frow(:)
82 integer(kind=kint),
allocatable :: ptype(:)
83 real(kind=
kreal),
allocatable :: dsub(:)
84 integer(kind=8),
allocatable :: bptr(:)
86 integer(kind=kint),
allocatable :: brank(:)
87 integer(kind=8),
allocatable :: bptru(:)
88 integer(kind=kint),
allocatable :: branku(:)
89 real(kind=
kreal),
allocatable :: lval(:)
90 real(kind=
kreal),
allocatable :: uval(:)
91 real(kind=
kreal),
allocatable :: cval(:)
92 real(kind=
kreal),
allocatable :: dval(:)
95 type hecmwst_mf_factor
96 integer(kind=kint) :: tile = 0
97 integer(kind=kint) :: nnode = 0
98 integer(kind=kint) :: nsuper = 0
99 integer(kind=kint) :: ndof_tot = 0
100 integer(kind=kint),
allocatable :: cdofptr(:)
101 integer(kind=kint),
allocatable :: pdof(:)
102 integer(kind=kint),
allocatable :: chead(:)
103 integer(kind=kint),
allocatable :: cnext(:)
104 integer(kind=kint),
allocatable :: cptr(:)
105 integer(kind=kint),
allocatable :: clist(:)
106 integer(kind=kint),
allocatable :: mirror(:)
107 integer(kind=kint),
allocatable :: mirroru(:)
111 integer(kind=kint) :: mode = 0
112 logical :: scan = .false.
113 logical :: blr = .false.
114 real(kind=
kreal) :: eps = 0.0d0
115 integer(kind=kint) :: blr_eta = 0
118 real(kind=
kreal) :: blr_beta = 1.0d0
119 logical :: blr_reuse = .false.
121 real(kind=
kreal) :: pivot_u = 0.0d0
122 real(kind=
kreal) :: pivot_zero = 0.0d0
124 logical :: lu = .false.
125 real(kind=
kreal) :: asym = 0.0d0
126 type(mf_snode),
allocatable :: sn(:)
128 integer(kind=8) :: factor_words = 0
129 integer(kind=8) :: factor_nnz = 0
130 integer(kind=8) :: stack_peak = 0
131 integer(kind=8) :: front_words = 0
132 integer(kind=kint) :: max_tiles = 0
133 integer(kind=kint) :: max_rows = 0
134 integer(kind=kint) :: max_tile_dim = 0
136 integer(kind=8) :: factor_words_act = 0
137 integer(kind=8) :: stack_peak_act = 0
138 integer(kind=8) :: front_words_act = 0
139 integer(kind=8) :: live_cb = 0
140 integer(kind=8) :: live_front = 0
141 integer(kind=8) :: front_peak = 0
142 integer(kind=kint) :: n_pos = 0
143 integer(kind=kint) :: n_neg = 0
144 integer(kind=kint) :: n_2x2 = 0
145 integer(kind=kint) :: n_swap = 0
146 integer(kind=kint) :: n_delay = 0
147 integer(kind=kint) :: max_growth = 0
148 integer(kind=8) :: blr_words_fr = 0
149 integer(kind=8) :: blr_tiles = 0
150 integer(kind=8) :: blr_tiles_lr = 0
151 integer(kind=8) :: blr_rank_sum = 0
152 integer(kind=kint) :: blr_rank_max = 0
153 integer(kind=8) :: blr_tiles_skip = 0
157 integer(kind=8),
allocatable :: bdptr(:)
158 integer(kind=kint),
allocatable :: bdist(:)
161 integer(kind=kint),
allocatable :: bdkc(:)
162 integer(kind=kint),
allocatable :: bdct(:)
163 logical :: factored = .false.
164 end type hecmwst_mf_factor
169 integer(kind=kint) :: nt = 0
170 integer(kind=kint) :: ntc = 0
171 integer(kind=kint) :: nrow = 0
172 integer(kind=kint) :: ncol = 0
173 integer(kind=kint),
allocatable :: tb(:)
174 integer(kind=8),
allocatable :: coloff(:)
175 integer(kind=kint),
allocatable :: dtile(:)
181 real(kind=
kreal),
allocatable :: fval(:)
182 real(kind=
kreal),
allocatable :: fvalu(:)
183 real(kind=
kreal),
allocatable :: pval(:)
184 real(kind=
kreal),
allocatable :: pvalu(:)
185 real(kind=
kreal),
allocatable :: wval(:)
186 real(kind=
kreal),
allocatable :: wk(:)
187 integer(kind=kint),
allocatable :: rowoff(:)
188 integer(kind=kint),
allocatable :: pos(:)
189 integer(kind=kint),
allocatable :: blk(:)
190 integer(kind=kint),
allocatable :: blkr(:)
191 integer(kind=kint),
allocatable :: cmapdof(:)
192 integer(kind=kint),
allocatable :: ctb(:)
193 integer(kind=8),
allocatable :: ccoloff(:)
194 real(kind=
kreal),
allocatable :: bval(:)
195 integer(kind=8),
allocatable :: boff(:)
196 integer(kind=kint),
allocatable :: brk(:)
197 real(kind=
kreal),
allocatable :: bvalu(:)
198 integer(kind=8),
allocatable :: boffu(:)
199 integer(kind=kint),
allocatable :: brku(:)
204 integer(kind=kint),
allocatable :: hdr(:)
205 integer(kind=kint),
allocatable :: tl(:)
206 real(kind=
kreal),
allocatable :: rv(:)
212 integer(kind=kint) :: n = 0
213 integer(kind=kint) :: nreq = 0
214 integer(kind=kint),
allocatable :: reqs(:)
215 type(mf_sendbox),
allocatable :: box(:)
221 function mf_pool_slot(pool)
result(ib)
223 type(mf_pool),
intent(inout) :: pool
224 integer(kind=kint) :: ib
225 type(mf_sendbox),
allocatable :: t(:)
226 integer(kind=kint) :: i, m
230 if (.not.
allocated(pool%box))
then
231 allocate(pool%box(64))
232 else if (pool%n >
size(pool%box))
then
235 do i = 1,
size(pool%box)
236 if (
allocated(pool%box(i)%hdr))
call move_alloc(pool%box(i)%hdr, t(i)%hdr)
237 if (
allocated(pool%box(i)%tl))
call move_alloc(pool%box(i)%tl, t(i)%tl)
238 if (
allocated(pool%box(i)%rv))
call move_alloc(pool%box(i)%rv, t(i)%rv)
240 call move_alloc(t, pool%box)
242 end function mf_pool_slot
245 subroutine mf_pool_req(pool, rq)
247 type(mf_pool),
intent(inout) :: pool
248 integer(kind=kint),
intent(in) :: rq
250 pool%nreq = pool%nreq + 1
251 call mf_grow_i(pool%reqs, max(pool%nreq, 64))
252 pool%reqs(pool%nreq) = rq
253 end subroutine mf_pool_req
256 subroutine mf_pool_wait(pool)
258 type(mf_pool),
intent(inout) :: pool
259 integer(kind=kint),
allocatable :: stats(:,:)
260 integer(kind=kint) :: i
262 if (pool%nreq > 0)
then
268 if (
allocated(pool%box(i)%hdr))
deallocate(pool%box(i)%hdr)
269 if (
allocated(pool%box(i)%tl))
deallocate(pool%box(i)%tl)
270 if (
allocated(pool%box(i)%rv))
deallocate(pool%box(i)%rv)
274 end subroutine mf_pool_wait
282 type(hecmwst_mf_graph),
intent(in) :: graph
283 type(hecmwst_mf_symbolic),
intent(in) :: sym
284 integer(kind=kint),
intent(in) :: tile
285 logical,
intent(in) :: blr
286 integer(kind=kint),
intent(in) :: eta
287 type(hecmwst_mf_factor),
intent(out) :: fct
289 integer(kind=8),
allocatable :: cbsize(:)
290 integer(kind=kint),
allocatable :: odofptr(:), wptr(:)
291 integer(kind=kint),
allocatable :: vst(:), vds(:), que(:), dmin(:), ctile(:)
292 integer(kind=kint) :: n, ns, s, k, d, c, j, hj, stamp
293 integer(kind=8) :: top
302 allocate(fct%cdofptr(n+1), odofptr(n+1))
306 fct%cdofptr(k+1) = fct%cdofptr(k) + sym%ndof(k)
307 odofptr(k+1) = odofptr(k) + sym%ndof(sym%invp(k))
309 fct%ndof_tot = fct%cdofptr(n+1) - 1
310 allocate(fct%pdof(fct%ndof_tot))
312 do d = 1, sym%ndof(k)
313 fct%pdof(fct%cdofptr(k)+d-1) = odofptr(sym%perm(k)) + d - 1
318 allocate(fct%chead(ns), fct%cnext(ns))
323 fct%cnext(s) = fct%chead(c)
326 allocate(fct%cptr(ns+1))
330 if (c /= 0) fct%cptr(c+1) = fct%cptr(c+1) + 1
334 fct%cptr(s+1) = fct%cptr(s) + fct%cptr(s+1)
336 allocate(fct%clist(max(fct%cptr(ns+1)-1, 1)), wptr(ns))
337 wptr(1:ns) = fct%cptr(1:ns)
341 fct%clist(wptr(c)) = s
342 wptr(c) = wptr(c) + 1
350 fct%blr_eta = max(eta, 0)
351 blradm = blr .and. eta > 0 .and. graph%nnode == n
353 allocate(fct%bdptr(ns+1), fct%bdkc(ns), fct%bdct(ns), vst(n), vds(n), que(n))
363 fct%max_tile_dim = tile
365 call mf_partition(sym, fct, s, 0, g)
366 if (blradm)
call mf_blr_admis(graph, sym, fct, s, g, eta, vst, vds, que, stamp, dmin, ctile)
367 fct%factor_words = fct%factor_words + g%coloff(g%ntc+1)
368 cbsize(s) = g%coloff(g%nt+1) - g%coloff(g%ntc+1)
369 fct%front_words = max(fct%front_words, g%coloff(g%nt+1))
370 fct%factor_nnz = fct%factor_nnz + g%coloff(g%ntc+1)
372 hj = g%tb(k) - g%tb(k-1)
373 fct%factor_nnz = fct%factor_nnz - int(hj, 8)*(hj-1)/2
375 fct%max_tiles = max(fct%max_tiles, g%nt)
376 fct%max_rows = max(fct%max_rows, g%nrow)
378 fct%max_tile_dim = max(fct%max_tile_dim, g%tb(j) - g%tb(j-1))
384 fct%stack_peak = max(fct%stack_peak, top)
387 top = top - cbsize(c)
390 top = top + cbsize(s)
391 fct%stack_peak = max(fct%stack_peak, top)
396 fct%factored = .false.
401 type(hecmwst_mf_factor),
intent(inout) :: fct
404 fct%factored = .false.
405 if (
allocated(fct%cdofptr))
deallocate(fct%cdofptr)
406 if (
allocated(fct%pdof))
deallocate(fct%pdof)
407 if (
allocated(fct%chead))
deallocate(fct%chead)
408 if (
allocated(fct%cnext))
deallocate(fct%cnext)
409 if (
allocated(fct%cptr))
deallocate(fct%cptr)
410 if (
allocated(fct%clist))
deallocate(fct%clist)
411 if (
allocated(fct%mirror))
deallocate(fct%mirror)
412 if (
allocated(fct%mirroru))
deallocate(fct%mirroru)
413 if (
allocated(fct%sn))
deallocate(fct%sn)
414 if (
allocated(fct%bdptr))
deallocate(fct%bdptr)
415 if (
allocated(fct%bdist))
deallocate(fct%bdist)
416 if (
allocated(fct%bdkc))
deallocate(fct%bdkc)
417 if (
allocated(fct%bdct))
deallocate(fct%bdct)
424 type(hecmwst_mf_factor),
intent(inout) :: fct
425 integer(kind=kint),
intent(in) :: root
426 integer(kind=kint) :: hdr(1)
427 integer(kind=kint) :: ns, len, comm, s
436 allocate(fct%bdptr(ns+1), fct%bdkc(ns), fct%bdct(ns), fct%bdist(max(len, 1)))
444 fct%bdptr(s+1) = fct%bdptr(s) + int(fct%bdkc(s), 8)*fct%bdct(s)
449 subroutine mf_grow_r(a, need)
451 real(kind=
kreal),
allocatable,
intent(inout) :: a(:)
452 integer(kind=8),
intent(in) :: need
453 real(kind=
kreal),
allocatable :: t(:)
456 if (
allocated(a))
then
457 if (
size(a, kind=8) >= need)
return
458 n = max(need,
size(a, kind=8) +
size(a, kind=8)/2)
460 t(1:
size(a, kind=8)) = a(:)
461 call move_alloc(t, a)
463 allocate(a(max(need, 1_8)))
465 end subroutine mf_grow_r
467 subroutine mf_grow_i(a, need)
469 integer(kind=kint),
allocatable,
intent(inout) :: a(:)
470 integer(kind=kint),
intent(in) :: need
471 integer(kind=kint),
allocatable :: t(:)
472 integer(kind=kint) :: n
474 if (
allocated(a))
then
475 if (
size(a) >= need)
return
476 n = max(need,
size(a) +
size(a)/2)
479 call move_alloc(t, a)
481 allocate(a(max(need, 1)))
483 end subroutine mf_grow_i
485 subroutine mf_grow_i8(a, need)
487 integer(kind=8),
allocatable,
intent(inout) :: a(:)
488 integer(kind=kint),
intent(in) :: need
489 integer(kind=8),
allocatable :: t(:)
490 integer(kind=kint) :: n
492 if (
allocated(a))
then
493 if (
size(a) >= need)
return
494 n = max(need,
size(a) +
size(a)/2)
497 call move_alloc(t, a)
499 allocate(a(max(need, 1)))
501 end subroutine mf_grow_i8
507 subroutine mf_partition(sym, fct, s, ndel, g)
509 type(hecmwst_mf_symbolic),
intent(in) :: sym
510 type(hecmwst_mf_factor),
intent(in) :: fct
511 integer(kind=kint),
intent(in) :: s, ndel
512 type(mf_grid),
intent(inout) :: g
513 integer(kind=kint) :: nown, nrow_nodes, ncol0, ncb, nt, i, j, nd, cur, acc, tile, pass
516 nown = sym%sptr(s+1) - sym%sptr(s)
517 nrow_nodes = sym%rptr(s+1) - sym%rptr(s)
518 ncol0 = fct%cdofptr(sym%sptr(s+1)) - fct%cdofptr(sym%sptr(s))
520 do i = nown + 1, nrow_nodes
521 ncb = ncb + sym%ndof(sym%rlist(sym%rptr(s)+i-1))
523 g%ncol = ncol0 + ndel
524 g%nrow = g%ncol + ncb
530 if (i == nown + 1)
then
534 if (pass == 2) g%tb(nt) = cur
538 if (pass == 2) g%tb(nt) = acc
542 nd = sym%ndof(sym%rlist(sym%rptr(s)+i-1))
543 if (acc > cur .and. acc + nd - cur > tile)
then
545 if (pass == 2) g%tb(nt) = acc
550 if (nown == nrow_nodes)
then
554 if (pass == 2) g%tb(nt) = cur
559 if (pass == 2) g%tb(nt) = acc
560 if (nown == nrow_nodes) g%ntc = nt
562 if (
allocated(g%tb))
then
563 if (
size(g%tb) < nt + 1)
deallocate(g%tb, g%coloff)
565 if (.not.
allocated(g%tb))
allocate(g%tb(0:nt), g%coloff(nt+1))
566 if (
allocated(g%dtile))
then
567 if (
size(g%dtile) < g%nrow)
deallocate(g%dtile)
569 if (.not.
allocated(g%dtile))
allocate(g%dtile(0:max(g%nrow, 1)-1))
576 g%coloff(j+1) = g%coloff(j) + int(g%nrow - g%tb(j-1), 8)*(g%tb(j) - g%tb(j-1))
577 g%dtile(g%tb(j-1):g%tb(j)-1) = j
579 end subroutine mf_partition
586 subroutine mf_blr_admis(graph, sym, fct, s, g, eta, vst, vds, que, stamp, dmin, ctile)
588 type(hecmwst_mf_graph),
intent(in) :: graph
589 type(hecmwst_mf_symbolic),
intent(in) :: sym
590 type(hecmwst_mf_factor),
intent(inout) :: fct
591 integer(kind=kint),
intent(in) :: s, eta
592 type(mf_grid),
intent(in) :: g
593 integer(kind=kint),
intent(inout) :: vst(:), vds(:), que(:), stamp
594 integer(kind=kint),
allocatable,
intent(inout) :: dmin(:), ctile(:)
595 integer(kind=kint) :: nown, nrow_nodes, nct, k, i, i0, j, p, nd, nb, d, head, tail, acc
596 integer(kind=8) :: base
598 nown = sym%sptr(s+1) - sym%sptr(s)
599 nrow_nodes = sym%rptr(s+1) - sym%rptr(s)
603 fct%bdptr(s+1) = fct%bdptr(s) + int(g%ntc, 8)*nct
605 call mf_grow_i(fct%bdist, int(fct%bdptr(s+1) - 1, kind=
kint))
606 call mf_grow_i(ctile, nrow_nodes - nown)
607 call mf_grow_i(dmin, nct)
609 do i = 1, nrow_nodes - nown
610 ctile(i) = g%dtile(acc) - g%ntc
611 acc = acc + sym%ndof(sym%rlist(sym%rptr(s)+nown+i-1))
613 base = fct%bdptr(s) - 1
619 do j = sym%sptr(s), sym%sptr(s+1) - 1
620 if (acc >= g%tb(k))
exit
621 if (acc >= g%tb(k-1))
then
628 acc = acc + sym%ndof(j)
630 do while (head < tail)
635 do p = graph%xadj(nd), graph%xadj(nd+1) - 1
637 if (vst(nb) /= stamp)
then
645 dmin(1:nct) = eta + 1
646 do i = 1, nrow_nodes - nown
647 nd = sym%perm(sym%rlist(sym%rptr(s)+nown+i-1))
648 if (vst(nd) == stamp)
then
650 dmin(i0) = min(dmin(i0), vds(nd))
654 fct%bdist(base + int(i0-1, 8)*g%ntc + k) = dmin(i0)
657 end subroutine mf_blr_admis
664 function mf_blr_skip(fct, s, ntc, k, i)
result(skip)
666 type(hecmwst_mf_factor),
intent(in) :: fct
667 integer(kind=kint),
intent(in) :: s, ntc, k, i
669 integer(kind=kint) :: i0, d
672 if (fct%blr_eta <= 0)
return
673 if (.not.
allocated(fct%bdptr))
return
674 if (k > fct%bdkc(s))
return
676 if (i0 < 1 .or. i0 > fct%bdct(s))
return
677 d = fct%bdist(fct%bdptr(s) - 1 + int(i0-1, 8)*fct%bdkc(s) + k)
678 skip = d >= 1 .and. d <= fct%blr_eta
679 end function mf_blr_skip
682 function mf_off(g, i, j)
result(off)
684 type(mf_grid),
intent(in) :: g
685 integer(kind=kint),
intent(in) :: i, j
686 integer(kind=8) :: off
688 off = g%coloff(j) + int(g%tb(i-1) - g%tb(j-1), 8)*(g%tb(j) - g%tb(j-1))
692 function mf_idx(g, r, c)
result(idx)
694 type(mf_grid),
intent(in) :: g
695 integer(kind=kint),
intent(in) :: r, c
696 integer(kind=8) :: idx
697 integer(kind=kint) :: ti, tj
701 idx = mf_off(g, ti, tj) + int(c - 1 - g%tb(tj-1), 8)*(g%tb(ti) - g%tb(ti-1)) + (r - 1 - g%tb(ti-1)) + 1
706 function mf_bidx(nt, k, i)
result(idx)
708 integer(kind=kint),
intent(in) :: nt, k, i
709 integer(kind=kint) :: idx
711 idx = (k-1)*(nt+1) - (k-1)*k/2 + (i - k + 1)
717 subroutine mf_update_ab(m, n, k, ra, fa, ua, va, rb, fb, ub, vb, c)
719 integer(kind=kint),
intent(in) :: m, n, k, ra, rb
720 real(kind=
kreal),
intent(in) :: fa(*), ua(*), va(*), fb(*), ub(*), vb(*)
721 real(kind=
kreal),
intent(inout) :: c(*)
722 real(kind=
kreal),
allocatable :: s(:), t(:)
724 if (ra == 0 .or. rb == 0)
return
725 if (ra < 0 .and. rb < 0)
then
727 else if (ra >= 0 .and. rb < 0)
then
728 allocate(t(int(n, 8)*ra))
732 else if (ra < 0)
then
733 allocate(t(int(m, 8)*rb))
738 allocate(s(int(ra, 8)*rb), t(int(m, 8)*rb))
744 end subroutine mf_update_ab
747 subroutine mf_col_copy(g, fval, x, i0, vec, put)
749 type(mf_grid),
intent(in) :: g
750 real(kind=
kreal),
intent(inout) :: fval(:)
751 integer(kind=kint),
intent(in) :: x, i0
752 real(kind=
kreal),
intent(inout) :: vec(*)
753 logical,
intent(in) :: put
754 integer(kind=kint) :: it, kx, r0, r1, h
755 integer(kind=8) :: base
757 if (i0 > g%nrow)
return
759 do it = g%dtile(i0-1), g%nt
760 h = g%tb(it) - g%tb(it-1)
761 r0 = max(i0, g%tb(it-1) + 1)
763 base = mf_off(g, it, kx) + int(x - 1 - g%tb(kx-1), 8)*h + (r0 - 1 - g%tb(it-1))
765 fval(base+1:base+r1-r0+1) = vec(r0-i0+1:r1-i0+1)
767 vec(r0-i0+1:r1-i0+1) = fval(base+1:base+r1-r0+1)
770 end subroutine mf_col_copy
774 subroutine mf_col_update(g, fval, ptype, dsub, x, q1, q2, vec, t)
776 type(mf_grid),
intent(in) :: g
777 real(kind=
kreal),
intent(in) :: fval(:)
778 integer(kind=kint),
intent(in) :: ptype(:)
779 real(kind=
kreal),
intent(in) :: dsub(:)
780 integer(kind=kint),
intent(in) :: x, q1, q2
781 real(kind=
kreal),
intent(inout) :: vec(*)
782 real(kind=
kreal),
intent(inout) :: t(*)
783 integer(kind=kint) :: q
784 real(kind=
kreal) :: l1, l2, d11, d21, d22
788 if (ptype(q) == 1)
then
789 t(q-q1+1) = fval(mf_idx(g, x, q)) * fval(mf_idx(g, q, q))
792 l1 = fval(mf_idx(g, x, q))
793 l2 = fval(mf_idx(g, x, q+1))
794 d11 = fval(mf_idx(g, q, q))
796 d22 = fval(mf_idx(g, q+1, q+1))
797 t(q-q1+1) = l1*d11 + l2*d21
798 t(q-q1+2) = l1*d21 + l2*d22
802 call mf_col_axpy(g, fval, x, q1, q2, t, vec)
803 end subroutine mf_col_update
807 subroutine mf_col_axpy(g, gval, i0, q1, q2, t, vec)
809 type(mf_grid),
intent(in) :: g
810 real(kind=
kreal),
intent(in) :: gval(:)
811 integer(kind=kint),
intent(in) :: i0, q1, q2
812 real(kind=
kreal),
intent(in) :: t(*)
813 real(kind=
kreal),
intent(inout) :: vec(*)
814 integer(kind=kint) :: q, it, kq, r0, r1, h
817 if (i0 > g%nrow)
return
820 do it = g%dtile(i0-1), g%nt
821 h = g%tb(it) - g%tb(it-1)
822 r0 = max(i0, g%tb(it-1) + 1)
824 b = mf_off(g, it, kq) + int(q - 1 - g%tb(kq-1), 8)*h + (r0 - 1 - g%tb(it-1))
825 vec(r0-i0+1:r1-i0+1) = vec(r0-i0+1:r1-i0+1) - gval(b+1:b+r1-r0+1) * t(q-q1+1)
828 end subroutine mf_col_axpy
831 subroutine mf_row_get(g, gval, x, q1, q2, t)
833 type(mf_grid),
intent(in) :: g
834 real(kind=
kreal),
intent(in) :: gval(:)
835 integer(kind=kint),
intent(in) :: x, q1, q2
836 real(kind=
kreal),
intent(out) :: t(*)
837 integer(kind=kint) :: q
840 t(q-q1+1) = gval(mf_idx(g, x, q))
842 end subroutine mf_row_get
845 subroutine mf_swap(g, fval, x, y)
847 type(mf_grid),
intent(in) :: g
848 real(kind=
kreal),
intent(inout) :: fval(:)
849 integer(kind=kint),
intent(in) :: x, y
850 integer(kind=kint) :: i
851 integer(kind=8) :: a, b
852 real(kind=
kreal) :: v
880 end subroutine mf_swap
883 subroutine mf_swap_row(g, fval, fvalu, x, y)
885 type(mf_grid),
intent(in) :: g
886 real(kind=
kreal),
intent(inout) :: fval(:), fvalu(:)
887 integer(kind=kint),
intent(in) :: x, y
888 integer(kind=kint) :: c
889 integer(kind=8) :: a, b
890 real(kind=
kreal) :: v
918 end subroutine mf_swap_row
921 subroutine mf_swap_lu(g, fval, fvalu, x, y)
923 type(mf_grid),
intent(in) :: g
924 real(kind=
kreal),
intent(inout) :: fval(:), fvalu(:)
925 integer(kind=kint),
intent(in) :: x, y
926 integer(kind=kint) :: i
927 integer(kind=8) :: a, b
928 real(kind=
kreal) :: v
969 end subroutine mf_swap_lu
973 subroutine mf_mirror(hecMAT, fct, ierr)
976 type(hecmwst_mf_factor),
intent(inout) :: fct
977 integer(kind=kint),
intent(out) :: ierr
978 integer(kind=kint),
allocatable :: ptr(:)
979 integer(kind=kint) :: i, j, kk, ku, l
982 if (
allocated(fct%mirror))
then
983 if (
size(fct%mirror) == hecmat%NPL .and.
size(fct%mirroru) == hecmat%NPU)
return
984 deallocate(fct%mirror, fct%mirroru)
986 if (hecmat%NPL /= hecmat%NPU)
then
990 allocate(fct%mirror(max(hecmat%NPL, 1)), fct%mirroru(max(hecmat%NPU, 1)), ptr(hecmat%NP))
992 ptr(1:hecmat%NP) = hecmat%indexU(0:hecmat%NP-1) + 1
994 do kk = hecmat%indexL(j-1)+1, hecmat%indexL(j)
997 if (ku <= hecmat%indexU(i))
then
998 if (hecmat%itemU(ku) /= j) ku = 0
1003 do l = hecmat%indexU(i-1)+1, hecmat%indexU(i)
1004 if (hecmat%itemU(l) == j) ku = l
1014 fct%mirroru(ku) = kk
1018 if (ierr == 0 .and. hecmat%NPU > 0)
then
1019 if (minval(fct%mirroru(1:hecmat%NPU)) == 0) ierr = -2
1022 if (ierr /= 0)
deallocate(fct%mirror, fct%mirroru)
1023 end subroutine mf_mirror
1027 subroutine mf_asymmetry(hecMAT, fct, amax)
1030 type(hecmwst_mf_factor),
intent(inout) :: fct
1031 real(kind=
kreal),
intent(in) :: amax
1032 integer(kind=kint) :: nd, i, kk, a, b
1033 integer(kind=8) :: base, base2
1034 real(kind=
kreal) :: asym
1036 fct%lu = (fct%mode == 2) .or. (fct%mode == 0 .and. .not. hecmat%symmetric)
1038 if (.not. fct%scan)
return
1042 base = int(i-1, 8)*nd*nd
1045 asym = max(asym, abs(hecmat%D(base + (a-1)*nd + b) - hecmat%D(base + (b-1)*nd + a)))
1049 do kk = 1, hecmat%NPL
1050 base = int(kk-1, 8)*nd*nd
1051 base2 = int(fct%mirror(kk)-1, 8)*nd*nd
1054 asym = max(asym, abs(hecmat%AL(base + (a-1)*nd + b) - hecmat%AU(base2 + (b-1)*nd + a)))
1058 if (amax > 0.0d0) fct%asym = asym / amax
1059 end subroutine mf_asymmetry
1064 subroutine mf_store_cb(g, gval, npiv, cndel, cnb, kbeg, ctb, ccoloff, sval, top, par)
1066 type(mf_grid),
intent(in) :: g
1067 real(kind=
kreal),
intent(in) :: gval(:)
1068 integer(kind=kint),
intent(in) :: npiv, cndel, cnb, kbeg
1069 integer(kind=kint),
intent(in) :: ctb(0:)
1070 integer(kind=8),
intent(in) :: ccoloff(:)
1071 real(kind=
kreal),
intent(inout) :: sval(:)
1072 integer(kind=8),
intent(in) :: top
1073 logical,
intent(in) :: par
1074 integer(kind=kint) :: i, j, hi, cc, r, m
1075 integer(kind=8) :: base, o
1077 do j = 1, min(cnb, kbeg)
1079 hi = ctb(i) - ctb(i-1)
1080 base = top + int(ctb(i-1), 8)*cndel
1082 do r = max(cc, ctb(i-1) + 1), ctb(i)
1083 if (r <= cndel)
then
1086 m = g%ncol + r - cndel
1088 sval(base + int(cc-1, 8)*hi + (r - ctb(i-1))) = gval(mf_idx(g, m, npiv + cc))
1094 do j = g%ntc + 1, g%nt
1097 hi = g%tb(i) - g%tb(i-1)
1098 m = g%tb(j) - g%tb(j-1)
1099 base = top + ccoloff(kbeg+j-g%ntc) + int(ctb(kbeg+i-g%ntc-1) - ctb(kbeg+j-g%ntc-1), 8)*m
1101 sval(base+1:base+int(hi, 8)*m) = gval(o+1:o+int(hi, 8)*m)
1106 end subroutine mf_store_cb
1115 type(hecmwst_mf_symbolic),
intent(in) :: sym
1116 type(hecmwst_mf_factor),
intent(inout) :: fct
1117 integer(kind=kint),
intent(out) :: ierr
1118 type(mf_work),
allocatable :: wrks(:)
1119 integer(kind=kint),
allocatable :: left(:)
1120 integer(kind=kint) :: s, k, nd, nthr
1121 real(kind=
kreal) :: amax, zero
1124 fct%factored = .false.
1126 if (sym%ndof(k) /= hecmat%NDOF)
then
1132 amax = maxval(abs(hecmat%D(1:hecmat%NP*nd*nd)))
1133 if (hecmat%NPL > 0) amax = max(amax, maxval(abs(hecmat%AL(1:hecmat%NPL*nd*nd))))
1134 if (hecmat%NPU > 0) amax = max(amax, maxval(abs(hecmat%AU(1:hecmat%NPU*nd*nd))))
1135 if (.not. (fct%pivot_u > 0.0d0)) fct%pivot_u = mf_pivot_u
1136 if (.not. (fct%pivot_zero > 0.0d0)) fct%pivot_zero = mf_pivot_zero
1139 if (.not. (fct%eps > 0.0d0)) fct%eps = mf_blr_eps
1141 zero = fct%pivot_zero * amax
1142 call mf_mirror(hecmat, fct, ierr)
1143 if (ierr /= 0)
return
1144 call mf_asymmetry(hecmat, fct, amax)
1152 fct%factor_words_act = 0
1153 fct%stack_peak_act = 0
1154 fct%front_words_act = 0
1158 fct%blr_words_fr = 0
1160 fct%blr_tiles_lr = 0
1161 fct%blr_rank_sum = 0
1162 fct%blr_rank_max = 0
1163 fct%blr_tiles_skip = 0
1167 allocate(wrks(0:nthr-1), left(fct%nsuper))
1168 do s = 1, fct%nsuper
1169 left(s) = fct%cptr(s+1) - fct%cptr(s)
1173 do s = 1, fct%nsuper
1176 if (fct%cptr(s+1) == fct%cptr(s))
then
1178 call mf_super_task(hecmat, sym, fct, wrks, left, s, zero, ierr)
1184 deallocate(wrks, left)
1185 if (ierr /= 0)
return
1186 fct%factored = .true.
1195 subroutine mf_super_task(hecMAT, sym, fct, wrks, left, s, zero, gierr)
1198 type(hecmwst_mf_symbolic),
intent(in) :: sym
1199 type(hecmwst_mf_factor),
intent(inout) :: fct
1200 type(mf_work),
intent(inout) :: wrks(0:)
1201 integer(kind=kint),
intent(inout) :: left(:)
1202 integer(kind=kint),
intent(in) :: s
1203 real(kind=
kreal),
intent(in) :: zero
1204 integer(kind=kint),
intent(inout) :: gierr
1205 integer(kind=kint) :: ierr, cur, tid, p, n, ss, c, l, nr, nteam
1219 nr = nr + fct%sn(c)%ncol - fct%sn(c)%npiv
1222 do l = sym%rptr(ss), sym%rptr(ss+1) - 1
1223 nr = nr + sym%ndof(sym%rlist(l))
1225 if (nr >= mf_par_rows .and. nteam > 1)
then
1229 type(mf_work) :: lwrk
1230 call mf_super_factor(hecmat, sym, fct, fct%sn(ss), lwrk, ss, zero, ierr)
1233 call mf_super_factor(hecmat, sym, fct, fct%sn(ss), wrks(tid), ss, zero, ierr)
1245 left(p) = left(p) - 1
1251 end subroutine mf_super_task
1269 type(hecmwst_matrix),
intent(in) :: hecmat
1270 type(hecmwst_mf_symbolic),
intent(in) :: sym
1271 type(hecmwst_mf_map),
intent(in) :: map
1272 type(hecmwst_mf_factor),
intent(inout) :: fct
1273 integer(kind=kint),
intent(out) :: ierr
1274 type(mf_work),
allocatable :: wrks(:)
1275 integer(kind=kint),
allocatable :: left(:)
1276 integer(kind=kint) :: s, k, nd, nthr, iu, gierr, iw(1)
1277 real(kind=kreal) :: amax, zero, rw(1)
1280 fct%factored = .false.
1282 if (sym%ndof(k) /= hecmat%NDOF)
then
1291 amax = maxval(abs(hecmat%D(1:hecmat%NP*nd*nd)))
1292 if (hecmat%NPL > 0) amax = max(amax, maxval(abs(hecmat%AL(1:hecmat%NPL*nd*nd))))
1293 if (hecmat%NPU > 0) amax = max(amax, maxval(abs(hecmat%AU(1:hecmat%NPU*nd*nd))))
1295 call hecmw_allreduce_r_comm(rw, 1, hecmw_max, map%comm)
1297 if (.not. (fct%pivot_u > 0.0d0)) fct%pivot_u = mf_pivot_u
1298 if (.not. (fct%pivot_zero > 0.0d0)) fct%pivot_zero = mf_pivot_zero
1300 if (.not. hecmw_mf_kernel_blr_available()) fct%blr = .false.
1301 if (.not. (fct%eps > 0.0d0)) fct%eps = mf_blr_eps
1303 zero = fct%pivot_zero * amax
1304 call mf_mirror(hecmat, fct, ierr)
1306 call hecmw_allreduce_i_comm(iw, 1, hecmw_min, map%comm)
1308 if (ierr /= 0)
return
1311 iw(1) = merge(1, 0, fct%scan)
1312 call hecmw_allreduce_i_comm(iw, 1, hecmw_max, map%comm)
1313 fct%scan = iw(1) /= 0
1314 call mf_asymmetry(hecmat, fct, amax)
1316 call hecmw_allreduce_r_comm(rw, 1, hecmw_max, map%comm)
1325 fct%factor_words_act = 0
1326 fct%stack_peak_act = 0
1327 fct%front_words_act = 0
1331 fct%blr_words_fr = 0
1333 fct%blr_tiles_lr = 0
1334 fct%blr_rank_sum = 0
1335 fct%blr_rank_max = 0
1336 fct%blr_tiles_skip = 0
1340 allocate(wrks(0:nthr-1), left(fct%nsuper))
1341 do s = 1, fct%nsuper
1342 left(s) = fct%cptr(s+1) - fct%cptr(s)
1347 do s = 1, fct%nsuper
1348 if (fct%cptr(s+1) == fct%cptr(s) .and. map%owner(s) == map%myrank .and. .not. map%upper(s))
then
1350 call mf_super_task_mpi(hecmat, sym, map, fct, wrks, left, s, zero, gierr)
1356 deallocate(wrks, left)
1363 call hecmw_allreduce_i_comm(iw, 1, hecmw_max, map%comm)
1365 if (gierr /= 0)
then
1369 do iu = 1, map%nupper
1371 if (map%myrank >= map%rbeg(s) .and. map%myrank < map%rbeg(s) + map%rcnt(s))
then
1372 call mf_super_factor_dist(hecmat, sym, map, fct, s, zero, gierr)
1376 call hecmw_allreduce_i_comm(iw, 1, hecmw_max, map%comm)
1379 if (ierr /= 0)
return
1380 fct%factored = .true.
1400 subroutine mf_super_factor_dist(hecMAT, sym, map, fct, s, zero, gierr)
1402 type(hecmwst_matrix),
intent(in) :: hecmat
1403 type(hecmwst_mf_symbolic),
intent(in) :: sym
1404 type(hecmwst_mf_map),
intent(in) :: map
1405 type(hecmwst_mf_factor),
intent(inout) :: fct
1406 integer(kind=kint),
intent(in) :: s
1407 real(kind=kreal),
intent(in) :: zero
1408 integer(kind=kint),
intent(inout) :: gierr
1410 type(mf_pool) :: pool
1412 integer(kind=kint),
allocatable :: ctb(:), towner(:), wrank(:), crow2t(:), tprow(:)
1413 integer(kind=8),
allocatable :: boff(:)
1414 real(kind=kreal),
allocatable,
target :: band(:), bandu(:), fsv(:), fsvu(:)
1415 integer(kind=8),
allocatable :: f2of(:)
1416 integer(kind=kint),
allocatable :: rowoff(:), pos(:), blk(:), blkr(:), mdel(:)
1417 integer(kind=kint),
allocatable :: perm(:), permr(:), itmp(:), swaps(:)
1418 integer(kind=kint),
allocatable :: cctb(:), cown(:), pmap(:), pair_i(:), pair_j(:)
1419 integer(kind=kint),
allocatable :: cctb2(:), cown2(:), ptv(:)
1420 integer(kind=kint),
allocatable :: plist(:), chk(:)
1421 logical,
allocatable :: pmem(:)
1422 real(kind=kreal),
allocatable :: pv(:), pvu(:), wk(:), tmat(:), tmatu(:), dgb(:), dgbu(:)
1423 real(kind=kreal),
allocatable :: pfull(:), pfullu(:), cmax(:), pband(:), pbandu(:)
1424 real(kind=kreal),
allocatable :: pfs(:), pfsu(:)
1425 real(kind=kreal),
allocatable,
target :: wscb(:), wfs(:), dum(:)
1426 real(kind=kreal),
allocatable :: tvec(:), tvecu(:), dsv(:), dgv(:), dd(:)
1427 real(kind=kreal),
allocatable :: sw1(:), sw2(:), swb(:)
1428 integer(kind=kint),
allocatable :: spr(:), spc(:), spf(:), sqr(:), sqc(:), sqf(:)
1429 real(kind=kreal),
allocatable,
target :: bval(:), bvalu(:)
1430 integer(kind=8),
allocatable :: bslot(:), bslotu(:), wfo(:), fxoff(:)
1431 integer(kind=kint),
allocatable :: brk(:), brku(:), xrk(:), xrku(:), xw(:), fxs(:)
1432 integer(kind=8),
allocatable :: xoff(:), xoffu(:)
1434 real(kind=kreal),
allocatable :: v(:)
1436 type(rbuf_t),
allocatable,
target :: xbuf(:), fxbuf(:)
1437 integer(kind=kint) :: ctli(8), stat(hecmw_status_size)
1438 integer(kind=kint) :: nd, ncol0, ndel, ncol, nrow, ncb, ncbt, nwk, nmine, myrows
1439 integer(kind=kint) :: master, me, i, j, k, c, t, m, nkc
1440 integer(kind=kint) :: npiv, nfs, p0, pa, pb, w, np, pk, pk0, mfs, info, nsw, n22, cndel
1441 integer(kind=kint) :: nown, nrow_nodes, hi, npair, req, cnb, kbegc, cnrow, wme, ferr
1442 integer(kind=kint) :: pr2, pc2, nfst, npl, nchk, rd, myfsr
1443 integer(kind=8) :: bandw, fw, btop, btopu, fsw
1444 real(kind=kreal) :: uinv
1445 logical :: ismaster, hasband, hasfs, fsok, acc, isroot, ruse_f
1446 integer(kind=kint) :: nswap_f, n2x2_f, npos_f, nneg_f, ntl_f, nlr_f, nsk_f, rmax_f
1447 integer(kind=8) :: rsum_f
1451 master = map%owner(s)
1452 ismaster = (master == me)
1453 nown = sym%sptr(s+1) - sym%sptr(s)
1454 nrow_nodes = sym%rptr(s+1) - sym%rptr(s)
1455 ncol0 = fct%cdofptr(sym%sptr(s+1)) - fct%cdofptr(sym%sptr(s))
1456 uinv = 1.0d0 / fct%pivot_u
1461 do j = fct%cptr(s), fct%cptr(s+1) - 1
1463 if (map%owner(c) == me)
call send_hdr(c)
1465 do j = fct%cptr(s), fct%cptr(s+1) - 1
1467 if (map%owner(c) /= me)
call recv_hdr(c)
1474 if (gierr == 0) gierr = ferr
1479 call mf_grow_i(fct%sn(s)%tbnd, 1)
1480 call mf_grow_i(fct%sn(s)%fsdof, 1)
1481 call mf_grow_i(fct%sn(s)%frow, 1)
1482 fct%sn(s)%tbnd(1) = 0
1483 call mf_pool_wait(pool)
1484 do j = fct%cptr(s), fct%cptr(s+1) - 1
1486 if (
allocated(fct%sn(c)%cval))
then
1487 fct%live_cb = fct%live_cb - fct%sn(c)%cbsize
1488 deallocate(fct%sn(c)%cval)
1497 allocate(mdel(max(fct%cptr(s+1) - fct%cptr(s), 1)))
1501 do j = fct%cptr(s), fct%cptr(s+1) - 1
1502 if (fct%clist(j) == c) mdel(j - fct%cptr(s) + 1) = m
1504 m = m + fct%sn(c)%ncol - fct%sn(c)%npiv
1508 if (ismaster) fct%max_growth = max(fct%max_growth, ndel)
1509 call mf_partition(sym, fct, s, ndel, g)
1513 call hecmw_mf_dist_rowtiles(sym, map, fct%tile, s, ndel, ncbt, ctb, towner)
1514 allocate(crow2t(max(ncb, 1)))
1516 crow2t(ctb(t-1)+1:ctb(t)) = t
1519 allocate(wrank(max(ncbt, 1)), tprow(0:max(ncbt, 1)), boff(max(ncbt, 1)))
1524 wrank(1) = towner(t)
1525 else if (towner(t) /= wrank(nwk))
then
1527 wrank(nwk) = towner(t)
1536 if (towner(t) == me)
then
1539 hi = ctb(t) - ctb(t-1)
1540 bandw = bandw + int(hi, 8)*(ncol + ctb(t))
1541 myrows = myrows + hi
1544 hasband = (nmine > 0)
1545 isroot = (ncbt == 0)
1547 call hecmw_mf_dist_fsgrid(map, s, g%ntc, pr2, pc2)
1548 nfst = mf_bidx(g%ntc, g%ntc, g%ntc)
1549 allocate(f2of(nfst))
1553 if (fso(i, j) == me)
then
1554 f2of(mf_bidx(g%ntc, j, i)) = fsw
1555 fsw = fsw + int(g%tb(i) - g%tb(i-1), 8)*(g%tb(j) - g%tb(j-1))
1563 allocate(plist(map%rcnt(s)), chk(map%rcnt(s)), pmem(map%rcnt(s)))
1564 pmem(1:map%rcnt(s)) = .false.
1565 pmem(master - map%rbeg(s) + 1) = .true.
1568 pmem(fso(i, j) - map%rbeg(s) + 1) = .true.
1572 pmem(wrank(i) - map%rbeg(s) + 1) = .true.
1575 do i = map%rbeg(s), map%rbeg(s) + map%rcnt(s) - 1
1576 if (pmem(i - map%rbeg(s) + 1))
then
1581 allocate(rowoff(nrow_nodes + 1), pos(sym%nnode))
1583 do i = 1, nrow_nodes
1584 k = sym%rlist(sym%rptr(s)+i-1)
1585 rowoff(i+1) = rowoff(i) + sym%ndof(k)
1586 if (i == nown) rowoff(i+1) = rowoff(i+1) + ndel
1592 if (hasband) fw = fw + bandw
1593 if (fct%lu) fw = 2*fw
1596 if (fct%lu)
allocate(fsvu(fsw))
1599 allocate(band(bandw))
1600 if (fct%lu)
allocate(bandu(bandw))
1602 fct%front_words_act = max(fct%front_words_act, fw)
1603 fct%live_front = fct%live_front + fw
1604 fct%front_peak = max(fct%front_peak, fct%live_front)
1608 if (hasfs) fsv(:) = 0.0d0
1610 if (hasfs .and. fct%lu) fsvu(:) = 0.0d0
1612 if (hasband) band(:) = 0.0d0
1614 if (hasband .and. fct%lu) bandu(:) = 0.0d0
1623 if (fct%blr .and. fct%blr_reuse .and.
allocated(fct%sn(s)%brank) .and.
allocated(fct%sn(s)%tbnd))
then
1624 if (fct%sn(s)%nt == g%nt .and. fct%sn(s)%ntc == g%ntc .and. fct%sn(s)%npiv == fct%sn(s)%ncol .and. &
1625 size(fct%sn(s)%tbnd) >= g%nt + 1)
then
1626 ruse_f = all(fct%sn(s)%tbnd(1:g%nt+1) == g%tb(0:g%nt))
1636 allocate(cctb(0:g%nt+2), cown(g%nt+2), pmap(max(nrow, 1)))
1637 do j = fct%cptr(s), fct%cptr(s+1) - 1
1639 call ea_setup(c, j - fct%cptr(s) + 1)
1644 do j = fct%cptr(s), fct%cptr(s+1) - 1
1645 if (fct%clist(j) == c)
call ea_setup(c, j - fct%cptr(s) + 1)
1647 call ea_recv_apply(c)
1662 if (ismaster .or. hasband .or. hasfs)
then
1664 if (gierr == 0)
call store_dist()
1666 if (ismaster .and. gierr == 0)
then
1667 fct%n_swap = fct%n_swap + nswap_f
1668 fct%n_2x2 = fct%n_2x2 + n2x2_f
1669 fct%n_delay = fct%n_delay + ncol - npiv
1671 fct%blr_tiles = fct%blr_tiles + ntl_f
1672 fct%blr_tiles_lr = fct%blr_tiles_lr + nlr_f
1673 fct%blr_rank_sum = fct%blr_rank_sum + rsum_f
1674 fct%blr_rank_max = max(fct%blr_rank_max, rmax_f)
1675 fct%blr_tiles_skip = fct%blr_tiles_skip + nsk_f
1676 fct%n_pos = fct%n_pos + npos_f
1677 fct%n_neg = fct%n_neg + nneg_f
1679 call mf_pool_wait(pool)
1681 do j = fct%cptr(s), fct%cptr(s+1) - 1
1683 if (
allocated(fct%sn(c)%cval))
then
1684 fct%live_cb = fct%live_cb - fct%sn(c)%cbsize
1685 deallocate(fct%sn(c)%cval)
1688 fct%live_front = fct%live_front - fw
1693 function fso(i0, j0)
result(r)
1694 integer(kind=kint),
intent(in) :: i0, j0
1695 integer(kind=kint) :: r
1696 r = hecmw_mf_dist_fsrank(map, s, pr2, pc2, i0, j0)
1700 function pslot(r0)
result(is0)
1701 integer(kind=kint),
intent(in) :: r0
1702 integer(kind=kint) :: is0, i0
1705 if (plist(i0) == r0) is0 = i0
1710 function fsbase(i0, j0)
result(off)
1711 integer(kind=kint),
intent(in) :: i0, j0
1712 integer(kind=8) :: off
1713 off = f2of(mf_bidx(g%ntc, j0, i0))
1717 function fsx(r, c)
result(ix)
1718 integer(kind=kint),
intent(in) :: r, c
1719 integer(kind=8) :: ix
1720 integer(kind=kint) :: ti, tj
1723 ix = f2of(mf_bidx(g%ntc, tj, ti)) + int(c - 1 - g%tb(tj-1), 8)*(g%tb(ti) - g%tb(ti-1)) + (r - g%tb(ti-1))
1727 subroutine send_hdr(c0)
1728 integer(kind=kint),
intent(in) :: c0
1729 integer(kind=kint) :: ib, dst, nt1, cnd, mm
1731 nt1 = fct%sn(c0)%nt + 1
1732 cnd = fct%sn(c0)%ncol - fct%sn(c0)%npiv
1734 ib = mf_pool_slot(pool)
1735 allocate(pool%box(ib)%hdr(5), pool%box(ib)%tl(mm))
1738 pool%box(ib)%hdr(1) = gierr
1739 if (gierr /= 0) ferr = gierr
1740 pool%box(ib)%hdr(2) = fct%sn(c0)%ncol
1741 pool%box(ib)%hdr(3) = fct%sn(c0)%npiv
1742 pool%box(ib)%hdr(4) = fct%sn(c0)%nt
1743 pool%box(ib)%hdr(5) = fct%sn(c0)%ntc
1744 pool%box(ib)%tl(1:nt1) = fct%sn(c0)%tbnd(1:nt1)
1745 pool%box(ib)%tl(nt1+1:nt1+cnd) = fct%sn(c0)%fsdof(fct%sn(c0)%npiv+1:fct%sn(c0)%ncol)
1746 pool%box(ib)%tl(nt1+cnd+1:nt1+2*cnd) = fct%sn(c0)%frow(fct%sn(c0)%npiv+1:fct%sn(c0)%ncol)
1747 do dst = map%rbeg(s), map%rbeg(s) + map%rcnt(s) - 1
1748 if (dst == me) cycle
1750 call hecmw_isend_int(pool%box(ib)%hdr, 5, dst, 16*c0, map%comm, req)
1751 call mf_pool_req(pool, req)
1752 call hecmw_isend_int(pool%box(ib)%tl, mm, dst, 16*c0+1, map%comm, req)
1753 call mf_pool_req(pool, req)
1755 end subroutine send_hdr
1758 subroutine recv_hdr(c0)
1759 integer(kind=kint),
intent(in) :: c0
1760 integer(kind=kint) :: hdr(5), nt1, cnd, mm
1761 integer(kind=kint),
allocatable :: tl(:)
1763 call hecmw_recv_int(hdr, 5, map%owner(c0), 16*c0, map%comm, stat)
1764 if (hdr(1) /= 0)
then
1768 fct%sn(c0)%ncol = hdr(2)
1769 fct%sn(c0)%npiv = hdr(3)
1770 fct%sn(c0)%nt = hdr(4)
1771 fct%sn(c0)%ntc = hdr(5)
1773 cnd = hdr(2) - hdr(3)
1774 call mf_grow_i(fct%sn(c0)%tbnd, nt1)
1775 call mf_grow_i(fct%sn(c0)%fsdof, max(hdr(2), 1))
1776 call mf_grow_i(fct%sn(c0)%frow, max(hdr(2), 1))
1779 call hecmw_recv_int(tl, mm, map%owner(c0), 16*c0+1, map%comm, stat)
1780 fct%sn(c0)%tbnd(1:nt1) = tl(1:nt1)
1781 fct%sn(c0)%fsdof(hdr(3)+1:hdr(2)) = tl(nt1+1:nt1+cnd)
1782 fct%sn(c0)%frow(hdr(3)+1:hdr(2)) = tl(nt1+cnd+1:nt1+2*cnd)
1784 end subroutine recv_hdr
1790 subroutine fct_sn_meta()
1791 integer(kind=kint) :: mm, kk, ii, cc, jj0
1793 allocate(blk(max(ncol, 1)), blkr(max(ncol, 1)))
1794 fct%sn(s)%ncol = ncol
1796 fct%sn(s)%ntc = g%ntc
1797 call mf_grow_i(fct%sn(s)%tbnd, g%nt + 1)
1798 fct%sn(s)%tbnd(1:g%nt+1) = g%tb(0:g%nt)
1799 call mf_grow_i(fct%sn(s)%fsdof, max(ncol, 1))
1800 call mf_grow_i(fct%sn(s)%frow, max(ncol, 1))
1801 call mf_grow_i(fct%sn(s)%ptype, max(ncol, 1))
1802 call mf_grow_r(fct%sn(s)%dsub, int(max(ncol, 1), 8))
1804 fct%sn(s)%ptype(1:ncol) = 1
1806 fct%sn(s)%ptype(1:ncol) = 0
1808 fct%sn(s)%dsub(1:ncol) = 0.0d0
1810 do kk = sym%sptr(s), sym%sptr(s+1) - 1
1811 do ii = fct%cdofptr(kk), fct%cdofptr(kk+1) - 1
1813 fct%sn(s)%fsdof(mm) = ii
1814 fct%sn(s)%frow(mm) = ii
1820 do jj0 = fct%sn(cc)%npiv + 1, fct%sn(cc)%ncol
1822 fct%sn(s)%fsdof(mm) = fct%sn(cc)%fsdof(jj0)
1823 fct%sn(s)%frow(mm) = fct%sn(cc)%frow(jj0)
1828 blkr(1:ncol) = blk(1:ncol)
1829 end subroutine fct_sn_meta
1835 function edest(pr, pc)
result(dr)
1836 integer(kind=kint),
intent(in) :: pr, pc
1837 integer(kind=kint) :: dr, rr, cc, gp
1840 if (rr <= ncol)
then
1841 if (pr2*pc2 == 1)
then
1844 gp = mod(g%dtile(rr-1) - 1, pr2)*pc2 + mod(g%dtile(cc-1) - 1, pc2)
1848 dr = map%rbeg(s) + gp - 1
1849 if (dr >= master) dr = dr + 1
1853 dr = towner(crow2t(rr - ncol))
1858 function bandix(t0, row, col)
result(ix)
1859 integer(kind=kint),
intent(in) :: t0, row, col
1860 integer(kind=8) :: ix
1861 ix = boff(t0) + int(col - 1, 8)*(ctb(t0) - ctb(t0-1)) + (row - ncol - ctb(t0-1))
1866 subroutine scatter_mine()
1867 integer(kind=kint) :: ii, kk, j0, coff, aa, bb, ki, roff, kk2
1868 integer(kind=8) :: bs, bs2
1872 kk = sym%sptr(s) + ii - 1
1875 bs = int(j0-1, 8)*nd*nd
1878 call put1(coff+aa, coff+bb, hecmat%D(bs + (aa-1)*nd + bb), .false.)
1882 call put1(coff+bb, coff+aa, hecmat%D(bs + (aa-1)*nd + bb), .true.)
1886 do kk2 = hecmat%indexL(j0-1)+1, hecmat%indexL(j0)
1887 ki = sym%invp(hecmat%itemL(kk2))
1889 roff = rowoff(pos(ki))
1890 bs = int(kk2-1, 8)*nd*nd
1891 if (.not. fct%lu)
then
1894 call put1(roff+aa, coff+bb, hecmat%AL(bs + (bb-1)*nd + aa), .false.)
1898 bs2 = int(fct%mirror(kk2)-1, 8)*nd*nd
1901 call put1(roff+bb, coff+aa, hecmat%AL(bs + (aa-1)*nd + bb), .true.)
1902 call put1(roff+bb, coff+aa, hecmat%AU(bs2 + (bb-1)*nd + aa), .false.)
1907 do kk2 = hecmat%indexU(j0-1)+1, hecmat%indexU(j0)
1908 ki = sym%invp(hecmat%itemU(kk2))
1910 roff = rowoff(pos(ki))
1911 bs = int(kk2-1, 8)*nd*nd
1912 if (.not. fct%lu)
then
1915 call put1(roff+aa, coff+bb, hecmat%AU(bs + (bb-1)*nd + aa), .false.)
1919 bs2 = int(fct%mirroru(kk2)-1, 8)*nd*nd
1922 call put1(roff+bb, coff+aa, hecmat%AU(bs + (aa-1)*nd + bb), .true.)
1923 call put1(roff+bb, coff+aa, hecmat%AL(bs2 + (bb-1)*nd + aa), .false.)
1930 end subroutine scatter_mine
1934 subroutine put1(pr, pc, v0, upface)
1935 integer(kind=kint),
intent(in) :: pr, pc
1936 real(kind=kreal),
intent(in) :: v0
1937 logical,
intent(in) :: upface
1938 integer(kind=kint) :: rr, cc, t0
1939 integer(kind=8) :: ix
1944 up = fct%lu .and. (upface .neqv. (pr < pc))
1945 if (rr <= ncol)
then
1946 if (edest(rr, cc) /= me)
return
1949 fsvu(ix) = fsvu(ix) + v0
1951 fsv(ix) = fsv(ix) + v0
1954 t0 = crow2t(rr - ncol)
1955 if (towner(t0) /= me)
return
1956 ix = bandix(t0, rr, cc)
1958 bandu(ix) = bandu(ix) + v0
1960 band(ix) = band(ix) + v0
1968 subroutine ea_setup(c0, jc)
1969 integer(kind=kint),
intent(in) :: c0, jc
1970 integer(kind=kint) :: jj0, rr, ll, kk, aa, ntc0, nde, nt2
1972 cndel = fct%sn(c0)%ncol - fct%sn(c0)%npiv
1978 cown(1) = map%owner(c0)
1981 ntc0 = fct%sn(c0)%ntc
1982 do jj0 = ntc0 + 1, fct%sn(c0)%nt
1984 cctb(cnb) = cndel + fct%sn(c0)%tbnd(jj0+1) - fct%sn(c0)%ncol
1986 if (map%upper(c0))
then
1987 nde = fct%sn(c0)%ncol - (fct%cdofptr(sym%sptr(c0+1)) - fct%cdofptr(sym%sptr(c0)))
1988 call hecmw_mf_dist_rowtiles(sym, map, fct%tile, c0, nde, nt2, cctb2, cown2)
1989 do jj0 = kbegc + 1, cnb
1990 cown(jj0) = cown2(jj0 - kbegc)
1993 do jj0 = kbegc + 1, cnb
1994 cown(jj0) = map%owner(c0)
2000 pmap(rr) = mdel(jc) + rr
2003 do ll = sym%rptr(c0) + (sym%sptr(c0+1) - sym%sptr(c0)), sym%rptr(c0+1) - 1
2004 kk = sym%cmap(sym%cmap_ptr(c0) + (ll - sym%rptr(c0) - (sym%sptr(c0+1) - sym%sptr(c0))))
2005 do aa = 1, sym%ndof(sym%rlist(ll))
2007 pmap(rr) = rowoff(kk) + aa
2010 end subroutine ea_setup
2014 subroutine ea_send(c0)
2015 integer(kind=kint),
intent(in) :: c0
2016 integer(kind=kint) :: dst, nlo, nup, ib
2018 if (.not.
allocated(fct%sn(c0)%cval))
return
2019 do dst = map%rbeg(s), map%rbeg(s) + map%rcnt(s) - 1
2020 if (dst == me) cycle
2021 call ea_count(c0, me, dst, nlo, nup)
2023 ib = mf_pool_slot(pool)
2024 allocate(pool%box(ib)%rv(nlo + nup))
2025 call ea_pack(c0, me, dst, pool%box(ib)%rv, nlo)
2027 call hecmw_isend_r(pool%box(ib)%rv, nlo + nup, dst, 16*c0+2, map%comm, req)
2028 call mf_pool_req(pool, req)
2030 end subroutine ea_send
2034 subroutine ea_recv_apply(c0)
2035 integer(kind=kint),
intent(in) :: c0
2036 integer(kind=kint) :: is0, src, nlo, nup, nsrc, jj0
2037 integer(kind=kint) :: srcs(cnb + 1)
2038 real(kind=kreal),
allocatable :: rb(:)
2044 if (srcs(is0) == cown(jj0))
exit
2046 if (is0 > nsrc)
then
2048 srcs(nsrc) = cown(jj0)
2053 call ea_count(c0, src, me, nlo, nup)
2056 call ea_apply_local(c0)
2058 allocate(rb(nlo + nup))
2059 call hecmw_recv_r(rb, nlo + nup, src, 16*c0+2, map%comm, stat)
2060 call ea_apply(c0, src, rb, nlo)
2064 end subroutine ea_recv_apply
2067 subroutine ea_count(c0, src, dst, nlo, nup)
2068 integer(kind=kint),
intent(in) :: c0, src, dst
2069 integer(kind=kint),
intent(out) :: nlo, nup
2070 integer(kind=kint) :: jj0, ii0, cc, rr, hr
2076 if (cown(ii0) /= src) cycle
2077 do cc = cctb(jj0-1) + 1, cctb(jj0)
2078 do rr = max(cc, cctb(ii0-1) + 1), cctb(ii0)
2079 hr = edest(pmap(rr), pmap(cc))
2080 if (hr /= dst) cycle
2082 if (fct%lu .and. rr > cc) nup = nup + 1
2087 end subroutine ea_count
2091 subroutine ea_pack(c0, src, dst, buf, nlo)
2092 integer(kind=kint),
intent(in) :: c0, src, dst, nlo
2093 real(kind=kreal),
intent(inout) :: buf(:)
2094 integer(kind=kint) :: jj0, ii0, cc, rr, hr, hh, il, iu
2095 integer(kind=8) :: sb, halfc, off
2097 halfc = ea_src_words(c0, src)
2103 if (cown(ii0) /= src) cycle
2104 hh = cctb(ii0) - cctb(ii0-1)
2105 do cc = cctb(jj0-1) + 1, cctb(jj0)
2106 do rr = max(cc, cctb(ii0-1) + 1), cctb(ii0)
2107 hr = edest(pmap(rr), pmap(cc))
2108 if (hr /= dst) cycle
2109 off = sb + int(cc - cctb(jj0-1) - 1, 8)*hh + (rr - cctb(ii0-1))
2111 buf(il) = fct%sn(c0)%cval(off)
2112 if (fct%lu .and. rr > cc)
then
2114 buf(nlo + iu) = fct%sn(c0)%cval(halfc + off)
2118 sb = sb + int(cctb(jj0) - cctb(jj0-1), 8)*hh
2121 end subroutine ea_pack
2124 subroutine ea_apply(c0, src, buf, nlo)
2125 integer(kind=kint),
intent(in) :: c0, src, nlo
2126 real(kind=kreal),
intent(in) :: buf(:)
2127 integer(kind=kint) :: jj0, ii0, cc, rr, hr, il, iu
2133 if (cown(ii0) /= src) cycle
2134 do cc = cctb(jj0-1) + 1, cctb(jj0)
2135 do rr = max(cc, cctb(ii0-1) + 1), cctb(ii0)
2136 hr = edest(pmap(rr), pmap(cc))
2139 call put1(pmap(rr), pmap(cc), buf(il), .false.)
2140 if (fct%lu .and. rr > cc)
then
2142 call put1(pmap(cc), pmap(rr), buf(nlo + iu), .false.)
2148 end subroutine ea_apply
2151 subroutine ea_apply_local(c0)
2152 integer(kind=kint),
intent(in) :: c0
2153 integer(kind=kint) :: jj0, ii0, cc, rr, hr, hh
2154 integer(kind=8) :: sb, halfc, off
2156 halfc = ea_src_words(c0, me)
2160 if (cown(ii0) /= me) cycle
2161 hh = cctb(ii0) - cctb(ii0-1)
2162 do cc = cctb(jj0-1) + 1, cctb(jj0)
2163 do rr = max(cc, cctb(ii0-1) + 1), cctb(ii0)
2164 hr = edest(pmap(rr), pmap(cc))
2166 off = sb + int(cc - cctb(jj0-1) - 1, 8)*hh + (rr - cctb(ii0-1))
2167 call put1(pmap(rr), pmap(cc), fct%sn(c0)%cval(off), .false.)
2168 if (fct%lu .and. rr > cc)
call put1(pmap(cc), pmap(rr), fct%sn(c0)%cval(halfc + off), .false.)
2171 sb = sb + int(cctb(jj0) - cctb(jj0-1), 8)*hh
2174 end subroutine ea_apply_local
2177 function ea_src_words(c0, src)
result(ww)
2178 integer(kind=kint),
intent(in) :: c0, src
2179 integer(kind=8) :: ww
2180 integer(kind=kint) :: jj0, ii0
2185 if (cown(ii0) /= src) cycle
2186 ww = ww + int(cctb(ii0) - cctb(ii0-1), 8)*(cctb(jj0) - cctb(jj0-1))
2189 end function ea_src_words
2193 subroutine factor_dist()
2194 integer(kind=kint) :: kk, iw2, mm, t2, h2, w2
2195 integer(kind=8) :: bb
2199 if (wrank(iw2) == me) wme = iw2
2202 mm = mf_bidx(g%nt, g%ntc, g%nt)
2203 allocate(brk(mm), bslot(mm))
2206 allocate(brku(mm), bslotu(mm))
2213 w2 = g%tb(kk) - g%tb(kk-1)
2215 if (towner(t2) /= me) cycle
2216 h2 = ctb(t2) - ctb(t2-1)
2218 bb = bb + int(h2, 8)*w2 + int(w2, 8)*min(h2, w2)
2220 bb = bb + int(h2, 8)*w2 + 2_8*int(w2, 8)*min(h2, w2)
2224 call mf_grow_r(bval, max(bb, 1_8))
2225 if (fct%lu)
call mf_grow_r(bvalu, max(bb, 1_8))
2229 allocate(perm(max(ncol, 1)), permr(max(ncol, 1)), itmp(5*max(ncol, 1)), swaps(2*max(ncol, 1)))
2230 allocate(cmax(max(ncol, 1)), ptv(max(ncol, 1)), dsv(max(ncol, 1)), dgv(max(ncol, 1)), dum(1))
2231 allocate(pair_i(g%nt*(g%nt+1)/2), pair_j(g%nt*(g%nt+1)/2))
2232 allocate(fxs(max(g%ntc, 1)), fxoff(max(g%ntc, 1)), wfo(max(g%ntc, 1)), fxbuf(max(npl, 1)))
2233 call mf_grow_r(wk, int(2*nrow, 8))
2234 call mf_grow_r(pv, int(max(ncol, 1), 8))
2235 if (.not. fct%lu)
then
2237 call mf_grow_r(pvu, 1_8)
2238 call mf_grow_r(pbandu, 1_8)
2239 call mf_grow_r(pfsu, 1_8)
2240 call mf_grow_r(pfullu, 1_8)
2241 call mf_grow_r(tmatu, 1_8)
2242 call mf_grow_r(dgbu, 1_8)
2243 call mf_grow_r(tvecu, 1_8)
2251 myfsr = fsrows_col(me, kk)
2254 do while (pa <= min(g%tb(kk), nfs))
2255 pb = min(g%tb(kk), nfs)
2258 call panel_attempt(kk)
2260 if (.not. fct%lu) fct%sn(s)%ptype(pa:pb) = 1
2264 call panel_piv_dist(kk)
2271 if (nfs <= g%tb(kk))
exit
2273 if (isroot .and. gierr == 0 .and. npiv < ncol)
call root_final()
2274 if (.not. fct%lu)
then
2275 if (ismaster) fct%sn(s)%frow(1:ncol) = fct%sn(s)%fsdof(1:ncol)
2276 if (gierr == 0)
call inertia_dist()
2278 fct%sn(s)%npiv = npiv
2279 end subroutine factor_dist
2282 function fsrows_col(r0, kk0)
result(nr0)
2283 integer(kind=kint),
intent(in) :: r0, kk0
2284 integer(kind=kint) :: nr0, i0
2287 do i0 = kk0 + 1, g%ntc
2288 if (fso(i0, kk0) == r0) nr0 = nr0 + g%tb(i0) - g%tb(i0-1)
2290 end function fsrows_col
2293 function isworker(r0)
result(yes)
2294 integer(kind=kint),
intent(in) :: r0
2296 integer(kind=kint) :: iw0
2300 if (wrank(iw0) == r0) yes = .true.
2302 end function isworker
2306 subroutine chk_build(kk)
2307 integer(kind=kint),
intent(in) :: kk
2308 integer(kind=kint) :: i0
2312 if (fsrows_col(plist(i0), kk) > 0 .or. isworker(plist(i0)))
then
2314 chk(nchk) = plist(i0)
2317 end subroutine chk_build
2320 function in_chk(r0)
result(yes)
2321 integer(kind=kint),
intent(in) :: r0
2323 integer(kind=kint) :: ic0
2327 if (chk(ic0) == r0) yes = .true.
2334 subroutine build_tmat()
2335 integer(kind=kint) :: jj2, xx, q
2336 real(kind=kreal) :: l1, l2, d11, d21, d22
2338 if (pk0 <= 0)
return
2339 call mf_grow_r(tmat, int(w, 8)*pk0)
2340 if (fct%lu)
call mf_grow_r(tmatu, int(w, 8)*pk0)
2343 if (.not. fct%lu)
then
2345 do while (q <= npiv)
2346 if (fct%sn(s)%ptype(q) == 1)
then
2347 tmat((jj2-1)*pk0 + (q-p0)) = fsv(fsx(xx, q)) * fsv(fsx(q, q))
2350 l1 = fsv(fsx(xx, q))
2351 l2 = fsv(fsx(xx, q+1))
2352 d11 = fsv(fsx(q, q))
2353 d21 = fct%sn(s)%dsub(q)
2354 d22 = fsv(fsx(q+1, q+1))
2355 tmat((jj2-1)*pk0 + (q-p0)) = l1*d11 + l2*d21
2356 tmat((jj2-1)*pk0 + (q-p0+1)) = l1*d21 + l2*d22
2362 tmat((jj2-1)*pk0 + (q-p0)) = fsvu(fsx(xx, q))
2363 tmatu((jj2-1)*pk0 + (q-p0)) = fsv(fsx(xx, q))
2367 end subroutine build_tmat
2372 subroutine fill_rd_panel(kk)
2373 integer(kind=kint),
intent(in) :: kk
2374 integer(kind=kint) :: jj2, xx, q, h2
2375 integer(kind=8) :: px, bx
2377 h2 = g%tb(kk) - g%tb(kk-1)
2378 call mf_grow_r(pv, int(mfs, 8)*w)
2379 if (fct%lu)
call mf_grow_r(pvu, int(mfs, 8)*w)
2382 px = int(jj2-1, 8)*mfs + jj2 - 1
2383 pv(px+1:px+mfs-jj2+1) = fsv(fsx(xx, xx):fsx(g%tb(kk), xx))
2386 pv(px+1:px+mfs-jj2+1) = pv(px+1:px+mfs-jj2+1) - fsv(bx+1:bx+mfs-jj2+1) * tmat((jj2-1)*pk0 + (q-p0))
2388 if (fct%lu .and. xx < g%tb(kk))
then
2389 pvu(px+2:px+mfs-jj2+1) = fsvu(fsx(xx+1, xx):fsx(g%tb(kk), xx))
2391 bx = fsx(xx+1, q) - 1
2392 pvu(px+2:px+mfs-jj2+1) = pvu(px+2:px+mfs-jj2+1) - fsvu(bx+1:bx+mfs-jj2) * tmatu((jj2-1)*pk0 + (q-p0))
2396 end subroutine fill_rd_panel
2399 subroutine rd_writeback(nc0)
2400 integer(kind=kint),
intent(in) :: nc0
2401 integer(kind=kint) :: jj2, xx
2402 integer(kind=8) :: px
2406 px = int(jj2-1, 8)*mfs + jj2 - 1
2407 fsv(fsx(xx, xx):fsx(g%tb(g%dtile(xx-1)), xx)) = pv(px+1:px+mfs-jj2+1)
2408 if (fct%lu .and. xx < g%tb(g%dtile(xx-1))) &
2409 fsvu(fsx(xx+1, xx):fsx(g%tb(g%dtile(xx-1)), xx)) = pvu(px+2:px+mfs-jj2+1)
2411 end subroutine rd_writeback
2416 subroutine fill_fs_panel(kk, solveit)
2417 integer(kind=kint),
intent(in) :: kk
2418 logical,
intent(in) :: solveit
2419 integer(kind=kint) :: jj2, t2, h2, r0, i2
2420 integer(kind=8) :: bx, px, ex
2422 call mf_grow_r(pfs, int(myfsr, 8)*w)
2423 if (fct%lu)
call mf_grow_r(pfsu, int(myfsr, 8)*w)
2426 do t2 = kk + 1, g%ntc
2427 if (fso(t2, kk) /= me) cycle
2428 h2 = g%tb(t2) - g%tb(t2-1)
2429 bx = fsbase(t2, kk) + int(pa + jj2 - 2 - g%tb(kk-1), 8)*h2
2430 ex = fsbase(t2, kk) + int(p0 - g%tb(kk-1), 8)*h2
2431 px = int(jj2-1, 8)*myfsr + r0
2432 pfs(px+1:px+h2) = fsv(bx+1:bx+h2)
2433 if (fct%lu) pfsu(px+1:px+h2) = fsvu(bx+1:bx+h2)
2435 call hecmw_mf_kernel_gemv(h2, pk0, h2, fsv(ex+1), tmat((jj2-1)*pk0+1), pfs(px+1))
2436 if (fct%lu)
call hecmw_mf_kernel_gemv(h2, pk0, h2, fsvu(ex+1), tmatu((jj2-1)*pk0+1), pfsu(px+1))
2442 if (.not. fct%lu)
then
2443 call hecmw_mf_kernel_trsm(myfsr, w, w, dgb, myfsr, pfs)
2445 call hecmw_mf_kernel_trsm_rt(myfsr, w, w, dgbu, .false., myfsr, pfs)
2446 call hecmw_mf_kernel_trsm_rt(myfsr, w, w, dgb, .true., myfsr, pfsu)
2449 px = int(jj2-1, 8)*myfsr
2451 if (.not. (abs(pfs(px+i2)) <= uinv)) cmax(jj2) = 1.0d0
2455 end subroutine fill_fs_panel
2458 subroutine fs_panel_writeback(kk, nc0)
2459 integer(kind=kint),
intent(in) :: kk, nc0
2460 integer(kind=kint) :: jj2, t2, h2, r0
2461 integer(kind=8) :: bx, px
2465 do t2 = kk + 1, g%ntc
2466 if (fso(t2, kk) /= me) cycle
2467 h2 = g%tb(t2) - g%tb(t2-1)
2468 bx = fsbase(t2, kk) + int(pa + jj2 - 2 - g%tb(kk-1), 8)*h2
2469 px = int(jj2-1, 8)*myfsr + r0
2470 fsv(bx+1:bx+h2) = pfs(px+1:px+h2)
2471 if (fct%lu) fsvu(bx+1:bx+h2) = pfsu(px+1:px+h2)
2475 end subroutine fs_panel_writeback
2480 subroutine fill_band_panel(solveit)
2481 logical,
intent(in) :: solveit
2482 integer(kind=kint) :: jj2, t2, h2, r0, i2
2483 integer(kind=8) :: bx, px
2485 call mf_grow_r(pband, int(myrows, 8)*w)
2486 if (fct%lu)
call mf_grow_r(pbandu, int(myrows, 8)*w)
2489 if (towner(t2) /= me) cycle
2490 h2 = ctb(t2) - ctb(t2-1)
2492 bx = boff(t2) + int(pa + jj2 - 2, 8)*h2
2493 px = int(jj2-1, 8)*myrows + r0
2494 pband(px+1:px+h2) = band(bx+1:bx+h2)
2495 if (fct%lu) pbandu(px+1:px+h2) = bandu(bx+1:bx+h2)
2497 call hecmw_mf_kernel_gemv(h2, pk0, h2, band(boff(t2) + int(p0, 8)*h2 + 1), &
2498 tmat((jj2-1)*pk0+1), pband(px+1))
2499 if (fct%lu)
call hecmw_mf_kernel_gemv(h2, pk0, h2, bandu(boff(t2) + int(p0, 8)*h2 + 1), &
2500 tmatu((jj2-1)*pk0+1), pbandu(px+1))
2505 if (.not. fct%lu)
then
2506 call hecmw_mf_kernel_trsm(myrows, w, w, dgb, myrows, pband)
2508 call hecmw_mf_kernel_trsm_rt(myrows, w, w, dgbu, .false., myrows, pband)
2509 call hecmw_mf_kernel_trsm_rt(myrows, w, w, dgb, .true., myrows, pbandu)
2512 px = int(jj2-1, 8)*myrows
2514 if (.not. (abs(pband(px+i2)) <= uinv)) cmax(jj2) = 1.0d0
2518 end subroutine fill_band_panel
2521 subroutine band_writeback()
2522 integer(kind=kint) :: jj2, t2, h2, r0
2523 integer(kind=8) :: bx, px
2527 if (towner(t2) /= me) cycle
2528 h2 = ctb(t2) - ctb(t2-1)
2530 bx = boff(t2) + int(pa + jj2 - 2, 8)*h2
2531 px = int(jj2-1, 8)*myrows + r0
2532 band(bx+1:bx+h2) = pband(px+1:px+h2)
2533 if (fct%lu) bandu(bx+1:bx+h2) = pbandu(px+1:px+h2)
2536 end subroutine band_writeback
2539 subroutine ctl_send_i(vals, n0)
2540 integer(kind=kint),
intent(in) :: vals(:)
2541 integer(kind=kint),
intent(in) :: n0
2542 integer(kind=kint) :: ib, i0
2544 ib = mf_pool_slot(pool)
2545 allocate(pool%box(ib)%hdr(n0))
2546 pool%box(ib)%hdr(1:n0) = vals(1:n0)
2548 if (plist(i0) == me) cycle
2550 call hecmw_isend_int(pool%box(ib)%hdr, n0, plist(i0), 16*s+3, map%comm, req)
2551 call mf_pool_req(pool, req)
2553 end subroutine ctl_send_i
2556 subroutine rd_send_i(vals, n0)
2557 integer(kind=kint),
intent(in) :: vals(:)
2558 integer(kind=kint),
intent(in) :: n0
2559 integer(kind=kint) :: ib, i0
2561 ib = mf_pool_slot(pool)
2562 allocate(pool%box(ib)%hdr(n0))
2563 pool%box(ib)%hdr(1:n0) = vals(1:n0)
2565 if (plist(i0) == me) cycle
2567 call hecmw_isend_int(pool%box(ib)%hdr, n0, plist(i0), 16*s+13, map%comm, req)
2568 call mf_pool_req(pool, req)
2570 end subroutine rd_send_i
2573 subroutine rd_send_box(ib, n0)
2574 integer(kind=kint),
intent(in) :: ib, n0
2575 integer(kind=kint) :: ic0
2578 if (chk(ic0) == me) cycle
2580 call hecmw_isend_r(pool%box(ib)%rv, n0, chk(ic0), 16*s+14, map%comm, req)
2581 call mf_pool_req(pool, req)
2583 end subroutine rd_send_box
2589 subroutine panel_attempt(kk)
2590 integer(kind=kint),
intent(in) :: kk
2591 real(kind=kreal),
allocatable :: rcm(:)
2592 integer(kind=kint) :: jj2, ic2, mm, ib, fac
2593 integer(kind=8) :: ox
2597 mfs = g%tb(kk) - pa + 1
2599 call mf_grow_r(dgb, int(w, 8)*w)
2600 if (fct%lu)
call mf_grow_r(dgbu, int(w, 8)*w)
2603 call fill_rd_panel(kk)
2604 if (.not. fct%lu)
then
2605 call hecmw_mf_kernel_panel_nopiv(mfs, w, mfs, pv, fct%pivot_u, zero, info)
2607 call hecmw_mf_kernel_panel_lu_nopiv(mfs, w, mfs, pv, mfs, pvu, fct%pivot_u, zero, info)
2608 if (info == 0 .and. mfs == w)
then
2610 pvu(int(jj2-1, 8)*mfs + jj2) = pv(int(jj2-1, 8)*mfs + jj2)
2617 dgb((jj2-1)*w+1:jj2*w) = pv(int(jj2-1, 8)*mfs + 1 : int(jj2-1, 8)*mfs + w)
2618 if (fct%lu) dgbu((jj2-1)*w+1:jj2*w) = pvu(int(jj2-1, 8)*mfs + 1 : int(jj2-1, 8)*mfs + w)
2626 ctli(5) = merge(1, 0, fsok)
2627 call rd_send_i(ctli, 8)
2629 if (fsok) mm = mm + fac*w*w
2630 if (mm > 0 .and. nchk > merge(1, 0, in_chk(me)))
then
2631 ib = mf_pool_slot(pool)
2632 allocate(pool%box(ib)%rv(mm))
2635 pool%box(ib)%rv(ox+1:ox+int(w, 8)*pk0) = tmat(1:int(w, 8)*pk0)
2636 ox = ox + int(w, 8)*pk0
2638 pool%box(ib)%rv(ox+1:ox+int(w, 8)*pk0) = tmatu(1:int(w, 8)*pk0)
2639 ox = ox + int(w, 8)*pk0
2643 pool%box(ib)%rv(ox+1:ox+int(w, 8)*w) = dgb(1:int(w, 8)*w)
2644 ox = ox + int(w, 8)*w
2646 pool%box(ib)%rv(ox+1:ox+int(w, 8)*w) = dgbu(1:int(w, 8)*w)
2649 call rd_send_box(ib, mm)
2652 call hecmw_recv_int(ctli, 8, rd, 16*s+13, map%comm, stat)
2657 fsok = (ctli(5) == 1)
2658 if (in_chk(me))
then
2660 if (fsok) mm = mm + fac*w*w
2663 call hecmw_recv_r(rcm, mm, rd, 16*s+14, map%comm, stat)
2666 call mf_grow_r(tmat, int(w, 8)*pk0)
2667 tmat(1:int(w, 8)*pk0) = rcm(ox+1:ox+int(w, 8)*pk0)
2668 ox = ox + int(w, 8)*pk0
2670 call mf_grow_r(tmatu, int(w, 8)*pk0)
2671 tmatu(1:int(w, 8)*pk0) = rcm(ox+1:ox+int(w, 8)*pk0)
2672 ox = ox + int(w, 8)*pk0
2676 dgb(1:int(w, 8)*w) = rcm(ox+1:ox+int(w, 8)*w)
2677 ox = ox + int(w, 8)*w
2678 if (fct%lu) dgbu(1:int(w, 8)*w) = rcm(ox+1:ox+int(w, 8)*w)
2687 if (in_chk(me))
then
2689 if (myfsr > 0)
call fill_fs_panel(kk, .true.)
2690 if (hasband)
call fill_band_panel(.true.)
2694 if (in_chk(me))
then
2696 if (cmax(jj2) /= 0.0d0) acc = .false.
2701 if (chk(ic2) == me) cycle
2702 call hecmw_recv_r(rcm, w, chk(ic2), 16*s+5, map%comm, stat)
2704 if (rcm(jj2) /= 0.0d0) acc = .false.
2710 ctli(2) = merge(1, 0, acc)
2711 call ctl_send_i(ctli, 8)
2713 if (in_chk(me))
then
2714 ib = mf_pool_slot(pool)
2715 allocate(pool%box(ib)%rv(w))
2716 pool%box(ib)%rv(1:w) = cmax(1:w)
2718 call hecmw_isend_r(pool%box(ib)%rv, w, master, 16*s+5, map%comm, req)
2719 call mf_pool_req(pool, req)
2721 call hecmw_recv_int(ctli, 8, master, 16*s+3, map%comm, stat)
2722 acc = (ctli(2) == 1)
2726 if (rd == me)
call rd_writeback(w)
2727 if (myfsr > 0)
call fs_panel_writeback(kk, w)
2728 if (hasband)
call band_writeback()
2730 end subroutine panel_attempt
2736 subroutine panel_piv_dist(kk)
2737 integer(kind=kint),
intent(in) :: kk
2738 real(kind=kreal),
allocatable :: rw(:)
2739 integer(kind=kint) :: ic2, x2, q2, mm, ib, mrows, fs1, i2, r0, nfs2
2740 integer(kind=8) :: ww
2742 mrows = nrow - pa + 1
2745 if (rd == me)
call fill_rd_panel(kk)
2746 if (myfsr > 0)
call fill_fs_panel(kk, .false.)
2747 if (hasband)
call fill_band_panel(.false.)
2749 call mf_grow_r(pfull, int(mrows, 8)*w)
2750 if (fct%lu)
call mf_grow_r(pfullu, int(mrows, 8)*w)
2751 ww = pmsg_words(kk, me, w)
2753 call mf_grow_r(swb, ww)
2754 call ppack_rows(kk, w, swb)
2755 call masterside(kk, w, me, swb, .true.)
2760 ww = pmsg_words(kk, r0, w)
2763 call hecmw_recv_r(rw, int(ww, kind=kint), r0, 16*s+5, map%comm, stat)
2764 call masterside(kk, w, r0, rw, .true.)
2767 if (.not. fct%lu)
then
2768 call hecmw_mf_kernel_panel_piv(mrows, w, mrows, pfull, blk(pa:pb), fct%pivot_u, mf_pivot_alpha, &
2769 zero, .true., wk, np, perm, fct%sn(s)%ptype(pa:pb), fct%sn(s)%dsub(pa:pb), nsw, n22, info)
2770 n2x2_f = n2x2_f + n22
2772 call hecmw_mf_kernel_panel_lu_piv(mrows, w, mrows, pfull, mrows, pfullu, blk(pa:pb), blkr(pa:pb), &
2773 fct%pivot_u, zero, .true., np, perm, permr, nsw, info)
2775 nswap_f = nswap_f + nsw
2784 swaps(2*npair-1) = x2
2785 swaps(2*npair) = nfs2
2799 call ctl_send_i(ctli, 8)
2800 mm = merge(2, 1, fct%lu)*w + 2*npair + w
2801 ib = mf_pool_slot(pool)
2802 allocate(pool%box(ib)%hdr(mm))
2803 pool%box(ib)%hdr(1:w) = perm(1:w)
2806 pool%box(ib)%hdr(q2+1:q2+w) = permr(1:w)
2809 pool%box(ib)%hdr(q2+1:q2+2*npair) = swaps(1:2*npair)
2811 pool%box(ib)%hdr(q2+1:q2+w) = fct%sn(s)%ptype(pa:pb)
2813 if (plist(ic2) == me) cycle
2815 call hecmw_isend_int(pool%box(ib)%hdr, mm, plist(ic2), 16*s+3, map%comm, req)
2816 call mf_pool_req(pool, req)
2820 if (.not. fct%lu .and. np > 0 .and. rd /= me)
then
2821 ib = mf_pool_slot(pool)
2822 allocate(pool%box(ib)%rv(np))
2823 pool%box(ib)%rv(1:np) = fct%sn(s)%dsub(pa:pa+np-1)
2825 call hecmw_isend_r(pool%box(ib)%rv, np, rd, 16*s+4, map%comm, req)
2826 call mf_pool_req(pool, req)
2830 ww = pmsg_words(kk, me, w)
2832 ib = mf_pool_slot(pool)
2833 allocate(pool%box(ib)%rv(ww))
2834 call ppack_rows(kk, w, pool%box(ib)%rv)
2836 call hecmw_isend_r(pool%box(ib)%rv, int(ww, kind=kint), master, 16*s+5, map%comm, req)
2837 call mf_pool_req(pool, req)
2839 call hecmw_recv_int(ctli, 8, master, 16*s+3, map%comm, stat)
2842 mm = merge(2, 1, fct%lu)*w + 2*npair + w
2843 call hecmw_recv_int(itmp, mm, master, 16*s+3, map%comm, stat)
2844 perm(1:w) = itmp(1:w)
2847 permr(1:w) = itmp(q2+1:q2+w)
2850 swaps(1:2*npair) = itmp(q2+1:q2+2*npair)
2852 fct%sn(s)%ptype(pa:pb) = itmp(q2+1:q2+w)
2853 if (.not. fct%lu .and. np > 0 .and. rd == me)
then
2854 call hecmw_recv_r(fct%sn(s)%dsub(pa:pa+np-1), np, master, 16*s+4, map%comm, stat)
2858 call perm_dance_apply()
2865 ww = pmsg_words(kk, r0, np)
2867 ib = mf_pool_slot(pool)
2868 allocate(pool%box(ib)%rv(ww))
2869 call masterside(kk, np, r0, pool%box(ib)%rv, .false.)
2871 call hecmw_isend_r(pool%box(ib)%rv, int(ww, kind=kint), r0, 16*s+4, map%comm, req)
2872 call mf_pool_req(pool, req)
2874 ww = pmsg_words(kk, me, np)
2876 call mf_grow_r(swb, ww)
2877 call masterside(kk, np, me, swb, .false.)
2878 call punpack_rows(kk, np, swb)
2881 ww = pmsg_words(kk, me, np)
2884 call hecmw_recv_r(rw, int(ww, kind=kint), master, 16*s+4, map%comm, stat)
2885 call punpack_rows(kk, np, rw)
2889 if (rd == me)
call rd_writeback(np)
2890 if (myfsr > 0)
call fs_panel_writeback(kk, np)
2891 if (hasband)
call band_scatter_writeback(np)
2897 if (.not. fct%lu)
then
2898 call dswap(1, x2, q2)
2900 call dswap(2, x2, q2)
2902 if (ismaster)
call swap_meta(x2, q2)
2904 if (.not. ismaster) nfs = ctli(4)
2905 end subroutine panel_piv_dist
2909 function pmsg_words(kk, r0, nc0)
result(ww)
2910 integer(kind=kint),
intent(in) :: kk, r0, nc0
2911 integer(kind=8) :: ww
2912 integer(kind=kint) :: nr0
2914 nr0 = fsrows_col(r0, kk) + wrows(r0)
2915 ww = int(nr0, 8)*nc0
2916 if (r0 == rd) ww = ww + int(mfs, 8)*nc0 - int(nc0, 8)*(nc0-1)/2
2919 if (r0 == rd) ww = ww - nc0
2921 end function pmsg_words
2925 subroutine ppack_rows(kk, nc0, buf)
2926 integer(kind=kint),
intent(in) :: kk, nc0
2927 real(kind=kreal),
intent(inout) :: buf(:)
2928 integer(kind=8) :: ox
2931 call pface_rows(kk, nc0, pv, pfs, pband, buf, ox, .true., .false.)
2932 if (fct%lu)
call pface_rows(kk, nc0, pvu, pfsu, pbandu, buf, ox, .true., .true.)
2933 end subroutine ppack_rows
2936 subroutine punpack_rows(kk, nc0, buf)
2937 integer(kind=kint),
intent(in) :: kk, nc0
2938 real(kind=kreal),
intent(inout) :: buf(:)
2939 integer(kind=8) :: ox
2942 call pface_rows(kk, nc0, pv, pfs, pband, buf, ox, .false., .false.)
2943 if (fct%lu)
call pface_rows(kk, nc0, pvu, pfsu, pbandu, buf, ox, .false., .true.)
2944 end subroutine punpack_rows
2947 subroutine pface_rows(kk, nc0, pva, pfsa, pbda, buf, ox, topack, uface)
2948 integer(kind=kint),
intent(in) :: kk, nc0
2949 real(kind=kreal),
intent(inout) :: pva(:), pfsa(:), pbda(:), buf(:)
2950 integer(kind=8),
intent(inout) :: ox
2951 logical,
intent(in) :: topack, uface
2952 integer(kind=kint) :: jj2, j0, ln
2953 integer(kind=8) :: px
2958 if (uface) j0 = jj2 + 1
2961 px = int(jj2-1, 8)*mfs + j0 - 1
2963 buf(ox+1:ox+ln) = pva(px+1:px+ln)
2965 pva(px+1:px+ln) = buf(ox+1:ox+ln)
2971 px = int(jj2-1, 8)*myfsr
2973 buf(ox+1:ox+myfsr) = pfsa(px+1:px+myfsr)
2975 pfsa(px+1:px+myfsr) = buf(ox+1:ox+myfsr)
2980 px = int(jj2-1, 8)*myrows
2982 buf(ox+1:ox+myrows) = pbda(px+1:px+myrows)
2984 pbda(px+1:px+myrows) = buf(ox+1:ox+myrows)
2989 end subroutine pface_rows
2994 subroutine masterside(kk, nc0, r0, buf, tofull)
2995 integer(kind=kint),
intent(in) :: kk, nc0, r0
2996 real(kind=kreal),
intent(inout) :: buf(:)
2997 logical,
intent(in) :: tofull
2998 integer(kind=8) :: ox
3001 call pface_master(kk, nc0, r0, pfull, buf, ox, tofull, .false.)
3002 if (fct%lu)
call pface_master(kk, nc0, r0, pfullu, buf, ox, tofull, .true.)
3003 end subroutine masterside
3006 subroutine pface_master(kk, nc0, r0, pf, buf, ox, tofull, uface)
3007 integer(kind=kint),
intent(in) :: kk, nc0, r0
3008 real(kind=kreal),
intent(inout) :: pf(:), buf(:)
3009 integer(kind=8),
intent(inout) :: ox
3010 logical,
intent(in) :: tofull, uface
3011 integer(kind=kint) :: jj2, j0, ln, t2, h2, r1, mrows2, fs2
3012 integer(kind=8) :: qx
3014 mrows2 = nrow - pa + 1
3017 qx = int(jj2-1, 8)*mrows2
3020 if (uface) j0 = jj2 + 1
3024 pf(qx+j0:qx+j0+ln-1) = buf(ox+1:ox+ln)
3026 buf(ox+1:ox+ln) = pf(qx+j0:qx+j0+ln-1)
3031 do t2 = kk + 1, g%ntc
3032 if (fso(t2, kk) /= r0) cycle
3033 h2 = g%tb(t2) - g%tb(t2-1)
3034 r1 = g%tb(t2-1) - pa + 2
3036 pf(qx+r1:qx+r1+h2-1) = buf(ox+1:ox+h2)
3038 buf(ox+1:ox+h2) = pf(qx+r1:qx+r1+h2-1)
3043 if (towner(t2) /= r0) cycle
3044 h2 = ctb(t2) - ctb(t2-1)
3045 r1 = fs2 + ctb(t2-1) + 1
3047 pf(qx+r1:qx+r1+h2-1) = buf(ox+1:ox+h2)
3049 buf(ox+1:ox+h2) = pf(qx+r1:qx+r1+h2-1)
3054 end subroutine pface_master
3057 function wrows(rk)
result(nr0)
3058 integer(kind=kint),
intent(in) :: rk
3059 integer(kind=kint) :: nr0, t2
3063 if (towner(t2) == rk) nr0 = nr0 + ctb(t2) - ctb(t2-1)
3070 subroutine perm_dance_apply()
3071 integer(kind=kint) :: jj2, t2, q2
3077 if (itmp(jj2) == perm(jj2)) cycle
3079 if (itmp(t2) == perm(jj2))
exit
3081 if (.not. fct%lu)
then
3082 call dswap(1, pa+jj2-1, pa+t2-1)
3084 q2 = fct%sn(s)%fsdof(pa+jj2-1)
3085 fct%sn(s)%fsdof(pa+jj2-1) = fct%sn(s)%fsdof(pa+t2-1)
3086 fct%sn(s)%fsdof(pa+t2-1) = q2
3089 call dswap(2, pa+jj2-1, pa+t2-1)
3090 if (ismaster)
call swap_meta(pa+jj2-1, pa+t2-1)
3092 itmp(t2) = itmp(jj2)
3093 itmp(jj2) = perm(jj2)
3096 itmp(1:w) = perm(1:w)
3098 if (itmp(jj2) == permr(jj2)) cycle
3100 if (itmp(t2) == permr(jj2))
exit
3102 call dswap(3, pa+jj2-1, pa+t2-1)
3104 q2 = fct%sn(s)%frow(pa+jj2-1)
3105 fct%sn(s)%frow(pa+jj2-1) = fct%sn(s)%frow(pa+t2-1)
3106 fct%sn(s)%frow(pa+t2-1) = q2
3108 blkr(pa+jj2-1) = blkr(pa+t2-1)
3111 itmp(t2) = itmp(jj2)
3112 itmp(jj2) = permr(jj2)
3115 end subroutine perm_dance_apply
3118 subroutine swap_meta(x0, y0)
3119 integer(kind=kint),
intent(in) :: x0, y0
3120 integer(kind=kint) :: q2
3122 q2 = fct%sn(s)%fsdof(x0)
3123 fct%sn(s)%fsdof(x0) = fct%sn(s)%fsdof(y0)
3124 fct%sn(s)%fsdof(y0) = q2
3129 q2 = fct%sn(s)%frow(x0)
3130 fct%sn(s)%frow(x0) = fct%sn(s)%frow(y0)
3131 fct%sn(s)%frow(y0) = q2
3136 end subroutine swap_meta
3139 subroutine band_scatter_writeback(nc0)
3140 integer(kind=kint),
intent(in) :: nc0
3141 integer(kind=kint) :: jj2, t2, h2
3142 integer(kind=8) :: px, bx
3146 if (towner(t2) /= me) cycle
3147 h2 = ctb(t2) - ctb(t2-1)
3148 px = int(jj2-1, 8)*myrows + tprow(t2)
3149 bx = boff(t2) + int(pa + jj2 - 2, 8)*h2
3150 band(bx+1:bx+h2) = pband(px+1:px+h2)
3151 if (fct%lu) bandu(bx+1:bx+h2) = pbandu(px+1:px+h2)
3154 end subroutine band_scatter_writeback
3162 subroutine dswap(kindx, x0, y0)
3163 integer(kind=kint),
intent(in) :: kindx, x0, y0
3164 integer(kind=kint) :: npr, p2, o1, o2, nsp, i0, r0, cnt, ib
3165 real(kind=kreal) :: v0
3167 call dswap_pairs(kindx, x0, y0, npr)
3170 o1 = edest(spr(p2), spc(p2))
3171 o2 = edest(sqr(p2), sqc(p2))
3174 v0 = pget(spr(p2), spc(p2), spf(p2))
3175 call pput(spr(p2), spc(p2), spf(p2), pget(sqr(p2), sqc(p2), sqf(p2)))
3176 call pput(sqr(p2), sqc(p2), sqf(p2), v0)
3182 if (nsp == 0)
return
3184 call mf_grow_r(sw1, int(npr, 8))
3185 call mf_grow_r(sw2, int(npr, 8))
3187 o1 = edest(spr(p2), spc(p2))
3188 o2 = edest(sqr(p2), sqc(p2))
3190 if (o1 == me) sw1(p2) = pget(spr(p2), spc(p2), spf(p2))
3191 if (o2 == me) sw2(p2) = pget(sqr(p2), sqc(p2), sqf(p2))
3196 cnt = split_count(npr, r0)
3198 call mf_grow_r(swb, int(cnt, 8))
3199 call hecmw_recv_r(swb, cnt, r0, 16*s+15, map%comm, stat)
3202 o1 = edest(spr(p2), spc(p2))
3203 o2 = edest(sqr(p2), sqc(p2))
3216 o1 = edest(spr(p2), spc(p2))
3217 o2 = edest(sqr(p2), sqc(p2))
3222 if (o1 == me)
call pput(spr(p2), spc(p2), spf(p2), sw1(p2))
3223 if (o2 == me)
call pput(sqr(p2), sqc(p2), sqf(p2), sw2(p2))
3228 cnt = split_count(npr, r0)
3230 ib = mf_pool_slot(pool)
3231 allocate(pool%box(ib)%rv(cnt))
3234 o1 = edest(spr(p2), spc(p2))
3235 o2 = edest(sqr(p2), sqc(p2))
3239 pool%box(ib)%rv(cnt) = sw1(p2)
3243 pool%box(ib)%rv(cnt) = sw2(p2)
3247 call hecmw_isend_r(pool%box(ib)%rv, cnt, r0, 16*s+15, map%comm, req)
3248 call mf_pool_req(pool, req)
3251 cnt = split_count(npr, me)
3253 ib = mf_pool_slot(pool)
3254 allocate(pool%box(ib)%rv(cnt))
3257 o1 = edest(spr(p2), spc(p2))
3258 o2 = edest(sqr(p2), sqc(p2))
3262 pool%box(ib)%rv(cnt) = pget(spr(p2), spc(p2), spf(p2))
3266 pool%box(ib)%rv(cnt) = pget(sqr(p2), sqc(p2), sqf(p2))
3270 call hecmw_isend_r(pool%box(ib)%rv, cnt, master, 16*s+15, map%comm, req)
3271 call mf_pool_req(pool, req)
3272 call mf_grow_r(swb, int(cnt, 8))
3273 call hecmw_recv_r(swb, cnt, master, 16*s+15, map%comm, stat)
3276 o1 = edest(spr(p2), spc(p2))
3277 o2 = edest(sqr(p2), sqc(p2))
3281 call pput(spr(p2), spc(p2), spf(p2), swb(cnt))
3285 call pput(sqr(p2), sqc(p2), sqf(p2), swb(cnt))
3290 end subroutine dswap
3295 subroutine dswap_pairs(kindx, x0, y0, npr)
3296 integer(kind=kint),
intent(in) :: kindx, x0, y0
3297 integer(kind=kint),
intent(out) :: npr
3298 integer(kind=kint) :: i0
3300 call mf_grow_i(spr, 2*nrow + 4)
3301 call mf_grow_i(spc, 2*nrow + 4)
3302 call mf_grow_i(spf, 2*nrow + 4)
3303 call mf_grow_i(sqr, 2*nrow + 4)
3304 call mf_grow_i(sqc, 2*nrow + 4)
3305 call mf_grow_i(sqf, 2*nrow + 4)
3307 if (kindx == 1)
then
3308 call addpair(npr, x0, x0, 0, y0, y0, 0)
3310 call addpair(npr, x0, i0, 0, y0, i0, 0)
3313 call addpair(npr, i0, x0, 0, y0, i0, 0)
3316 call addpair(npr, i0, x0, 0, i0, y0, 0)
3318 else if (kindx == 2)
then
3319 call addpair(npr, x0, x0, 0, y0, y0, 0)
3321 call addpair(npr, x0, i0, 0, y0, i0, 0)
3322 call addpair(npr, x0, i0, 1, y0, i0, 1)
3325 call addpair(npr, i0, x0, 0, y0, i0, 1)
3326 call addpair(npr, y0, i0, 0, i0, x0, 1)
3328 call addpair(npr, y0, x0, 0, y0, x0, 1)
3330 call addpair(npr, i0, x0, 0, i0, y0, 0)
3331 call addpair(npr, i0, x0, 1, i0, y0, 1)
3335 call addpair(npr, x0, i0, 0, y0, i0, 0)
3338 call addpair(npr, i0, x0, 1, y0, i0, 0)
3340 call addpair(npr, y0, x0, 1, y0, y0, 0)
3342 call addpair(npr, i0, x0, 1, i0, y0, 1)
3345 end subroutine dswap_pairs
3347 subroutine addpair(npr, r1, c1, f1, r2, c2, f2)
3348 integer(kind=kint),
intent(inout) :: npr
3349 integer(kind=kint),
intent(in) :: r1, c1, f1, r2, c2, f2
3358 end subroutine addpair
3361 function split_count(npr, r0)
result(cnt)
3362 integer(kind=kint),
intent(in) :: npr, r0
3363 integer(kind=kint) :: cnt, p2, o1, o2
3367 o1 = edest(spr(p2), spc(p2))
3368 o2 = edest(sqr(p2), sqc(p2))
3370 if (o1 == r0) cnt = cnt + 1
3371 if (o2 == r0) cnt = cnt + 1
3373 end function split_count
3376 function pget(r, c, f0)
result(v0)
3377 integer(kind=kint),
intent(in) :: r, c, f0
3378 real(kind=kreal) :: v0
3379 integer(kind=kint) :: t0
3380 integer(kind=8) :: ix
3390 t0 = crow2t(r - ncol)
3391 ix = bandix(t0, r, c)
3400 subroutine pput(r, c, f0, v0)
3401 integer(kind=kint),
intent(in) :: r, c, f0
3402 real(kind=kreal),
intent(in) :: v0
3403 integer(kind=kint) :: t0
3404 integer(kind=8) :: ix
3414 t0 = crow2t(r - ncol)
3415 ix = bandix(t0, r, c)
3429 subroutine postcol(kk)
3430 integer(kind=kint),
intent(in) :: kk
3431 real(kind=kreal),
allocatable :: rcm(:)
3432 integer(kind=kint) :: q2, x2, mm, ib, ndl, fac
3433 integer(kind=8) :: ox
3437 ndl = g%tb(kk) - npiv
3440 if (.not. fct%lu)
then
3441 do q2 = p0 + 1, npiv
3442 ptv(q2-p0) = fct%sn(s)%ptype(q2)
3443 dsv(q2-p0) = fct%sn(s)%dsub(q2)
3444 dgv(q2-p0) = fsv(fsx(q2, q2))
3447 call mf_grow_r(tvec, max(int(ndl, 8)*pk, 1_8))
3448 if (fct%lu)
call mf_grow_r(tvecu, max(int(ndl, 8)*pk, 1_8))
3449 do x2 = npiv + 1, g%tb(kk)
3450 if (.not. fct%lu)
then
3451 call build_tcol(x2, tvec(int(x2-npiv-1, 8)*pk + 1))
3453 do q2 = p0 + 1, npiv
3454 tvec(int(x2-npiv-1, 8)*pk + (q2-p0)) = fsvu(fsx(x2, q2))
3455 tvecu(int(x2-npiv-1, 8)*pk + (q2-p0)) = fsv(fsx(x2, q2))
3459 if (.not. fct%lu)
call rd_send_i(ptv, pk)
3460 mm = merge(0, 2*pk, fct%lu) + fac*ndl*pk
3461 if (mm > 0 .and. npl > 1)
then
3462 ib = mf_pool_slot(pool)
3463 allocate(pool%box(ib)%rv(mm))
3465 if (.not. fct%lu)
then
3466 pool%box(ib)%rv(1:pk) = dsv(1:pk)
3467 pool%box(ib)%rv(pk+1:2*pk) = dgv(1:pk)
3470 pool%box(ib)%rv(ox+1:ox+int(ndl, 8)*pk) = tvec(1:int(ndl, 8)*pk)
3471 ox = ox + int(ndl, 8)*pk
3472 if (fct%lu) pool%box(ib)%rv(ox+1:ox+int(ndl, 8)*pk) = tvecu(1:int(ndl, 8)*pk)
3473 call rd_send_box_all(ib, mm)
3476 if (.not. fct%lu)
call hecmw_recv_int(ptv, pk, rd, 16*s+13, map%comm, stat)
3477 mm = merge(0, 2*pk, fct%lu) + fac*ndl*pk
3480 call hecmw_recv_r(rcm, mm, rd, 16*s+14, map%comm, stat)
3482 if (.not. fct%lu)
then
3483 dsv(1:pk) = rcm(1:pk)
3484 dgv(1:pk) = rcm(pk+1:2*pk)
3487 call mf_grow_r(tvec, max(int(ndl, 8)*pk, 1_8))
3488 tvec(1:int(ndl, 8)*pk) = rcm(ox+1:ox+int(ndl, 8)*pk)
3489 ox = ox + int(ndl, 8)*pk
3491 call mf_grow_r(tvecu, max(int(ndl, 8)*pk, 1_8))
3492 tvecu(1:int(ndl, 8)*pk) = rcm(ox+1:ox+int(ndl, 8)*pk)
3497 if (fct%blr .and. hasband)
call compress_band(kk)
3498 if (hasband)
call exchange_band(kk)
3499 call fs_elim_exchange(kk)
3501 end subroutine postcol
3505 subroutine rd_send_box_all(ib, n0)
3506 integer(kind=kint),
intent(in) :: ib, n0
3507 integer(kind=kint) :: i0
3510 if (plist(i0) == me) cycle
3512 call hecmw_isend_r(pool%box(ib)%rv, n0, plist(i0), 16*s+14, map%comm, req)
3513 call mf_pool_req(pool, req)
3515 end subroutine rd_send_box_all
3520 function fs_needs(r0, i0, kk0)
result(yes)
3521 integer(kind=kint),
intent(in) :: r0, i0, kk0
3523 integer(kind=kint) :: j0
3528 if (fso(i0, j0) == r0) yes = .true.
3532 if (fso(j0, i0) == r0) yes = .true.
3534 end function fs_needs
3539 subroutine fs_elim_exchange(kk)
3540 integer(kind=kint),
intent(in) :: kk
3541 integer(kind=kint) :: i0, i2, t2, h2, dst, ib, is2
3542 integer(kind=8) :: sz, ox, bx
3546 if (dst == me) cycle
3548 do t2 = kk + 1, g%ntc
3549 if (fso(t2, kk) /= me) cycle
3550 if (.not. fs_needs(dst, t2, kk)) cycle
3551 sz = sz + merge(2, 1, fct%lu)*int(g%tb(t2) - g%tb(t2-1), 8)*pk
3554 ib = mf_pool_slot(pool)
3555 allocate(pool%box(ib)%rv(sz))
3557 do t2 = kk + 1, g%ntc
3558 if (fso(t2, kk) /= me) cycle
3559 if (.not. fs_needs(dst, t2, kk)) cycle
3560 h2 = g%tb(t2) - g%tb(t2-1)
3561 bx = fsbase(t2, kk) + int(p0 - g%tb(kk-1), 8)*h2
3562 pool%box(ib)%rv(ox+1:ox+int(h2, 8)*pk) = fsv(bx+1:bx+int(h2, 8)*pk)
3563 ox = ox + int(h2, 8)*pk
3565 pool%box(ib)%rv(ox+1:ox+int(h2, 8)*pk) = fsvu(bx+1:bx+int(h2, 8)*pk)
3566 ox = ox + int(h2, 8)*pk
3570 call hecmw_isend_r(pool%box(ib)%rv, int(sz, kind=kint), dst, 16*s+14, map%comm, req)
3571 call mf_pool_req(pool, req)
3575 do i2 = kk + 1, g%ntc
3576 if (fso(i2, kk) == me)
then
3580 if (.not. fs_needs(me, i2, kk)) cycle
3581 if (fxs(i2) >= 0) cycle
3582 is2 = pslot(fso(i2, kk))
3584 do t2 = kk + 1, g%ntc
3585 if (fso(t2, kk) /= fso(i2, kk)) cycle
3586 if (.not. fs_needs(me, t2, kk)) cycle
3587 h2 = g%tb(t2) - g%tb(t2-1)
3590 sz = sz + merge(2, 1, fct%lu)*int(h2, 8)*pk
3592 if (
allocated(fxbuf(is2)%v))
then
3593 if (
size(fxbuf(is2)%v, kind=8) < sz)
deallocate(fxbuf(is2)%v)
3595 if (.not.
allocated(fxbuf(is2)%v))
allocate(fxbuf(is2)%v(max(sz, 1_8)))
3596 call hecmw_recv_r(fxbuf(is2)%v, int(sz, kind=kint), fso(i2, kk), 16*s+14, map%comm, stat)
3598 end subroutine fs_elim_exchange
3602 subroutine fs_elim_blk(t2, kk, uface, fb)
3603 integer(kind=kint),
intent(in) :: t2, kk
3604 logical,
intent(in) :: uface
3605 real(kind=kreal),
pointer,
contiguous,
intent(out) :: fb(:)
3606 integer(kind=kint) :: h2
3607 integer(kind=8) :: bx
3609 h2 = g%tb(t2) - g%tb(t2-1)
3610 if (fso(t2, kk) == me)
then
3611 bx = fsbase(t2, kk) + int(p0 - g%tb(kk-1), 8)*h2
3613 fb => fsvu(bx+1 : bx+int(h2, 8)*pk)
3615 fb => fsv(bx+1 : bx+int(h2, 8)*pk)
3619 if (uface) bx = bx + int(h2, 8)*pk
3620 fb => fxbuf(fxs(t2))%v(bx+1 : bx+int(h2, 8)*pk)
3622 end subroutine fs_elim_blk
3626 subroutine build_tcol(x2, tv)
3627 integer(kind=kint),
intent(in) :: x2
3628 real(kind=kreal),
intent(out) :: tv(*)
3629 integer(kind=kint) :: q2
3630 real(kind=kreal) :: l1, l2, d11, d21, d22
3633 do while (q2 <= npiv)
3634 if (fct%sn(s)%ptype(q2) == 1)
then
3635 tv(q2-p0) = fsv(fsx(x2, q2)) * fsv(fsx(q2, q2))
3638 l1 = fsv(fsx(x2, q2))
3639 l2 = fsv(fsx(x2, q2+1))
3640 d11 = fsv(fsx(q2, q2))
3641 d21 = fct%sn(s)%dsub(q2)
3642 d22 = fsv(fsx(q2+1, q2+1))
3643 tv(q2-p0) = l1*d11 + l2*d21
3644 tv(q2-p0+1) = l1*d21 + l2*d22
3648 end subroutine build_tcol
3652 subroutine compress_band(kk)
3653 integer(kind=kint),
intent(in) :: kk
3654 integer(kind=kint) :: t2, h2, r2, idx2, rcp
3655 logical,
allocatable :: skpl(:), skpu(:)
3659 allocate(skpl(ncbt), skpu(ncbt))
3661 if (towner(t2) /= me) cycle
3662 h2 = ctb(t2) - ctb(t2-1)
3663 idx2 = mf_bidx(g%nt, kk, g%ntc + t2)
3664 skpl(t2) = mf_blr_skip(fct, s, g%ntc, kk, g%ntc + t2)
3667 if (.not. skpl(t2)) skpl(t2) = fct%sn(s)%brank(idx2) < 0
3669 if (.not. skpu(t2)) skpu(t2) = fct%sn(s)%branku(idx2) < 0
3672 if (.not. skpl(t2))
then
3675 btop = btop + int(h2, 8)*pk + int(pk, 8)*min(h2, pk)
3677 btop = btop + int(h2, 8)*pk + 2_8*int(pk, 8)*min(h2, pk)
3680 if (fct%lu .and. .not. skpu(t2))
then
3681 bslotu(idx2) = btopu
3682 btopu = btopu + int(h2, 8)*pk + int(pk, 8)*min(h2, pk)
3685 call mf_grow_r(bval, max(btop, 1_8))
3686 if (fct%lu)
call mf_grow_r(bvalu, max(btopu, 1_8))
3689 if (towner(t2) /= me) cycle
3690 h2 = ctb(t2) - ctb(t2-1)
3691 idx2 = mf_bidx(g%nt, kk, g%ntc + t2)
3692 rcp = int(fct%blr_beta*real((h2*pk - 1)/(h2 + pk), kind=kreal))
3693 if (.not. skpl(t2))
then
3694 bval(bslot(idx2)+1:bslot(idx2)+int(h2, 8)*pk) = &
3695 band(boff(t2)+int(p0, 8)*h2+1:boff(t2)+int(p0, 8)*h2+int(h2, 8)*pk)
3696 call hecmw_mf_kernel_compress(h2, pk, h2, bval(bslot(idx2)+1), fct%eps, pk, &
3697 bval(bslot(idx2)+int(h2, 8)*pk+1), r2, rcp)
3700 if (fct%lu .and. .not. skpu(t2))
then
3701 bvalu(bslotu(idx2)+1:bslotu(idx2)+int(h2, 8)*pk) = &
3702 bandu(boff(t2)+int(p0, 8)*h2+1:boff(t2)+int(p0, 8)*h2+int(h2, 8)*pk)
3703 call hecmw_mf_kernel_compress(h2, pk, h2, bvalu(bslotu(idx2)+1), fct%eps, pk, &
3704 bvalu(bslotu(idx2)+int(h2, 8)*pk+1), r2, rcp)
3710 if (towner(t2) /= me) cycle
3711 idx2 = mf_bidx(g%nt, kk, g%ntc + t2)
3717 if (brk(idx2) >= 0)
then
3719 rsum_f = rsum_f + brk(idx2)
3720 rmax_f = max(rmax_f, brk(idx2))
3728 if (brku(idx2) >= 0)
then
3730 rsum_f = rsum_f + brku(idx2)
3731 rmax_f = max(rmax_f, brku(idx2))
3735 end subroutine compress_band
3740 subroutine exchange_band(kk)
3741 integer(kind=kint),
intent(in) :: kk
3742 integer(kind=kint),
allocatable :: irk(:)
3743 integer(kind=kint) :: t2, h2, r2, ru2, iw2, ib, nti
3744 integer(kind=8) :: sz, ox
3746 if (.not.
allocated(xw))
then
3747 allocate(xw(max(ncbt, 1)), xoff(max(ncbt, 1)), xrk(max(ncbt, 1)))
3748 allocate(xoffu(max(ncbt, 1)), xrku(max(ncbt, 1)), xbuf(max(nwk, 1)))
3755 if (towner(t2) /= me) cycle
3756 h2 = ctb(t2) - ctb(t2-1)
3758 sz = sz + xtile_words(t2, kk, .false.)
3759 if (fct%lu) sz = sz + xtile_words(t2, kk, .true.)
3761 ib = mf_pool_slot(pool)
3762 allocate(pool%box(ib)%rv(max(sz, 1_8)))
3763 if (fct%blr)
allocate(pool%box(ib)%hdr(max(merge(2, 1, fct%lu)*nti, 1)))
3767 if (towner(t2) /= me) cycle
3768 call xtile_pack(t2, kk, .false., pool%box(ib)%rv, ox)
3771 pool%box(ib)%hdr(merge(2, 1, fct%lu)*(nti-1)+1) = xrank_of(t2, kk, .false.)
3772 if (fct%lu) pool%box(ib)%hdr(merge(2, 1, fct%lu)*(nti-1)+2) = xrank_of(t2, kk, .true.)
3774 if (fct%lu)
call xtile_pack(t2, kk, .true., pool%box(ib)%rv, ox)
3776 do iw2 = wme + 1, nwk
3779 call hecmw_isend_int(pool%box(ib)%hdr, merge(2, 1, fct%lu)*max(nti, 1), wrank(iw2), 16*s+6, map%comm, req)
3780 call mf_pool_req(pool, req)
3783 call hecmw_isend_r(pool%box(ib)%rv, int(max(sz, 1_8), kind=kint), wrank(iw2), 16*s+7, map%comm, req)
3784 call mf_pool_req(pool, req)
3791 if (towner(t2) == wrank(iw2)) nti = nti + 1
3793 allocate(irk(merge(2, 1, fct%lu)*max(nti, 1)))
3795 call hecmw_recv_int(irk, merge(2, 1, fct%lu)*max(nti, 1), wrank(iw2), 16*s+6, map%comm, stat)
3802 if (towner(t2) /= wrank(iw2)) cycle
3803 h2 = ctb(t2) - ctb(t2-1)
3805 r2 = irk(merge(2, 1, fct%lu)*(nti-1)+1)
3810 sz = sz + int(h2, 8)*pk
3812 sz = sz + int(h2, 8)*r2 + int(pk, 8)*r2
3815 ru2 = irk(merge(2, 1, fct%lu)*(nti-1)+2)
3819 sz = sz + int(h2, 8)*pk
3821 sz = sz + int(h2, 8)*ru2 + int(pk, 8)*ru2
3826 if (
allocated(xbuf(iw2)%v))
then
3827 if (
size(xbuf(iw2)%v, kind=8) < sz)
deallocate(xbuf(iw2)%v)
3829 if (.not.
allocated(xbuf(iw2)%v))
allocate(xbuf(iw2)%v(max(sz, 1_8)))
3830 call hecmw_recv_r(xbuf(iw2)%v, int(max(sz, 1_8), kind=kint), wrank(iw2), 16*s+7, map%comm, stat)
3832 end subroutine exchange_band
3835 function xrank_of(t2, kk, uface)
result(r2)
3836 integer(kind=kint),
intent(in) :: t2, kk
3837 logical,
intent(in) :: uface
3838 integer(kind=kint) :: r2
3843 r2 = brku(mf_bidx(g%nt, kk, g%ntc + t2))
3845 r2 = brk(mf_bidx(g%nt, kk, g%ntc + t2))
3848 end function xrank_of
3851 function xtile_words(t2, kk, uface)
result(ww)
3852 integer(kind=kint),
intent(in) :: t2, kk
3853 logical,
intent(in) :: uface
3854 integer(kind=8) :: ww
3855 integer(kind=kint) :: h2, r2
3857 h2 = ctb(t2) - ctb(t2-1)
3858 r2 = xrank_of(t2, kk, uface)
3862 ww = int(h2, 8)*r2 + int(pk, 8)*r2
3864 end function xtile_words
3867 subroutine xtile_pack(t2, kk, uface, buf, ox)
3868 integer(kind=kint),
intent(in) :: t2, kk
3869 logical,
intent(in) :: uface
3870 real(kind=kreal),
intent(inout) :: buf(:)
3871 integer(kind=8),
intent(inout) :: ox
3872 integer(kind=kint) :: h2, r2, idx2
3873 integer(kind=8) :: bx
3875 h2 = ctb(t2) - ctb(t2-1)
3876 r2 = xrank_of(t2, kk, uface)
3878 bx = boff(t2) + int(p0, 8)*h2
3880 buf(ox+1:ox+int(h2, 8)*pk) = bandu(bx+1:bx+int(h2, 8)*pk)
3882 buf(ox+1:ox+int(h2, 8)*pk) = band(bx+1:bx+int(h2, 8)*pk)
3884 ox = ox + int(h2, 8)*pk
3886 idx2 = mf_bidx(g%nt, kk, g%ntc + t2)
3889 buf(ox+1:ox+int(h2, 8)*r2) = bvalu(bx+1:bx+int(h2, 8)*r2)
3890 buf(ox+int(h2, 8)*r2+1:ox+int(h2, 8)*r2+int(pk, 8)*r2) = &
3891 bvalu(bx+int(h2, 8)*pk+1:bx+int(h2, 8)*pk+int(pk, 8)*r2)
3894 buf(ox+1:ox+int(h2, 8)*r2) = bval(bx+1:bx+int(h2, 8)*r2)
3895 buf(ox+int(h2, 8)*r2+1:ox+int(h2, 8)*r2+int(pk, 8)*r2) = &
3896 bval(bx+int(h2, 8)*pk+1:bx+int(h2, 8)*pk+int(pk, 8)*r2)
3898 ox = ox + int(h2, 8)*r2 + int(pk, 8)*r2
3900 end subroutine xtile_pack
3905 subroutine update_col(kk)
3906 integer(kind=kint),
intent(in) :: kk
3907 integer(kind=kint) :: t2, h2, i2, j2, l2, q2, x2, mn2, r2
3908 integer(kind=8) :: ox
3909 real(kind=kreal),
pointer,
contiguous :: fb(:)
3913 if (.not. fct%lu)
then
3914 call mf_grow_r(dd, int(pk, 8)*pk)
3915 dd(1:int(pk, 8)*pk) = 0.0d0
3917 dd(int(q2-1, 8)*pk + q2) = dgv(q2)
3921 do i2 = kk + 1, g%ntc
3923 if (fso(i2, j2) == me)
then
3924 if (wfo(i2) < 0)
then
3926 ox = ox + int(g%tb(i2) - g%tb(i2-1), 8)*pk
3931 call mf_grow_r(wfs, max(ox, 1_8))
3933 do i2 = kk + 1, g%ntc
3934 if (wfo(i2) < 0) cycle
3935 h2 = g%tb(i2) - g%tb(i2-1)
3936 call fs_elim_blk(i2, kk, .false., fb)
3937 call hecmw_mf_kernel_scale(h2, pk, pk, dd, ptv, dsv, h2, fb, wfs(wfo(i2)+1))
3941 call mf_grow_r(wscb, max(int(myrows, 8)*pk, 1_8))
3944 if (towner(t2) /= me) cycle
3945 h2 = ctb(t2) - ctb(t2-1)
3946 r2 = xrank_of(t2, kk, .false.)
3948 call hecmw_mf_kernel_scale(h2, pk, pk, dd, ptv, dsv, h2, &
3949 band(boff(t2) + int(p0, 8)*h2 + 1), wscb(int(tprow(t2), 8)*pk + 1))
3950 else if (r2 > 0)
then
3952 ox = bslot(mf_bidx(g%nt, kk, g%ntc + t2))
3953 call hecmw_mf_kernel_scale_rows(pk, r2, pk, dd, ptv, dsv, pk, &
3954 bval(ox + int(h2, 8)*pk + 1), pk, bval(ox + int(h2, 8)*pk + int(pk, 8)*mn2 + 1))
3962 do i2 = kk + 1, g%ntc
3964 if (fso(i2, j2) /= me) cycle
3971 if (towner(t2) /= me) cycle
3972 do j2 = kk + 1, g%ntc + t2
3974 pair_i(npair) = g%ntc + t2
3980 call update_pair(kk, pair_i(l2), pair_j(l2))
3987 h2 = g%tb(kk) - g%tb(kk-1)
3988 do x2 = npiv + 1, g%tb(kk)
3989 ox = fsbase(kk, kk) + int(p0 - g%tb(kk-1), 8)*h2 + (x2 - 1 - g%tb(kk-1))
3990 call hecmw_mf_kernel_gemv(g%tb(kk) - x2 + 1, pk, h2, fsv(ox+1), &
3991 tvec(int(x2-npiv-1, 8)*pk + 1), fsv(fsx(x2, x2)))
3993 if (x2 < g%tb(kk))
call hecmw_mf_kernel_gemv(g%tb(kk) - x2, pk, h2, fsvu(ox+2), &
3994 tvecu(int(x2-npiv-1, 8)*pk + 1), fsvu(fsx(x2+1, x2)))
3999 do t2 = kk + 1, g%ntc
4000 if (fso(t2, kk) /= me) cycle
4001 h2 = g%tb(t2) - g%tb(t2-1)
4003 do x2 = npiv + 1, g%tb(kk)
4004 call hecmw_mf_kernel_gemv(h2, pk, h2, fsv(ox + int(p0 - g%tb(kk-1), 8)*h2 + 1), &
4005 tvec(int(x2-npiv-1, 8)*pk + 1), fsv(ox + int(x2 - 1 - g%tb(kk-1), 8)*h2 + 1))
4006 if (fct%lu)
call hecmw_mf_kernel_gemv(h2, pk, h2, fsvu(ox + int(p0 - g%tb(kk-1), 8)*h2 + 1), &
4007 tvecu(int(x2-npiv-1, 8)*pk + 1), fsvu(ox + int(x2 - 1 - g%tb(kk-1), 8)*h2 + 1))
4014 if (towner(t2) /= me) cycle
4015 h2 = ctb(t2) - ctb(t2-1)
4016 do x2 = npiv + 1, g%tb(kk)
4017 call hecmw_mf_kernel_gemv(h2, pk, h2, band(boff(t2) + int(p0, 8)*h2 + 1), &
4018 tvec(int(x2-npiv-1, 8)*pk + 1), band(boff(t2) + int(x2-1, 8)*h2 + 1))
4019 if (fct%lu)
call hecmw_mf_kernel_gemv(h2, pk, h2, bandu(boff(t2) + int(p0, 8)*h2 + 1), &
4020 tvecu(int(x2-npiv-1, 8)*pk + 1), bandu(boff(t2) + int(x2-1, 8)*h2 + 1))
4025 end subroutine update_col
4031 subroutine update_pair(kk, i2, j2)
4032 integer(kind=kint),
intent(in) :: kk, i2, j2
4033 real(kind=kreal),
pointer,
contiguous :: fa(:), ua(:), va(:), fb(:), ub(:), vb(:)
4034 integer(kind=kint) :: hi2, wj2, ra2, rb2, t2, rua2, rlb2
4035 integer(kind=8) :: ox, oc
4037 hi2 = g%tb(i2) - g%tb(i2-1)
4038 wj2 = g%tb(j2) - g%tb(j2-1)
4039 if (i2 <= g%ntc)
then
4041 if (.not. fct%lu)
then
4042 call fs_elim_blk(j2, kk, .false., fb)
4043 call hecmw_mf_kernel_gemm(hi2, wj2, pk, hi2, wfs(wfo(i2)+1), wj2, fb, hi2, fsv(oc+1))
4045 call fs_elim_blk(i2, kk, .false., fa)
4046 call fs_elim_blk(j2, kk, .true., fb)
4047 call hecmw_mf_kernel_gemm(hi2, wj2, pk, hi2, fa, wj2, fb, hi2, fsv(oc+1))
4048 call fs_elim_blk(i2, kk, .true., fa)
4049 call fs_elim_blk(j2, kk, .false., fb)
4050 call hecmw_mf_kernel_gemm(hi2, wj2, pk, hi2, fa, wj2, fb, hi2, fsvu(oc+1))
4055 oc = boff(t2) + int(hi2, 8)*g%tb(j2-1)
4056 if (.not. fct%lu)
then
4057 ra2 = xrank_of(t2, kk, .false.)
4059 fa => wscb(int(tprow(t2), 8)*pk + 1 : int(tprow(t2), 8)*pk + int(hi2, 8)*pk)
4063 ox = bslot(mf_bidx(g%nt, kk, i2))
4065 ua => bval(ox+1 : ox+int(hi2, 8)*max(ra2, 1))
4066 va => bval(ox + int(hi2, 8)*pk + int(pk, 8)*min(hi2, pk) + 1 : &
4067 ox + int(hi2, 8)*pk + int(pk, 8)*min(hi2, pk) + int(pk, 8)*max(ra2, 1))
4069 call bside(kk, j2, .false., rb2, fb, ub, vb)
4070 call mf_update_ab(hi2, wj2, pk, ra2, fa, ua, va, rb2, fb, ub, vb, band(oc+1:oc+int(hi2, 8)*wj2))
4072 ra2 = xrank_of(t2, kk, .false.)
4073 rua2 = xrank_of(t2, kk, .true.)
4074 call aside_lu(kk, t2, hi2, .false., ra2, fa, ua, va)
4075 call bside(kk, j2, .true., rb2, fb, ub, vb)
4076 call mf_update_ab(hi2, wj2, pk, ra2, fa, ua, va, rb2, fb, ub, vb, band(oc+1:oc+int(hi2, 8)*wj2))
4077 call aside_lu(kk, t2, hi2, .true., rua2, fa, ua, va)
4078 call bside(kk, j2, .false., rlb2, fb, ub, vb)
4079 call mf_update_ab(hi2, wj2, pk, rua2, fa, ua, va, rlb2, fb, ub, vb, bandu(oc+1:oc+int(hi2, 8)*wj2))
4081 end subroutine update_pair
4084 subroutine aside_lu(kk, t2, hi2, uface, ra2, fa, ua, va)
4085 integer(kind=kint),
intent(in) :: kk, t2, hi2
4086 logical,
intent(in) :: uface
4087 integer(kind=kint),
intent(in) :: ra2
4088 real(kind=kreal),
pointer,
contiguous,
intent(out) :: fa(:), ua(:), va(:)
4089 integer(kind=8) :: ox
4093 fa => bandu(boff(t2) + int(p0, 8)*hi2 + 1 : boff(t2) + int(p0, 8)*hi2 + int(hi2, 8)*pk)
4095 fa => band(boff(t2) + int(p0, 8)*hi2 + 1 : boff(t2) + int(p0, 8)*hi2 + int(hi2, 8)*pk)
4102 ox = bslotu(mf_bidx(g%nt, kk, g%ntc + t2))
4103 ua => bvalu(ox+1 : ox+int(hi2, 8)*max(ra2, 1))
4104 va => bvalu(ox + int(hi2, 8)*pk + 1 : ox + int(hi2, 8)*pk + int(pk, 8)*max(ra2, 1))
4106 ox = bslot(mf_bidx(g%nt, kk, g%ntc + t2))
4107 ua => bval(ox+1 : ox+int(hi2, 8)*max(ra2, 1))
4108 va => bval(ox + int(hi2, 8)*pk + 1 : ox + int(hi2, 8)*pk + int(pk, 8)*max(ra2, 1))
4111 end subroutine aside_lu
4116 subroutine bside(kk, j2, uface, rb2, fb, ub, vb)
4117 integer(kind=kint),
intent(in) :: kk, j2
4118 logical,
intent(in) :: uface
4119 integer(kind=kint),
intent(out) :: rb2
4120 real(kind=kreal),
pointer,
contiguous,
intent(out) :: fb(:), ub(:), vb(:)
4121 integer(kind=kint) :: tj2, hj2, iwx
4122 integer(kind=8) :: ox
4124 hj2 = g%tb(j2) - g%tb(j2-1)
4129 if (j2 <= g%ntc)
then
4130 call fs_elim_blk(j2, kk, uface, fb)
4134 if (towner(tj2) == me)
then
4135 rb2 = xrank_of(tj2, kk, uface)
4138 fb => bandu(boff(tj2) + int(p0, 8)*hj2 + 1 : boff(tj2) + int(p0, 8)*hj2 + int(hj2, 8)*pk)
4140 fb => band(boff(tj2) + int(p0, 8)*hj2 + 1 : boff(tj2) + int(p0, 8)*hj2 + int(hj2, 8)*pk)
4144 ox = bslotu(mf_bidx(g%nt, kk, j2))
4145 ub => bvalu(ox+1 : ox+int(hj2, 8)*max(rb2, 1))
4146 vb => bvalu(ox + int(hj2, 8)*pk + 1 : ox + int(hj2, 8)*pk + int(pk, 8)*max(rb2, 1))
4148 ox = bslot(mf_bidx(g%nt, kk, j2))
4149 ub => bval(ox+1 : ox+int(hj2, 8)*max(rb2, 1))
4150 vb => bval(ox + int(hj2, 8)*pk + 1 : ox + int(hj2, 8)*pk + int(pk, 8)*max(rb2, 1))
4163 fb => xbuf(iwx)%v(ox+1 : ox+int(hj2, 8)*pk)
4165 ub => xbuf(iwx)%v(ox+1 : ox+int(hj2, 8)*max(rb2, 1))
4166 vb => xbuf(iwx)%v(ox + int(hj2, 8)*max(rb2, 0) + 1 : ox + int(hj2, 8)*max(rb2, 0) + int(pk, 8)*max(rb2, 1))
4169 end subroutine bside
4173 subroutine inertia_dist()
4174 integer(kind=kint) :: kk2, x2
4175 real(kind=kreal) :: d11, d21, d22, det
4178 if (fso(kk2, kk2) /= me) cycle
4179 x2 = g%tb(kk2-1) + 1
4180 do while (x2 <= min(g%tb(kk2), npiv))
4181 if (fct%sn(s)%ptype(x2) == 1)
then
4182 if (fsv(fsx(x2, x2)) > 0.0d0)
then
4189 d11 = fsv(fsx(x2, x2))
4190 d21 = fct%sn(s)%dsub(x2)
4191 d22 = fsv(fsx(x2+1, x2+1))
4192 det = d11*d22 - d21*d21
4193 if (det < 0.0d0)
then
4196 else if (d11 + d22 > 0.0d0)
then
4205 end subroutine inertia_dist
4210 function rfp_count(r0)
result(cnt)
4211 integer(kind=kint),
intent(in) :: r0
4212 integer(kind=kint) :: cnt, c0, r1
4217 if (fso(g%dtile(r1-1), g%dtile(c0-1)) == r0)
then
4218 cnt = cnt + merge(2, 1, fct%lu)
4219 if (fct%lu .and. r1 == c0) cnt = cnt - 1
4223 end function rfp_count
4227 subroutine rfp_move(r0, buf, topanel)
4228 integer(kind=kint),
intent(in) :: r0
4229 real(kind=kreal),
intent(inout) :: buf(:)
4230 logical,
intent(in) :: topanel
4231 integer(kind=kint) :: c0, r1, jj2, m0, cnt
4232 integer(kind=8) :: px
4238 px = int(jj2-1, 8)*m0
4240 if (fso(g%dtile(r1-1), g%dtile(c0-1)) /= r0) cycle
4243 pv(px + (r1 - pa + 1)) = buf(cnt)
4245 buf(cnt) = pv(px + (r1 - pa + 1))
4249 do r1 = c0 + 1, ncol
4250 if (fso(g%dtile(r1-1), g%dtile(c0-1)) /= r0) cycle
4253 pvu(px + (r1 - pa + 1)) = buf(cnt)
4255 buf(cnt) = pvu(px + (r1 - pa + 1))
4260 end subroutine rfp_move
4264 subroutine rfp_mine(buf, topack)
4265 real(kind=kreal),
intent(inout) :: buf(:)
4266 logical,
intent(in) :: topack
4267 integer(kind=kint) :: c0, r1, cnt
4272 if (fso(g%dtile(r1-1), g%dtile(c0-1)) /= me) cycle
4275 buf(cnt) = fsv(fsx(r1, c0))
4277 fsv(fsx(r1, c0)) = buf(cnt)
4281 do r1 = c0 + 1, ncol
4282 if (fso(g%dtile(r1-1), g%dtile(c0-1)) /= me) cycle
4285 buf(cnt) = fsvu(fsx(r1, c0))
4287 fsvu(fsx(r1, c0)) = buf(cnt)
4292 end subroutine rfp_mine
4299 subroutine root_final()
4300 real(kind=kreal),
allocatable :: rw(:)
4301 integer(kind=kint) :: ic2, q2, mm, ib, r0, cnt, m0
4309 call mf_grow_r(pv, int(m0, 8)*w)
4310 if (fct%lu)
call mf_grow_r(pvu, int(m0, 8)*w)
4311 call mf_grow_r(wk, int(2*m0, 8))
4314 call mf_grow_r(swb, int(cnt, 8))
4315 call rfp_mine(swb, .true.)
4316 call rfp_move(me, swb, .true.)
4324 call hecmw_recv_r(rw, cnt, r0, 16*s+5, map%comm, stat)
4325 call rfp_move(r0, rw, .true.)
4328 if (.not. fct%lu)
then
4329 call hecmw_mf_kernel_panel_piv(m0, w, m0, pv, blk(pa:pb), fct%pivot_u, mf_pivot_alpha, &
4330 zero, .false., wk, np, perm, fct%sn(s)%ptype(pa:pb), fct%sn(s)%dsub(pa:pb), nsw, n22, info)
4331 n2x2_f = n2x2_f + n22
4333 call hecmw_mf_kernel_panel_lu_piv(m0, w, m0, pv, m0, pvu, blk(pa:pb), blkr(pa:pb), &
4334 fct%pivot_u, zero, .false., np, perm, permr, nsw, info)
4336 nswap_f = nswap_f + nsw
4339 if (info /= 0) ctli(2) = fct%sn(s)%fsdof(pa + info - 1)
4340 call ctl_send_i(ctli, 8)
4342 mm = merge(2, 1, fct%lu)*w + w
4343 ib = mf_pool_slot(pool)
4344 allocate(pool%box(ib)%hdr(mm))
4345 pool%box(ib)%hdr(1:w) = perm(1:w)
4348 pool%box(ib)%hdr(q2+1:q2+w) = permr(1:w)
4351 pool%box(ib)%hdr(q2+1:q2+w) = fct%sn(s)%ptype(pa:pb)
4353 if (plist(ic2) == me) cycle
4355 call hecmw_isend_int(pool%box(ib)%hdr, mm, plist(ic2), 16*s+3, map%comm, req)
4356 call mf_pool_req(pool, req)
4358 if (.not. fct%lu)
then
4359 ib = mf_pool_slot(pool)
4360 allocate(pool%box(ib)%rv(w))
4361 pool%box(ib)%rv(1:w) = fct%sn(s)%dsub(pa:pb)
4363 if (plist(ic2) == me) cycle
4365 call hecmw_isend_r(pool%box(ib)%rv, w, plist(ic2), 16*s+4, map%comm, req)
4366 call mf_pool_req(pool, req)
4373 ib = mf_pool_slot(pool)
4374 allocate(pool%box(ib)%rv(cnt))
4375 call rfp_mine(pool%box(ib)%rv, .true.)
4377 call hecmw_isend_r(pool%box(ib)%rv, cnt, master, 16*s+5, map%comm, req)
4378 call mf_pool_req(pool, req)
4380 call hecmw_recv_int(ctli, 8, master, 16*s+3, map%comm, stat)
4381 if (ctli(2) /= 0) info = 1
4383 mm = merge(2, 1, fct%lu)*w + w
4384 call hecmw_recv_int(itmp, mm, master, 16*s+3, map%comm, stat)
4385 perm(1:w) = itmp(1:w)
4388 permr(1:w) = itmp(q2+1:q2+w)
4391 fct%sn(s)%ptype(pa:pb) = itmp(q2+1:q2+w)
4392 if (.not. fct%lu)
call hecmw_recv_r(fct%sn(s)%dsub(pa:pb), w, master, 16*s+4, map%comm, stat)
4396 if (gierr == 0) gierr = ctli(2)
4401 call perm_dance_apply()
4408 ib = mf_pool_slot(pool)
4409 allocate(pool%box(ib)%rv(cnt))
4410 call rfp_move(r0, pool%box(ib)%rv, .false.)
4412 call hecmw_isend_r(pool%box(ib)%rv, cnt, r0, 16*s+4, map%comm, req)
4413 call mf_pool_req(pool, req)
4417 call mf_grow_r(swb, int(cnt, 8))
4418 call rfp_move(me, swb, .false.)
4419 call rfp_mine(swb, .false.)
4425 call hecmw_recv_r(rw, cnt, master, 16*s+4, map%comm, stat)
4426 call rfp_mine(rw, .false.)
4431 end subroutine root_final
4438 subroutine store_dist()
4439 integer(kind=kint) :: kk, t2, h2, mw, idx2, r2, i2, j2, cc2, rr2, colx, cnbp, kbegp, cndelp
4440 integer(kind=kint) :: ic2, r0, cnt
4441 integer(kind=8) :: off, offu, pw0, bx, sb, halfp
4442 integer(kind=kint),
allocatable :: ctbp(:)
4443 real(kind=kreal),
allocatable :: dsq(:), rw(:)
4447 if (min(g%tb(kk), npiv) - g%tb(kk-1) <= 0)
exit
4451 idx2 = mf_bidx(g%nt, nkc, g%nt)
4452 call mf_grow_i(fct%sn(s)%brank, idx2)
4453 call mf_grow_i8(fct%sn(s)%bptr, idx2 + 1)
4455 call mf_grow_i(fct%sn(s)%branku, idx2)
4456 call mf_grow_i8(fct%sn(s)%bptru, idx2 + 1)
4463 mw = min(g%tb(kk), npiv) - g%tb(kk-1)
4465 if (.not. tile_is_mine(kk, i2)) cycle
4466 h2 = g%tb(i2) - g%tb(i2-1)
4467 idx2 = mf_bidx(g%nt, kk, i2)
4468 pw0 = pw0 + int(h2, 8)*mw
4470 if (fct%blr .and. i2 > g%ntc) r2 = brk(idx2)
4471 fct%sn(s)%brank(idx2) = r2
4472 fct%sn(s)%bptr(idx2) = off
4474 off = off + int(h2, 8)*mw
4476 off = off + int(r2, 8)*(h2 + mw)
4480 if (fct%blr .and. i2 > g%ntc) r2 = brku(idx2)
4481 fct%sn(s)%branku(idx2) = r2
4482 fct%sn(s)%bptru(idx2) = offu
4484 offu = offu + int(h2, 8)*mw
4486 offu = offu + int(r2, 8)*(h2 + mw)
4491 call mf_grow_r(fct%sn(s)%lval, max(off, 1_8))
4492 if (fct%lu)
call mf_grow_r(fct%sn(s)%uval, max(offu, 1_8))
4494 mw = min(g%tb(kk), npiv) - g%tb(kk-1)
4496 if (.not. tile_is_mine(kk, i2)) cycle
4497 call store_tile(kk, i2, mw)
4500 fct%sn(s)%pwords = off + offu
4503 fct%factor_words_act = fct%factor_words_act + off
4504 if (fct%lu .and. fct%blr) fct%factor_words_act = fct%factor_words_act + offu
4506 fct%blr_words_fr = fct%blr_words_fr + merge(2, 1, fct%lu)*pw0
4511 cndelp = ncol - npiv
4512 if (cndelp > 0)
then
4513 allocate(dsq(merge(2, 1, fct%lu)*int(cndelp, 8)*cndelp))
4515 call dsq_mine(dsq, cndelp)
4519 cnt = dsq_count(r0, cndelp)
4522 call hecmw_recv_r(rw, cnt, r0, 16*s+15, map%comm, stat)
4523 call dsq_unpack(r0, cndelp, rw, dsq)
4527 cnt = dsq_count(me, cndelp)
4529 i2 = mf_pool_slot(pool)
4530 allocate(pool%box(i2)%rv(cnt))
4531 call dsq_pack(cndelp, pool%box(i2)%rv)
4533 call hecmw_isend_r(pool%box(i2)%rv, cnt, master, 16*s+15, map%comm, req)
4534 call mf_pool_req(pool, req)
4540 allocate(ctbp(0:g%nt - g%ntc + 2))
4543 if (cndelp > 0)
then
4548 do j2 = g%ntc + 1, g%nt
4550 ctbp(cnbp) = cndelp + g%tb(j2) - ncol
4555 if (.not. cbblk_is_mine(i2, kbegp)) cycle
4556 halfp = halfp + int(ctbp(i2) - ctbp(i2-1), 8)*(ctbp(j2) - ctbp(j2-1))
4559 fct%sn(s)%cbsize = merge(2, 1, fct%lu)*halfp
4560 if (
allocated(fct%sn(s)%cval))
deallocate(fct%sn(s)%cval)
4561 if (fct%sn(s)%cbsize > 0)
then
4562 allocate(fct%sn(s)%cval(fct%sn(s)%cbsize))
4566 if (.not. cbblk_is_mine(i2, kbegp)) cycle
4567 h2 = ctbp(i2) - ctbp(i2-1)
4568 do cc2 = ctbp(j2-1) + 1, ctbp(j2)
4569 colx = cb_col(cc2, cndelp)
4570 if (i2 <= kbegp)
then
4572 do rr2 = max(cc2, ctbp(i2-1) + 1), ctbp(i2)
4573 fct%sn(s)%cval(sb + int(cc2 - ctbp(j2-1) - 1, 8)*h2 + (rr2 - ctbp(i2-1))) = &
4574 dsq(int(cc2-1, 8)*cndelp + rr2)
4575 if (fct%lu) fct%sn(s)%cval(halfp + sb + int(cc2 - ctbp(j2-1) - 1, 8)*h2 + (rr2 - ctbp(i2-1))) = &
4576 dsq(int(cndelp, 8)*cndelp + int(cc2-1, 8)*cndelp + rr2)
4580 bx = boff(t2) + int(colx - 1, 8)*h2
4581 fct%sn(s)%cval(sb + int(cc2 - ctbp(j2-1) - 1, 8)*h2 + 1 : &
4582 sb + int(cc2 - ctbp(j2-1) - 1, 8)*h2 + h2) = band(bx+1:bx+h2)
4583 if (fct%lu) fct%sn(s)%cval(halfp + sb + int(cc2 - ctbp(j2-1) - 1, 8)*h2 + 1 : &
4584 halfp + sb + int(cc2 - ctbp(j2-1) - 1, 8)*h2 + h2) = bandu(bx+1:bx+h2)
4587 sb = sb + int(ctbp(j2) - ctbp(j2-1), 8)*h2
4590 fct%live_cb = fct%live_cb + fct%sn(s)%cbsize
4591 fct%stack_peak_act = max(fct%stack_peak_act, fct%live_cb)
4594 if (
allocated(dsq))
deallocate(dsq)
4595 end subroutine store_dist
4598 function dsq_count(r0, cndelP)
result(cnt)
4599 integer(kind=kint),
intent(in) :: r0, cndelp
4600 integer(kind=kint) :: cnt, cc2, rr2
4604 do rr2 = cc2, cndelp
4605 if (fso(g%dtile(npiv+rr2-1), g%dtile(npiv+cc2-1)) == r0) cnt = cnt + merge(2, 1, fct%lu)
4608 end function dsq_count
4611 subroutine dsq_pack(cndelP, buf)
4612 integer(kind=kint),
intent(in) :: cndelp
4613 real(kind=kreal),
intent(inout) :: buf(:)
4614 integer(kind=kint) :: cc2, rr2, cnt
4618 do rr2 = cc2, cndelp
4619 if (fso(g%dtile(npiv+rr2-1), g%dtile(npiv+cc2-1)) /= me) cycle
4621 buf(cnt) = fsv(fsx(npiv+rr2, npiv+cc2))
4624 buf(cnt) = fsvu(fsx(npiv+rr2, npiv+cc2))
4628 end subroutine dsq_pack
4632 subroutine dsq_unpack(r0, cndelP, buf, dsq)
4633 integer(kind=kint),
intent(in) :: r0, cndelp
4634 real(kind=kreal),
intent(in) :: buf(:)
4635 real(kind=kreal),
intent(inout) :: dsq(:)
4636 integer(kind=kint) :: cc2, rr2, cnt
4640 do rr2 = cc2, cndelp
4641 if (fso(g%dtile(npiv+rr2-1), g%dtile(npiv+cc2-1)) /= r0) cycle
4643 dsq(int(cc2-1, 8)*cndelp + rr2) = buf(cnt)
4646 dsq(int(cndelp, 8)*cndelp + int(cc2-1, 8)*cndelp + rr2) = buf(cnt)
4650 end subroutine dsq_unpack
4653 subroutine dsq_mine(dsq, cndelP)
4654 real(kind=kreal),
intent(inout) :: dsq(:)
4655 integer(kind=kint),
intent(in) :: cndelp
4656 integer(kind=kint) :: cc2, rr2
4659 do rr2 = cc2, cndelp
4660 if (fso(g%dtile(npiv+rr2-1), g%dtile(npiv+cc2-1)) /= me) cycle
4661 dsq(int(cc2-1, 8)*cndelp + rr2) = fsv(fsx(npiv+rr2, npiv+cc2))
4662 if (fct%lu) dsq(int(cndelp, 8)*cndelp + int(cc2-1, 8)*cndelp + rr2) = fsvu(fsx(npiv+rr2, npiv+cc2))
4665 end subroutine dsq_mine
4668 function tile_is_mine(kk, i2)
result(mine0)
4669 integer(kind=kint),
intent(in) :: kk, i2
4672 if (i2 <= g%ntc)
then
4673 mine0 = (fso(i2, kk) == me)
4675 mine0 = (towner(i2 - g%ntc) == me)
4677 end function tile_is_mine
4680 function cbblk_is_mine(i2, kbegP)
result(mine0)
4681 integer(kind=kint),
intent(in) :: i2, kbegp
4684 if (i2 <= kbegp)
then
4687 mine0 = (towner(i2 - kbegp) == me)
4689 end function cbblk_is_mine
4692 function cb_col(cc2, cndelP)
result(colx)
4693 integer(kind=kint),
intent(in) :: cc2, cndelp
4694 integer(kind=kint) :: colx
4696 if (cc2 <= cndelp)
then
4699 colx = ncol + (cc2 - cndelp)
4705 subroutine store_tile(kk, i2, mw)
4706 integer(kind=kint),
intent(in) :: kk, i2, mw
4707 integer(kind=kint) :: h2, r2, t2, idx2
4708 integer(kind=8) :: bx, dx
4710 h2 = g%tb(i2) - g%tb(i2-1)
4711 idx2 = mf_bidx(g%nt, kk, i2)
4712 dx = fct%sn(s)%bptr(idx2)
4713 if (i2 <= g%ntc)
then
4715 fct%sn(s)%lval(dx+1:dx+int(h2, 8)*mw) = fsv(bx+1:bx+int(h2, 8)*mw)
4717 dx = fct%sn(s)%bptru(idx2)
4718 fct%sn(s)%uval(dx+1:dx+int(h2, 8)*mw) = fsvu(bx+1:bx+int(h2, 8)*mw)
4723 r2 = fct%sn(s)%brank(idx2)
4725 bx = boff(t2) + int(g%tb(kk-1), 8)*h2
4726 fct%sn(s)%lval(dx+1:dx+int(h2, 8)*mw) = band(bx+1:bx+int(h2, 8)*mw)
4729 fct%sn(s)%lval(dx+1:dx+int(h2, 8)*r2) = bval(bx+1:bx+int(h2, 8)*r2)
4730 fct%sn(s)%lval(dx+int(h2, 8)*r2+1:dx+int(h2, 8)*r2+int(mw, 8)*r2) = &
4731 bval(bx+int(h2, 8)*mw+1:bx+int(h2, 8)*mw+int(mw, 8)*r2)
4734 dx = fct%sn(s)%bptru(idx2)
4735 r2 = fct%sn(s)%branku(idx2)
4737 bx = boff(t2) + int(g%tb(kk-1), 8)*h2
4738 fct%sn(s)%uval(dx+1:dx+int(h2, 8)*mw) = bandu(bx+1:bx+int(h2, 8)*mw)
4741 fct%sn(s)%uval(dx+1:dx+int(h2, 8)*r2) = bvalu(bx+1:bx+int(h2, 8)*r2)
4742 fct%sn(s)%uval(dx+int(h2, 8)*r2+1:dx+int(h2, 8)*r2+int(mw, 8)*r2) = &
4743 bvalu(bx+int(h2, 8)*mw+1:bx+int(h2, 8)*mw+int(mw, 8)*r2)
4746 end subroutine store_tile
4748 end subroutine mf_super_factor_dist
4752 subroutine mf_super_task_mpi(hecMAT, sym, map, fct, wrks, left, s, zero, gierr)
4754 type(hecmwst_matrix),
intent(in) :: hecmat
4755 type(hecmwst_mf_symbolic),
intent(in) :: sym
4756 type(hecmwst_mf_map),
intent(in) :: map
4757 type(hecmwst_mf_factor),
intent(inout) :: fct
4758 type(mf_work),
intent(inout) :: wrks(0:)
4759 integer(kind=kint),
intent(inout) :: left(:)
4760 integer(kind=kint),
intent(in) :: s
4761 real(kind=kreal),
intent(in) :: zero
4762 integer(kind=kint),
intent(inout) :: gierr
4763 integer(kind=kint) :: ierr, cur, tid, p, n, ss, c, l, nr, nteam
4777 nr = nr + fct%sn(c)%ncol - fct%sn(c)%npiv
4780 do l = sym%rptr(ss), sym%rptr(ss+1) - 1
4781 nr = nr + sym%ndof(sym%rlist(l))
4783 if (nr >= mf_par_rows .and. nteam > 1)
then
4787 type(mf_work) :: lwrk
4788 call mf_super_factor(hecmat, sym, fct, fct%sn(ss), lwrk, ss, zero, ierr)
4791 call mf_super_factor(hecmat, sym, fct, fct%sn(ss), wrks(tid), ss, zero, ierr)
4800 if (map%upper(p))
return
4804 left(p) = left(p) - 1
4810 end subroutine mf_super_task_mpi
4815 subroutine mf_super_factor(hecMAT, sym, fct, sn, wrk, s, zero, ierr)
4817 type(hecmwst_matrix),
intent(in) :: hecmat
4818 type(hecmwst_mf_symbolic),
intent(in) :: sym
4819 type(hecmwst_mf_factor),
intent(inout) :: fct
4820 type(mf_snode),
intent(inout) :: sn
4821 type(mf_work),
intent(inout) :: wrk
4822 integer(kind=kint),
intent(in) :: s
4823 real(kind=kreal),
intent(in) :: zero
4824 integer(kind=kint),
intent(out) :: ierr
4825 integer(kind=kint) :: nown, nrow_nodes, ncol0, ndel, c, i, i0, j, k, l, nd, a, b, m, r, cc, pr, pc
4826 integer(kind=kint) :: j0, ki, kk, coff, roff, cnown, cnb, cndel, cnrow, hi, cnt, cntc, kbeg, nkc
4827 integer(kind=kint) :: nswap, n2x2, npos, nneg, ntl, nlr, nsk, rmax
4828 integer(kind=8) :: o, o2, base, base2, half, pw, fw, rsum, pwa, pwu
4829 integer(kind=8),
allocatable :: pcb(:)
4830 real(kind=kreal) :: v
4831 logical :: par, ruse
4835 nown = sym%sptr(s+1) - sym%sptr(s)
4836 nrow_nodes = sym%rptr(s+1) - sym%rptr(s)
4837 ncol0 = fct%cdofptr(sym%sptr(s+1)) - fct%cdofptr(sym%sptr(s))
4841 ndel = ndel + fct%sn(c)%ncol - fct%sn(c)%npiv
4845 fct%max_growth = max(fct%max_growth, ndel)
4847 call mf_partition(sym, fct, s, ndel, wrk%g)
4848 call mf_grow_i(wrk%rowoff, nrow_nodes + 1)
4849 call mf_grow_i(wrk%pos, sym%nnode)
4850 if (
allocated(wrk%ccoloff))
then
4851 if (
size(wrk%ccoloff) < wrk%g%nt + 2)
deallocate(wrk%ccoloff, wrk%ctb)
4853 if (.not.
allocated(wrk%ccoloff))
allocate(wrk%ccoloff(wrk%g%nt+2), wrk%ctb(0:wrk%g%nt+1))
4858 if (fct%blr .and. fct%blr_reuse .and.
allocated(sn%brank) .and.
allocated(sn%tbnd))
then
4859 if (sn%nt == wrk%g%nt .and. sn%ntc == wrk%g%ntc .and. sn%npiv == sn%ncol .and. &
4860 size(sn%tbnd) >= wrk%g%nt + 1)
then
4861 ruse = all(sn%tbnd(1:wrk%g%nt+1) == wrk%g%tb(0:wrk%g%nt))
4864 sn%ncol = wrk%g%ncol
4867 call mf_grow_i(sn%tbnd, wrk%g%nt + 1)
4868 sn%tbnd(1:wrk%g%nt+1) = wrk%g%tb(0:wrk%g%nt)
4869 call mf_grow_i(sn%fsdof, wrk%g%ncol)
4870 call mf_grow_i(sn%frow, wrk%g%ncol)
4871 call mf_grow_i(sn%ptype, wrk%g%ncol)
4872 call mf_grow_r(sn%dsub, int(wrk%g%ncol, 8))
4873 call mf_grow_i(wrk%blk, wrk%g%ncol)
4874 call mf_grow_i(wrk%blkr, wrk%g%ncol)
4876 do k = sym%sptr(s), sym%sptr(s+1) - 1
4877 do i = fct%cdofptr(k), fct%cdofptr(k+1) - 1
4886 do j = fct%sn(c)%npiv + 1, fct%sn(c)%ncol
4888 sn%fsdof(m) = fct%sn(c)%fsdof(j)
4889 sn%frow(m) = fct%sn(c)%frow(j)
4894 wrk%blkr(1:wrk%g%ncol) = wrk%blk(1:wrk%g%ncol)
4899 do i = 1, nrow_nodes
4900 k = sym%rlist(sym%rptr(s)+i-1)
4901 wrk%rowoff(i+1) = wrk%rowoff(i) + sym%ndof(k)
4902 if (i == nown) wrk%rowoff(i+1) = wrk%rowoff(i+1) + ndel
4907 call mf_grow_r(wrk%fval, wrk%g%coloff(wrk%g%nt+1))
4908 fw = wrk%g%coloff(wrk%g%nt+1)
4909 if (fct%lu) fw = 2*fw
4911 fct%front_words_act = max(fct%front_words_act, fw)
4912 fct%live_front = fct%live_front + fw
4913 fct%front_peak = max(fct%front_peak, fct%live_front)
4915 if (fct%lu)
call mf_grow_r(wrk%fvalu, wrk%g%coloff(wrk%g%nt+1))
4919 wrk%fval(wrk%g%coloff(j)+1:wrk%g%coloff(j+1)) = 0.0d0
4920 if (fct%lu) wrk%fvalu(wrk%g%coloff(j)+1:wrk%g%coloff(j+1)) = 0.0d0
4931 do i = i0, min(i0+7, nown)
4932 k = sym%sptr(s) + i - 1
4934 coff = wrk%rowoff(i)
4935 base = int(j0-1, 8)*nd*nd
4938 o = mf_idx(wrk%g, coff+a, coff+b)
4939 wrk%fval(o) = wrk%fval(o) + hecmat%D(base + (a-1)*nd + b)
4943 o = mf_idx(wrk%g, coff+b, coff+a)
4944 wrk%fvalu(o) = wrk%fvalu(o) + hecmat%D(base + (a-1)*nd + b)
4948 do kk = hecmat%indexL(j0-1)+1, hecmat%indexL(j0)
4949 ki = sym%invp(hecmat%itemL(kk))
4951 roff = wrk%rowoff(wrk%pos(ki))
4952 base = int(kk-1, 8)*nd*nd
4953 if (.not. fct%lu)
then
4956 o = mf_idx(wrk%g, roff+a, coff+b)
4957 wrk%fval(o) = wrk%fval(o) + hecmat%AL(base + (b-1)*nd + a)
4961 base2 = int(fct%mirror(kk)-1, 8)*nd*nd
4964 o = mf_idx(wrk%g, roff+b, coff+a)
4965 wrk%fvalu(o) = wrk%fvalu(o) + hecmat%AL(base + (a-1)*nd + b)
4966 wrk%fval(o) = wrk%fval(o) + hecmat%AU(base2 + (b-1)*nd + a)
4971 do kk = hecmat%indexU(j0-1)+1, hecmat%indexU(j0)
4972 ki = sym%invp(hecmat%itemU(kk))
4974 roff = wrk%rowoff(wrk%pos(ki))
4975 base = int(kk-1, 8)*nd*nd
4976 if (.not. fct%lu)
then
4979 o = mf_idx(wrk%g, roff+a, coff+b)
4980 wrk%fval(o) = wrk%fval(o) + hecmat%AU(base + (b-1)*nd + a)
4984 base2 = int(fct%mirroru(kk)-1, 8)*nd*nd
4987 o = mf_idx(wrk%g, roff+b, coff+a)
4988 wrk%fvalu(o) = wrk%fvalu(o) + hecmat%AU(base + (a-1)*nd + b)
4989 wrk%fval(o) = wrk%fval(o) + hecmat%AL(base2 + (b-1)*nd + a)
5004 cnown = sym%sptr(c+1) - sym%sptr(c)
5005 cndel = fct%sn(c)%ncol - fct%sn(c)%npiv
5007 cntc = fct%sn(c)%ntc
5008 if (
size(wrk%ccoloff) < cnt + 2)
then
5009 deallocate(wrk%ccoloff, wrk%ctb)
5010 allocate(wrk%ccoloff(cnt+2), wrk%ctb(0:cnt+1))
5018 do j = cntc + 1, cnt
5020 wrk%ctb(cnb) = cndel + fct%sn(c)%tbnd(j+1) - fct%sn(c)%ncol
5022 cnrow = wrk%ctb(cnb)
5025 wrk%ccoloff(j+1) = wrk%ccoloff(j) + int(cnrow - wrk%ctb(j-1), 8)*(wrk%ctb(j) - wrk%ctb(j-1))
5027 half = wrk%ccoloff(cnb+1)
5028 call mf_grow_i(wrk%cmapdof, cnrow)
5030 wrk%cmapdof(r) = m + r
5034 do l = sym%rptr(c) + cnown, sym%rptr(c+1) - 1
5035 k = sym%cmap(sym%cmap_ptr(c) + (l - sym%rptr(c) - cnown))
5036 do a = 1, sym%ndof(sym%rlist(l))
5038 wrk%cmapdof(r) = wrk%rowoff(k) + a
5047 base = wrk%ccoloff(j) + int(wrk%ctb(i-1) - wrk%ctb(j-1), 8)*(wrk%ctb(j) - wrk%ctb(j-1))
5048 hi = wrk%ctb(i) - wrk%ctb(i-1)
5049 do cc = wrk%ctb(j-1) + 1, wrk%ctb(j)
5050 do r = max(cc, wrk%ctb(i-1) + 1), wrk%ctb(i)
5051 o = base + int(cc - 1 - wrk%ctb(j-1), 8)*hi + (r - 1 - wrk%ctb(i-1)) + 1
5052 v = fct%sn(c)%cval(o)
5053 if (.not. fct%lu)
then
5054 o = mf_idx(wrk%g, max(wrk%cmapdof(r), wrk%cmapdof(cc)), min(wrk%cmapdof(r), wrk%cmapdof(cc)))
5055 wrk%fval(o) = wrk%fval(o) + v
5058 pc = wrk%cmapdof(cc)
5060 o2 = mf_idx(wrk%g, pr, pc)
5061 wrk%fval(o2) = wrk%fval(o2) + v
5063 o2 = mf_idx(wrk%g, pc, pr)
5064 wrk%fvalu(o2) = wrk%fvalu(o2) + v
5068 o2 = mf_idx(wrk%g, pc, pr)
5069 wrk%fval(o2) = wrk%fval(o2) + fct%sn(c)%cval(o + half)
5071 o2 = mf_idx(wrk%g, pr, pc)
5072 wrk%fvalu(o2) = wrk%fvalu(o2) + fct%sn(c)%cval(o + half)
5083 fct%live_cb = fct%live_cb - fct%sn(c)%cbsize
5085 deallocate(fct%sn(c)%cval)
5099 call mf_factor_front_lu(sn, wrk%g, wrk%fval, wrk%fvalu, wrk%pval, wrk%pvalu, wrk%wk, wrk%blk, wrk%blkr, &
5100 sym%sparent(s) == 0, fct%pivot_u, zero, fct%blr, fct%eps, wrk%bval, wrk%boff, wrk%brk, &
5101 wrk%bvalu, wrk%boffu, wrk%brku, fct, s, ruse, nswap, ntl, nlr, nsk, rsum, rmax, ierr)
5103 call mf_factor_front(sn, wrk%g, wrk%fval, wrk%pval, wrk%wval, wrk%wk, wrk%blk, &
5104 sym%sparent(s) == 0, fct%pivot_u, zero, fct%blr, fct%eps, wrk%bval, wrk%boff, wrk%brk, &
5105 fct, s, ruse, nswap, n2x2, npos, nneg, ntl, nlr, nsk, rsum, rmax, ierr)
5106 sn%frow(1:wrk%g%ncol) = sn%fsdof(1:wrk%g%ncol)
5110 fct%live_front = fct%live_front - fw
5115 fct%n_swap = fct%n_swap + nswap
5116 fct%n_2x2 = fct%n_2x2 + n2x2
5117 fct%n_pos = fct%n_pos + npos
5118 fct%n_neg = fct%n_neg + nneg
5119 fct%n_delay = fct%n_delay + wrk%g%ncol - sn%npiv
5120 fct%blr_tiles = fct%blr_tiles + ntl
5121 fct%blr_tiles_lr = fct%blr_tiles_lr + nlr
5122 fct%blr_rank_sum = fct%blr_rank_sum + rsum
5123 fct%blr_rank_max = max(fct%blr_rank_max, rmax)
5124 fct%blr_tiles_skip = fct%blr_tiles_skip + nsk
5130 m = min(wrk%g%tb(k), sn%npiv) - wrk%g%tb(k-1)
5132 pw = pw + int(wrk%g%nrow - wrk%g%tb(k-1), 8)*m
5134 allocate(pcb(0:wrk%g%ntc))
5138 m = min(wrk%g%tb(k), sn%npiv) - wrk%g%tb(k-1)
5141 pcb(k) = pcb(k-1) + int(wrk%g%nrow - wrk%g%tb(k-1), 8)*m
5143 if (.not. fct%blr)
then
5144 call mf_grow_r(sn%lval, pw)
5145 if (fct%lu)
call mf_grow_r(sn%uval, pw)
5149 m = min(wrk%g%tb(k), sn%npiv) - wrk%g%tb(k-1)
5152 hi = wrk%g%tb(i) - wrk%g%tb(i-1)
5153 o = mf_off(wrk%g, i, k)
5154 sn%lval(base:base+int(hi, 8)*m-1) = wrk%fval(o+1:o+int(hi, 8)*m)
5155 if (fct%lu) sn%uval(base:base+int(hi, 8)*m-1) = wrk%fvalu(o+1:o+int(hi, 8)*m)
5156 base = base + int(hi, 8)*m
5162 if (fct%lu) sn%pwords = 2*pw
5164 fct%factor_words_act = fct%factor_words_act + sn%pwords
5167 call mf_store_blr(wrk%g, sn%npiv, nkc, wrk%fval, wrk%bval, wrk%boff, wrk%brk, &
5168 sn%brank, sn%bptr, sn%lval, pwa, par)
5170 call mf_store_blr(wrk%g, sn%npiv, nkc, wrk%fvalu, wrk%bvalu, wrk%boffu, wrk%brku, &
5171 sn%branku, sn%bptru, sn%uval, pwu, par)
5176 fct%factor_words_act = fct%factor_words_act + pwa
5177 fct%blr_words_fr = fct%blr_words_fr + pw
5178 if (fct%lu) fct%blr_words_fr = fct%blr_words_fr + pw
5184 cndel = wrk%g%ncol - sn%npiv
5192 do j = wrk%g%ntc + 1, wrk%g%nt
5194 wrk%ctb(cnb) = cndel + wrk%g%tb(j) - wrk%g%ncol
5196 cnrow = wrk%ctb(cnb)
5199 wrk%ccoloff(j+1) = wrk%ccoloff(j) + int(cnrow - wrk%ctb(j-1), 8)*(wrk%ctb(j) - wrk%ctb(j-1))
5201 half = wrk%ccoloff(cnb+1)
5203 if (fct%lu) sn%cbsize = 2*half
5204 if (
allocated(sn%cval))
deallocate(sn%cval)
5205 if (sn%cbsize > 0)
then
5206 allocate(sn%cval(sn%cbsize))
5207 call mf_store_cb(wrk%g, wrk%fval, sn%npiv, cndel, cnb, kbeg, wrk%ctb, wrk%ccoloff, sn%cval, 0_8, par)
5208 if (fct%lu)
call mf_store_cb(wrk%g, wrk%fvalu, sn%npiv, cndel, cnb, kbeg, wrk%ctb, wrk%ccoloff, sn%cval, half, par)
5210 fct%live_cb = fct%live_cb + sn%cbsize
5211 fct%stack_peak_act = max(fct%stack_peak_act, fct%live_cb)
5215 fct%live_front = fct%live_front - fw
5218 end subroutine mf_super_factor
5223 subroutine mf_store_blr(g, npiv, nkc, gval, bval, boff, brk, brank, bptr, dst, pwa, par)
5225 type(mf_grid),
intent(in) :: g
5226 integer(kind=kint),
intent(in) :: npiv, nkc
5227 real(kind=kreal),
intent(in) :: gval(:), bval(:)
5228 integer(kind=8),
intent(in) :: boff(:)
5229 integer(kind=kint),
intent(in) :: brk(:)
5230 integer(kind=kint),
allocatable,
intent(inout) :: brank(:)
5231 integer(kind=8),
allocatable,
intent(inout) :: bptr(:)
5232 real(kind=kreal),
allocatable,
intent(inout) :: dst(:)
5233 integer(kind=8),
intent(out) :: pwa
5234 logical,
intent(in) :: par
5235 integer(kind=kint) :: ntile, k, i, m, hi, r, idx
5236 integer(kind=8) :: off, base, ob, o
5240 ntile = ntile + g%nt - k + 1
5242 call mf_grow_i(brank, ntile)
5243 call mf_grow_i8(bptr, ntile + 1)
5246 m = min(g%tb(k), npiv) - g%tb(k-1)
5248 idx = mf_bidx(g%nt, k, i)
5249 hi = g%tb(i) - g%tb(i-1)
5251 if (i > g%ntc) r = brk(idx)
5255 off = off + int(hi, 8)*m
5257 off = off + int(r, 8)*(hi + m)
5263 call mf_grow_r(dst, pwa)
5267 m = min(g%tb(k), npiv) - g%tb(k-1)
5269 idx = mf_bidx(g%nt, k, i)
5270 hi = g%tb(i) - g%tb(i-1)
5275 dst(base+1:base+int(hi, 8)*m) = gval(o+1:o+int(hi, 8)*m)
5276 else if (r > 0)
then
5278 dst(base+1:base+int(hi, 8)*r) = bval(ob+1:ob+int(hi, 8)*r)
5279 dst(base+int(hi, 8)*r+1:base+int(hi, 8)*r+int(m, 8)*r) = &
5280 bval(ob+int(hi, 8)*m+1:ob+int(hi, 8)*m+int(m, 8)*r)
5286 end subroutine mf_store_blr
5295 subroutine mf_factor_front(sn, g, fval, pval, wval, wk, blk, isroot, u, zero, blr, eps, &
5296 bval, boff, brk, fct, s, ruse, nswap, n2x2, npos, nneg, ntl, nlr, nsk, rsum, rmax, ierr)
5298 type(mf_snode),
intent(inout) :: sn
5299 type(mf_grid),
intent(in) :: g
5300 real(kind=kreal),
allocatable,
intent(inout) :: fval(:)
5301 real(kind=kreal),
allocatable,
intent(inout) :: pval(:), wval(:), wk(:)
5302 integer(kind=kint),
intent(inout) :: blk(:)
5303 logical,
intent(in) :: isroot, blr
5304 real(kind=kreal),
intent(in) :: u, zero, eps
5305 real(kind=kreal),
allocatable,
intent(inout) :: bval(:)
5306 integer(kind=8),
allocatable,
intent(inout) :: boff(:)
5307 integer(kind=kint),
allocatable,
intent(inout) :: brk(:)
5308 type(hecmwst_mf_factor),
intent(in) :: fct
5309 integer(kind=kint),
intent(in) :: s
5310 logical,
intent(in) :: ruse
5311 integer(kind=kint),
intent(out) :: nswap, n2x2, npos, nneg, ntl, nlr, nsk, rmax, ierr
5312 integer(kind=8),
intent(out) :: rsum
5313 integer(kind=kint),
allocatable :: perm(:), itmp(:), ip(:), jp(:)
5314 integer(kind=kint) :: npiv, nfs, k, p0, pa, pb, w, m, jj, np, nsw, n22, info, x, i, j, hi, wj, pk, pt
5315 integer(kind=kint) :: np2, l2, mn, r, ra, rb, rcp
5316 integer(kind=8) :: oik, ojk, oij, okk, btop, ob, oa, oav, ob2, obv, iw
5317 real(kind=kreal) :: d11, d21, d22, det
5318 logical,
allocatable :: skp(:)
5333 call mf_grow_i(brk, mf_bidx(g%nt, g%ntc, g%nt))
5334 call mf_grow_i8(boff, mf_bidx(g%nt, g%ntc, g%nt))
5339 w = g%tb(k) - g%tb(k-1)
5340 do i = g%ntc+1, g%nt
5341 hi = g%tb(i) - g%tb(i-1)
5342 ob = ob + int(hi, 8)*w + 2_8*int(w, 8)*min(hi, w)
5345 call mf_grow_r(bval, max(ob, 1_8))
5347 brk(1:mf_bidx(g%nt, g%ntc, g%nt)) = -1
5351 sn%ptype(1:g%ncol) = 0
5352 sn%dsub(1:g%ncol) = 0.0d0
5355 allocate(perm(g%ncol), itmp(g%ncol), ip(g%nt*(g%nt+1)/2), jp(g%nt*(g%nt+1)/2))
5356 call mf_grow_r(pval, int(g%nrow, 8))
5357 call mf_grow_r(wk, int(2*g%nrow, 8))
5362 do while (pa <= min(g%tb(k), nfs))
5363 pb = min(g%tb(k), nfs)
5366 call mf_grow_r(pval, int(m, 8)*w)
5367 call mf_grow_r(wk, int(2*m, 8))
5368 call fill_panel(npiv > p0)
5369 call hecmw_mf_kernel_panel_nopiv(m, w, m, pval, u, zero, info)
5373 call mf_col_copy(g, fval, pa+jj-1, pa+jj-1, pval(int(jj-1, 8)*m + jj), .true.)
5379 call fill_panel(npiv > p0)
5380 call hecmw_mf_kernel_panel_piv(m, w, m, pval, blk(pa:pb), u, mf_pivot_alpha, zero, &
5381 .true., wk, np, perm, sn%ptype(pa:pb), sn%dsub(pa:pb), nsw, n22, info)
5382 call permute_rows(pa, w)
5384 call mf_col_copy(g, fval, pa+jj-1, pa+jj-1, pval(int(jj-1, 8)*m + jj), .true.)
5392 call mf_swap(g, fval, x, nfs)
5393 call swap_pos(x, nfs)
5408 if (blr .and. pk > 0)
then
5414 do i = g%ntc+1, g%nt
5415 skp(i) = mf_blr_skip(fct, s, g%ntc, k, i)
5416 if (.not. skp(i) .and. ruse) skp(i) = sn%brank(mf_bidx(g%nt, k, i)) < 0
5418 hi = g%tb(i) - g%tb(i-1)
5419 boff(mf_bidx(g%nt, k, i)) = btop
5420 btop = btop + int(hi, 8)*pk + 2_8*int(pk, 8)*min(hi, pk)
5422 call mf_grow_r(bval, btop)
5424 do i = g%ntc+1, g%nt
5427 hi = g%tb(i) - g%tb(i-1)
5428 ob = boff(mf_bidx(g%nt, k, i))
5429 oik = mf_off(g, i, k)
5430 bval(ob+1:ob+int(hi, 8)*pk) = fval(oik+1:oik+int(hi, 8)*pk)
5431 rcp = int(fct%blr_beta*real((hi*pk - 1)/(hi + pk), kind=kreal))
5432 call hecmw_mf_kernel_compress(hi, pk, hi, bval(ob+1), eps, pk, bval(ob+int(hi, 8)*pk+1), &
5434 brk(mf_bidx(g%nt, k, i)) = r
5438 do i = g%ntc+1, g%nt
5444 r = brk(mf_bidx(g%nt, k, i))
5453 okk = mf_off(g, k, k)
5457 call mf_grow_r(wval, int(g%nrow - g%tb(k), 8)*pk)
5461 hi = g%tb(i) - g%tb(i-1)
5462 oik = mf_off(g, i, k)
5464 if (i > g%ntc) r = brk(mf_bidx(g%nt, k, i))
5466 call hecmw_mf_kernel_scale(hi, pk, g%tb(k) - g%tb(k-1), fval(okk+1), sn%ptype(p0+1:npiv), &
5467 sn%dsub(p0+1:npiv), hi, fval(oik+1), wval(int(g%tb(i-1) - g%tb(k), 8)*pk + 1))
5468 else if (r > 0)
then
5470 ob = boff(mf_bidx(g%nt, k, i))
5471 call hecmw_mf_kernel_scale_rows(pk, r, g%tb(k) - g%tb(k-1), fval(okk+1), sn%ptype(p0+1:npiv), &
5472 sn%dsub(p0+1:npiv), pk, bval(ob+int(hi, 8)*pk+1), pk, bval(ob+int(hi, 8)*pk+int(pk, 8)*mn+1))
5490 hi = g%tb(i) - g%tb(i-1)
5491 wj = g%tb(j) - g%tb(j-1)
5492 oij = mf_off(g, i, j)
5493 iw = int(g%tb(i-1) - g%tb(k), 8)*pk
5497 if (i > g%ntc) ra = brk(mf_bidx(g%nt, k, i))
5499 oa = boff(mf_bidx(g%nt, k, i))
5500 oav = oa + int(hi, 8)*pk + int(pk, 8)*min(hi, pk)
5505 if (j > g%ntc) rb = brk(mf_bidx(g%nt, k, j))
5507 ob2 = boff(mf_bidx(g%nt, k, j))
5508 obv = ob2 + int(wj, 8)*pk
5510 call mf_update_ab(hi, wj, pk, ra, wval(iw+1), bval(oa+1), bval(oav+1), &
5511 rb, fval(mf_off(g, j, k)+1), bval(ob2+1), bval(obv+1), fval(oij+1))
5516 call mf_grow_r(wval, int(g%nrow - g%tb(k), 8)*pk)
5520 hi = g%tb(i) - g%tb(i-1)
5521 oik = mf_off(g, i, k)
5522 call hecmw_mf_kernel_scale(hi, pk, g%tb(k) - g%tb(k-1), fval(okk+1), sn%ptype(p0+1:npiv), &
5523 sn%dsub(p0+1:npiv), hi, fval(oik+1), wval(int(g%tb(i-1) - g%tb(k), 8)*pk + 1))
5540 hi = g%tb(i) - g%tb(i-1)
5541 wj = g%tb(j) - g%tb(j-1)
5542 ojk = mf_off(g, j, k)
5543 oij = mf_off(g, i, j)
5544 call hecmw_mf_kernel_gemm(hi, wj, pk, hi, wval(int(g%tb(i-1) - g%tb(k), 8)*pk + 1), wj, &
5545 fval(ojk+1), hi, fval(oij+1))
5551 hi = g%tb(i) - g%tb(i-1)
5552 oik = mf_off(g, i, k)
5553 call mf_grow_r(wval, int(hi, 8)*pk)
5554 call hecmw_mf_kernel_scale(hi, pk, g%tb(k) - g%tb(k-1), fval(okk+1), sn%ptype(p0+1:npiv), &
5555 sn%dsub(p0+1:npiv), hi, fval(oik+1), wval)
5557 wj = g%tb(j) - g%tb(j-1)
5558 ojk = mf_off(g, j, k)
5559 oij = mf_off(g, i, j)
5560 call hecmw_mf_kernel_gemm(hi, wj, pk, hi, wval, wj, fval(ojk+1), hi, fval(oij+1))
5564 do x = npiv + 1, g%tb(k)
5565 call mf_col_copy(g, fval, x, x, pval, .false.)
5566 call mf_col_update(g, fval, sn%ptype, sn%dsub, x, p0+1, npiv, pval, wk)
5567 call mf_col_copy(g, fval, x, x, pval, .true.)
5570 if (nfs <= g%tb(k))
exit
5573 if (isroot .and. npiv < g%ncol)
then
5578 call mf_grow_r(pval, int(m, 8)*w)
5579 call mf_grow_r(wk, int(2*m, 8))
5580 call fill_panel(.false.)
5581 call hecmw_mf_kernel_panel_piv(m, w, m, pval, blk(pa:pb), u, mf_pivot_alpha, zero, &
5582 .false., wk, np, perm, sn%ptype(pa:pb), sn%dsub(pa:pb), nsw, n22, info)
5583 call permute_rows(pa, w)
5587 ierr = sn%fsdof(pa + info - 1)
5588 deallocate(perm, itmp, ip, jp)
5592 call mf_col_copy(g, fval, pa+jj-1, pa+jj-1, pval(int(jj-1, 8)*m + jj), .true.)
5599 do while (x <= npiv)
5602 if (fval(mf_idx(g, x, x)) > 0.0d0)
then
5609 d11 = fval(mf_idx(g, x, x))
5611 d22 = fval(mf_idx(g, x+1, x+1))
5612 det = d11*d22 - d21*d21
5613 if (det < 0.0d0)
then
5616 else if (d11 + d22 > 0.0d0)
then
5624 deallocate(perm, itmp, ip, jp)
5630 subroutine fill_panel(upd)
5631 logical,
intent(in) :: upd
5632 integer(kind=kint) :: jj
5635 call mf_col_copy(g, fval, pa+jj-1, pa+jj-1, pval(int(jj-1, 8)*m + jj), .false.)
5637 call mf_col_update(g, fval, sn%ptype, sn%dsub, pa+jj-1, p0+1, npiv, &
5638 pval(int(jj-1, 8)*m + jj), wk)
5640 end subroutine fill_panel
5644 subroutine permute_rows(pa, w)
5645 integer(kind=kint),
intent(in) :: pa, w
5646 integer(kind=kint) :: jj, t, q
5652 if (itmp(jj) == perm(jj)) cycle
5654 if (itmp(t) == perm(jj))
exit
5656 call mf_swap(g, fval, pa+jj-1, pa+t-1)
5657 q = sn%fsdof(pa+jj-1)
5658 sn%fsdof(pa+jj-1) = sn%fsdof(pa+t-1)
5659 sn%fsdof(pa+t-1) = q
5663 end subroutine permute_rows
5665 subroutine swap_pos(x, y)
5666 integer(kind=kint),
intent(in) :: x, y
5667 integer(kind=kint) :: t
5670 sn%fsdof(x) = sn%fsdof(y)
5675 end subroutine swap_pos
5677 end subroutine mf_factor_front
5684 subroutine mf_factor_front_lu(sn, g, fval, fvalu, pval, pvalu, wk, blk, blkr, isroot, u, zero, &
5685 blr, eps, bval, boff, brk, bvalu, boffu, brku, fct, s, ruse, nswap, ntl, nlr, nsk, &
5688 type(mf_snode),
intent(inout) :: sn
5689 type(mf_grid),
intent(in) :: g
5690 real(kind=kreal),
allocatable,
intent(inout) :: fval(:), fvalu(:)
5691 real(kind=kreal),
allocatable,
intent(inout) :: pval(:), pvalu(:), wk(:)
5692 integer(kind=kint),
intent(inout) :: blk(:), blkr(:)
5693 logical,
intent(in) :: isroot, blr
5694 real(kind=kreal),
intent(in) :: u, zero, eps
5695 real(kind=kreal),
allocatable,
intent(inout) :: bval(:), bvalu(:)
5696 integer(kind=8),
allocatable,
intent(inout) :: boff(:), boffu(:)
5697 integer(kind=kint),
allocatable,
intent(inout) :: brk(:), brku(:)
5698 type(hecmwst_mf_factor),
intent(in) :: fct
5699 integer(kind=kint),
intent(in) :: s
5700 logical,
intent(in) :: ruse
5701 integer(kind=kint),
intent(out) :: nswap, ntl, nlr, nsk, rmax, ierr
5702 integer(kind=8),
intent(out) :: rsum
5703 integer(kind=kint),
allocatable :: permc(:), permr(:), itmp(:), ip(:), jp(:)
5704 integer(kind=kint) :: npiv, nfs, k, p0, pa, pb, w, m, np, nsw, info, x, i, j, hi, wj, pk
5705 integer(kind=kint) :: np2, l2, r, rla, rlb, rua, rub, rcp
5706 integer(kind=8) :: oik, ojk, oij, btop, btopu, ob, ola, olav, olb, olbv, oua, ouav, oub, oubv
5707 logical,
allocatable :: skpl(:), skpu(:)
5720 call mf_grow_i(brk, mf_bidx(g%nt, g%ntc, g%nt))
5721 call mf_grow_i(brku, mf_bidx(g%nt, g%ntc, g%nt))
5722 call mf_grow_i8(boff, mf_bidx(g%nt, g%ntc, g%nt))
5723 call mf_grow_i8(boffu, mf_bidx(g%nt, g%ntc, g%nt))
5728 w = g%tb(k) - g%tb(k-1)
5729 do i = g%ntc+1, g%nt
5730 hi = g%tb(i) - g%tb(i-1)
5731 ob = ob + int(hi, 8)*w + int(w, 8)*min(hi, w)
5734 call mf_grow_r(bval, max(ob, 1_8))
5735 call mf_grow_r(bvalu, max(ob, 1_8))
5736 allocate(skpl(g%nt), skpu(g%nt))
5737 brk(1:mf_bidx(g%nt, g%ntc, g%nt)) = -1
5738 brku(1:mf_bidx(g%nt, g%ntc, g%nt)) = -1
5742 sn%ptype(1:g%ncol) = 1
5743 sn%dsub(1:g%ncol) = 0.0d0
5746 allocate(permc(g%ncol), permr(g%ncol), itmp(g%ncol), ip(g%nt*(g%nt+1)/2), jp(g%nt*(g%nt+1)/2))
5747 call mf_grow_r(pval, int(g%nrow, 8))
5748 call mf_grow_r(pvalu, int(g%nrow, 8))
5749 call mf_grow_r(wk, int(2*g%nrow, 8))
5754 do while (pa <= min(g%tb(k), nfs))
5755 pb = min(g%tb(k), nfs)
5758 call mf_grow_r(pval, int(m, 8)*w)
5759 call mf_grow_r(pvalu, int(m, 8)*w)
5760 call fill_panel(npiv > p0)
5761 call hecmw_mf_kernel_panel_lu_nopiv(m, w, m, pval, m, pvalu, u, zero, info)
5768 call fill_panel(npiv > p0)
5769 call hecmw_mf_kernel_panel_lu_piv(m, w, m, pval, m, pvalu, blk(pa:pb), blkr(pa:pb), u, &
5770 zero, .true., np, permc, permr, nsw, info)
5771 call permute_front(pa, w)
5778 call mf_swap_lu(g, fval, fvalu, x, nfs)
5779 call swap_pos(x, nfs)
5793 if (blr .and. pk > 0)
then
5796 do i = g%ntc+1, g%nt
5797 skpl(i) = mf_blr_skip(fct, s, g%ntc, k, i)
5800 if (.not. skpl(i)) skpl(i) = sn%brank(mf_bidx(g%nt, k, i)) < 0
5801 if (.not. skpu(i)) skpu(i) = sn%branku(mf_bidx(g%nt, k, i)) < 0
5803 hi = g%tb(i) - g%tb(i-1)
5804 if (.not. skpl(i))
then
5805 boff(mf_bidx(g%nt, k, i)) = btop
5806 btop = btop + int(hi, 8)*pk + int(pk, 8)*min(hi, pk)
5808 if (.not. skpu(i))
then
5809 boffu(mf_bidx(g%nt, k, i)) = btopu
5810 btopu = btopu + int(hi, 8)*pk + int(pk, 8)*min(hi, pk)
5813 call mf_grow_r(bval, btop)
5814 call mf_grow_r(bvalu, btopu)
5816 do i = g%ntc+1, g%nt
5817 if (skpl(i) .and. skpu(i)) cycle
5819 hi = g%tb(i) - g%tb(i-1)
5820 oik = mf_off(g, i, k)
5821 rcp = int(fct%blr_beta*real((hi*pk - 1)/(hi + pk), kind=kreal))
5822 if (.not. skpl(i))
then
5823 ob = boff(mf_bidx(g%nt, k, i))
5824 bval(ob+1:ob+int(hi, 8)*pk) = fval(oik+1:oik+int(hi, 8)*pk)
5825 call hecmw_mf_kernel_compress(hi, pk, hi, bval(ob+1), eps, pk, bval(ob+int(hi, 8)*pk+1), &
5827 brk(mf_bidx(g%nt, k, i)) = r
5829 if (.not. skpu(i))
then
5830 ob = boffu(mf_bidx(g%nt, k, i))
5831 bvalu(ob+1:ob+int(hi, 8)*pk) = fvalu(oik+1:oik+int(hi, 8)*pk)
5832 call hecmw_mf_kernel_compress(hi, pk, hi, bvalu(ob+1), eps, pk, bvalu(ob+int(hi, 8)*pk+1), &
5834 brku(mf_bidx(g%nt, k, i)) = r
5839 do i = g%ntc+1, g%nt
5850 r = brk(mf_bidx(g%nt, k, i))
5856 r = brku(mf_bidx(g%nt, k, i))
5880 hi = g%tb(i) - g%tb(i-1)
5881 wj = g%tb(j) - g%tb(j-1)
5882 oik = mf_off(g, i, k)
5883 ojk = mf_off(g, j, k)
5884 oij = mf_off(g, i, j)
5892 rla = brk(mf_bidx(g%nt, k, i))
5893 rua = brku(mf_bidx(g%nt, k, i))
5895 ola = boff(mf_bidx(g%nt, k, i))
5896 olav = ola + int(hi, 8)*pk
5899 oua = boffu(mf_bidx(g%nt, k, i))
5900 ouav = oua + int(hi, 8)*pk
5910 rlb = brk(mf_bidx(g%nt, k, j))
5911 rub = brku(mf_bidx(g%nt, k, j))
5913 olb = boff(mf_bidx(g%nt, k, j))
5914 olbv = olb + int(wj, 8)*pk
5917 oub = boffu(mf_bidx(g%nt, k, j))
5918 oubv = oub + int(wj, 8)*pk
5921 call mf_update_ab(hi, wj, pk, rla, fval(oik+1), bval(ola+1), bval(olav+1), &
5922 rub, fvalu(ojk+1), bvalu(oub+1), bvalu(oubv+1), fval(oij+1))
5923 call mf_update_ab(hi, wj, pk, rua, fvalu(oik+1), bvalu(oua+1), bvalu(ouav+1), &
5924 rlb, fval(ojk+1), bval(olb+1), bval(olbv+1), fvalu(oij+1))
5934 hi = g%tb(i) - g%tb(i-1)
5935 wj = g%tb(j) - g%tb(j-1)
5936 oik = mf_off(g, i, k)
5937 ojk = mf_off(g, j, k)
5938 oij = mf_off(g, i, j)
5939 call hecmw_mf_kernel_gemm(hi, wj, pk, hi, fval(oik+1), wj, fvalu(ojk+1), hi, fval(oij+1))
5940 call hecmw_mf_kernel_gemm(hi, wj, pk, hi, fvalu(oik+1), wj, fval(ojk+1), hi, fvalu(oij+1))
5945 do x = npiv + 1, g%tb(k)
5946 call update_pos(x, p0 + 1, npiv)
5949 if (nfs <= g%tb(k))
exit
5952 if (isroot .and. npiv < g%ncol)
then
5957 call mf_grow_r(pval, int(m, 8)*w)
5958 call mf_grow_r(pvalu, int(m, 8)*w)
5959 call fill_panel(.false.)
5960 call hecmw_mf_kernel_panel_lu_piv(m, w, m, pval, m, pvalu, blk(pa:pb), blkr(pa:pb), u, &
5961 zero, .false., np, permc, permr, nsw, info)
5962 call permute_front(pa, w)
5965 ierr = sn%fsdof(pa + info - 1)
5966 deallocate(permc, permr, itmp, ip, jp)
5973 deallocate(permc, permr, itmp, ip, jp)
5980 subroutine fill_panel(upd)
5981 logical,
intent(in) :: upd
5982 integer(kind=kint) :: jj, xx
5986 call mf_col_copy(g, fval, xx, xx, pval(int(jj-1, 8)*m + jj), .false.)
5988 call mf_row_get(g, fvalu, xx, p0+1, npiv, wk)
5989 call mf_col_axpy(g, fval, xx, p0+1, npiv, wk, pval(int(jj-1, 8)*m + jj))
5991 if (xx < g%nrow)
then
5992 call mf_col_copy(g, fvalu, xx, xx+1, pvalu(int(jj-1, 8)*m + jj + 1), .false.)
5994 call mf_row_get(g, fval, xx, p0+1, npiv, wk)
5995 call mf_col_axpy(g, fvalu, xx+1, p0+1, npiv, wk, pvalu(int(jj-1, 8)*m + jj + 1))
5999 end subroutine fill_panel
6002 subroutine write_back(n)
6003 integer(kind=kint),
intent(in) :: n
6004 integer(kind=kint) :: jj, xx
6008 call mf_col_copy(g, fval, xx, xx, pval(int(jj-1, 8)*m + jj), .true.)
6009 if (xx < g%nrow)
call mf_col_copy(g, fvalu, xx, xx+1, pvalu(int(jj-1, 8)*m + jj + 1), .true.)
6011 end subroutine write_back
6014 subroutine update_pos(x, q1, q2)
6015 integer(kind=kint),
intent(in) :: x, q1, q2
6017 call mf_col_copy(g, fval, x, x, pval, .false.)
6018 call mf_row_get(g, fvalu, x, q1, q2, wk)
6019 call mf_col_axpy(g, fval, x, q1, q2, wk, pval)
6020 call mf_col_copy(g, fval, x, x, pval, .true.)
6021 if (x < g%nrow)
then
6022 call mf_col_copy(g, fvalu, x, x+1, pvalu, .false.)
6023 call mf_row_get(g, fval, x, q1, q2, wk)
6024 call mf_col_axpy(g, fvalu, x+1, q1, q2, wk, pvalu)
6025 call mf_col_copy(g, fvalu, x, x+1, pvalu, .true.)
6027 end subroutine update_pos
6031 subroutine permute_front(pa, w)
6032 integer(kind=kint),
intent(in) :: pa, w
6033 integer(kind=kint) :: jj, t, q
6039 if (itmp(jj) == permc(jj)) cycle
6041 if (itmp(t) == permc(jj))
exit
6043 call mf_swap_lu(g, fval, fvalu, pa+jj-1, pa+t-1)
6044 call swap_pos(pa+jj-1, pa+t-1)
6046 itmp(jj) = permc(jj)
6048 itmp(1:w) = permc(1:w)
6050 if (itmp(jj) == permr(jj)) cycle
6052 if (itmp(t) == permr(jj))
exit
6054 call mf_swap_row(g, fval, fvalu, pa+jj-1, pa+t-1)
6055 q = sn%frow(pa+jj-1)
6056 sn%frow(pa+jj-1) = sn%frow(pa+t-1)
6059 blkr(pa+jj-1) = blkr(pa+t-1)
6062 itmp(jj) = permr(jj)
6064 end subroutine permute_front
6066 subroutine swap_pos(x, y)
6067 integer(kind=kint),
intent(in) :: x, y
6068 integer(kind=kint) :: t
6071 sn%fsdof(x) = sn%fsdof(y)
6074 sn%frow(x) = sn%frow(y)
6082 end subroutine swap_pos
6084 end subroutine mf_factor_front_lu
6087 subroutine mf_stored_grid(sn, g, tw)
6089 type(mf_snode),
intent(in) :: sn
6090 type(mf_grid),
intent(inout) :: g
6091 integer(kind=kint),
intent(out) :: tw(:)
6092 integer(kind=kint) :: j
6097 if (
allocated(g%tb))
then
6098 if (
size(g%tb) < g%nt + 1)
deallocate(g%tb, g%coloff)
6100 if (.not.
allocated(g%tb))
allocate(g%tb(0:g%nt), g%coloff(g%nt+1))
6101 g%tb(0:g%nt) = sn%tbnd(1:g%nt+1)
6105 tw(j) = max(min(g%tb(j), sn%npiv) - g%tb(j-1), 0)
6106 g%coloff(j+1) = g%coloff(j) + int(g%nrow - g%tb(j-1), 8)*tw(j)
6108 end subroutine mf_stored_grid
6116 type(hecmwst_mf_symbolic),
intent(in) :: sym
6117 type(hecmwst_mf_factor),
intent(inout) :: fct
6118 real(kind=kreal),
intent(in) :: b(:)
6119 real(kind=kreal),
intent(out) :: x(:)
6120 real(kind=kreal),
allocatable :: z(:), xx(:)
6121 integer(kind=kint),
allocatable :: map(:), left(:)
6122 integer(kind=kint) :: ns, s, i
6125 allocate(z(fct%ndof_tot), xx(fct%ndof_tot), map(fct%ndof_tot), left(ns))
6127 do i = 1, fct%ndof_tot
6128 z(i) = b(fct%pdof(i))
6133 left(s) = fct%cptr(s+1) - fct%cptr(s)
6138 if (fct%cptr(s+1) == fct%cptr(s))
then
6140 call mf_fwd_task(sym, fct, left, s, z, map)
6150 if (sym%sparent(s) == 0)
then
6160 do i = 1, fct%ndof_tot
6161 x(fct%pdof(i)) = xx(i)
6164 deallocate(z, xx, map, left)
6171 recursive subroutine bwd_task(s0)
6172 integer(kind=kint),
intent(in) :: s0
6173 integer(kind=kint) :: j, c
6175 call mf_super_bwd(sym, fct, fct%sn(s0), s0, z, xx)
6176 do j = fct%cptr(s0), fct%cptr(s0+1) - 1
6182 end subroutine bwd_task
6188 subroutine mf_fwd_task(sym, fct, left, s, z, map)
6190 type(hecmwst_mf_symbolic),
intent(in) :: sym
6191 type(hecmwst_mf_factor),
intent(inout) :: fct
6192 integer(kind=kint),
intent(inout) :: left(:)
6193 integer(kind=kint),
intent(in) :: s
6194 real(kind=kreal),
intent(inout) :: z(:)
6195 integer(kind=kint),
intent(inout) :: map(:)
6196 integer(kind=kint) :: p, n, ss
6200 call mf_super_fwd(sym, fct, fct%sn(ss), ss, z, map)
6206 left(p) = left(p) - 1
6212 end subroutine mf_fwd_task
6226 type(hecmwst_mf_symbolic),
intent(in) :: sym
6227 type(hecmwst_mf_map),
intent(in) :: map
6228 type(hecmwst_mf_factor),
intent(inout) :: fct
6229 real(kind=kreal),
intent(in) :: b(:)
6230 real(kind=kreal),
intent(out) :: x(:)
6231 type(mf_pool) :: pool, gpool
6232 real(kind=kreal),
allocatable :: z(:), xx(:)
6233 real(kind=kreal),
allocatable :: v(:), vb(:), tok(:), tbuf(:), vsg(:), xv(:)
6234 integer(kind=kint),
allocatable :: dmap(:), left(:)
6235 integer(kind=kint),
allocatable :: ctb(:), towner(:), crow2t(:), tprow(:), wrank(:), rdof(:), rowoff(:)
6236 integer(kind=kint),
allocatable :: cctb2(:), cown2(:), pmapr(:), pdofm(:), crown(:), mdel(:)
6237 integer(kind=kint),
allocatable :: cur5(:), stw(:), sta(:), stb(:)
6238 logical,
allocatable :: held(:)
6239 integer(kind=kint) :: ns, s, i, iu, j, c, p, stat(hecmw_status_size), req
6240 integer(kind=kint) :: ncol, npv, ntc0, nt0, nrow0, ncb0, ncbt, nwk, wme, myrows, master, ndel0
6241 integer(kind=kint) :: cnrow0, cndel0, cncbt
6242 integer(kind=kint) :: pr2, pc2, kend, sg, myfst
6243 logical :: ismaster, hasband
6246 allocate(z(fct%ndof_tot), xx(fct%ndof_tot), dmap(fct%ndof_tot), left(ns))
6248 do i = 1, fct%ndof_tot
6249 z(i) = b(fct%pdof(i))
6254 left(s) = fct%cptr(s+1) - fct%cptr(s)
6259 if (fct%cptr(s+1) == fct%cptr(s) .and. map%owner(s) == map%myrank .and. .not. map%upper(s))
then
6261 call mf_fwd_task_mpi(sym, map, fct, left, s, z, dmap)
6268 do iu = 1, map%nupper
6270 if (map%myrank >= map%rbeg(s) .and. map%myrank < map%rbeg(s) + map%rcnt(s))
call fwd_dist(s)
6273 do iu = map%nupper, 1, -1
6275 if (map%myrank >= map%rbeg(s) .and. map%myrank < map%rbeg(s) + map%rcnt(s))
call bwd_dist(s)
6278 if (map%owner(s) /= map%myrank .or. map%upper(s) .or. sym%sparent(s) == 0) cycle
6279 if (map%upper(sym%sparent(s)) .and. map%owner(sym%sparent(s)) /= map%myrank)
call recv_xx(s)
6284 if (map%owner(s) /= map%myrank .or. map%upper(s)) cycle
6287 if (.not. map%upper(p)) cycle
6290 call bwd_task_mpi(s)
6295 call mf_pool_wait(gpool)
6297 z(1:fct%ndof_tot) = 0.0d0
6299 if (map%owner(s) /= map%myrank) cycle
6300 do i = 1, fct%sn(s)%npiv
6301 z(fct%pdof(fct%sn(s)%fsdof(i))) = xx(fct%sn(s)%fsdof(i))
6304 x(1:fct%ndof_tot) = z(1:fct%ndof_tot)
6305 call hecmw_allreduce_r_comm(x, fct%ndof_tot, hecmw_sum, map%comm)
6306 deallocate(z, xx, dmap, left)
6311 subroutine geo_1d(s0)
6312 integer(kind=kint),
intent(in) :: s0
6313 integer(kind=kint) :: t2, i2, k2, d2, nown2, a2
6315 ncol = fct%sn(s0)%ncol
6316 npv = fct%sn(s0)%npiv
6317 ntc0 = fct%sn(s0)%ntc
6319 nrow0 = fct%sn(s0)%tbnd(nt0+1)
6321 master = map%owner(s0)
6322 ismaster = (master == map%myrank)
6323 ndel0 = ncol - (fct%cdofptr(sym%sptr(s0+1)) - fct%cdofptr(sym%sptr(s0)))
6324 call hecmw_mf_dist_rowtiles(sym, map, fct%tile, s0, ndel0, ncbt, ctb, towner)
6325 if (
allocated(crow2t))
then
6326 if (
size(crow2t) < ncb0)
deallocate(crow2t)
6328 if (.not.
allocated(crow2t))
allocate(crow2t(max(ncb0, 1)))
6330 crow2t(ctb(t2-1)+1:ctb(t2)) = t2
6332 if (
allocated(wrank))
then
6333 if (
size(wrank) < ncbt + 1)
deallocate(wrank, tprow)
6335 if (.not.
allocated(wrank))
allocate(wrank(max(ncbt, 1)), tprow(0:max(ncbt, 1)))
6343 wrank(1) = towner(t2)
6344 else if (towner(t2) /= wrank(nwk))
then
6346 wrank(nwk) = towner(t2)
6348 if (wrank(nwk) == map%myrank) wme = nwk
6350 if (towner(t2) == map%myrank) myrows = myrows + ctb(t2) - ctb(t2-1)
6352 hasband = (myrows > 0)
6353 nown2 = sym%sptr(s0+1) - sym%sptr(s0)
6354 if (
allocated(rowoff))
then
6355 if (
size(rowoff) < sym%rptr(s0+1) - sym%rptr(s0) + 1)
deallocate(rowoff)
6357 if (.not.
allocated(rowoff))
allocate(rowoff(sym%rptr(s0+1) - sym%rptr(s0) + 1))
6359 do i2 = 1, sym%rptr(s0+1) - sym%rptr(s0)
6360 k2 = sym%rlist(sym%rptr(s0)+i2-1)
6361 rowoff(i2+1) = rowoff(i2) + sym%ndof(k2)
6362 if (i2 == nown2) rowoff(i2+1) = rowoff(i2+1) + ndel0
6365 if (
allocated(rdof))
then
6366 if (
size(rdof) < ncb0)
deallocate(rdof)
6368 if (.not.
allocated(rdof))
allocate(rdof(max(ncb0, 1)))
6370 do k2 = sym%rptr(s0) + nown2, sym%rptr(s0+1) - 1
6371 do a2 = 1, sym%ndof(sym%rlist(k2))
6373 rdof(i2) = fct%cdofptr(sym%rlist(k2)) + a2 - 1
6377 if (
allocated(mdel))
then
6378 if (
size(mdel) < fct%cptr(s0+1) - fct%cptr(s0))
deallocate(mdel)
6380 if (.not.
allocated(mdel))
allocate(mdel(max(fct%cptr(s0+1) - fct%cptr(s0), 1)))
6381 d2 = fct%cdofptr(sym%sptr(s0+1)) - fct%cdofptr(sym%sptr(s0))
6384 do j = fct%cptr(s0), fct%cptr(s0+1) - 1
6385 if (fct%clist(j) == c) mdel(j - fct%cptr(s0) + 1) = d2
6387 d2 = d2 + fct%sn(c)%ncol - fct%sn(c)%npiv
6392 call hecmw_mf_dist_fsgrid(map, s0, ntc0, pr2, pc2)
6395 if (min(fct%sn(s0)%tbnd(t2+1), npv) - fct%sn(s0)%tbnd(t2) > 0) kend = t2
6400 if (fso_s(i2, t2) == map%myrank) myfst = myfst + 1
6403 if (
allocated(vsg))
then
6404 if (
size(vsg) < ncol)
deallocate(vsg)
6406 if (.not.
allocated(vsg))
allocate(vsg(max(ncol, 1)))
6407 if (
allocated(held))
then
6408 if (
size(held) < ntc0)
deallocate(held)
6410 if (.not.
allocated(held))
allocate(held(max(ntc0, 1)))
6411 held(1:ntc0) = .false.
6412 if (
allocated(cur5))
then
6413 if (
size(cur5) < map%rcnt(s0))
deallocate(cur5)
6415 if (.not.
allocated(cur5))
allocate(cur5(max(map%rcnt(s0), 1)))
6416 if (
allocated(stw))
then
6417 if (
size(stw) < ntc0 + nwk + 2)
deallocate(stw, sta, stb)
6419 if (.not.
allocated(stw))
allocate(stw(ntc0 + nwk + 2), sta(ntc0 + nwk + 2), stb(ntc0 + nwk + 2))
6420 call mf_grow_r(xv, int(max(fct%tile, ncol, 1), 8))
6421 end subroutine geo_1d
6424 function fso_s(i0, j0)
result(r)
6425 integer(kind=kint),
intent(in) :: i0, j0
6426 integer(kind=kint) :: r
6428 r = hecmw_mf_dist_fsrank(map, sg, pr2, pc2, i0, j0)
6432 function drank_s(row)
result(dr)
6433 integer(kind=kint),
intent(in) :: row
6434 integer(kind=kint) :: dr
6436 if (row <= ncol)
then
6439 dr = towner(crow2t(row - ncol))
6441 end function drank_s
6446 subroutine crows_setup(s0, c0, jc)
6447 integer(kind=kint),
intent(in) :: s0, c0, jc
6448 integer(kind=kint) :: rr, ll, kk2, aa, nde, nown2, cnown2
6450 cndel0 = fct%sn(c0)%ncol - fct%sn(c0)%npiv
6451 cnrow0 = cndel0 + fct%sn(c0)%tbnd(fct%sn(c0)%nt+1) - fct%sn(c0)%ncol
6452 if (
allocated(pmapr))
then
6453 if (
size(pmapr) < cnrow0)
deallocate(pmapr, pdofm, crown)
6455 if (.not.
allocated(pmapr))
allocate(pmapr(max(cnrow0, 1)), pdofm(max(cnrow0, 1)), crown(max(cnrow0, 1)))
6457 if (map%upper(c0))
then
6458 nde = fct%sn(c0)%ncol - (fct%cdofptr(sym%sptr(c0+1)) - fct%cdofptr(sym%sptr(c0)))
6459 call hecmw_mf_dist_rowtiles(sym, map, fct%tile, c0, nde, cncbt, cctb2, cown2)
6460 do rr = 1, cnrow0 - cndel0
6462 if (rr <= cctb2(ll))
exit
6464 crown(cndel0 + rr) = cown2(ll)
6467 crown(cndel0+1:cnrow0) = map%owner(c0)
6469 crown(1:cndel0) = map%owner(c0)
6470 nown2 = sym%sptr(s0+1) - sym%sptr(s0)
6471 cnown2 = sym%sptr(c0+1) - sym%sptr(c0)
6473 pmapr(rr) = mdel(jc) + rr
6474 pdofm(rr) = fct%sn(c0)%frow(fct%sn(c0)%npiv + rr)
6477 do ll = sym%rptr(c0) + cnown2, sym%rptr(c0+1) - 1
6478 kk2 = sym%cmap(sym%cmap_ptr(c0) + (ll - sym%rptr(c0) - cnown2))
6479 do aa = 1, sym%ndof(sym%rlist(ll))
6481 pmapr(rr) = rowoff(kk2) + aa
6482 if (kk2 <= nown2)
then
6483 pdofm(rr) = fct%cdofptr(sym%rlist(sym%rptr(s0)+kk2-1)) + aa - 1
6489 end subroutine crows_setup
6493 subroutine dv_send(c0)
6494 integer(kind=kint),
intent(in) :: c0
6495 integer(kind=kint) :: dst, rr, nsl, il, ib
6497 if (.not.
allocated(fct%sn(c0)%dval))
return
6498 do dst = map%rbeg(s), map%rbeg(s) + map%rcnt(s) - 1
6499 if (dst == map%myrank) cycle
6502 if (crown(rr) == map%myrank .and. drank_s(pmapr(rr)) == dst) nsl = nsl + 1
6505 ib = mf_pool_slot(pool)
6506 allocate(pool%box(ib)%rv(nsl))
6509 if (crown(rr) /= map%myrank) cycle
6510 if (drank_s(pmapr(rr)) /= dst) cycle
6512 pool%box(ib)%rv(il) = fct%sn(c0)%dval(dv_local(c0, rr))
6515 call hecmw_isend_r(pool%box(ib)%rv, nsl, dst, 16*c0+8, map%comm, req)
6516 call mf_pool_req(pool, req)
6518 end subroutine dv_send
6521 function dv_local(c0, rr)
result(il)
6522 integer(kind=kint),
intent(in) :: c0, rr
6523 integer(kind=kint) :: il, t2, ll
6525 if (.not. map%upper(c0))
then
6527 else if (rr <= cndel0)
then
6533 if (rr - cndel0 <= cctb2(t2))
exit
6536 if (map%owner(c0) == map%myrank) il = cndel0
6538 if (cown2(ll) == map%myrank) il = il + cctb2(ll) - cctb2(ll-1)
6540 il = il + (rr - cndel0 - cctb2(t2-1))
6542 end function dv_local
6546 subroutine dv_apply(s0, c0)
6547 integer(kind=kint),
intent(in) :: s0, c0
6548 real(kind=kreal),
allocatable :: rb(:)
6549 integer(kind=kint) :: is0, src, nsl, il, rr, nsrc
6550 integer(kind=kint) :: srcs(cnrow0 + 1)
6555 if (srcs(is0) == crown(rr))
exit
6557 if (is0 > nsrc)
then
6559 srcs(nsrc) = crown(rr)
6566 if (crown(rr) == src .and. drank_s(pmapr(rr)) == map%myrank) nsl = nsl + 1
6569 if (src == map%myrank)
then
6571 if (crown(rr) /= map%myrank) cycle
6572 if (drank_s(pmapr(rr)) /= map%myrank) cycle
6573 call dv_add(rr, fct%sn(c0)%dval(dv_local(c0, rr)))
6577 call hecmw_recv_r(rb, nsl, src, 16*c0+8, map%comm, stat)
6580 if (crown(rr) /= src) cycle
6581 if (drank_s(pmapr(rr)) /= map%myrank) cycle
6583 call dv_add(rr, rb(il))
6588 end subroutine dv_apply
6591 subroutine dv_add(rr, val)
6592 integer(kind=kint),
intent(in) :: rr
6593 real(kind=kreal),
intent(in) :: val
6594 integer(kind=kint) :: pr, t2
6597 if (pr <= ncol)
then
6598 v(dmap(pdofm(rr))) = v(dmap(pdofm(rr))) + val
6600 t2 = crow2t(pr - ncol)
6601 vb(tprow(t2) + (pr - ncol - ctb(t2-1))) = vb(tprow(t2) + (pr - ncol - ctb(t2-1))) + val
6603 end subroutine dv_add
6612 subroutine fwd_dist(s0)
6613 integer(kind=kint),
intent(in) :: s0
6614 integer(kind=kint) :: d2, d0, d1, kk, tw2, hk2, t2, h2, r2, i2, idx2, jc, ib, ndv, ofs
6615 integer(kind=kint) :: rdk, nxt, dst
6616 integer(kind=8) :: okk2, oik2
6620 call mf_grow_r(v, int(max(ncol, 1), 8))
6621 call mf_grow_r(vb, int(max(myrows, 1), 8))
6622 call mf_grow_r(tbuf, int(max(ncol, nrow0 - ncol, 1), 8))
6624 d0 = fct%cdofptr(sym%sptr(s0))
6625 d1 = fct%cdofptr(sym%sptr(s0+1)) - 1
6627 dmap(fct%sn(s0)%frow(d2)) = d2
6628 if (fct%sn(s0)%frow(d2) >= d0 .and. fct%sn(s0)%frow(d2) <= d1)
then
6629 v(d2) = z(fct%sn(s0)%frow(d2))
6635 vb(1:max(myrows, 1)) = 0.0d0
6636 do j = fct%cptr(s0), fct%cptr(s0+1) - 1
6638 call crows_setup(s0, c, j - fct%cptr(s0) + 1)
6644 do j = fct%cptr(s0), fct%cptr(s0+1) - 1
6645 if (fct%clist(j) == c) jc = j - fct%cptr(s0) + 1
6647 call crows_setup(s0, c, jc)
6648 call dv_apply(s0, c)
6651 partic = ismaster .or. hasband .or. myfst > 0
6656 if (dst == me())
then
6657 vsg(fct%sn(s0)%tbnd(i2)+1:fct%sn(s0)%tbnd(i2+1)) = v(fct%sn(s0)%tbnd(i2)+1:fct%sn(s0)%tbnd(i2+1))
6660 ib = mf_pool_slot(pool)
6661 h2 = fct%sn(s0)%tbnd(i2+1) - fct%sn(s0)%tbnd(i2)
6662 allocate(pool%box(ib)%rv(h2))
6663 pool%box(ib)%rv(1:h2) = v(fct%sn(s0)%tbnd(i2)+1:fct%sn(s0)%tbnd(i2+1))
6665 call hecmw_isend_r(pool%box(ib)%rv, h2, dst, 16*s0+14, map%comm, req)
6666 call mf_pool_req(pool, req)
6669 else if (partic)
then
6671 if (fwd_fh(i2) /= me()) cycle
6672 h2 = fct%sn(s0)%tbnd(i2+1) - fct%sn(s0)%tbnd(i2)
6673 call hecmw_recv_r(vsg(fct%sn(s0)%tbnd(i2)+1:), h2, master, 16*s0+14, map%comm, stat)
6680 tw2 = min(fct%sn(s0)%tbnd(kk+1), npv) - fct%sn(s0)%tbnd(kk)
6681 hk2 = fct%sn(s0)%tbnd(kk+1) - fct%sn(s0)%tbnd(kk)
6683 if (rdk == me())
then
6684 if (.not. held(kk))
then
6685 call hecmw_recv_r(vsg(fct%sn(s0)%tbnd(kk)+1:), hk2, fwd_ph(kk, kk), 16*s0+13, map%comm, stat)
6688 okk2 = 1 + fct%sn(s0)%bptr(mf_bidx(nt0, kk, kk))
6689 call hecmw_mf_kernel_trsv(tw2, hk2, fct%sn(s0)%lval(okk2), vsg(fct%sn(s0)%tbnd(kk)+1))
6690 if (hk2 > tw2)
call hecmw_mf_kernel_gemv(hk2 - tw2, tw2, hk2, fct%sn(s0)%lval(okk2 + tw2), &
6691 vsg(fct%sn(s0)%tbnd(kk)+1), vsg(fct%sn(s0)%tbnd(kk)+tw2+1))
6692 xv(1:tw2) = vsg(fct%sn(s0)%tbnd(kk)+1:fct%sn(s0)%tbnd(kk)+tw2)
6694 ib = mf_pool_slot(pool)
6695 allocate(pool%box(ib)%rv(tw2))
6696 pool%box(ib)%rv(1:tw2) = xv(1:tw2)
6697 do dst = map%rbeg(s0), map%rbeg(s0) + map%rcnt(s0) - 1
6698 if (dst == me()) cycle
6699 if (.not. fwd_xdest(dst, kk)) cycle
6701 call hecmw_isend_r(pool%box(ib)%rv, tw2, dst, 16*s0+9, map%comm, req)
6702 call mf_pool_req(pool, req)
6704 if (.not. fct%lu)
call hecmw_mf_kernel_dsolve(tw2, hk2, fct%sn(s0)%lval(okk2), &
6705 fct%sn(s0)%ptype(fct%sn(s0)%tbnd(kk)+1:), fct%sn(s0)%dsub(fct%sn(s0)%tbnd(kk)+1:), &
6706 vsg(fct%sn(s0)%tbnd(kk)+1))
6708 call fwd_fin(s0, kk)
6709 else if (fwd_xdest(me(), kk))
then
6710 call hecmw_recv_r(xv, tw2, rdk, 16*s0+9, map%comm, stat)
6713 do i2 = kk + 1, ntc0
6714 if (fso_s(i2, kk) /= me()) cycle
6715 if (.not. held(i2))
then
6716 call hecmw_recv_r(vsg(fct%sn(s0)%tbnd(i2)+1:), fct%sn(s0)%tbnd(i2+1) - fct%sn(s0)%tbnd(i2), &
6717 fwd_ph(i2, kk), 16*s0+13, map%comm, stat)
6720 h2 = fct%sn(s0)%tbnd(i2+1) - fct%sn(s0)%tbnd(i2)
6721 oik2 = 1 + fct%sn(s0)%bptr(mf_bidx(nt0, kk, i2))
6722 call hecmw_mf_kernel_gemv(h2, tw2, h2, fct%sn(s0)%lval(oik2), xv, vsg(fct%sn(s0)%tbnd(i2)+1))
6725 if (kk < kend .and. kk + 1 < i2)
then
6726 nxt = fso_s(i2, kk+1)
6727 else if (i2 <= kend)
then
6730 if (nxt >= 0 .and. nxt /= me())
then
6731 ib = mf_pool_slot(pool)
6732 allocate(pool%box(ib)%rv(h2))
6733 pool%box(ib)%rv(1:h2) = vsg(fct%sn(s0)%tbnd(i2)+1:fct%sn(s0)%tbnd(i2+1))
6735 call hecmw_isend_r(pool%box(ib)%rv, h2, nxt, 16*s0+13, map%comm, req)
6736 call mf_pool_req(pool, req)
6743 if (towner(t2) /= me()) cycle
6744 h2 = ctb(t2) - ctb(t2-1)
6745 idx2 = mf_bidx(nt0, kk, ntc0 + t2)
6746 r2 = fct%sn(s0)%brank(idx2)
6747 oik2 = 1 + fct%sn(s0)%bptr(idx2)
6750 call hecmw_mf_kernel_gemv(h2, tw2, h2, fct%sn(s0)%lval(oik2), xv, vb(tprow(t2)+1))
6752 call hecmw_mf_kernel_mult_tv(tw2, r2, tw2, fct%sn(s0)%lval(oik2 + int(h2, 8)*r2), xv, tbuf)
6753 call hecmw_mf_kernel_gemv(h2, r2, h2, fct%sn(s0)%lval(oik2), tbuf, vb(tprow(t2)+1))
6759 do i2 = kend + 1, ntc0
6760 if (held(i2))
call fwd_fin(s0, i2)
6765 if (fwd_lh(i2) == me()) cycle
6766 call hecmw_recv_r(v(fct%sn(s0)%tbnd(i2)+1:), fct%sn(s0)%tbnd(i2+1) - fct%sn(s0)%tbnd(i2), &
6767 fwd_lh(i2), 16*s0+7, map%comm, stat)
6770 z(fct%sn(s0)%frow(d2)) = v(d2)
6776 if (ismaster) ndv = ndv + (ncol - npv)
6777 if (
allocated(fct%sn(s0)%dval))
deallocate(fct%sn(s0)%dval)
6779 allocate(fct%sn(s0)%dval(ndv))
6782 fct%sn(s0)%dval(1:ncol-npv) = v(npv+1:ncol)
6785 if (myrows > 0) fct%sn(s0)%dval(ofs+1:ofs+myrows) = vb(1:myrows)
6787 call mf_pool_wait(pool)
6788 do j = fct%cptr(s0), fct%cptr(s0+1) - 1
6790 if (
allocated(fct%sn(c)%dval))
deallocate(fct%sn(c)%dval)
6792 end subroutine fwd_dist
6795 function me()
result(r)
6796 integer(kind=kint) :: r
6803 function fwd_fh(i2)
result(r)
6804 integer(kind=kint),
intent(in) :: i2
6805 integer(kind=kint) :: r
6807 if (kend >= 1 .and. i2 > 1)
then
6809 else if (i2 <= kend)
then
6817 function fwd_ph(i2, kk)
result(r)
6818 integer(kind=kint),
intent(in) :: i2, kk
6819 integer(kind=kint) :: r
6829 function fwd_lh(i2)
result(r)
6830 integer(kind=kint),
intent(in) :: i2
6831 integer(kind=kint) :: r
6833 if (i2 <= kend)
then
6835 else if (kend >= 1)
then
6844 function fwd_xdest(r0, kk)
result(yes)
6845 integer(kind=kint),
intent(in) :: r0, kk
6847 integer(kind=kint) :: i0
6850 do i0 = kk + 1, ntc0
6851 if (fso_s(i0, kk) == r0) yes = .true.
6855 if (wrank(i0) == r0) yes = .true.
6857 end function fwd_xdest
6860 subroutine fwd_fin(s0, i2)
6861 integer(kind=kint),
intent(in) :: s0, i2
6862 integer(kind=kint) :: h2, ib
6864 h2 = fct%sn(s0)%tbnd(i2+1) - fct%sn(s0)%tbnd(i2)
6866 v(fct%sn(s0)%tbnd(i2)+1:fct%sn(s0)%tbnd(i2+1)) = vsg(fct%sn(s0)%tbnd(i2)+1:fct%sn(s0)%tbnd(i2+1))
6868 ib = mf_pool_slot(pool)
6869 allocate(pool%box(ib)%rv(h2))
6870 pool%box(ib)%rv(1:h2) = vsg(fct%sn(s0)%tbnd(i2)+1:fct%sn(s0)%tbnd(i2+1))
6872 call hecmw_isend_r(pool%box(ib)%rv, h2, master, 16*s0+7, map%comm, req)
6873 call mf_pool_req(pool, req)
6876 end subroutine fwd_fin
6884 subroutine bwd_dist(s0)
6885 integer(kind=kint),
intent(in) :: s0
6886 integer(kind=kint) :: d2, kk, tw2, hk2, t2, h2, i2, j2, iw2, ib, nr2, il, r0
6887 integer(kind=kint) :: rdk, nst, st, holder, dst
6888 integer(kind=8) :: okk2, okku2, oik2
6892 if (.not. (ismaster .or. hasband .or. myfst > 0))
return
6893 call mf_grow_r(v, int(max(ncol, 1), 8))
6894 call mf_grow_r(vb, int(max(myrows, 1), 8))
6895 call mf_grow_r(tok, int(max(ncol, 1), 8))
6896 call mf_grow_r(tbuf, int(max(ncol, nrow0 - ncol, 1), 8))
6897 cur5(1:map%rcnt(s0)) = ntc0
6901 if (map%owner(p) /= me())
call recv_xx(s0)
6904 if (wrank(iw2) == me()) cycle
6907 if (towner(t2) == wrank(iw2)) nr2 = nr2 + ctb(t2) - ctb(t2-1)
6910 ib = mf_pool_slot(pool)
6911 allocate(pool%box(ib)%rv(nr2))
6914 if (towner(t2) /= wrank(iw2)) cycle
6915 do d2 = ctb(t2-1) + 1, ctb(t2)
6917 pool%box(ib)%rv(il) = xx(rdof(d2))
6921 call hecmw_isend_r(pool%box(ib)%rv, nr2, wrank(iw2), 16*s0+10, map%comm, req)
6922 call mf_pool_req(pool, req)
6925 v(d2) = z(fct%sn(s0)%frow(d2))
6927 do d2 = npv + 1, ncol
6928 v(d2) = xx(fct%sn(s0)%fsdof(d2))
6933 h2 = fct%sn(s0)%tbnd(i2+1) - fct%sn(s0)%tbnd(i2)
6934 do dst = map%rbeg(s0), map%rbeg(s0) + map%rcnt(s0) - 1
6935 if (.not. bwd_ides(dst, i2)) cycle
6936 if (dst == me())
then
6937 vsg(fct%sn(s0)%tbnd(i2)+1:fct%sn(s0)%tbnd(i2+1)) = v(fct%sn(s0)%tbnd(i2)+1:fct%sn(s0)%tbnd(i2+1))
6940 ib = mf_pool_slot(pool)
6941 allocate(pool%box(ib)%rv(h2))
6942 pool%box(ib)%rv(1:h2) = v(fct%sn(s0)%tbnd(i2)+1:fct%sn(s0)%tbnd(i2+1))
6944 call hecmw_isend_r(pool%box(ib)%rv, h2, dst, 16*s0+6, map%comm, req)
6945 call mf_pool_req(pool, req)
6951 if (.not. bwd_ides(me(), i2)) cycle
6952 call hecmw_recv_r(vsg(fct%sn(s0)%tbnd(i2)+1:), fct%sn(s0)%tbnd(i2+1) - fct%sn(s0)%tbnd(i2), &
6953 master, 16*s0+6, map%comm, stat)
6960 if (towner(t2) /= me()) cycle
6961 do d2 = ctb(t2-1) + 1, ctb(t2)
6962 vb(tprow(t2) + (d2 - ctb(t2-1))) = xx(rdof(d2))
6966 call hecmw_recv_r(vb, myrows, master, 16*s0+10, map%comm, stat)
6971 tw2 = min(fct%sn(s0)%tbnd(kk+1), npv) - fct%sn(s0)%tbnd(kk)
6972 hk2 = fct%sn(s0)%tbnd(kk+1) - fct%sn(s0)%tbnd(kk)
6974 if (.not. (rdk == me() .or. hasband .or. fwd_xdest(me(), kk))) cycle
6978 do while (i2 <= ntc0)
6981 do while (j2 + 1 <= ntc0)
6982 if (fso_s(j2+1, kk) /= r0)
exit
6993 stw(nst) = wrank(iw2)
6998 if (rdk == me()) tok(1:tw2) = vsg(fct%sn(s0)%tbnd(kk)+1:fct%sn(s0)%tbnd(kk)+tw2)
7000 if (stw(st) == me())
then
7001 if (holder /= me())
call hecmw_recv_r(tok, tw2, holder, 16*s0+11, map%comm, stat)
7002 if (sta(st) > 0)
then
7003 do i2 = sta(st), stb(st)
7004 if (.not. held(i2))
call bwd_drain(s0, fso_s(i2, i2), i2)
7005 h2 = fct%sn(s0)%tbnd(i2+1) - fct%sn(s0)%tbnd(i2)
7007 oik2 = 1 + fct%sn(s0)%bptru(mf_bidx(nt0, kk, i2))
7008 call hecmw_mf_kernel_gemv_t(h2, tw2, h2, fct%sn(s0)%uval(oik2), &
7009 vsg(fct%sn(s0)%tbnd(i2)+1), tok)
7011 oik2 = 1 + fct%sn(s0)%bptr(mf_bidx(nt0, kk, i2))
7012 call hecmw_mf_kernel_gemv_t(h2, tw2, h2, fct%sn(s0)%lval(oik2), &
7013 vsg(fct%sn(s0)%tbnd(i2)+1), tok)
7017 call chain_apply(s0, kk, tw2)
7021 if (holder == me())
then
7022 ib = mf_pool_slot(pool)
7023 allocate(pool%box(ib)%rv(tw2))
7024 pool%box(ib)%rv(1:tw2) = tok(1:tw2)
7026 call hecmw_isend_r(pool%box(ib)%rv, tw2, stw(st), 16*s0+11, map%comm, req)
7027 call mf_pool_req(pool, req)
7032 if (rdk /= me())
then
7034 if (holder == me())
then
7035 ib = mf_pool_slot(pool)
7036 allocate(pool%box(ib)%rv(tw2))
7037 pool%box(ib)%rv(1:tw2) = tok(1:tw2)
7039 call hecmw_isend_r(pool%box(ib)%rv, tw2, rdk, 16*s0+11, map%comm, req)
7040 call mf_pool_req(pool, req)
7043 if (holder /= me())
call hecmw_recv_r(tok, tw2, holder, 16*s0+11, map%comm, stat)
7044 okk2 = 1 + fct%sn(s0)%bptr(mf_bidx(nt0, kk, kk))
7046 okku2 = 1 + fct%sn(s0)%bptru(mf_bidx(nt0, kk, kk))
7047 if (hk2 > tw2)
call hecmw_mf_kernel_gemv_t(hk2 - tw2, tw2, hk2, fct%sn(s0)%uval(okku2 + tw2), &
7048 vsg(fct%sn(s0)%tbnd(kk)+tw2+1), tok)
7049 call hecmw_mf_kernel_usolve(tw2, hk2, fct%sn(s0)%lval(okk2), hk2, fct%sn(s0)%uval(okku2), tok)
7051 if (hk2 > tw2)
call hecmw_mf_kernel_gemv_t(hk2 - tw2, tw2, hk2, fct%sn(s0)%lval(okk2 + tw2), &
7052 vsg(fct%sn(s0)%tbnd(kk)+tw2+1), tok)
7053 call hecmw_mf_kernel_trsv_t(tw2, hk2, fct%sn(s0)%lval(okk2), tok)
7055 vsg(fct%sn(s0)%tbnd(kk)+1:fct%sn(s0)%tbnd(kk)+tw2) = tok(1:tw2)
7058 do dst = map%rbeg(s0), map%rbeg(s0) + map%rcnt(s0) - 1
7059 if (dst == me()) cycle
7062 if (fso_s(kk, j2) == dst) sent = .true.
7064 if (.not. sent) cycle
7066 ib = mf_pool_slot(pool)
7067 allocate(pool%box(ib)%rv(hk2))
7068 pool%box(ib)%rv(1:hk2) = vsg(fct%sn(s0)%tbnd(kk)+1:fct%sn(s0)%tbnd(kk+1))
7071 call hecmw_isend_r(pool%box(ib)%rv, hk2, dst, 16*s0+5, map%comm, req)
7072 call mf_pool_req(pool, req)
7075 v(fct%sn(s0)%tbnd(kk)+1:fct%sn(s0)%tbnd(kk+1)) = vsg(fct%sn(s0)%tbnd(kk)+1:fct%sn(s0)%tbnd(kk+1))
7077 ib = mf_pool_slot(pool)
7078 allocate(pool%box(ib)%rv(hk2))
7079 pool%box(ib)%rv(1:hk2) = vsg(fct%sn(s0)%tbnd(kk)+1:fct%sn(s0)%tbnd(kk+1))
7081 call hecmw_isend_r(pool%box(ib)%rv, hk2, master, 16*s0+7, map%comm, req)
7082 call mf_pool_req(pool, req)
7088 if (fso_s(kk, kk) == me()) cycle
7089 call hecmw_recv_r(v(fct%sn(s0)%tbnd(kk)+1:), fct%sn(s0)%tbnd(kk+1) - fct%sn(s0)%tbnd(kk), &
7090 fso_s(kk, kk), 16*s0+7, map%comm, stat)
7093 xx(fct%sn(s0)%fsdof(d2)) = v(d2)
7095 do j = fct%cptr(s0), fct%cptr(s0+1) - 1
7097 if (map%owner(c) /= me())
call send_xx(c)
7100 call mf_pool_wait(pool)
7101 end subroutine bwd_dist
7105 function bwd_ides(r0, i2)
result(yes)
7106 integer(kind=kint),
intent(in) :: r0, i2
7108 integer(kind=kint) :: kk2
7110 if (i2 <= kend)
then
7111 yes = (fso_s(i2, i2) == r0)
7114 do kk2 = 1, min(kend, i2 - 1)
7115 if (fso_s(i2, kk2) == r0) yes = .true.
7118 end function bwd_ides
7122 subroutine bwd_drain(s0, src, ineed)
7123 integer(kind=kint),
intent(in) :: s0, src, ineed
7124 integer(kind=kint) :: slot, i0, kk2
7127 slot = src - map%rbeg(s0) + 1
7132 if (i0 <= kend .and. fso_s(i0, i0) == src .and. src /= me())
then
7134 if (fso_s(i0, kk2) == me()) mine = .true.
7140 call hecmw_recv_r(vsg(fct%sn(s0)%tbnd(i0)+1:), fct%sn(s0)%tbnd(i0+1) - fct%sn(s0)%tbnd(i0), &
7141 src, 16*s0+5, map%comm, stat)
7144 if (i0 == ineed)
exit
7146 end subroutine bwd_drain
7149 subroutine chain_apply(s0, kk, tw2)
7150 integer(kind=kint),
intent(in) :: s0, kk, tw2
7151 integer(kind=kint) :: t2, h2, r2, idx2
7152 integer(kind=8) :: oik2
7155 if (towner(t2) /= map%myrank) cycle
7156 h2 = ctb(t2) - ctb(t2-1)
7157 idx2 = mf_bidx(nt0, kk, ntc0 + t2)
7159 r2 = fct%sn(s0)%branku(idx2)
7160 oik2 = 1 + fct%sn(s0)%bptru(idx2)
7163 call hecmw_mf_kernel_gemv_t(h2, tw2, h2, fct%sn(s0)%uval(oik2), vb(tprow(t2)+1), tok)
7165 call hecmw_mf_kernel_mult_tv(h2, r2, h2, fct%sn(s0)%uval(oik2), vb(tprow(t2)+1), tbuf)
7166 call hecmw_mf_kernel_gemv(tw2, r2, tw2, fct%sn(s0)%uval(oik2 + int(h2, 8)*r2), tbuf, tok)
7169 r2 = fct%sn(s0)%brank(idx2)
7170 oik2 = 1 + fct%sn(s0)%bptr(idx2)
7173 call hecmw_mf_kernel_gemv_t(h2, tw2, h2, fct%sn(s0)%lval(oik2), vb(tprow(t2)+1), tok)
7175 call hecmw_mf_kernel_mult_tv(h2, r2, h2, fct%sn(s0)%lval(oik2), vb(tprow(t2)+1), tbuf)
7176 call hecmw_mf_kernel_gemv(tw2, r2, tw2, fct%sn(s0)%lval(oik2 + int(h2, 8)*r2), tbuf, tok)
7180 end subroutine chain_apply
7184 subroutine xx_dofs(c0, mm, list)
7185 integer(kind=kint),
intent(in) :: c0
7186 integer(kind=kint),
intent(out) :: mm
7187 integer(kind=kint),
intent(out),
allocatable :: list(:)
7188 integer(kind=kint) :: cnown, d, k, l, idx
7190 cnown = sym%sptr(c0+1) - sym%sptr(c0)
7191 mm = fct%sn(c0)%ncol - fct%sn(c0)%npiv
7192 do l = sym%rptr(c0) + cnown, sym%rptr(c0+1) - 1
7193 mm = mm + sym%ndof(sym%rlist(l))
7195 allocate(list(max(mm, 1)))
7197 do d = fct%sn(c0)%npiv + 1, fct%sn(c0)%ncol
7199 list(idx) = fct%sn(c0)%fsdof(d)
7201 do l = sym%rptr(c0) + cnown, sym%rptr(c0+1) - 1
7203 do d = fct%cdofptr(k), fct%cdofptr(k+1) - 1
7208 end subroutine xx_dofs
7212 subroutine send_xx(c0)
7213 integer(kind=kint),
intent(in) :: c0
7214 integer(kind=kint),
allocatable :: list(:)
7215 integer(kind=kint) :: mm, i0, ib
7217 call xx_dofs(c0, mm, list)
7219 ib = mf_pool_slot(gpool)
7220 allocate(gpool%box(ib)%rv(mm))
7222 gpool%box(ib)%rv(i0) = xx(list(i0))
7226 call hecmw_isend_r(gpool%box(ib)%rv, mm, map%owner(c0), 16*c0+12, map%comm, req)
7227 call mf_pool_req(gpool, req)
7228 end subroutine send_xx
7231 subroutine recv_xx(s0)
7232 integer(kind=kint),
intent(in) :: s0
7233 integer(kind=kint),
allocatable :: list(:)
7234 real(kind=kreal),
allocatable :: rv(:)
7235 integer(kind=kint) :: mm, i0
7237 call xx_dofs(s0, mm, list)
7240 call hecmw_recv_r(rv, mm, map%owner(sym%sparent(s0)), 16*s0+12, map%comm, stat)
7242 xx(list(i0)) = rv(i0)
7244 deallocate(rv, list)
7245 end subroutine recv_xx
7248 recursive subroutine bwd_task_mpi(s0)
7249 integer(kind=kint),
intent(in) :: s0
7250 integer(kind=kint) :: jj, cc
7252 call mf_super_bwd(sym, fct, fct%sn(s0), s0, z, xx)
7253 do jj = fct%cptr(s0), fct%cptr(s0+1) - 1
7256 call bwd_task_mpi(cc)
7259 end subroutine bwd_task_mpi
7265 subroutine mf_fwd_task_mpi(sym, map, fct, left, s, z, dmap)
7267 type(hecmwst_mf_symbolic),
intent(in) :: sym
7268 type(hecmwst_mf_map),
intent(in) :: map
7269 type(hecmwst_mf_factor),
intent(inout) :: fct
7270 integer(kind=kint),
intent(inout) :: left(:)
7271 integer(kind=kint),
intent(in) :: s
7272 real(kind=kreal),
intent(inout) :: z(:)
7273 integer(kind=kint),
intent(inout) :: dmap(:)
7274 integer(kind=kint) :: p, n, ss
7278 call mf_super_fwd(sym, fct, fct%sn(ss), ss, z, dmap)
7281 if (map%upper(p))
return
7285 left(p) = left(p) - 1
7291 end subroutine mf_fwd_task_mpi
7300 subroutine mf_super_fwd(sym, fct, sn, s, z, map)
7302 type(hecmwst_mf_symbolic),
intent(in) :: sym
7303 type(hecmwst_mf_factor),
intent(inout) :: fct
7304 type(mf_snode),
intent(inout) :: sn
7305 integer(kind=kint),
intent(in) :: s
7306 real(kind=kreal),
intent(inout) :: z(:)
7307 integer(kind=kint),
intent(inout) :: map(:)
7309 real(kind=kreal),
allocatable :: v(:), t(:)
7310 integer(kind=kint),
allocatable :: tw(:), rowoff(:)
7311 integer(kind=kint) :: nown, nrow_nodes, ncol0, ndel, d0, d1, d, c, i, j, k, l, a, r, rr
7312 integer(kind=kint) :: cnown, cndel, hk, hi
7313 integer(kind=8) :: okk, oik
7316 call mf_stored_grid(sn, g, tw)
7318 d0 = fct%cdofptr(sym%sptr(s))
7319 d1 = fct%cdofptr(sym%sptr(s+1)) - 1
7322 if (sn%frow(d) >= d0 .and. sn%frow(d) <= d1)
then
7323 v(d) = z(sn%frow(d))
7328 v(sn%ncol+1:g%nrow) = 0.0d0
7330 nown = sym%sptr(s+1) - sym%sptr(s)
7331 nrow_nodes = sym%rptr(s+1) - sym%rptr(s)
7333 ndel = sn%ncol - ncol0
7334 allocate(rowoff(nrow_nodes+1))
7336 do i = 1, nrow_nodes
7337 k = sym%rlist(sym%rptr(s)+i-1)
7338 rowoff(i+1) = rowoff(i) + sym%ndof(k)
7339 if (i == nown) rowoff(i+1) = rowoff(i+1) + ndel
7343 cnown = sym%sptr(c+1) - sym%sptr(c)
7344 cndel = fct%sn(c)%ncol - fct%sn(c)%npiv
7346 d = fct%sn(c)%frow(fct%sn(c)%npiv + j)
7347 v(map(d)) = v(map(d)) + fct%sn(c)%dval(j)
7350 do l = sym%rptr(c) + cnown, sym%rptr(c+1) - 1
7351 k = sym%cmap(sym%cmap_ptr(c) + (l - sym%rptr(c) - cnown))
7352 do a = 1, sym%ndof(sym%rlist(l))
7355 d = fct%cdofptr(sym%rlist(sym%rptr(s)+k-1)) + a - 1
7356 v(map(d)) = v(map(d)) + fct%sn(c)%dval(r)
7358 v(rowoff(k)+a) = v(rowoff(k)+a) + fct%sn(c)%dval(r)
7362 deallocate(fct%sn(c)%dval)
7367 if (fct%blr)
allocate(t(maxval(g%tb(1:g%nt) - g%tb(0:g%nt-1))))
7369 if (tw(k) == 0)
exit
7370 hk = g%tb(k) - g%tb(k-1)
7371 okk = 1 + g%coloff(k)
7372 if (fct%blr) okk = 1 + sn%bptr(mf_bidx(g%nt, k, k))
7373 call hecmw_mf_kernel_trsv(tw(k), hk, sn%lval(okk), v(g%tb(k-1)+1))
7374 if (hk > tw(k))
call hecmw_mf_kernel_gemv(hk - tw(k), tw(k), hk, sn%lval(okk + tw(k)), v(g%tb(k-1)+1), &
7375 v(g%tb(k-1)+tw(k)+1))
7377 hi = g%tb(i) - g%tb(i-1)
7380 rr = sn%brank(mf_bidx(g%nt, k, i))
7381 oik = 1 + sn%bptr(mf_bidx(g%nt, k, i))
7383 oik = 1 + g%coloff(k) + int(g%tb(i-1) - g%tb(k-1), 8)*tw(k)
7387 call hecmw_mf_kernel_gemv(hi, tw(k), hi, sn%lval(oik), v(g%tb(k-1)+1), v(g%tb(i-1)+1))
7389 call hecmw_mf_kernel_mult_tv(tw(k), rr, tw(k), sn%lval(oik + int(hi, 8)*rr), v(g%tb(k-1)+1), t)
7390 call hecmw_mf_kernel_gemv(hi, rr, hi, sn%lval(oik), t, v(g%tb(i-1)+1))
7393 if (.not. fct%lu)
call hecmw_mf_kernel_dsolve(tw(k), hk, sn%lval(okk), sn%ptype(g%tb(k-1)+1:), &
7394 sn%dsub(g%tb(k-1)+1:), v(g%tb(k-1)+1))
7396 if (
allocated(t))
deallocate(t)
7399 z(sn%frow(d)) = v(d)
7401 if (g%nrow > sn%npiv)
then
7402 allocate(sn%dval(g%nrow - sn%npiv))
7403 sn%dval(1:g%nrow-sn%npiv) = v(sn%npiv+1:g%nrow)
7406 end subroutine mf_super_fwd
7410 subroutine mf_super_bwd(sym, fct, sn, s, z, xx)
7412 type(hecmwst_mf_symbolic),
intent(in) :: sym
7413 type(hecmwst_mf_factor),
intent(in) :: fct
7414 type(mf_snode),
intent(in) :: sn
7415 integer(kind=kint),
intent(in) :: s
7416 real(kind=kreal),
intent(in) :: z(:)
7417 real(kind=kreal),
intent(inout) :: xx(:)
7419 real(kind=kreal),
allocatable :: v(:), t(:)
7420 integer(kind=kint),
allocatable :: tw(:)
7421 integer(kind=kint) :: nown, d, i, k, l, m, hk, hi, rr
7422 integer(kind=8) :: okk, okku, oik
7425 call mf_stored_grid(sn, g, tw)
7428 v(d) = z(sn%frow(d))
7430 do d = sn%npiv + 1, sn%ncol
7431 v(d) = xx(sn%fsdof(d))
7433 nown = sym%sptr(s+1) - sym%sptr(s)
7435 do l = sym%rptr(s) + nown, sym%rptr(s+1) - 1
7437 do d = fct%cdofptr(k), fct%cdofptr(k+1) - 1
7443 if (fct%blr)
allocate(t(maxval(g%tb(1:g%nt) - g%tb(0:g%nt-1))))
7445 if (tw(k) == 0) cycle
7446 hk = g%tb(k) - g%tb(k-1)
7447 okk = 1 + g%coloff(k)
7450 okk = 1 + sn%bptr(mf_bidx(g%nt, k, k))
7451 if (fct%lu) okku = 1 + sn%bptru(mf_bidx(g%nt, k, k))
7455 hi = g%tb(i) - g%tb(i-1)
7458 rr = sn%branku(mf_bidx(g%nt, k, i))
7459 oik = 1 + sn%bptru(mf_bidx(g%nt, k, i))
7461 oik = 1 + g%coloff(k) + int(g%tb(i-1) - g%tb(k-1), 8)*tw(k)
7465 call hecmw_mf_kernel_gemv_t(hi, tw(k), hi, sn%uval(oik), v(g%tb(i-1)+1), v(g%tb(k-1)+1))
7467 call hecmw_mf_kernel_mult_tv(hi, rr, hi, sn%uval(oik), v(g%tb(i-1)+1), t)
7468 call hecmw_mf_kernel_gemv(tw(k), rr, tw(k), sn%uval(oik + int(hi, 8)*rr), t, v(g%tb(k-1)+1))
7471 if (hk > tw(k))
call hecmw_mf_kernel_gemv_t(hk - tw(k), tw(k), hk, sn%uval(okku + tw(k)), &
7472 v(g%tb(k-1)+tw(k)+1), v(g%tb(k-1)+1))
7473 call hecmw_mf_kernel_usolve(tw(k), hk, sn%lval(okk), hk, sn%uval(okku), v(g%tb(k-1)+1))
7476 hi = g%tb(i) - g%tb(i-1)
7479 rr = sn%brank(mf_bidx(g%nt, k, i))
7480 oik = 1 + sn%bptr(mf_bidx(g%nt, k, i))
7482 oik = 1 + g%coloff(k) + int(g%tb(i-1) - g%tb(k-1), 8)*tw(k)
7486 call hecmw_mf_kernel_gemv_t(hi, tw(k), hi, sn%lval(oik), v(g%tb(i-1)+1), v(g%tb(k-1)+1))
7488 call hecmw_mf_kernel_mult_tv(hi, rr, hi, sn%lval(oik), v(g%tb(i-1)+1), t)
7489 call hecmw_mf_kernel_gemv(tw(k), rr, tw(k), sn%lval(oik + int(hi, 8)*rr), t, v(g%tb(k-1)+1))
7492 if (hk > tw(k))
call hecmw_mf_kernel_gemv_t(hk - tw(k), tw(k), hk, sn%lval(okk + tw(k)), &
7493 v(g%tb(k-1)+tw(k)+1), v(g%tb(k-1)+1))
7494 call hecmw_mf_kernel_trsv_t(tw(k), hk, sn%lval(okk), v(g%tb(k-1)+1))
7497 if (
allocated(t))
deallocate(t)
7500 xx(sn%fsdof(d)) = v(d)
7503 end subroutine mf_super_bwd
7508 type(hecmwst_mf_factor),
intent(in) :: fct
7509 integer(kind=kint),
intent(in) :: s
7510 integer(kind=8) :: w
7512 w = fct%sn(s)%pwords
7517 type(hecmwst_mf_factor),
intent(in) :: fct
7519 write(*,
'(a,i0,a,i0,a,i0)')
'[DIRECTmf]: tile size = ', fct%tile,
', tiles (max per front) = ', fct%max_tiles, &
7520 ', tile dim (max) = ', fct%max_tile_dim
7521 write(*,
'(a,i0,a,f10.3,a)')
'[DIRECTmf]: factor words (estimate) = ', fct%factor_words, &
7522 ' (', real(fct%factor_words, kind=kreal)*8.0d0/1024.0d0**3,
' GB)'
7523 write(*,
'(a,i0,a,i0,a,f10.3,a)')
'[DIRECTmf]: stack peak words = ', fct%stack_peak,
', front words = ', &
7524 fct%front_words,
' (', real(fct%stack_peak + fct%front_words, kind=kreal)*8.0d0/1024.0d0**3,
' GB)'
7525 write(*,
'(a,i0,a,i0,a,i0,a,f10.3,a)')
'[DIRECTmf]: LU mode estimate: factor words = ', 2*fct%factor_words, &
7526 ', stack peak words = ', 2*fct%stack_peak,
', front words = ', 2*fct%front_words,
' (', &
7527 real(2*(fct%factor_words + fct%stack_peak + fct%front_words), kind=kreal)*8.0d0/1024.0d0**3,
' GB)'
MPI distribution layer of the multifrontal solver: the supernode-to-rank mapping and the broadcast of...
Compressed (node level) graph of a block sparse matrix for the multifrontal solver.
Dense kernels of the multifrontal solver, symmetric (LDLt) and unsymmetric (LU) mode.
subroutine, public hecmw_mf_kernel_gemm(m, n, k, ldw, w, ldb, b, ldc, c)
c(m,n) <- c - w(m,k) * b(n,k)^T
logical function, public hecmw_mf_kernel_blr_available()
Whether the BLR compression is available (LAPACK build).
subroutine, public hecmw_mf_kernel_mult_nn(m, n, k, lda, a, ldb, b, ldc, c)
c(m,n) <- a(m,k) * b(k,n)
subroutine, public hecmw_mf_kernel_mult_tn(m, n, k, lda, a, ldb, b, ldc, c)
c(m,n) <- a(k,m)^T * b(k,n)
Numeric factorization (tiled multifrontal LDLt with threshold pivoting and delayed pivots) and the tr...
subroutine, public hecmw_mf_numeric_solve_mpi(sym, map, fct, b, x)
x = A^-1 b on the supernodes distributed by map; b and x are the replicated global vectors in the ori...
subroutine, public hecmw_mf_numeric_finalize(fct)
subroutine, public hecmw_mf_numeric_print(fct)
integer(kind=kint), parameter, public mf_blr_eta
BLR admissibility: default cluster distance threshold, a contribution tile within it of the pivot clu...
subroutine, public hecmw_mf_numeric_factor(hecMAT, sym, fct, ierr)
Factor the matrix of hecMAT in the order of sym, in LDLt mode (lower part referenced) or in LU mode w...
subroutine, public hecmw_mf_numeric_admis_bcast(fct, root)
Broadcast the admissibility cache from the rank that built it (the distributed driver builds the grap...
subroutine, public hecmw_mf_numeric_solve(sym, fct, b, x)
x = A^-1 b with the factor; b and x are in the original DOF numbering. The forward solve works on the...
subroutine, public hecmw_mf_numeric_init(graph, sym, tile, blr, eta, fct)
Estimates of the factor and stack sizes from the symbolic structure and the work space; tile is the t...
subroutine, public hecmw_mf_numeric_factor_mpi(hecMAT, sym, map, fct, ierr)
Factor the matrix with the supernodes distributed by map: every rank first factors its subtrees with ...
integer(kind=8) function, public hecmw_mf_numeric_front_words(fct, s)
Stored panel words of supernode s held on this rank (test support).
Symbolic factorization (elimination tree, column counts, supernodes, row structure) for the multifron...
integer(kind=kint) function hecmw_comm_get_size()
integer(kind=4), parameter kint
integer(kind=kint) function hecmw_comm_get_comm()
integer(kind=4), parameter kreal
integer(kind=kint), parameter hecmw_status_size
integer(kind=kint) function hecmw_comm_get_rank()
subroutine hecmw_bcast_i_comm(val, n, nbase, comm)
subroutine hecmw_waitall(cnt, reqs, stats)