16 logical,
save :: INITIALIZED = .false.
17 integer,
save :: SymType = 0
19 integer,
parameter :: DEBUG = 0
25 type (hecmwst_local_mesh),
intent(in) :: hecmesh
26 type (hecmwst_matrix ),
intent(inout) :: hecmat
28 logical,
intent(in) :: is_sym
48 type (hecmwst_local_mesh),
intent(in) :: hecmesh
49 type (hecmwst_matrix ),
intent(inout) :: hecmat
51 integer(kind=kint),
intent(out) :: istat
52 type (hecmwst_matrix),
intent(in),
optional :: conmat
53 integer :: method_org, precond_org
54 logical :: fg_eliminate
56 integer(kind=kint) :: num_lagrange_global
62 precond_org = hecmw_mat_get_precond(hecmat)
63 if (precond_org >= 30 .and. precond_org <= 32)
then
64 fg_eliminate = .false.
65 call hecmw_mat_set_precond(hecmat, precond_org - 20)
68 if (precond_org == 5) fg_amg = .true.
71 num_lagrange_global = fstrmat%num_lagrange
72 call hecmw_allreduce_i1(hecmesh, num_lagrange_global, hecmw_sum)
74 if (num_lagrange_global == 0)
then
75 if ((debug >= 1 .and.
myrank==0) .or. debug >= 2)
write(0,*)
'DEBUG: no contact'
77 method_org = hecmw_mat_get_method(hecmat)
78 call hecmw_mat_set_method(hecmat, 1)
84 call hecmw_mat_set_method(hecmat, method_org)
87 if ((debug >= 1 .and.
myrank==0) .or. debug >= 2)
write(0,*)
'DEBUG: with contact'
88 if (fg_eliminate)
then
90 call solve_eliminate(hecmesh, hecmat, fstrmat, conmat)
92 call solve_eliminate(hecmesh, hecmat, fstrmat)
96 write(0,*)
'ERROR: solve_no_eliminate not implemented for ParaCon'
97 call hecmw_abort( hecmw_comm_get_comm())
99 if (fstrmat%num_lagrange > 0)
then
100 call solve_no_eliminate(hecmesh, hecmat, fstrmat)
107 if (precond_org >= 30)
then
108 call hecmw_mat_set_precond(hecmat, precond_org)
111 istat = hecmw_mat_get_flag_diverged(hecmat)
118 subroutine solve_eliminate(hecMESH,hecMAT,fstrMAT,conMAT)
120 type (hecmwst_local_mesh),
intent(in),
target :: hecmesh
121 type (hecmwst_matrix ),
intent(inout) :: hecmat
123 type (hecmwst_matrix),
intent(in),
optional :: conmat
124 integer(kind=kint) :: ndof
125 integer(kind=kint),
allocatable :: iw2(:), iws(:)
126 real(kind=kreal),
allocatable :: wsl(:), wsu(:)
127 type(hecmwst_local_matrix),
target :: btmat
128 type(hecmwst_local_matrix) :: bttmat
129 type(hecmwst_matrix) :: hectkt
130 type(hecmwst_local_mesh),
pointer :: hecmeshtmp
131 type (hecmwst_local_matrix),
pointer :: bt_all
132 real(kind=kreal) :: t1
133 integer(kind=kint) :: method_org, i, ilag
134 logical :: fg_paracon
135 type (hecmwst_local_matrix),
pointer :: bkmat, bttkmat, bttktmat, bktmp
136 integer(kind=kint),
allocatable :: mark(:)
137 type(fstrst_contact_comm) :: concomm
138 integer(kind=kint) :: n_contact_dof, n_slave
139 integer(kind=kint),
allocatable :: contact_dofs(:), slaves(:)
140 real(kind=kreal),
allocatable :: btot(:)
143 if ((debug >= 1 .and.
myrank==0) .or. debug >= 2)
write(0,*)
'DEBUG: solve_eliminate start', hecmw_wtime()-t1
149 allocate(iw2(hecmat%NP*ndof))
151 allocate(iw2(hecmat%N*ndof))
153 if (debug >= 2)
write(0,*)
' DEBUG2[',
myrank,
']: num_lagrange',fstrmat%num_lagrange
154 allocate(iws(fstrmat%num_lagrange), wsl(fstrmat%num_lagrange), &
155 wsu(fstrmat%num_lagrange))
158 call choose_slaves(hecmat, fstrmat, iw2, iws, wsl, wsu, fg_paracon)
159 if (debug >= 2)
write(0,*)
' DEBUG2: slave DOFs chosen', hecmw_wtime()-t1
162 call make_btmat(hecmat, fstrmat, iw2, wsl, btmat, fg_paracon)
163 if (debug >= 2)
write(0,*)
' DEBUG2: make T done', hecmw_wtime()-t1
165 write(1000+
myrank,*)
'BTmat (local)'
166 call hecmw_localmat_write(btmat, 1000+
myrank)
169 call make_bttmat(hecmat, fstrmat, iw2, iws, wsu, bttmat, fg_paracon)
170 if (debug >= 2)
write(0,*)
' DEBUG2: make Tt done', hecmw_wtime()-t1
172 write(1000+
myrank,*)
'BTtmat (local)'
173 call hecmw_localmat_write(bttmat, 1000+
myrank)
178 call make_contact_dof_list(hecmat, fstrmat, n_contact_dof, contact_dofs)
182 call fstr_contact_comm_init(concomm, hecmesh, ndof, n_contact_dof, contact_dofs)
183 if (debug >= 2)
write(0,*)
' DEBUG2: make contact comm_table done', hecmw_wtime()-t1
187 call copy_mesh(hecmesh, hecmeshtmp, fg_paracon)
190 call hecmw_localmat_assemble(btmat, hecmesh, hecmeshtmp)
191 if (debug >= 2)
write(0,*)
' DEBUG2: assemble T done', hecmw_wtime()-t1
192 if (btmat%nc /= hecmesh%n_node)
then
193 if (debug >= 2)
write(0,*)
' DEBUG2[',
myrank,
']: node migrated with T',btmat%nc-hecmesh%n_node
196 write(1000+
myrank,*)
'BTmat (assembled)'
197 call hecmw_localmat_write(btmat, 1000+
myrank)
201 call hecmw_localmat_assemble(bttmat, hecmesh, hecmeshtmp)
202 if (debug >= 2)
write(0,*)
' DEBUG2: assemble Tt done', hecmw_wtime()-t1
203 if (bttmat%nc /= btmat%nc)
then
204 if (debug >= 2)
write(0,*)
' DEBUG2[',
myrank,
']: node migrated with BTtmat',bttmat%nc-btmat%nc
208 write(1000+
myrank,*)
'BTtmat (assembled)'
209 call hecmw_localmat_write(bttmat, 1000+
myrank)
213 allocate(mark(btmat%nr * btmat%ndof))
214 call mark_slave_dof(btmat, mark, n_slave, slaves)
215 call place_one_on_diag_of_unmarked_dof(btmat, mark)
216 call place_one_on_diag_of_unmarked_dof(bttmat, mark)
217 if (debug >= 2)
write(0,*)
' DEBUG2: place 1 on diag of T and Tt done', hecmw_wtime()-t1
219 write(1000+
myrank,*)
'BTmat (1s on non-slave diag)'
220 call hecmw_localmat_write(btmat, 1000+
myrank)
221 write(1000+
myrank,*)
'BTtmat (1s on non-slave diag)'
222 call hecmw_localmat_write(bttmat, 1000+
myrank)
227 call hecmw_localmat_init_with_hecmat(bkmat, conmat, fstrmat%num_lagrange)
229 write(1000+
myrank,*)
'BKmat (conMAT local)'
230 call hecmw_localmat_write(bkmat, 1000+
myrank)
234 call hecmw_localmat_assemble(bkmat, hecmesh, hecmeshtmp)
235 if (debug >= 2)
write(0,*)
' DEBUG2: assemble K (conMAT) done', hecmw_wtime()-t1
236 if (bkmat%nc /= bttmat%nc)
then
237 if (debug >= 2)
write(0,*)
' DEBUG2[',
myrank,
']: node migrated with BKmat',bkmat%nc-bttmat%nc
242 write(1000+
myrank,*)
'BKmat (conMAT assembled)'
243 call hecmw_localmat_write(bkmat, 1000+
myrank)
247 call hecmw_localmat_add_hecmat(bkmat, hecmat)
248 if (debug >= 2)
write(0,*)
' DEBUG2: add hecMAT to K done', hecmw_wtime()-t1
250 write(1000+
myrank,*)
'BKmat (hecMAT added)'
251 call hecmw_localmat_write(bkmat, 1000+
myrank)
256 call hecmw_localmat_multmat(bttmat, bkmat, hecmeshtmp, bttkmat)
257 if (debug >= 2)
write(0,*)
' DEBUG2: multiply Tt and K done', hecmw_wtime()-t1
259 write(1000+
myrank,*)
'BTtKmat'
260 call hecmw_localmat_write(bttkmat, 1000+
myrank)
265 call hecmw_localmat_multmat(bttkmat, btmat, hecmeshtmp, bttktmat)
266 if (debug >= 2)
write(0,*)
' DEBUG2: multiply TtK and T done', hecmw_wtime()-t1
268 write(1000+
myrank,*)
'BTtKTmat'
269 call hecmw_localmat_write(bttktmat, 1000+
myrank)
271 call hecmw_localmat_free(bttkmat)
277 call place_num_on_diag_of_marked_dof(bttktmat, 1.0d0, mark)
279 write(1000+
myrank,*)
'BTtKTmat (place 1.0 on slave diag)'
280 call hecmw_localmat_write(bttktmat, 1000+
myrank)
282 call hecmw_mat_init(hectkt)
283 call hecmw_localmat_make_hecmat(hecmat, bttktmat, hectkt)
284 if (debug >= 2)
write(0,*)
' DEBUG2: convert TtKT to hecTKT done', hecmw_wtime()-t1
285 call hecmw_localmat_free(bttktmat)
290 if (hecmesh%n_neighbor_pe > 0)
then
292 allocate(hecmeshtmp, bt_all)
293 call update_comm_table(hecmesh, btmat, hecmeshtmp, bt_all)
294 if (debug >= 2)
write(0,*)
' DEBUG2: Updating communication table done', hecmw_wtime()-t1
295 call hecmw_localmat_free(btmat)
298 hecmeshtmp => hecmesh
303 call hecmw_mat_init(hectkt)
304 call hecmw_trimatmul_ttkt(hecmeshtmp, bttmat, hecmat, bt_all, iws, fstrmat%num_lagrange, hectkt)
306 if ((debug >= 1 .and.
myrank==0) .or. debug >= 2)
write(0,*)
'DEBUG: calculated TtKT', hecmw_wtime()-t1
310 write(1000+
myrank,*)
'======================================================================='
311 write(1000+
myrank,*)
'RHS(original)----------------------------------------------------------'
312 write(1000+
myrank,*)
'size of hecMAT%B',
size(hecmat%B)
313 write(1000+
myrank,*)
'hecMAT%B: 1-',hecmat%N*ndof
314 write(1000+
myrank,*) hecmat%B(1:hecmat%N*ndof)
315 if (fstrmat%num_lagrange > 0)
then
316 write(1000+
myrank,*)
'hecMAT%B(lag):',hecmat%NP*ndof+1,
'-',hecmat%NP*ndof+fstrmat%num_lagrange
317 write(1000+
myrank,*) hecmat%B(hecmat%NP*ndof+1:hecmat%NP*ndof+fstrmat%num_lagrange)
319 if (n_slave > 0)
then
320 write(1000+
myrank,*)
'hecMAT%B(slave):',slaves(:)
321 write(1000+
myrank,*) hecmat%B(slaves(:))
328 write(1000+
myrank,*)
'RHS(conMAT)------------------------------------------------------------'
329 write(1000+
myrank,*)
'size of conMAT%B',
size(conmat%B)
330 write(1000+
myrank,*)
'conMAT%B: 1-',conmat%N*ndof
331 write(1000+
myrank,*) conmat%B(1:conmat%N*ndof)
332 write(1000+
myrank,*)
'conMAT%B(external): ',conmat%N*ndof+1,
'-',conmat%NP*ndof
333 write(1000+
myrank,*) conmat%B(conmat%N*ndof+1:conmat%NP*ndof)
334 if (fstrmat%num_lagrange > 0)
then
335 write(1000+
myrank,*)
'conMAT%B(lag):',conmat%NP*ndof+1,
'-',conmat%NP*ndof+fstrmat%num_lagrange
336 write(1000+
myrank,*) conmat%B(conmat%NP*ndof+1:conmat%NP*ndof+fstrmat%num_lagrange)
338 if (n_slave > 0)
then
339 write(1000+
myrank,*)
'conMAT%B(slave):',slaves(:)
340 write(1000+
myrank,*) conmat%B(slaves(:))
344 allocate(btot(hecmat%NP*ndof+fstrmat%num_lagrange))
345 call assemble_b_paracon(hecmesh, hecmat, conmat, fstrmat%num_lagrange, concomm, btot)
346 if (debug >= 2)
write(0,*)
' DEBUG2: assemble conMAT%B done', hecmw_wtime()-t1
348 write(1000+
myrank,*)
'RHS(conMAT assembled)--------------------------------------------------'
349 write(1000+
myrank,*)
'size of Btot',
size(btot)
350 write(1000+
myrank,*)
'Btot: 1-',conmat%N*ndof
351 write(1000+
myrank,*) btot(1:conmat%N*ndof)
352 if (fstrmat%num_lagrange > 0)
then
353 write(1000+
myrank,*)
'Btot(lag):',conmat%NP*ndof+1,
'-',conmat%NP*ndof+fstrmat%num_lagrange
354 write(1000+
myrank,*) btot(conmat%NP*ndof+1:conmat%NP*ndof+fstrmat%num_lagrange)
356 if (n_slave > 0)
then
357 write(1000+
myrank,*)
'Btot(slave):',slaves(:)
358 write(1000+
myrank,*) btot(slaves(:))
363 btot(i)=btot(i)+hecmat%B(i)
367 write(1000+
myrank,*)
'RHS(total)-------------------------------------------------------------'
368 write(1000+
myrank,*)
'size of Btot',
size(btot)
369 write(1000+
myrank,*)
'Btot: 1-',conmat%N*ndof
370 write(1000+
myrank,*) btot(1:conmat%N*ndof)
371 if (fstrmat%num_lagrange > 0)
then
372 write(1000+
myrank,*)
'Btot(lag):',conmat%NP*ndof+1,
'-',conmat%NP*ndof+fstrmat%num_lagrange
373 write(1000+
myrank,*) btot(conmat%NP*ndof+1:conmat%NP*ndof+fstrmat%num_lagrange)
375 if (n_slave > 0)
then
376 write(1000+
myrank,*)
'Btot(slave):',slaves(:)
377 write(1000+
myrank,*) btot(slaves(:))
382 allocate(hectkt%B(hectkt%NP*ndof + fstrmat%num_lagrange))
383 allocate(hectkt%X(hectkt%NP*ndof + fstrmat%num_lagrange))
385 do i=1, hectkt%N*ndof
386 hectkt%B(i) = btot(i)
389 do i=1, hectkt%N*ndof
390 hectkt%B(i) = hecmat%B(i)
393 do ilag=1, fstrmat%num_lagrange
394 hectkt%B(hectkt%NP*ndof + ilag) = hecmat%B(hecmat%NP*ndof + ilag)
396 do i=1, hectkt%N*ndof
397 hectkt%X(i) = hecmat%X(i)
399 do ilag=1,fstrmat%num_lagrange
400 hectkt%X(iws(ilag)) = 0.d0
405 call make_new_b_paracon(hecmesh, hecmat, conmat, btot, hecmeshtmp, hectkt, bttmat, bkmat, &
406 iws, wsl, fstrmat%num_lagrange, concomm, hectkt%B)
408 call make_new_b(hecmesh, hecmat, bttmat, iws, wsl, &
409 fstrmat%num_lagrange, hectkt%B)
411 if ((debug >= 1 .and.
myrank==0) .or. debug >= 2)
write(0,*)
'DEBUG: converted RHS', hecmw_wtime()-t1
413 write(1000+
myrank,*)
'RHS(converted)---------------------------------------------------------'
414 write(1000+
myrank,*)
'size of hecTKT%B',
size(hectkt%B)
415 write(1000+
myrank,*)
'hecTKT%B: 1-',hectkt%N*ndof
416 write(1000+
myrank,*) hectkt%B(1:hectkt%N*ndof)
424 if ((debug >= 1 .and.
myrank==0) .or. debug >= 2)
write(0,*)
'DEBUG: linear solver finished', hecmw_wtime()-t1
428 write(1000+
myrank,*)
'Solution(converted)----------------------------------------------------'
429 write(1000+
myrank,*)
'size of hecTKT%X',
size(hectkt%X)
430 write(1000+
myrank,*)
'hecTKT%X: 1-',hectkt%N*ndof
431 write(1000+
myrank,*) hectkt%X(1:hectkt%N*ndof)
434 hecmat%Iarray=hectkt%Iarray
435 hecmat%Rarray=hectkt%Rarray
439 call comp_x_slave_paracon(hecmesh, hecmat, btot, hecmeshtmp, hectkt, btmat, &
440 fstrmat%num_lagrange, iws, wsl, concomm, n_slave, slaves)
442 call hecmw_localmat_mulvec(bt_all, hectkt%X, hecmat%X)
443 call subst_blag(hecmat, iws, wsl, fstrmat%num_lagrange)
445 if ((debug >= 1 .and.
myrank==0) .or. debug >= 2)
write(0,*)
'DEBUG: recovered slave disp', hecmw_wtime()-t1
449 call comp_lag_paracon(hecmesh, hecmat, btot, hecmeshtmp, hectkt, bkmat, &
450 n_slave, slaves, fstrmat%num_lagrange, iws, wsu, concomm)
453 call comp_lag(hecmat, iws, wsu, fstrmat%num_lagrange)
455 if ((debug >= 1 .and.
myrank==0) .or. debug >= 2)
write(0,*)
'DEBUG: calculated lag', hecmw_wtime()-t1
462 write(1000+
myrank,*)
'Solution(original)-----------------------------------------------------'
463 write(1000+
myrank,*)
'size of hecMAT%X',
size(hecmat%X)
464 write(1000+
myrank,*)
'hecMAT%X: 1-',hecmat%N*ndof
465 write(1000+
myrank,*) hecmat%X(1:hecmat%N*ndof)
466 if (fstrmat%num_lagrange > 0)
then
467 write(1000+
myrank,*)
'hecMAT%X(lag):',hecmat%NP*ndof+1,
'-',hecmat%NP*ndof+fstrmat%num_lagrange
468 write(1000+
myrank,*) hecmat%X(hecmat%NP*ndof+1:hecmat%NP*ndof+fstrmat%num_lagrange)
470 if (n_slave > 0)
then
471 write(1000+
myrank,*)
'hecMAT%X(slave):',slaves(:)
472 write(1000+
myrank,*) hecmat%X(slaves(:))
477 if (fg_paracon .and. debug >= 2)
then
480 call check_solution(hecmesh, hecmat, fstrmat, btot, hecmeshtmp, hectkt, bkmat, &
481 concomm, n_slave, slaves)
485 call hecmw_localmat_free(bt_all)
486 call hecmw_localmat_free(bttmat)
487 call hecmw_mat_finalize(hectkt)
489 call hecmw_localmat_free(bkmat)
491 call fstr_contact_comm_finalize(concomm)
492 call free_mesh(hecmeshtmp, fg_paracon)
493 deallocate(hecmeshtmp)
495 if (hecmesh%n_neighbor_pe > 0)
then
496 call free_mesh(hecmeshtmp)
497 deallocate(hecmeshtmp)
502 if ((debug >= 1 .and.
myrank==0) .or. debug >= 2)
write(0,*)
'DEBUG: solve_eliminate end', hecmw_wtime()-t1
503 end subroutine solve_eliminate
505 subroutine choose_slaves(hecMAT, fstrMAT, iw2, iwS, wSL, wSU, fg_paracon)
507 type (hecmwst_matrix ),
intent(in) :: hecmat
509 integer,
intent(out) :: iw2(:), iws(:)
510 real(kind=kreal),
intent(out) :: wsl(:)
511 real(kind=kreal),
intent(out) :: wsu(:)
512 logical,
intent(in) :: fg_paracon
513 integer :: ndof, n, i, j, idof, jdof, l, ls, le, idx, imax
514 real(kind=kreal) :: val, vmax
515 integer,
allocatable :: iw1l(:), iw1u(:)
516 integer(kind=kint) :: n_slave_in,n_slave_out
520 if (fstrmat%num_lagrange == 0)
return
531 allocate(iw1l(n*ndof))
532 allocate(iw1u(n*ndof))
538 do i=1,fstrmat%num_lagrange
539 ls=fstrmat%indexL_lagrange(i-1)+1
540 le=fstrmat%indexL_lagrange(i)
542 j=fstrmat%itemL_lagrange(l)
545 iw1l(idx)=iw1l(idx)+1
551 ls=fstrmat%indexU_lagrange(i-1)+1
552 le=fstrmat%indexU_lagrange(i)
554 j=fstrmat%itemU_lagrange(l)
557 iw1u(idx)=iw1u(idx)+1
569 do i=1,fstrmat%num_lagrange
570 ls=fstrmat%indexL_lagrange(i-1)+1
571 le=fstrmat%indexL_lagrange(i)
575 j=fstrmat%itemL_lagrange(l)
578 val=fstrmat%AL_lagrange((l-1)*ndof+jdof)
579 if (iw1l(idx) == 1 .and. iw1u(idx) == 1 .and. abs(val) > abs(vmax))
then
585 if (imax == -1) stop
"ERROR: iterative solver for contact failed"
587 iws(i)=imax; wsl(i)=-1.d0/vmax
597 if (iw2(i) > 0) n_slave_in = n_slave_in + 1
600 do i=hecmat%N*ndof,n*ndof
601 if (iw2(i) > 0) n_slave_out = n_slave_out + 1
603 if (debug >= 2)
write(0,*)
' DEBUG2[',
myrank,
']: n_slave(in,out,tot)',n_slave_in,n_slave_out,n_slave_in+n_slave_out
605 deallocate(iw1l, iw1u)
607 call make_wsu(fstrmat, n, ndof, iw2, wsu)
608 end subroutine choose_slaves
610 subroutine make_wsu(fstrMAT, n, ndof, iw2, wSU)
613 integer(kind=kint),
intent(in) :: n, ndof
614 integer(kind=kint),
intent(in) :: iw2(:)
615 real(kind=kreal),
intent(out) :: wsu(:)
616 integer(kind=kint) :: i, idof, idx, js, je, j, k
618 if (fstrmat%num_lagrange == 0)
return
624 if (iw2(idx) > 0)
then
625 js=fstrmat%indexU_lagrange(i-1)+1
626 je=fstrmat%indexU_lagrange(i)
628 k=fstrmat%itemU_lagrange(j)
629 if (k==iw2(idx))
then
630 wsu(iw2(idx)) = -1.0/fstrmat%AU_lagrange((j-1)*ndof+idof)
637 end subroutine make_wsu
639 subroutine make_btmat(hecMAT, fstrMAT, iw2, wSL, BTmat, fg_paracon)
641 type (hecmwst_matrix ),
intent(inout) :: hecmat
643 integer,
intent(in) :: iw2(:)
644 real(kind=kreal),
intent(in) :: wsl(:)
645 type (hecmwst_local_matrix),
intent(out) :: btmat
646 logical,
intent(in) :: fg_paracon
647 type (hecmwst_local_matrix) :: tmat
648 integer :: ndof, n, i, nnz, l, js, je, j, k, jdof, kk, jj
649 real(kind=kreal) :: factor
660 tmat%nnz=fstrmat%numL_lagrange*ndof-fstrmat%num_lagrange
662 tmat%nnz=tmat%nr+fstrmat%numL_lagrange*ndof-2*fstrmat%num_lagrange
666 allocate(tmat%index(0:tmat%nr))
667 allocate(tmat%item(tmat%nnz), tmat%A(tmat%nnz))
672 nnz=ndof*(fstrmat%indexL_lagrange(iw2(i))-fstrmat%indexL_lagrange(iw2(i)-1))-1
680 tmat%index(i)=tmat%index(i-1)+nnz
682 if (tmat%nnz /= tmat%index(tmat%nr))
then
683 write(0,*) tmat%nnz, tmat%index(tmat%nr)
684 stop
'ERROR: Tmat%nnz wrong'
690 js=fstrmat%indexL_lagrange(iw2(i)-1)+1
691 je=fstrmat%indexL_lagrange(iw2(i))
694 k=fstrmat%itemL_lagrange(j)
700 tmat%A(l)=fstrmat%AL_lagrange(jj)*factor
705 if (.not. fg_paracon)
then
711 if (l /= tmat%index(i)+1)
then
712 write(0,*) l, tmat%index(i)+1
713 stop
'ERROR: Tmat%index wrong'
718 call hecmw_localmat_blocking(tmat, ndof, btmat)
719 call hecmw_localmat_free(tmat)
720 end subroutine make_btmat
722 subroutine make_bttmat(hecMAT, fstrMAT, iw2, iwS, wSU, BTtmat, fg_paracon)
724 type (hecmwst_matrix ),
intent(inout) :: hecmat
726 integer,
intent(in) :: iw2(:), iws(:)
727 real(kind=kreal),
intent(in) :: wsu(:)
728 type (hecmwst_local_matrix),
intent(out) :: bttmat
729 logical,
intent(in) :: fg_paracon
730 type (hecmwst_local_matrix) :: ttmat
731 integer :: ndof, n, i, nnz, l, js, je, j, k, idof, jdof, idx
742 ttmat%nnz=fstrmat%numU_lagrange*ndof-fstrmat%num_lagrange
744 ttmat%nnz=ttmat%nr+fstrmat%numU_lagrange*ndof-2*fstrmat%num_lagrange
748 allocate(ttmat%index(0:ttmat%nr))
749 allocate(ttmat%item(ttmat%nnz), ttmat%A(ttmat%nnz))
755 if (iw2(idx) <= 0)
then
756 if (fstrmat%num_lagrange == 0)
then
764 nnz=fstrmat%indexU_lagrange(i)-fstrmat%indexU_lagrange(i-1)
766 nnz=fstrmat%indexU_lagrange(i)-fstrmat%indexU_lagrange(i-1)+1
772 ttmat%index(idx)=ttmat%index(idx-1)+nnz
775 if (ttmat%nnz /= ttmat%index(ttmat%nr))
then
776 write(0,*) ttmat%nnz, ttmat%index(ttmat%nr)
777 stop
'ERROR: Ttmat%nnz wrong'
783 l=ttmat%index(idx-1)+1
784 if (iw2(idx) <= 0)
then
785 if (.not. fg_paracon)
then
791 if (fstrmat%num_lagrange > 0)
then
793 js=fstrmat%indexU_lagrange(i-1)+1
794 je=fstrmat%indexU_lagrange(i)
796 k=fstrmat%itemU_lagrange(j)
798 ttmat%A(l)=fstrmat%AU_lagrange((j-1)*ndof+idof)*wsu(k)
805 if (l /= ttmat%index(idx)+1)
then
806 write(0,*) l, ttmat%index(idx)+1
807 stop
'ERROR: Ttmat%index wrong'
813 call hecmw_localmat_blocking(ttmat, ndof, bttmat)
814 call hecmw_localmat_free(ttmat)
815 end subroutine make_bttmat
817 subroutine make_contact_dof_list(hecMAT, fstrMAT, n_contact_dof, contact_dofs)
819 type (hecmwst_matrix),
intent(in) :: hecmat
821 integer(kind=kint),
intent(out) :: n_contact_dof
822 integer(kind=kint),
allocatable,
intent(out) :: contact_dofs(:)
823 integer(kind=kint) :: ndof, icnt, i, ls, le, l, j, jdof, idx, k, idof
824 integer(kind=kint),
allocatable :: iw(:)
825 real(kind=kreal) :: val
826 if (fstrmat%num_lagrange == 0)
then
832 allocate(iw(hecmat%NP * ndof))
834 do i = 1, fstrmat%num_lagrange
835 ls = fstrmat%indexL_lagrange(i-1)+1
836 le = fstrmat%indexL_lagrange(i)
838 j = fstrmat%itemL_lagrange(l)
839 ljdof1:
do jdof = 1, ndof
840 val = fstrmat%AL_lagrange((l-1)*ndof+jdof)
842 if (val == 0.0d0) cycle
843 idx = (j-1)*ndof+jdof
845 if (iw(k) == idx) cycle ljdof1
855 ls = fstrmat%indexU_lagrange(i-1)+1
856 le = fstrmat%indexU_lagrange(i)
858 j = fstrmat%itemU_lagrange(l)
859 lidof1:
do idof = 1, ndof
860 val = fstrmat%AU_lagrange((l-1)*ndof+idof)
862 if (val == 0.0d0) cycle
863 idx = (i-1)*ndof+idof
865 if (iw(k) == idx) cycle lidof1
873 call quick_sort(iw, 1, icnt)
874 allocate(contact_dofs(icnt))
875 contact_dofs(1:icnt) = iw(1:icnt)
878 if (debug >= 2)
write(0,*)
' DEBUG2: contact_dofs',contact_dofs(:)
879 end subroutine make_contact_dof_list
881 subroutine make_new_b(hecMESH, hecMAT, BTtmat, iwS, wSL, num_lagrange, Bnew)
883 type(hecmwst_local_mesh),
intent(in) :: hecmesh
884 type(hecmwst_matrix),
intent(in) :: hecmat
885 type(hecmwst_local_matrix),
intent(in) :: bttmat
886 integer(kind=kint),
intent(in) :: iws(:)
887 real(kind=kreal),
intent(in) :: wsl(:)
888 integer(kind=kint),
intent(in) :: num_lagrange
889 real(kind=kreal),
intent(out) :: bnew(:)
890 real(kind=kreal),
allocatable :: btmp(:)
891 integer(kind=kint) :: npndof, nndof, i
893 npndof=hecmat%NP*hecmat%NDOF
894 nndof=hecmat%N*hecmat%NDOF
896 allocate(btmp(npndof))
903 bnew(iws(i))=wsl(i)*hecmat%B(npndof+i)
906 call hecmw_matvec(hecmesh, hecmat, bnew, btmp)
908 btmp(i)=hecmat%B(i)+btmp(i)
911 call hecmw_localmat_mulvec(bttmat, btmp, bnew)
914 end subroutine make_new_b
916 subroutine assemble_b_paracon(hecMESH, hecMAT, conMAT, num_lagrange, conCOMM, Btot)
918 type (hecmwst_local_mesh),
intent(in) :: hecmesh
919 type (hecmwst_matrix),
intent(in) :: hecmat, conmat
920 integer(kind=kint),
intent(in) :: num_lagrange
921 type (fstrst_contact_comm),
intent(in) :: concomm
922 real(kind=kreal),
intent(out) :: btot(:)
923 integer(kind=kint) :: ndof, nndof, npndof, i
925 nndof = hecmat%N * ndof
926 npndof = hecmat%NP * ndof
927 do i=1,npndof+num_lagrange
928 btot(i) = conmat%B(i)
930 call fstr_contact_comm_reduce_r(concomm, btot, hecmw_sum)
931 end subroutine assemble_b_paracon
933 subroutine make_new_b_paracon(hecMESH, hecMAT, conMAT, Btot, hecMESHtmp, hecTKT, BTtmat, BKmat, &
934 iwS, wSL, num_lagrange, conCOMM, Bnew)
936 type(hecmwst_local_mesh),
intent(in) :: hecmesh, hecmeshtmp
937 type(hecmwst_matrix),
intent(in) :: hecmat, conmat, hectkt
938 real(kind=kreal),
intent(in) :: btot(:)
939 type(hecmwst_local_matrix),
intent(in) :: bttmat, bkmat
940 integer(kind=kint),
intent(in) :: iws(:)
941 real(kind=kreal),
intent(in) :: wsl(:)
942 integer(kind=kint),
intent(in) :: num_lagrange
943 type (fstrst_contact_comm),
intent(in) :: concomm
944 real(kind=kreal),
intent(out) :: bnew(:)
945 real(kind=kreal),
allocatable :: btmp(:)
946 integer(kind=kint) :: npndof, nndof, npndof_new, i
951 npndof = hecmat%NP*hecmat%NDOF
952 nndof = hecmat%N *hecmat%NDOF
953 npndof_new = hectkt%NP*hectkt%NDOF
954 allocate(btmp(npndof_new))
962 bnew(iws(i))=wsl(i)*btot(npndof+i)
965 call fstr_contact_comm_reduce_r(concomm, bnew, hecmw_sum)
967 call hecmw_update_3_r(hecmeshtmp, bnew, hectkt%NP)
968 call hecmw_localmat_mulvec(bkmat, bnew, btmp)
970 btmp(i)=btot(i)+btmp(i)
973 call hecmw_update_3_r(hecmeshtmp, btmp, hectkt%NP)
974 call hecmw_localmat_mulvec(bttmat, btmp, bnew)
976 end subroutine make_new_b_paracon
978 subroutine subst_blag(hecMAT, iwS, wSL, num_lagrange)
980 type(hecmwst_matrix),
intent(inout) :: hecmat
981 integer(kind=kint),
intent(in) :: iws(:)
982 real(kind=kreal),
intent(in) :: wsl(:)
983 integer(kind=kint),
intent(in) :: num_lagrange
984 integer(kind=kint) :: npndof, i, ilag
986 npndof=hecmat%NP*hecmat%NDOF
989 hecmat%X(ilag)=hecmat%X(ilag)-wsl(i)*hecmat%B(npndof+i)
991 end subroutine subst_blag
993 subroutine comp_x_slave_paracon(hecMESH, hecMAT, Btot, hecMESHtmp, hecTKT, BTmat, &
994 num_lagrange, iwS, wSL, conCOMM, n_slave, slaves)
996 type (hecmwst_local_mesh),
intent(in) :: hecmesh, hecmeshtmp
997 type (hecmwst_matrix),
intent(inout) :: hecmat
998 type (hecmwst_matrix),
intent(in) :: hectkt
999 real(kind=kreal),
intent(in) :: btot(:)
1000 type (hecmwst_local_matrix),
intent(in) :: btmat
1001 integer(kind=kint),
intent(in) :: num_lagrange
1002 integer(kind=kint),
intent(in) :: iws(:)
1003 real(kind=kreal),
intent(in) :: wsl(:)
1004 type (fstrst_contact_comm),
intent(in) :: concomm
1005 integer(kind=kint),
intent(in) :: n_slave
1006 integer(kind=kint),
intent(in) :: slaves(:)
1007 integer(kind=kint) :: ndof, ndof2, npndof, nndof, ilag, islave, i
1008 real(kind=kreal),
allocatable :: xtmp(:)
1011 npndof = hecmat%NP * ndof
1012 nndof = hecmat%N * ndof
1017 call hecmw_update_3_r(hecmeshtmp, hectkt%X, hectkt%NP)
1018 call hecmw_localmat_mulvec(btmat, hectkt%X, hecmat%X)
1021 allocate(xtmp(npndof))
1023 do ilag = 1, num_lagrange
1026 xtmp(islave) = wsl(ilag) * btot(npndof + ilag)
1030 call fstr_contact_comm_reduce_r(concomm, xtmp, hecmw_sum)
1035 hecmat%X(islave) = hecmat%X(islave) - xtmp(islave)
1038 end subroutine comp_x_slave_paracon
1040 subroutine comp_lag(hecMAT, iwS, wSU, num_lagrange)
1042 type(hecmwst_matrix),
intent(inout) :: hecmat
1043 integer(kind=kint),
intent(in) :: iws(:)
1044 real(kind=kreal),
intent(in) :: wsu(:)
1045 integer(kind=kint),
intent(in) :: num_lagrange
1046 integer(kind=kint) :: ndof, ndof2, ilag, is, i, idof
1047 integer(kind=kint) :: js, je, j, jj, ij0, j0, jdof
1048 real(kind=kreal),
pointer :: xlag(:)
1053 xlag=>hecmat%X(hecmat%NP*ndof+1:hecmat%NP*ndof+num_lagrange)
1055 do ilag=1,num_lagrange
1058 idof=mod(is-1, ndof)+1
1059 xlag(ilag)=hecmat%B(is)
1061 js=hecmat%indexL(i-1)+1
1065 ij0=(j-1)*ndof2+(idof-1)*ndof
1068 xlag(ilag)=xlag(ilag)-hecmat%AL(ij0+jdof)*hecmat%X(j0+jdof)
1072 ij0=(i-1)*ndof2+(idof-1)*ndof
1075 xlag(ilag)=xlag(ilag)-hecmat%D(ij0+jdof)*hecmat%X(j0+jdof)
1078 js=hecmat%indexU(i-1)+1
1082 ij0=(j-1)*ndof2+(idof-1)*ndof
1085 xlag(ilag)=xlag(ilag)-hecmat%AU(ij0+jdof)*hecmat%X(j0+jdof)
1088 xlag(ilag)=-wsu(ilag)*xlag(ilag)
1090 end subroutine comp_lag
1092 subroutine comp_lag_paracon(hecMESH, hecMAT, Btot, hecMESHtmp, hecTKT, BKmat, &
1093 n_slave, slaves, num_lagrange, iwS, wSU, conCOMM)
1095 type (hecmwst_local_mesh),
intent(in) :: hecmesh, hecmeshtmp
1096 type (hecmwst_matrix),
intent(inout) :: hecmat, hectkt
1097 real(kind=kreal),
intent(in) :: btot(:)
1098 type (hecmwst_local_matrix),
intent(in) :: bkmat
1099 integer(kind=kint),
intent(in) :: n_slave, num_lagrange
1100 integer(kind=kint),
intent(in) :: slaves(:), iws(:)
1101 real(kind=kreal),
intent(in) :: wsu(:)
1102 type (fstrst_contact_comm),
intent(in) :: concomm
1103 integer(kind=kint) :: ndof, npndof, nndof, npndof_new, i, ilag, islave
1104 real(kind=kreal),
allocatable :: btmp(:)
1105 real(kind=kreal),
pointer :: xlag(:)
1106 integer(kind=kint),
allocatable :: slaves_ndup(:)
1108 npndof = hecmat%NP * ndof
1109 nndof = hecmat%N * ndof
1110 npndof_new = hectkt%NP * ndof
1115 hectkt%X(1:nndof) = hecmat%X(1:nndof)
1116 call hecmw_update_3_r(hecmeshtmp, hectkt%X, hectkt%NP)
1117 allocate(btmp(npndof))
1118 call hecmw_localmat_mulvec(bkmat, hectkt%X, btmp)
1126 btmp(islave) = btot(islave) - btmp(islave)
1132 call fstr_contact_comm_bcast_r(concomm, btmp)
1135 xlag => hecmat%X(npndof+1:npndof+num_lagrange)
1136 do ilag = 1, num_lagrange
1139 xlag(ilag)=-wsu(ilag)*btmp(islave)
1142 end subroutine comp_lag_paracon
1144 subroutine comp_lag_paracon2(hecMESH, hecMAT, conMAT, num_lagrange, iwS, wSU, conCOMM)
1146 type (hecmwst_local_mesh),
intent(in) :: hecmesh
1147 type (hecmwst_matrix),
intent(inout) :: hecmat
1148 type (hecmwst_matrix),
intent(in) :: conmat
1149 integer(kind=kint),
intent(in) :: num_lagrange
1150 integer(kind=kint),
intent(in) :: iws(:)
1151 real(kind=kreal),
intent(in) :: wsu(:)
1152 type (fstrst_contact_comm),
intent(in) :: concomm
1153 integer(kind=kint) :: ndof, ndof2, npndof, nndof, i, ilag, irow, idof, js, je, j, jcol, jdof, islave
1154 real(kind=kreal),
allocatable :: btmp(:)
1155 real(kind=kreal),
pointer :: xlag(:)
1158 npndof = hecmat%NP * ndof
1159 nndof = hecmat%N * ndof
1168 allocate(btmp(npndof))
1169 call hecmw_matresid(hecmesh, hecmat, hecmat%X, hecmat%B, btmp)
1170 call fstr_contact_comm_bcast_r(concomm, btmp)
1173 do ilag = 1, num_lagrange
1175 btmp(i) = btmp(i) + conmat%B(i)
1176 irow = (i + ndof - 1) / ndof
1177 idof = i - ndof*(irow-1)
1179 js = conmat%indexL(irow-1)+1
1180 je = conmat%indexL(irow)
1182 jcol = conmat%itemL(j)
1184 btmp(i) = btmp(i) - conmat%AL(ndof2*(j-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(jcol-1)+jdof)
1189 btmp(i) = btmp(i) - conmat%D(ndof2*(irow-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(irow-1)+jdof)
1192 js = conmat%indexU(irow-1)+1
1193 je = conmat%indexU(irow)
1195 jcol = conmat%itemU(j)
1197 btmp(i) = btmp(i) - conmat%AU(ndof2*(j-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(jcol-1)+jdof)
1203 xlag => hecmat%X(npndof+1:npndof+num_lagrange)
1204 do ilag = 1, num_lagrange
1207 xlag(ilag)=-wsu(ilag)*btmp(islave)
1210 end subroutine comp_lag_paracon2
1212 subroutine check_solution(hecMESH, hecMAT, fstrMAT, Btot, hecMESHtmp, hecTKT, BKmat, &
1213 conCOMM, n_slave, slaves)
1215 type (hecmwst_local_mesh),
intent(in) :: hecmesh, hecmeshtmp
1216 type (hecmwst_matrix),
intent(inout) :: hecmat, hectkt
1218 real(kind=kreal),
target,
intent(in) :: btot(:)
1219 type (hecmwst_local_matrix),
intent(in) :: bkmat
1220 type (fstrst_contact_comm),
intent(in) :: concomm
1221 integer(kind=kint),
intent(in) :: n_slave
1222 integer(kind=kint),
intent(in) :: slaves(:)
1223 integer(kind=kint) :: ndof, nndof, npndof, num_lagrange, i, ls, le, l, j, idof, jdof
1224 real(kind=kreal),
allocatable,
target :: r(:)
1225 real(kind=kreal),
allocatable :: btmp(:)
1226 real(kind=kreal),
pointer :: rlag(:), blag(:), xlag(:)
1227 real(kind=kreal) :: rnrm2, rlagnrm2
1228 real(kind=kreal) :: bnrm2, blagnrm2
1230 nndof = hecmat%N * ndof
1231 npndof = hecmat%NP * ndof
1232 num_lagrange = fstrmat%num_lagrange
1234 allocate(r(npndof + num_lagrange))
1236 allocate(btmp(npndof))
1238 rlag => r(npndof+1:npndof+num_lagrange)
1239 blag => btot(npndof+1:npndof+num_lagrange)
1240 xlag => hecmat%X(npndof+1:npndof+num_lagrange)
1246 hectkt%X(1:nndof) = hecmat%X(1:nndof)
1247 call hecmw_update_3_r(hecmeshtmp, hectkt%X, hectkt%NP)
1248 call hecmw_localmat_mulvec(bkmat, hectkt%X, btmp)
1249 r(1:nndof) = btot(1:nndof) - btmp(1:nndof)
1253 if (fstrmat%num_lagrange > 0)
then
1255 ls = fstrmat%indexU_lagrange(i-1)+1
1256 le = fstrmat%indexU_lagrange(i)
1258 j = fstrmat%itemU_lagrange(l)
1260 btmp(ndof*(i-1)+idof) = btmp(ndof*(i-1)+idof) + fstrmat%AU_lagrange(ndof*(l-1)+idof) * xlag(j)
1265 call fstr_contact_comm_reduce_r(concomm, btmp, hecmw_sum)
1266 r(1:nndof) = r(1:nndof) - btmp(1:nndof)
1269 call hecmw_update_3_r(hecmesh, hecmat%X, hecmat%NP)
1270 do i = 1, num_lagrange
1272 ls = fstrmat%indexL_lagrange(i-1)+1
1273 le = fstrmat%indexL_lagrange(i)
1275 j = fstrmat%itemL_lagrange(l)
1277 rlag(i) = rlag(i) - fstrmat%AL_lagrange(ndof*(l-1)+jdof) * hecmat%X(ndof*(j-1)+jdof)
1283 if (debug >= 3)
then
1284 write(1000+
myrank,*)
'Residual---------------------------------------------------------------'
1285 write(1000+
myrank,*)
'size of R',
size(r)
1286 write(1000+
myrank,*)
'R: 1-',hecmat%N*ndof
1287 write(1000+
myrank,*) r(1:hecmat%N*ndof)
1288 if (fstrmat%num_lagrange > 0)
then
1289 write(1000+
myrank,*)
'R(lag):',hecmat%NP*ndof+1,
'-',hecmat%NP*ndof+fstrmat%num_lagrange
1290 write(1000+
myrank,*) r(hecmat%NP*ndof+1:hecmat%NP*ndof+fstrmat%num_lagrange)
1292 if (n_slave > 0)
then
1293 write(1000+
myrank,*)
'R(slave):',slaves(:)
1294 write(1000+
myrank,*) r(slaves(:))
1298 call hecmw_innerproduct_r(hecmesh, ndof, r, r, rnrm2)
1299 call hecmw_innerproduct_r(hecmesh, ndof, btot, btot, bnrm2)
1301 do i = 1, num_lagrange
1302 rlagnrm2 = rlagnrm2 + rlag(i)*rlag(i)
1304 call hecmw_allreduce_r1(hecmesh, rlagnrm2, hecmw_sum)
1306 do i = 1, num_lagrange
1307 blagnrm2 = blagnrm2 + blag(i)*blag(i)
1309 call hecmw_allreduce_r1(hecmesh, blagnrm2, hecmw_sum)
1312 write(0,*)
'INFO: resid(x,lag,tot)',sqrt(rnrm2),sqrt(rlagnrm2),sqrt(rnrm2+rlagnrm2)
1313 write(0,*)
'INFO: rhs (x,lag,tot)',sqrt(bnrm2),sqrt(blagnrm2),sqrt(bnrm2+blagnrm2)
1315 end subroutine check_solution
1317 subroutine check_solution2(hecMESH, hecMAT, conMAT, fstrMAT, n_contact_dof, contact_dofs, &
1318 conCOMM, n_slave, slaves)
1320 type (hecmwst_local_mesh),
intent(in) :: hecmesh
1321 type (hecmwst_matrix),
intent(inout) :: hecmat
1322 type (hecmwst_matrix),
intent(in) :: conmat
1324 integer(kind=kint),
intent(in) :: n_contact_dof
1325 integer(kind=kint),
intent(in) :: contact_dofs(:)
1326 type (fstrst_contact_comm),
intent(in) :: concomm
1327 integer(kind=kint),
intent(in) :: n_slave
1328 integer(kind=kint),
intent(in) :: slaves(:)
1329 integer(kind=kint) :: ndof, ndof2, nndof, npndof, num_lagrange
1330 integer(kind=kint) :: icon, i, irow, idof, js, je, j, jcol, jdof, ls, le, l
1331 real(kind=kreal),
allocatable,
target :: r(:)
1332 real(kind=kreal),
allocatable :: r_con(:)
1333 real(kind=kreal),
pointer :: rlag(:), blag(:), xlag(:)
1334 real(kind=kreal) :: rnrm2, rlagnrm2
1337 nndof = hecmat%N * ndof
1338 npndof = hecmat%NP * ndof
1339 num_lagrange = fstrmat%num_lagrange
1341 allocate(r(npndof + num_lagrange))
1343 allocate(r_con(npndof))
1346 rlag => r(npndof+1:npndof+num_lagrange)
1347 blag => conmat%B(npndof+1:npndof+num_lagrange)
1348 xlag => hecmat%X(npndof+1:npndof+num_lagrange)
1360 call hecmw_matresid(hecmesh, hecmat, hecmat%X, hecmat%B, r)
1362 if (debug >= 3)
then
1363 write(1000+
myrank,*)
'Residual(original)-----------------------------------------------------'
1364 write(1000+
myrank,*)
'size of R',
size(r)
1365 write(1000+
myrank,*)
'R: 1-',hecmat%N*ndof
1366 write(1000+
myrank,*) r(1:hecmat%N*ndof)
1367 if (n_slave > 0)
then
1368 write(1000+
myrank,*)
'R(slave):',slaves(:)
1369 write(1000+
myrank,*) r(slaves(:))
1377 r_con(i) = conmat%B(i)
1379 do icon = 1, n_contact_dof
1380 i = contact_dofs(icon)
1381 irow = (i + ndof - 1) / ndof
1382 idof = i - ndof*(irow-1)
1384 js = conmat%indexL(irow-1)+1
1385 je = conmat%indexL(irow)
1387 jcol = conmat%itemL(j)
1389 r_con(i) = r_con(i) - conmat%AL(ndof2*(j-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(jcol-1)+jdof)
1394 r_con(i) = r_con(i) - conmat%D(ndof2*(irow-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(irow-1)+jdof)
1397 js = conmat%indexU(irow-1)+1
1398 je = conmat%indexU(irow)
1400 jcol = conmat%itemU(j)
1402 r_con(i) = r_con(i) - conmat%AU(ndof2*(j-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(jcol-1)+jdof)
1408 if (fstrmat%num_lagrange > 0)
then
1410 ls = fstrmat%indexU_lagrange(i-1)+1
1411 le = fstrmat%indexU_lagrange(i)
1413 j = fstrmat%itemU_lagrange(l)
1415 r_con(ndof*(i-1)+idof) = r_con(ndof*(i-1)+idof) - fstrmat%AU_lagrange(ndof*(l-1)+idof) * xlag(j)
1421 if (debug >= 3)
then
1422 write(1000+
myrank,*)
'Residual(contact,local)------------------------------------------------'
1423 write(1000+
myrank,*)
'size of R_con',
size(r_con)
1424 write(1000+
myrank,*)
'R_con: 1-',hecmat%N*ndof
1425 write(1000+
myrank,*) r_con(1:hecmat%N*ndof)
1426 write(1000+
myrank,*)
'R_con(external): ',hecmat%N*ndof+1,
'-',hecmat%NP*ndof
1427 write(1000+
myrank,*) r_con(hecmat%N*ndof+1:hecmat%NP*ndof)
1428 if (n_slave > 0)
then
1429 write(1000+
myrank,*)
'R_con(slave):',slaves(:)
1430 write(1000+
myrank,*) r_con(slaves(:))
1435 call fstr_contact_comm_reduce_r(concomm, r_con, hecmw_sum)
1437 if (debug >= 3)
then
1438 write(1000+
myrank,*)
'Residual(contact,assembled)--------------------------------------------'
1439 write(1000+
myrank,*)
'size of R_con',
size(r_con)
1440 write(1000+
myrank,*)
'R_con: 1-',hecmat%N*ndof
1441 write(1000+
myrank,*) r_con(1:hecmat%N*ndof)
1442 if (n_slave > 0)
then
1443 write(1000+
myrank,*)
'R_con(slave):',slaves(:)
1444 write(1000+
myrank,*) r_con(slaves(:))
1450 r(i) = r(i) + r_con(i)
1453 if (debug >= 3)
then
1454 write(1000+
myrank,*)
'Residual(total)--------------------------------------------------------'
1455 write(1000+
myrank,*)
'size of R',
size(r)
1456 write(1000+
myrank,*)
'R: 1-',hecmat%N*ndof
1457 write(1000+
myrank,*) r(1:hecmat%N*ndof)
1458 if (n_slave > 0)
then
1459 write(1000+
myrank,*)
'R(slave):',slaves(:)
1460 write(1000+
myrank,*) r(slaves(:))
1466 do i = 1, num_lagrange
1468 ls = fstrmat%indexL_lagrange(i-1)+1
1469 le = fstrmat%indexL_lagrange(i)
1471 j = fstrmat%itemL_lagrange(l)
1473 rlag(i) = rlag(i) - fstrmat%AL_lagrange(ndof*(l-1)+jdof) * hecmat%X(ndof*(j-1)+jdof)
1478 if (debug >= 3)
then
1479 write(1000+
myrank,*)
'Residual(lagrange)-----------------------------------------------------'
1480 if (fstrmat%num_lagrange > 0)
then
1481 write(1000+
myrank,*)
'R(lag):',hecmat%NP*ndof+1,
'-',hecmat%NP*ndof+fstrmat%num_lagrange
1482 write(1000+
myrank,*) r(hecmat%NP*ndof+1:hecmat%NP*ndof+fstrmat%num_lagrange)
1486 call hecmw_innerproduct_r(hecmesh, ndof, r, r, rnrm2)
1488 do i = 1, num_lagrange
1489 rlagnrm2 = rlagnrm2 + rlag(i)*rlag(i)
1491 call hecmw_allreduce_r1(hecmesh, rlagnrm2, hecmw_sum)
1493 if (
myrank == 0)
write(0,*)
'INFO: resid(x,lag,tot)',sqrt(rnrm2),sqrt(rlagnrm2),sqrt(rnrm2+rlagnrm2)
1494 end subroutine check_solution2
1496 subroutine mark_slave_dof(BTmat, mark, n_slave, slaves)
1498 type (hecmwst_local_matrix),
intent(in) :: btmat
1499 integer(kind=kint),
intent(out) :: mark(:)
1500 integer(kind=kint),
intent(out) :: n_slave
1501 integer(kind=kint),
allocatable,
intent(out) :: slaves(:)
1502 integer(kind=kint) :: ndof, ndof2, irow, js, je, j, jcol, idof, jdof, i
1506 do irow = 1, btmat%nr
1507 js = btmat%index(irow-1)+1
1508 je = btmat%index(irow)
1510 jcol = btmat%item(j)
1512 if (mark(ndof*(irow-1)+idof) == 1) cycle
1514 if (irow == jcol .and. idof == jdof) cycle
1515 if (btmat%A(ndof2*(j-1)+ndof*(idof-1)+jdof) /= 0.0d0)
then
1516 mark(ndof*(irow-1)+idof) = 1
1524 do i = 1, btmat%nr * ndof
1525 if (mark(i) /= 0) n_slave = n_slave + 1
1527 if (debug >= 2)
write(0,*)
' DEBUG2: n_slave',n_slave
1528 allocate(slaves(n_slave))
1530 do i = 1, btmat%nr * ndof
1531 if (mark(i) /= 0)
then
1532 n_slave = n_slave + 1
1536 if (debug >= 2)
write(0,*)
' DEBUG2: slaves',slaves(:)
1537 end subroutine mark_slave_dof
1539 subroutine place_one_on_diag_of_unmarked_dof(BTmat, mark)
1541 type (hecmwst_local_matrix),
intent(inout) :: btmat
1542 integer(kind=kint),
intent(in) :: mark(:)
1543 type (hecmwst_local_matrix) :: imat, wmat
1544 integer(kind=kint) :: ndof, ndof2, i, irow, js, je, j, jcol, idof, jdof, n_slave, n_other
1552 allocate(imat%index(0:imat%nr))
1553 allocate(imat%item(imat%nnz))
1559 allocate(imat%A(ndof2 * imat%nnz))
1563 do irow = 1, imat%nr
1565 if (mark(ndof*(irow-1)+idof) == 0)
then
1566 imat%A(ndof2*(irow-1)+ndof*(idof-1)+idof) = 1.0d0
1567 n_other = n_other + 1
1569 n_slave = n_slave + 1
1573 if (debug >= 2)
write(0,*)
' DEBUG2: n_slave,n_other',n_slave,n_other
1574 call hecmw_localmat_add(btmat, imat, wmat)
1575 call hecmw_localmat_free(btmat)
1576 call hecmw_localmat_free(imat)
1579 btmat%nnz = wmat%nnz
1580 btmat%ndof = wmat%ndof
1581 btmat%index => wmat%index
1582 btmat%item => wmat%item
1584 end subroutine place_one_on_diag_of_unmarked_dof
1586 subroutine place_num_on_diag_of_marked_dof(BTtKTmat, num, mark)
1588 type (hecmwst_local_matrix),
intent(inout) :: bttktmat
1589 real(kind=kreal),
intent(in) :: num
1590 integer(kind=kint),
intent(in) :: mark(:)
1591 integer(kind=kint) :: ndof, ndof2, irow, js, je, j, jcol, idof, jdof
1592 integer(kind=kint) :: n_error = 0
1593 ndof = bttktmat%ndof
1595 do irow = 1, bttktmat%nr
1596 js = bttktmat%index(irow-1)+1
1597 je = bttktmat%index(irow)
1599 jcol = bttktmat%item(j)
1600 if (irow /= jcol) cycle
1602 if (mark(ndof*(irow-1)+idof) == 1)
then
1603 if (bttktmat%A(ndof2*(j-1)+ndof*(idof-1)+idof) /= 0.0d0) &
1604 stop
'ERROR: nonzero diag on slave dof of BTtKTmat'
1605 bttktmat%A(ndof2*(j-1)+ndof*(idof-1)+idof) = num
1607 if (bttktmat%A(ndof2*(j-1)+ndof*(idof-1)+idof) == 0.0d0)
then
1615 if (n_error > 0)
then
1616 write(0,*)
'n_error',n_error
1617 stop
'ERROR: zero diag on non-slave dof of BTtKTmat'
1619 end subroutine place_num_on_diag_of_marked_dof
1621 subroutine update_comm_table(hecMESH, BTmat, hecMESHtmp, BT_all)
1623 type (hecmwst_local_mesh),
intent(in),
target :: hecmesh
1624 type(hecmwst_local_matrix),
intent(in) :: btmat
1625 type(hecmwst_local_mesh),
intent(inout),
target :: hecmeshtmp
1626 type (hecmwst_local_matrix),
intent(out) :: bt_all
1627 type(hecmwst_local_matrix),
allocatable :: bt_exp(:)
1628 integer(kind=kint) :: n_send, idom, irank, n_curexp, n_oldexp, n_orgexp
1629 integer(kind=kint) :: idx_0, idx_n, k, knod, n_newexp, j, jnod
1630 integer(kind=kint),
pointer :: cur_export(:), org_export(:)
1631 integer(kind=kint),
pointer :: old_export(:)
1632 integer(kind=kint),
allocatable,
target :: old_export_item(:)
1633 integer(kind=kint),
allocatable :: new_export(:)
1634 integer(kind=kint) :: sendbuf(2), recvbuf(2)
1635 integer(kind=kint) :: n_oldimp, n_newimp, n_orgimp, i0, n_curimp
1636 integer(kind=kint),
allocatable :: old_import(:)
1637 integer(kind=kint),
pointer :: org_import(:), cur_import(:)
1638 integer(kind=kint) :: tag
1639 type (hecmwst_local_matrix) :: bt_imp
1640 integer(kind=kint) :: nnz
1641 integer(kind=kint),
allocatable :: nnz_imp(:)
1642 integer(kind=kint),
allocatable :: index_imp(:), item_imp(:)
1643 real(kind=kreal),
allocatable :: val_imp(:)
1644 integer(kind=kint),
allocatable :: requests(:)
1645 integer(kind=kint),
allocatable :: statuses(:,:)
1646 integer(kind=kint) :: nr_imp, jj, ndof2, idx_0_tmp, idx_n_tmp
1647 integer(kind=kint) :: cnt, ks, ke, iimp, i, ii, ierror
1649 call copy_mesh(hecmesh, hecmeshtmp)
1650 allocate(bt_exp(hecmesh%n_neighbor_pe))
1651 call extract_bt_exp(btmat, hecmesh, bt_exp)
1654 allocate(statuses(hecmw_status_size,2*hecmesh%n_neighbor_pe))
1655 allocate(requests(2*hecmesh%n_neighbor_pe))
1657 allocate(old_export_item(hecmesh%export_index(hecmesh%n_neighbor_pe)))
1660 do idom = 1,hecmesh%n_neighbor_pe
1661 irank = hecmesh%neighbor_pe(idom)
1662 allocate(cur_export(bt_exp(idom)%nnz))
1663 call extract_cols(bt_exp(idom), cur_export, n_curexp)
1664 if (debug >= 1)
write(0,*)
'DEBUG: extract_cols done'
1666 idx_0 = hecmesh%export_index(idom-1)
1667 idx_n = hecmesh%export_index(idom)
1668 n_orgexp = idx_n - idx_0
1669 org_export => hecmesh%export_item(idx_0+1:idx_n)
1671 old_export => old_export_item(idx_0+1:idx_n)
1673 knod = org_export(k)
1674 if (.not. is_included(cur_export, n_curexp, knod))
then
1675 n_oldexp = n_oldexp + 1
1676 old_export(n_oldexp) = k
1679 if (debug >= 1)
write(0,*)
'DEBUG: making old_export done'
1681 call reorder_current_export(cur_export, n_curexp, org_export, n_orgexp, n_newexp, hecmesh%nn_internal)
1682 if (debug >= 1)
write(0,*)
'DEBUG: reorder_current_export done'
1684 if (n_curexp /= n_orgexp - n_oldexp + n_newexp) &
1685 stop
'ERROR: unknown error(num of export nodes)'
1687 call convert_bt_exp_col_id(bt_exp(idom), cur_export, n_curexp)
1688 if (debug >= 1)
write(0,*)
'DEBUG: convert_BT_expx_col_id done'
1690 call append_commtable(hecmeshtmp%n_neighbor_pe, hecmeshtmp%export_index, &
1691 hecmeshtmp%export_item, idom, cur_export, n_curexp)
1692 if (debug >= 1)
write(0,*)
'DEBUG: append_commtable (export) done'
1693 deallocate(cur_export)
1694 cur_export => hecmeshtmp%export_item(hecmeshtmp%export_index(idom-1)+1:hecmeshtmp%export_index(idom))
1696 sendbuf(1) = n_oldexp
1697 sendbuf(2) = n_newexp
1699 call hecmw_isend_int(sendbuf, 2, irank, tag, &
1700 hecmesh%MPI_COMM, requests(idom))
1701 if (n_oldexp > 0)
then
1704 call hecmw_isend_int(old_export, n_oldexp, irank, tag, &
1705 hecmesh%MPI_COMM, requests(hecmesh%n_neighbor_pe+n_send))
1708 if (debug >= 1)
write(0,*)
'DEBUG: isend n_oldexp, n_newexp, old_export done'
1709 do idom = 1,hecmesh%n_neighbor_pe
1710 irank = hecmesh%neighbor_pe(idom)
1713 call hecmw_recv_int(recvbuf, 2, irank, tag, &
1714 hecmesh%MPI_COMM, statuses(:,1))
1715 n_oldimp = recvbuf(1)
1716 n_newimp = recvbuf(2)
1717 if (n_oldimp > 0)
then
1718 allocate(old_import(n_oldimp))
1720 call hecmw_recv_int(old_import, n_oldimp, irank, tag, &
1721 hecmesh%MPI_COMM, statuses(:,1))
1724 idx_0 = hecmesh%import_index(idom-1)
1725 idx_n = hecmesh%import_index(idom)
1726 n_orgimp = idx_n - idx_0
1727 org_import => hecmesh%import_item(idx_0+1:idx_n)
1728 call append_nodes(hecmeshtmp, n_newimp, i0)
1729 if (debug >= 1)
write(0,*)
'DEBUG: append_nodes done'
1730 n_curimp = n_orgimp - n_oldimp + n_newimp
1731 allocate(cur_import(n_curimp))
1732 call make_cur_import(org_import, n_orgimp, old_import, n_oldimp, &
1733 n_newimp, i0, cur_import)
1734 if (n_oldimp > 0)
deallocate(old_import)
1735 if (debug >= 1)
write(0,*)
'DEBUG: make_cur_import done'
1736 call append_commtable(hecmeshtmp%n_neighbor_pe, hecmeshtmp%import_index, &
1737 hecmeshtmp%import_item, idom, cur_import, n_curimp)
1738 if (debug >= 1)
write(0,*)
'DEBUG: append_commtable (import) done'
1739 deallocate(cur_import)
1742 if (debug >= 1)
write(0,*)
'DEBUG: recv n_oldimp, n_newimp, old_import done'
1743 call hecmw_waitall(hecmesh%n_neighbor_pe + n_send, requests, statuses)
1744 deallocate(old_export_item)
1747 do idom = 1,hecmesh%n_neighbor_pe
1748 irank = hecmesh%neighbor_pe(idom)
1749 sendbuf(1) = bt_exp(idom)%nr
1750 sendbuf(2) = bt_exp(idom)%nnz
1752 call hecmw_isend_int(sendbuf, 2, irank, tag, &
1753 hecmesh%MPI_COMM, requests(2*idom-1))
1755 call hecmw_isend_int(bt_exp(idom)%index(0:bt_exp(idom)%nr), bt_exp(idom)%nr+1, &
1756 irank, tag, hecmesh%MPI_COMM, requests(2*idom))
1758 if (debug >= 1)
write(0,*)
'DEBUG: isend BT_exp (nnz and index) done'
1760 bt_imp%nc = hecmeshtmp%n_node - hecmeshtmp%nn_internal
1762 allocate(bt_imp%index(0:hecmesh%import_index(hecmesh%n_neighbor_pe)))
1764 allocate(nnz_imp(hecmesh%n_neighbor_pe))
1765 do idom = 1,hecmesh%n_neighbor_pe
1766 irank = hecmesh%neighbor_pe(idom)
1768 call hecmw_recv_int(recvbuf, 2, irank, tag, &
1769 hecmesh%MPI_COMM, statuses(:,1))
1771 nnz_imp(idom) = recvbuf(2)
1772 idx_0 = hecmesh%import_index(idom-1)
1773 idx_n = hecmesh%import_index(idom)
1774 if (nr_imp /= idx_n - idx_0) &
1775 stop
'ERROR: num of rows of BT_imp incorrect'
1776 bt_imp%nr = bt_imp%nr + nr_imp
1777 bt_imp%nnz = bt_imp%nnz + nnz_imp(idom)
1778 allocate(index_imp(0:nr_imp))
1780 call hecmw_recv_int(index_imp(0), nr_imp+1, irank, tag, &
1781 hecmesh%MPI_COMM, statuses(:,1))
1782 if (index_imp(nr_imp) /= nnz_imp(idom))
then
1783 if (debug >= 1)
write(0,*)
'ERROR: num of nonzero of BT_imp incorrect'
1784 if (debug >= 1)
write(0,*)
'nr_imp, index_imp(nr_imp), nnz_imp', &
1785 nr_imp, index_imp(nr_imp), nnz_imp(idom)
1789 jj = hecmesh%import_item(idx_0+j) - hecmesh%nn_internal
1790 bt_imp%index(jj) = index_imp(j) - index_imp(j-1)
1792 deallocate(index_imp)
1794 if (debug >= 1)
write(0,*)
'DEBUG: recv BT_imp (nnz and index) done'
1795 do j = 1, hecmesh%import_index(hecmesh%n_neighbor_pe)
1796 bt_imp%index(j) = bt_imp%index(j-1) + bt_imp%index(j)
1798 if (bt_imp%index(hecmesh%import_index(hecmesh%n_neighbor_pe)) /= bt_imp%nnz) &
1799 stop
'ERROR: total num of nonzero of BT_imp incorrect'
1800 ndof2 = btmat%ndof ** 2
1801 allocate(bt_imp%item(bt_imp%nnz),bt_imp%A(bt_imp%nnz * ndof2))
1802 call hecmw_waitall(hecmesh%n_neighbor_pe * 2, requests, statuses)
1805 do idom = 1,hecmesh%n_neighbor_pe
1806 irank = hecmesh%neighbor_pe(idom)
1808 call hecmw_isend_int(bt_exp(idom)%item, bt_exp(idom)%nnz, &
1809 irank, tag, hecmesh%MPI_COMM, requests(2*idom-1))
1811 call hecmw_isend_r(bt_exp(idom)%A, bt_exp(idom)%nnz * ndof2, &
1812 irank, tag, hecmesh%MPI_COMM, requests(2*idom))
1814 if (debug >= 1)
write(0,*)
'DEBUG: isend BT_exp (item and val) done'
1815 do idom = 1,hecmesh%n_neighbor_pe
1816 irank = hecmesh%neighbor_pe(idom)
1817 idx_0 = hecmesh%import_index(idom-1)
1818 idx_n = hecmesh%import_index(idom)
1819 allocate(item_imp(nnz_imp(idom)))
1821 call hecmw_recv_int(item_imp, nnz_imp(idom), &
1822 irank, tag, hecmesh%MPI_COMM, statuses(:,1))
1823 allocate(val_imp(nnz_imp(idom) * ndof2))
1825 call hecmw_recv_r(val_imp, nnz_imp(idom) * ndof2, &
1826 irank, tag, hecmesh%MPI_COMM, statuses(:,1))
1829 idx_0_tmp = hecmeshtmp%import_index(idom-1)
1830 idx_n_tmp = hecmeshtmp%import_index(idom)
1831 cur_import => hecmeshtmp%import_item(idx_0_tmp+1:idx_n_tmp)
1832 n_curimp = idx_n_tmp - idx_0_tmp
1833 n_orgimp = idx_n - idx_0
1836 jj = hecmesh%import_item(idx_0+j) - hecmesh%nn_internal
1837 ks = bt_imp%index(jj-1)
1838 ke = bt_imp%index(jj)
1841 iimp = item_imp(cnt)
1842 if (iimp <= 0 .or. n_curimp < iimp) &
1843 stop
'ERROR: received column id out of range'
1844 bt_imp%item(k) = cur_import(iimp)
1845 bt_imp%A((k-1)*ndof2+1:k*ndof2) = val_imp((cnt-1)*ndof2+1:cnt*ndof2)
1848 deallocate(item_imp, val_imp)
1851 if (debug >= 1)
write(0,*)
'DEBUG: recv BT_imp (item and val) done'
1852 call hecmw_waitall(hecmesh%n_neighbor_pe * 2, requests, statuses)
1854 deallocate(statuses)
1855 deallocate(requests)
1858 bt_all%nr = btmat%nr + bt_imp%nr
1859 bt_all%nc = btmat%nc + bt_imp%nc
1860 bt_all%nnz = btmat%nnz + bt_imp%nnz
1861 bt_all%ndof = btmat%ndof
1862 allocate(bt_all%index(0:bt_all%nr))
1863 allocate(bt_all%item(bt_all%nnz))
1864 allocate(bt_all%A(bt_all%nnz * ndof2))
1867 bt_all%index(i) = btmat%index(i)
1870 bt_all%index(btmat%nr+i) = bt_all%index(btmat%nr+i-1) + &
1871 bt_imp%index(i) - bt_imp%index(i-1)
1874 bt_all%item(i) = btmat%item(i)
1875 bt_all%A((i-1)*ndof2+1:i*ndof2) = btmat%A((i-1)*ndof2+1:i*ndof2)
1877 do i = 1, bt_imp%nnz
1879 bt_all%item(ii) = bt_imp%item(i)
1880 bt_all%A((ii-1)*ndof2+1:ii*ndof2) = bt_imp%A((i-1)*ndof2+1:i*ndof2)
1882 if (debug >= 1)
write(0,*)
'DEBUG: making BT_all done'
1885 do idom=1,hecmesh%n_neighbor_pe
1886 call hecmw_localmat_free(bt_exp(idom))
1889 end subroutine update_comm_table
1891 subroutine copy_mesh(src, dst, fg_paracon)
1893 type (hecmwst_local_mesh),
intent(in) :: src
1894 type (hecmwst_local_mesh),
intent(out) :: dst
1895 logical,
intent(in),
optional :: fg_paracon
1897 dst%MPI_COMM = src%MPI_COMM
1898 dst%PETOT = src%PETOT
1899 dst%PEsmpTOT = src%PEsmpTOT
1900 dst%my_rank = src%my_rank
1901 dst%n_subdomain = src%n_subdomain
1902 dst%n_node = src%n_node
1903 dst%nn_internal = src%nn_internal
1904 dst%n_elem = src%n_elem
1905 dst%ne_internal = src%ne_internal
1906 dst%n_elem_type = src%n_elem_type
1907 dst%n_dof = src%n_dof
1908 dst%n_neighbor_pe = src%n_neighbor_pe
1909 allocate(dst%neighbor_pe(dst%n_neighbor_pe))
1910 dst%neighbor_pe(:) = src%neighbor_pe(:)
1911 allocate(dst%import_index(0:dst%n_neighbor_pe))
1912 allocate(dst%export_index(0:dst%n_neighbor_pe))
1913 if (
present(fg_paracon) .and. fg_paracon)
then
1914 dst%import_index(:)= src%import_index(:)
1915 dst%export_index(:)= src%export_index(:)
1916 allocate(dst%import_item(dst%import_index(dst%n_neighbor_pe)))
1917 dst%import_item(:) = src%import_item(:)
1918 allocate(dst%export_item(dst%export_index(dst%n_neighbor_pe)))
1919 dst%export_item(:) = src%export_item(:)
1920 allocate(dst%global_node_ID(dst%n_node))
1921 dst%global_node_ID(1:dst%n_node) = src%global_node_ID(1:dst%n_node)
1923 dst%import_index(:)= 0
1924 dst%export_index(:)= 0
1925 dst%import_item => null()
1926 dst%export_item => null()
1928 allocate(dst%node_ID(2*dst%n_node))
1929 dst%node_ID(1:2*dst%n_node) = src%node_ID(1:2*dst%n_node)
1930 allocate(dst%elem_type_item(dst%n_elem_type))
1931 dst%elem_type_item(:) = src%elem_type_item(:)
1935 dst%node => src%node
1936 end subroutine copy_mesh
1938 subroutine free_mesh(hecMESH, fg_paracon)
1940 type (hecmwst_local_mesh),
intent(inout) :: hecmesh
1941 logical,
intent(in),
optional :: fg_paracon
1942 deallocate(hecmesh%neighbor_pe)
1943 deallocate(hecmesh%import_index)
1944 deallocate(hecmesh%export_index)
1945 if (
present(fg_paracon) .and. fg_paracon)
then
1946 deallocate(hecmesh%import_item)
1947 deallocate(hecmesh%export_item)
1948 deallocate(hecmesh%global_node_ID)
1950 if (
associated(hecmesh%import_item))
deallocate(hecmesh%import_item)
1951 if (
associated(hecmesh%export_item))
deallocate(hecmesh%export_item)
1953 deallocate(hecmesh%node_ID)
1954 deallocate(hecmesh%elem_type_item)
1956 end subroutine free_mesh
1958 subroutine extract_bt_exp(BTmat, hecMESH, BT_exp)
1960 type(hecmwst_local_matrix),
intent(in) :: btmat
1961 type(hecmwst_local_mesh),
intent(in) :: hecmesh
1962 type(hecmwst_local_matrix),
intent(out) :: bt_exp(:)
1963 integer(kind=kint) :: i, j, k, n, idx_0, idx_n, jrow, ndof2
1964 ndof2 = btmat%ndof ** 2
1965 do i = 1,hecmesh%n_neighbor_pe
1966 idx_0 = hecmesh%export_index(i-1)
1967 idx_n = hecmesh%export_index(i)
1968 bt_exp(i)%nr = idx_n - idx_0
1969 bt_exp(i)%nc = btmat%nc
1971 bt_exp(i)%ndof = btmat%ndof
1972 allocate(bt_exp(i)%index(0:bt_exp(i)%nr))
1973 bt_exp(i)%index(0) = 0
1974 do j = 1,bt_exp(i)%nr
1975 jrow = hecmesh%export_item(j + idx_0)
1976 n = btmat%index(jrow) - btmat%index(jrow-1)
1977 bt_exp(i)%nnz = bt_exp(i)%nnz + n
1978 bt_exp(i)%index(j) = bt_exp(i)%index(j-1) + n
1980 allocate(bt_exp(i)%item(bt_exp(i)%nnz))
1981 allocate(bt_exp(i)%A(bt_exp(i)%nnz * ndof2))
1983 do j = 1,bt_exp(i)%nr
1984 jrow = hecmesh%export_item(j + idx_0)
1985 do k = btmat%index(jrow-1)+1,btmat%index(jrow)
1988 bt_exp(i)%item(n) = btmat%item(k)
1989 bt_exp(i)%A(ndof2*(n-1)+1:ndof2*n) = btmat%A(ndof2*(k-1)+1:ndof2*k)
1993 end subroutine extract_bt_exp
1995 subroutine extract_cols(BT_exp, cur_export, n_curexp)
1997 type(hecmwst_local_matrix),
intent(in) :: bt_exp
1998 integer(kind=kint),
intent(out) :: cur_export(:)
1999 integer(kind=kint),
intent(out) :: n_curexp
2002 cur_export(1:bt_exp%nnz) = bt_exp%item(1:bt_exp%nnz)
2003 call quick_sort(cur_export, 1, bt_exp%nnz)
2004 call unique(cur_export, bt_exp%nnz, n_curexp)
2007 end subroutine extract_cols
2009 subroutine reorder_current_export(cur_export, n_curexp, org_export, n_orgexp, n_newexp, nn_internal)
2011 integer(kind=kint),
intent(inout) :: cur_export(:)
2012 integer(kind=kint),
intent(in) :: n_curexp
2013 integer(kind=kint),
intent(in) :: org_export(:)
2014 integer(kind=kint),
intent(in) :: n_orgexp
2015 integer(kind=kint),
intent(out) :: n_newexp
2016 integer(kind=kint),
intent(in) :: nn_internal
2017 integer(kind=kint),
allocatable :: new_export(:)
2018 integer(kind=kint) :: j, jnod
2020 allocate(new_export(n_curexp))
2022 jnod = cur_export(j)
2023 if (jnod > nn_internal) &
2024 stop
'ERROR: unknown error (jnod)'
2025 if (.not. is_included(org_export, n_orgexp, jnod))
then
2026 n_newexp = n_newexp + 1
2027 new_export(n_newexp) = jnod
2029 else if (n_newexp > 0)
then
2030 cur_export(j - n_newexp) = jnod
2034 cur_export(n_curexp - n_newexp + j) = new_export(j)
2036 deallocate(new_export)
2039 end subroutine reorder_current_export
2041 subroutine convert_bt_exp_col_id(BT_exp, cur_export, n_curexp)
2043 type(hecmwst_local_matrix),
intent(inout) :: bt_exp
2044 integer(kind=kint),
intent(in) :: cur_export(:)
2045 integer(kind=kint),
intent(in) :: n_curexp
2046 integer(kind=kint) :: i, icol, j
2049 do i = 1, bt_exp%nnz
2050 icol = bt_exp%item(i)
2053 if (icol == cur_export(j))
then
2059 if (.not. found)
then
2061 stop
'ERROR: unknown error (item not found in cur_export)'
2064 end subroutine convert_bt_exp_col_id
2066 subroutine append_commtable(n, index, item, idom, cur, ncur)
2068 integer(kind=kint),
intent(in) :: n, idom, ncur
2069 integer(kind=kint),
pointer :: index(:), item(:)
2070 integer(kind=kint),
pointer :: cur(:)
2071 integer(kind=kint),
allocatable :: tmp_index(:), tmp_item(:)
2072 integer(kind=kint) :: norg, j
2073 allocate(tmp_index(0:n))
2074 tmp_index(:) = index(:)
2076 allocate(tmp_item(norg))
2078 tmp_item(:) = item(:)
2079 if (
associated(item))
deallocate(item)
2081 allocate(item(norg + ncur))
2083 index(j) = index(j) + ncur
2085 do j = 1,tmp_index(idom)
2086 item(j) = tmp_item(j)
2089 item(tmp_index(idom)+j) = cur(j)
2091 do j = tmp_index(idom)+1,tmp_index(n)
2092 item(j+ncur) = tmp_item(j)
2094 deallocate(tmp_index, tmp_item)
2095 end subroutine append_commtable
2097 subroutine append_nodes(hecMESHtmp, n_newimp, i0)
2099 type(hecmwst_local_mesh),
intent(inout) :: hecmeshtmp
2100 integer(kind=kint),
intent(in) :: n_newimp
2101 integer(kind=kint),
intent(out) :: i0
2102 i0 = hecmeshtmp%n_node
2103 hecmeshtmp%n_node = hecmeshtmp%n_node + n_newimp
2104 end subroutine append_nodes
2106 subroutine make_cur_import(org_import, n_orgimp, old_import, n_oldimp, &
2107 n_newimp, i0, cur_import)
2109 integer(kind=kint),
intent(in) :: org_import(:), old_import(:)
2110 integer(kind=kint),
intent(in) :: n_orgimp, n_oldimp, n_newimp, i0
2111 integer(kind=kint),
intent(out) :: cur_import(:)
2113 integer(kind=kint) :: ndel, i, j
2116 do while (i <= n_orgimp .and. ndel < n_oldimp)
2117 if (org_import(i) == old_import(ndel+1))
then
2120 cur_import(i-ndel) = org_import(i)
2124 if (ndel /= n_oldimp) stop
'ERROR: unknown error (ndel)'
2126 cur_import(j-ndel) = org_import(j)
2130 cur_import(i + j) = i0+j
2132 end subroutine make_cur_import
2134 recursive subroutine quick_sort(array, id1, id2)
2136 integer(kind=kint),
intent(inout) :: array(:)
2137 integer(kind=kint),
intent(in) :: id1, id2
2138 integer(kind=kint) :: pivot, center, left, right, tmp
2139 if (id1 >= id2)
return
2140 center = (id1 + id2) / 2
2141 pivot = array(center)
2145 do while (array(left) < pivot)
2148 do while (pivot < array(right))
2151 if (left >= right)
exit
2153 array(left) = array(right)
2158 if (id1 < left-1)
call quick_sort(array, id1, left-1)
2159 if (right+1 < id2)
call quick_sort(array, right+1, id2)
2161 end subroutine quick_sort
2163 subroutine unique(array, len, newlen)
2165 integer(kind=kint),
intent(inout) :: array(:)
2166 integer(kind=kint),
intent(in) :: len
2167 integer(kind=kint),
intent(out) :: newlen
2168 integer(kind=kint) :: i, ndup
2171 if (array(i) == array(i - 1 - ndup))
then
2173 else if (ndup > 0)
then
2174 array(i - ndup) = array(i)
2178 end subroutine unique
2180 function is_included(array, len, ival)
2182 logical :: is_included
2183 integer(kind=kint),
intent(in) :: array(:)
2184 integer(kind=kint),
intent(in) :: len
2185 integer(kind=kint),
intent(in) :: ival
2186 integer(kind=kint) :: i
2187 is_included = .false.
2189 if (array(i) == ival)
then
2190 is_included = .true.
2194 end function is_included
2200 subroutine solve_no_eliminate(hecMESH,hecMAT,fstrMAT)
2202 type (hecmwst_local_mesh),
intent(in) :: hecmesh
2203 type (hecmwst_matrix ),
intent(inout) :: hecmat
2205 integer :: ndof, ndof2, nb_lag, ndofextra
2206 integer :: i, ls, le, l, j, jb_lag, ib_lag, idof, jdof, ilag, k
2207 integer :: idx, idx_lag_s, idx_lag_e, ll
2208 integer,
allocatable :: iwur(:), iwuc(:), iwlr(:), iwlc(:)
2209 type(hecmwst_matrix) :: hecmatlag
2210 real(kind=kreal) :: t1
2213 if (debug >= 1)
write(0,*)
'DEBUG: solve_no_eliminate, start', hecmw_wtime()-t1
2215 call hecmw_mat_init(hecmatlag)
2219 nb_lag = (fstrmat%num_lagrange + 2)/3
2220 hecmatlag%NDOF = ndof
2221 hecmatlag%N = hecmat%N + nb_lag
2222 hecmatlag%NP = hecmat%NP + nb_lag
2226 ndofextra = hecmatlag%N*ndof - hecmat%N*ndof - fstrmat%num_lagrange
2227 if (debug >= 1)
write(0,*)
'DEBUG: num_lagrange,nb_lag,ndofextra=',fstrmat%num_lagrange,nb_lag,ndofextra
2230 allocate(iwur(hecmat%N))
2231 allocate(iwuc(nb_lag))
2235 ls=fstrmat%indexU_lagrange(i-1)+1
2236 le=fstrmat%indexU_lagrange(i)
2238 j=fstrmat%itemU_lagrange(l)
2243 if (iwuc(j) > 0) iwur(i) = iwur(i) + 1
2249 allocate(iwlr(nb_lag))
2250 allocate(iwlc(hecmat%N))
2252 do ib_lag = 1, nb_lag
2254 i = hecmat%N + ib_lag
2256 ilag = (ib_lag-1)*ndof + idof
2257 if (ilag > fstrmat%num_lagrange)
exit
2258 ls=fstrmat%indexL_lagrange(ilag-1)+1
2259 le=fstrmat%indexL_lagrange(ilag)
2261 j=fstrmat%itemL_lagrange(l)
2266 if (iwlc(j) > 0) iwlr(ib_lag) = iwlr(ib_lag) + 1
2272 allocate(hecmatlag%indexU(0:hecmatlag%NP))
2273 hecmatlag%indexU(0) = 0
2275 hecmatlag%indexU(i) = hecmatlag%indexU(i-1) + &
2276 (hecmat%indexU(i) - hecmat%indexU(i-1)) + iwur(i)
2278 do i = hecmat%N+1, hecmatlag%N
2279 hecmatlag%indexU(i) = hecmatlag%indexU(i-1)
2281 do i = hecmatlag%N+1, hecmatlag%NP
2282 hecmatlag%indexU(i) = hecmatlag%indexU(i-1) + &
2283 (hecmat%indexU(i-nb_lag) - hecmat%indexU(i-1-nb_lag))
2285 hecmatlag%NPU = hecmatlag%indexU(hecmatlag%NP)
2289 allocate(hecmatlag%indexL(0:hecmatlag%NP))
2291 hecmatlag%indexL(i) = hecmat%indexL(i)
2293 do i = hecmat%N+1, hecmatlag%N
2294 hecmatlag%indexL(i) = hecmatlag%indexL(i-1) + iwlr(i-hecmat%N)
2296 do i = hecmatlag%N+1, hecmatlag%NP
2297 hecmatlag%indexL(i) = hecmatlag%indexL(i-1) + &
2298 (hecmat%indexL(i-nb_lag) - hecmat%indexL(i-1-nb_lag))
2300 hecmatlag%NPL = hecmatlag%indexL(hecmatlag%NP)
2304 allocate(hecmatlag%itemU(hecmatlag%NPU))
2305 allocate(hecmatlag%AU(hecmatlag%NPU*ndof2))
2308 idx = hecmatlag%indexU(i-1)+1
2310 ls = hecmat%indexU(i-1)+1
2311 le = hecmat%indexU(i)
2313 ll = hecmat%itemU(l)
2314 if (ll > hecmat%N) cycle
2315 hecmatlag%itemU(idx) = ll
2316 hecmatlag%AU((idx-1)*ndof2+1:idx*ndof2)=hecmat%AU((l-1)*ndof2+1:l*ndof2)
2321 ls=fstrmat%indexU_lagrange(i-1)+1
2322 le=fstrmat%indexU_lagrange(i)
2324 j=fstrmat%itemU_lagrange(l)
2330 if (iwuc(j) > 0)
then
2331 hecmatlag%itemU(idx) = hecmat%N + j
2337 ls=fstrmat%indexU_lagrange(i-1)+1
2338 le=fstrmat%indexU_lagrange(i)
2340 j=fstrmat%itemU_lagrange(l)
2342 jdof = j - (jb_lag - 1)*ndof
2343 do k = idx_lag_s, idx_lag_e
2344 if (hecmatlag%itemU(k) < hecmat%N + jb_lag) cycle
2345 if (hecmatlag%itemU(k) > hecmat%N + jb_lag) cycle
2348 hecmatlag%AU((k-1)*ndof2+(idof-1)*ndof+jdof) = &
2349 fstrmat%AU_lagrange((l-1)*ndof+idof)
2354 ls = hecmat%indexU(i-1)+1
2355 le = hecmat%indexU(i)
2357 ll = hecmat%itemU(l)
2358 if (ll <= hecmat%N) cycle
2359 hecmatlag%itemU(idx) = ll + nb_lag
2360 hecmatlag%AU((idx-1)*ndof2+1:idx*ndof2)=hecmat%AU((l-1)*ndof2+1:l*ndof2)
2363 if (idx /= hecmatlag%indexU(i)+1) stop
'ERROR idx indexU'
2365 do i = hecmat%N, hecmat%NP
2366 idx = hecmatlag%indexU(i+nb_lag-1)+1
2367 ls=hecmat%indexU(i-1)+1
2370 ll = hecmat%itemU(l)
2371 hecmatlag%itemU(idx) = ll + nb_lag
2372 hecmatlag%AU((idx-1)*ndof2+1:idx*ndof2)=hecmat%AU((l-1)*ndof2+1:l*ndof2)
2378 allocate(hecmatlag%itemL(hecmatlag%NPL))
2379 allocate(hecmatlag%AL(hecmatlag%NPL*ndof2))
2381 do i = 1, hecmat%indexL(hecmat%N)
2382 hecmatlag%itemL(i) = hecmat%itemL(i)
2384 do i = 1, hecmat%indexL(hecmat%N)*ndof2
2385 hecmatlag%AL(i) = hecmat%AL(i)
2388 idx = hecmat%indexL(hecmat%N) + 1
2389 do ib_lag = 1, nb_lag
2391 i = hecmat%N + ib_lag
2393 ilag = (ib_lag-1)*ndof + idof
2394 if (ilag > fstrmat%num_lagrange)
exit
2395 ls=fstrmat%indexL_lagrange(ilag-1)+1
2396 le=fstrmat%indexL_lagrange(ilag)
2398 j=fstrmat%itemL_lagrange(l)
2404 if (iwlc(j) > 0)
then
2405 hecmatlag%itemL(idx) = j
2410 if (idx /= hecmatlag%indexL(hecmat%N + ib_lag)+1)
then
2411 stop
'ERROR idx indexL'
2415 ilag = (ib_lag-1)*ndof + idof
2416 if (ilag > fstrmat%num_lagrange)
exit
2417 ls=fstrmat%indexL_lagrange(ilag-1)+1
2418 le=fstrmat%indexL_lagrange(ilag)
2420 j=fstrmat%itemL_lagrange(l)
2421 do k = idx_lag_s, idx_lag_e
2422 if (hecmatlag%itemL(k) < j) cycle
2423 if (hecmatlag%itemL(k) > j) cycle
2426 hecmatlag%AL((k-1)*ndof2+(idof-1)*ndof+jdof) = &
2427 fstrmat%AL_lagrange((l-1)*ndof+jdof)
2433 do i = hecmat%N+1, hecmat%NP
2434 idx = hecmatlag%indexL(i+nb_lag-1)+1
2435 ls=hecmat%indexL(i-1)+1
2438 ll = hecmat%itemL(l)
2439 if (ll <= hecmat%N)
then
2440 hecmatlag%itemL(idx) = ll
2442 hecmatlag%itemL(idx) = ll + nb_lag
2444 hecmatlag%AL((idx-1)*ndof2+1:idx*ndof2)=hecmat%AL((l-1)*ndof2+1:l*ndof2)
2449 deallocate(iwur, iwuc, iwlr, iwlc)
2451 allocate(hecmatlag%D(hecmatlag%NP*ndof2))
2454 do i = 1, hecmat%N*ndof2
2455 hecmatlag%D(i) = hecmat%D(i)
2458 do idof = ndof - ndofextra + 1, ndof
2459 hecmatlag%D((hecmatlag%N-1)*ndof2 + (idof-1)*ndof + idof) = 1.d0
2462 do i = hecmat%N*ndof2+1, hecmat%NP*ndof2
2463 hecmatlag%D(i + nb_lag*ndof2) = hecmat%D(i)
2466 allocate(hecmatlag%B(hecmatlag%NP*ndof))
2467 allocate(hecmatlag%X(hecmatlag%NP*ndof))
2472 do i = 1, hecmat%N*ndof
2473 hecmatlag%B(i) = hecmat%B(i)
2476 do i = 1, fstrmat%num_lagrange
2477 hecmatlag%B(hecmat%N*ndof + i) = hecmat%B(hecmat%NP*ndof + i)
2480 do i = hecmat%N*ndof+1, hecmat%NP*ndof
2481 hecmatlag%B(i + nb_lag*ndof) = hecmat%B(i)
2484 hecmatlag%Iarray=hecmat%Iarray
2485 hecmatlag%Rarray=hecmat%Rarray
2487 if (debug >= 1)
write(0,*)
'DEBUG: made hecMATLag', hecmw_wtime()-t1
2489 if (hecmesh%n_neighbor_pe > 0)
then
2490 do i = 1, hecmesh%import_index(hecmesh%n_neighbor_pe)
2491 hecmesh%import_item(i) = hecmesh%import_item(i) + nb_lag
2497 hecmat%Iarray=hecmatlag%Iarray
2498 hecmat%Rarray=hecmatlag%Rarray
2500 if (hecmesh%n_neighbor_pe > 0)
then
2501 do i = 1, hecmesh%import_index(hecmesh%n_neighbor_pe)
2502 hecmesh%import_item(i) = hecmesh%import_item(i) - nb_lag
2508 do i = 1, hecmat%N*ndof
2509 hecmat%X(i) = hecmatlag%X(i)
2512 do i = hecmat%N*ndof+1, hecmat%NP*ndof
2513 hecmat%X(i) = hecmatlag%X(i + nb_lag*ndof)
2516 do i = 1, fstrmat%num_lagrange
2517 hecmat%X(hecmat%NP*ndof + i) = hecmatlag%X(hecmat%N*ndof + i)
2520 call hecmw_mat_finalize(hecmatlag)
2522 if (debug >= 1)
write(0,*)
'DEBUG: solve_no_eliminate end', hecmw_wtime()-t1
2523 end subroutine solve_no_eliminate
subroutine hecmw_solve(hecmesh, hecmat)
This module defined coomon data and basic structures for analysis.
integer(kind=kint) myrank
PARALLEL EXECUTION.
logical paracontactflag
PARALLEL CONTACT FLAG.