FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
solve_LINEQ_iter_contact.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
8  use m_fstr
10  use hecmw_solver
11 
12  private
14  public :: solve_lineq_iter_contact
15 
16  logical, save :: INITIALIZED = .false.
17  integer, save :: SymType = 0
18 
19  integer, parameter :: DEBUG = 0 ! 0: no message, 1: some messages, 2: more messages, 3: vector output, 4: matrix output
20 
21 contains
22 
23  subroutine solve_lineq_iter_contact_init(hecMESH,hecMAT,fstrMAT,is_sym)
24  implicit none
25  type (hecmwst_local_mesh), intent(in) :: hecmesh
26  type (hecmwst_matrix ), intent(inout) :: hecmat
27  type (fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
28  logical, intent(in) :: is_sym
29 
30  if (initialized) then
31  initialized = .false.
32  endif
33 
34  hecmat%Iarray(98) = 1
35  hecmat%Iarray(97) = 1
36 
37  if (is_sym) then
38  symtype = 1
39  else
40  symtype = 0
41  endif
42 
43  initialized = .true.
44  end subroutine solve_lineq_iter_contact_init
45 
46  subroutine solve_lineq_iter_contact(hecMESH,hecMAT,fstrMAT,istat,conMAT)
47  implicit none
48  type (hecmwst_local_mesh), intent(in) :: hecmesh
49  type (hecmwst_matrix ), intent(inout) :: hecmat
50  type (fstrst_matrix_contact_lagrange), intent(inout) :: fstrmat
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
55  logical :: fg_amg
56  integer(kind=kint) :: num_lagrange_global
57 
58  hecmat%Iarray(97) = 1
59 
60  ! set if use eliminate version or not
61  fg_amg = .false.
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)
66  else
67  fg_eliminate = .true.
68  if (precond_org == 5) fg_amg = .true.
69  endif
70 
71  num_lagrange_global = fstrmat%num_lagrange
72  call hecmw_allreduce_i1(hecmesh, num_lagrange_global, hecmw_sum)
73 
74  if (num_lagrange_global == 0) then
75  if ((debug >= 1 .and. myrank==0) .or. debug >= 2) write(0,*) 'DEBUG: no contact'
76  ! use CG because the matrix is symmetric
77  method_org = hecmw_mat_get_method(hecmat)
78  call hecmw_mat_set_method(hecmat, 1)
79  ! avoid ML when no contact
80  !if (fg_amg) call hecmw_mat_set_precond(hecMAT, 3) ! set diag-scaling
81  ! solve
82  call hecmw_solve(hecmesh, hecmat)
83  ! restore solver setting
84  call hecmw_mat_set_method(hecmat, method_org)
85  !if (fg_amg) call hecmw_mat_set_precond(hecMAT, 5)
86  else
87  if ((debug >= 1 .and. myrank==0) .or. debug >= 2) write(0,*) 'DEBUG: with contact'
88  if (fg_eliminate) then
89  if(paracontactflag.and.present(conmat)) then
90  call solve_eliminate(hecmesh, hecmat, fstrmat, conmat)
91  else
92  call solve_eliminate(hecmesh, hecmat, fstrmat)
93  endif
94  else
95  if(paracontactflag.and.present(conmat)) then
96  write(0,*) 'ERROR: solve_no_eliminate not implemented for ParaCon'
97  call hecmw_abort( hecmw_comm_get_comm())
98  endif
99  if (fstrmat%num_lagrange > 0) then
100  call solve_no_eliminate(hecmesh, hecmat, fstrmat)
101  else
102  call hecmw_solve(hecmesh, hecmat)
103  endif
104  endif
105  endif
106 
107  if (precond_org >= 30) then
108  call hecmw_mat_set_precond(hecmat, precond_org)
109  endif
110 
111  istat = hecmw_mat_get_flag_diverged(hecmat)
112  end subroutine solve_lineq_iter_contact
113 
114  !!
115  !! Solve with elimination of Lagrange-multipliers
116  !!
117 
118  subroutine solve_eliminate(hecMESH,hecMAT,fstrMAT,conMAT)
119  implicit none
120  type (hecmwst_local_mesh), intent(in), target :: hecmesh
121  type (hecmwst_matrix ), intent(inout) :: hecmat
122  type (fstrst_matrix_contact_lagrange), intent(inout) :: fstrmat
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(:)
141 
142  t1 = hecmw_wtime()
143  if ((debug >= 1 .and. myrank==0) .or. debug >= 2) write(0,*) 'DEBUG: solve_eliminate start', hecmw_wtime()-t1
144 
145  fg_paracon = (paracontactflag.and.present(conmat))
146 
147  ndof=hecmat%NDOF
148  if(fg_paracon) then
149  allocate(iw2(hecmat%NP*ndof))
150  else
151  allocate(iw2(hecmat%N*ndof))
152  endif
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))
156 
157  ! choose slave DOFs to be eliminated with Lag. DOFs
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
160 
161  ! make transformation matrix and its transpose
162  call make_btmat(hecmat, fstrmat, iw2, wsl, btmat, fg_paracon)
163  if (debug >= 2) write(0,*) ' DEBUG2: make T done', hecmw_wtime()-t1
164  if (debug >= 4) then
165  write(1000+myrank,*) 'BTmat (local)'
166  call hecmw_localmat_write(btmat, 1000+myrank)
167  endif
168 
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
171  if (debug >= 4) then
172  write(1000+myrank,*) 'BTtmat (local)'
173  call hecmw_localmat_write(bttmat, 1000+myrank)
174  endif
175 
176  if(fg_paracon) then
177  ! make contact dof list
178  call make_contact_dof_list(hecmat, fstrmat, n_contact_dof, contact_dofs)
179 
180  ! make comm_table for contact dof
181  ! call fstr_contact_comm_init(conCOMM, hecMESH, ndof, fstrMAT%num_lagrange, iwS)
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
184 
185  ! copy hecMESH to hecMESHtmp
186  allocate(hecmeshtmp)
187  call copy_mesh(hecmesh, hecmeshtmp, fg_paracon)
188 
189  ! communicate and complete BTmat (update hecMESHtmp)
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
194  endif
195  if (debug >= 4) then
196  write(1000+myrank,*) 'BTmat (assembled)'
197  call hecmw_localmat_write(btmat, 1000+myrank)
198  endif
199 
200  ! communicate and complete BTtmat (update hecMESHtmp)
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
205  btmat%nc = bttmat%nc
206  endif
207  if (debug >= 4) then
208  write(1000+myrank,*) 'BTtmat (assembled)'
209  call hecmw_localmat_write(bttmat, 1000+myrank)
210  endif
211 
212  ! place 1 on diag of non-slave dofs of BTmat and BTtmat
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
218  if (debug >= 4) then
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)
223  endif
224 
225  ! init BKmat and substitute conMAT
226  allocate(bkmat)
227  call hecmw_localmat_init_with_hecmat(bkmat, conmat, fstrmat%num_lagrange)
228  if (debug >= 4) then
229  write(1000+myrank,*) 'BKmat (conMAT local)'
230  call hecmw_localmat_write(bkmat, 1000+myrank)
231  endif
232 
233  ! communicate and complete BKmat (update hecMESHtmp)
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
238  btmat%nc = bkmat%nc
239  bttmat%nc = bkmat%nc
240  endif
241  if (debug >= 4) then
242  write(1000+myrank,*) 'BKmat (conMAT assembled)'
243  call hecmw_localmat_write(bkmat, 1000+myrank)
244  endif
245 
246  ! add hecMAT to BKmat
247  call hecmw_localmat_add_hecmat(bkmat, hecmat)
248  if (debug >= 2) write(0,*) ' DEBUG2: add hecMAT to K done', hecmw_wtime()-t1
249  if (debug >= 4) then
250  write(1000+myrank,*) 'BKmat (hecMAT added)'
251  call hecmw_localmat_write(bkmat, 1000+myrank)
252  endif
253 
254  ! compute BTtKmat = BTtmat * BKmat (update hecMESHtmp)
255  allocate(bttkmat)
256  call hecmw_localmat_multmat(bttmat, bkmat, hecmeshtmp, bttkmat)
257  if (debug >= 2) write(0,*) ' DEBUG2: multiply Tt and K done', hecmw_wtime()-t1
258  if (debug >= 4) then
259  write(1000+myrank,*) 'BTtKmat'
260  call hecmw_localmat_write(bttkmat, 1000+myrank)
261  endif
262 
263  ! compute BTtKTmat = BTtKmat * BTmat (update hecMESHtmp)
264  allocate(bttktmat)
265  call hecmw_localmat_multmat(bttkmat, btmat, hecmeshtmp, bttktmat)
266  if (debug >= 2) write(0,*) ' DEBUG2: multiply TtK and T done', hecmw_wtime()-t1
267  if (debug >= 4) then
268  write(1000+myrank,*) 'BTtKTmat'
269  call hecmw_localmat_write(bttktmat, 1000+myrank)
270  endif
271  call hecmw_localmat_free(bttkmat)
272  deallocate(bttkmat)
273 
274  ! shrink comm_table
275  ! call hecmw_localmat_shrink_comm_table(BTtKTmat, hecMESHtmp)
276 
277  call place_num_on_diag_of_marked_dof(bttktmat, 1.0d0, mark)
278  if (debug >= 4) then
279  write(1000+myrank,*) 'BTtKTmat (place 1.0 on slave diag)'
280  call hecmw_localmat_write(bttktmat, 1000+myrank)
281  endif
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)
286  deallocate(bttktmat)
287  !
288  bt_all => btmat
289  else
290  if (hecmesh%n_neighbor_pe > 0) then
291  ! update communication table
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)
296  else
297  ! in serial computation
298  hecmeshtmp => hecmesh
299  bt_all => btmat
300  end if
301 
302  ! calc trimatmul in hecmwST_matrix data structure
303  call hecmw_mat_init(hectkt)
304  call hecmw_trimatmul_ttkt(hecmeshtmp, bttmat, hecmat, bt_all, iws, fstrmat%num_lagrange, hectkt)
305  endif
306  if ((debug >= 1 .and. myrank==0) .or. debug >= 2) write(0,*) 'DEBUG: calculated TtKT', hecmw_wtime()-t1
307 
308  ! original RHS
309  if (debug >= 3) then
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)
318  endif
319  if (n_slave > 0) then
320  write(1000+myrank,*) 'hecMAT%B(slave):',slaves(:)
321  write(1000+myrank,*) hecmat%B(slaves(:))
322  endif
323  endif
324 
325  if (fg_paracon) then
326  ! contact RHS
327  if (debug >= 3) then
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)
337  endif
338  if (n_slave > 0) then
339  write(1000+myrank,*) 'conMAT%B(slave):',slaves(:)
340  write(1000+myrank,*) conmat%B(slaves(:))
341  endif
342  endif
343 
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
347  if (debug >= 3) then
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)
355  endif
356  if (n_slave > 0) then
357  write(1000+myrank,*) 'Btot(slave):',slaves(:)
358  write(1000+myrank,*) btot(slaves(:))
359  endif
360  endif
361 
362  do i=1,hecmat%N*ndof
363  btot(i)=btot(i)+hecmat%B(i)
364  enddo
365  ! assembled RHS
366  if (debug >= 3) then
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)
374  endif
375  if (n_slave > 0) then
376  write(1000+myrank,*) 'Btot(slave):',slaves(:)
377  write(1000+myrank,*) btot(slaves(:))
378  endif
379  endif
380  endif
381 
382  allocate(hectkt%B(hectkt%NP*ndof + fstrmat%num_lagrange))
383  allocate(hectkt%X(hectkt%NP*ndof + fstrmat%num_lagrange))
384  if (fg_paracon) then
385  do i=1, hectkt%N*ndof
386  hectkt%B(i) = btot(i)
387  enddo
388  else
389  do i=1, hectkt%N*ndof
390  hectkt%B(i) = hecmat%B(i)
391  enddo
392  endif
393  do ilag=1, fstrmat%num_lagrange
394  hectkt%B(hectkt%NP*ndof + ilag) = hecmat%B(hecmat%NP*ndof + ilag)
395  enddo
396  do i=1, hectkt%N*ndof
397  hectkt%X(i) = hecmat%X(i)
398  enddo
399  do ilag=1,fstrmat%num_lagrange
400  hectkt%X(iws(ilag)) = 0.d0
401  enddo
402 
403  ! make new RHS
404  if (fg_paracon) then
405  call make_new_b_paracon(hecmesh, hecmat, conmat, btot, hecmeshtmp, hectkt, bttmat, bkmat, &
406  iws, wsl, fstrmat%num_lagrange, concomm, hectkt%B)
407  else
408  call make_new_b(hecmesh, hecmat, bttmat, iws, wsl, &
409  fstrmat%num_lagrange, hectkt%B)
410  endif
411  if ((debug >= 1 .and. myrank==0) .or. debug >= 2) write(0,*) 'DEBUG: converted RHS', hecmw_wtime()-t1
412  if (debug >= 3) then
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)
417  endif
418 
419  ! ! use CG when the matrix is symmetric
420  ! if (SymType == 1) call hecmw_mat_set_method(hecTKT, 1)
421 
422  ! solve
423  call hecmw_solve(hecmeshtmp, hectkt)
424  if ((debug >= 1 .and. myrank==0) .or. debug >= 2) write(0,*) 'DEBUG: linear solver finished', hecmw_wtime()-t1
425 
426  ! solution in converged system
427  if (debug >= 3) then
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)
432  endif
433 
434  hecmat%Iarray=hectkt%Iarray
435  hecmat%Rarray=hectkt%Rarray
436 
437  ! calc u_s
438  if (fg_paracon) then
439  call comp_x_slave_paracon(hecmesh, hecmat, btot, hecmeshtmp, hectkt, btmat, &
440  fstrmat%num_lagrange, iws, wsl, concomm, n_slave, slaves)
441  else
442  call hecmw_localmat_mulvec(bt_all, hectkt%X, hecmat%X) !!!<== maybe, BT_all should be BTmat ???
443  call subst_blag(hecmat, iws, wsl, fstrmat%num_lagrange)
444  endif
445  if ((debug >= 1 .and. myrank==0) .or. debug >= 2) write(0,*) 'DEBUG: recovered slave disp', hecmw_wtime()-t1
446 
447  ! calc lambda
448  if (fg_paracon) then
449  call comp_lag_paracon(hecmesh, hecmat, btot, hecmeshtmp, hectkt, bkmat, &
450  n_slave, slaves, fstrmat%num_lagrange, iws, wsu, concomm)
451  ! call comp_lag_paracon2(hecMESH, hecMAT, conMAT, fstrMAT%num_lagrange, iwS, wSU, conCOMM)
452  else
453  call comp_lag(hecmat, iws, wsu, fstrmat%num_lagrange)
454  endif
455  if ((debug >= 1 .and. myrank==0) .or. debug >= 2) write(0,*) 'DEBUG: calculated lag', hecmw_wtime()-t1
456 
457  ! write(0,*) 'size of conMAT%X',size(conMAT%X)
458  ! conMAT%X(:) = hecMAT%X(:)
459 
460  ! solution in original system
461  if (debug >= 3) then
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)
469  endif
470  if (n_slave > 0) then
471  write(1000+myrank,*) 'hecMAT%X(slave):',slaves(:)
472  write(1000+myrank,*) hecmat%X(slaves(:))
473  endif
474  endif
475 
476  ! check solution in the original system
477  if (fg_paracon .and. debug >= 2) then
478  ! call check_solution2(hecMESH, hecMAT, conMAT, fstrMAT, n_contact_dof, contact_dofs, &
479  ! conCOMM, n_slave, slaves)
480  call check_solution(hecmesh, hecmat, fstrmat, btot, hecmeshtmp, hectkt, bkmat, &
481  concomm, n_slave, slaves)
482  endif
483 
484  ! free matrices
485  call hecmw_localmat_free(bt_all)
486  call hecmw_localmat_free(bttmat)
487  call hecmw_mat_finalize(hectkt)
488  if (fg_paracon) then
489  call hecmw_localmat_free(bkmat)
490  deallocate(bkmat)
491  call fstr_contact_comm_finalize(concomm)
492  call free_mesh(hecmeshtmp, fg_paracon)
493  deallocate(hecmeshtmp)
494  else
495  if (hecmesh%n_neighbor_pe > 0) then
496  call free_mesh(hecmeshtmp)
497  deallocate(hecmeshtmp)
498  deallocate(bt_all)
499  endif
500  endif
501  deallocate(iw2, iws)
502  if ((debug >= 1 .and. myrank==0) .or. debug >= 2) write(0,*) 'DEBUG: solve_eliminate end', hecmw_wtime()-t1
503  end subroutine solve_eliminate
504 
505  subroutine choose_slaves(hecMAT, fstrMAT, iw2, iwS, wSL, wSU, fg_paracon)
506  implicit none
507  type (hecmwst_matrix ), intent(in) :: hecmat
508  type (fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
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
517 
518  iw2=-1
519 
520  if (fstrmat%num_lagrange == 0) return
521 
522  ndof=hecmat%NDOF
523  iws=0
524 
525  if (fg_paracon) then
526  n = hecmat%NP
527  else
528  n = hecmat%N
529  endif
530 
531  allocate(iw1l(n*ndof))
532  allocate(iw1u(n*ndof))
533  iw1l=0
534  iw1u=0
535 
536  ! Count how many times each dof appear in Lagrange matrix
537  ! lower
538  do i=1,fstrmat%num_lagrange
539  ls=fstrmat%indexL_lagrange(i-1)+1
540  le=fstrmat%indexL_lagrange(i)
541  do l=ls,le
542  j=fstrmat%itemL_lagrange(l)
543  do jdof=1,ndof
544  idx=(j-1)*ndof+jdof
545  iw1l(idx)=iw1l(idx)+1
546  enddo
547  enddo
548  enddo
549  ! upper
550  do i=1,n
551  ls=fstrmat%indexU_lagrange(i-1)+1
552  le=fstrmat%indexU_lagrange(i)
553  do l=ls,le
554  j=fstrmat%itemU_lagrange(l)
555  do idof=1,ndof
556  idx=(i-1)*ndof+idof
557  iw1u(idx)=iw1u(idx)+1
558  enddo
559  enddo
560  enddo
561  !!$ write(0,*) 'iw1L, iw1U:'
562  !!$ do i=1,n*ndof
563  !!$ if (iw1L(i) > 0 .or. iw1U(i) > 0) write(0,*) i, iw1L(i), iw1U(i)
564  !!$ enddo
565 
566  ! Choose dofs that
567  ! - appear only onece in both lower and upper Lag. and
568  ! - has greatest coefficient among them (in lower Lag.)
569  do i=1,fstrmat%num_lagrange
570  ls=fstrmat%indexL_lagrange(i-1)+1
571  le=fstrmat%indexL_lagrange(i)
572  vmax = 0.d0
573  imax = -1
574  do l=ls,le
575  j=fstrmat%itemL_lagrange(l)
576  do jdof=1,ndof
577  idx=(j-1)*ndof+jdof
578  val=fstrmat%AL_lagrange((l-1)*ndof+jdof)
579  if (iw1l(idx) == 1 .and. iw1u(idx) == 1 .and. abs(val) > abs(vmax)) then
580  imax=idx
581  vmax=val
582  endif
583  enddo
584  enddo
585  if (imax == -1) stop "ERROR: iterative solver for contact failed"
586  iw2(imax)=i
587  iws(i)=imax; wsl(i)=-1.d0/vmax
588  enddo
589  !!$ write(0,*) 'iw2:'
590  !!$ do i=1,n*ndof
591  !!$ if (iw2(i) > 0) write(0,*) i, iw2(i), iw1L(i), iw1U(i)
592  !!$ enddo
593  !!$ write(0,*) 'iwS:'
594  !!$ write(0,*) iwS(:)
595  n_slave_in = 0
596  do i=1,hecmat%N*ndof
597  if (iw2(i) > 0) n_slave_in = n_slave_in + 1
598  enddo
599  n_slave_out = 0
600  do i=hecmat%N*ndof,n*ndof
601  if (iw2(i) > 0) n_slave_out = n_slave_out + 1
602  enddo
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
604 
605  deallocate(iw1l, iw1u)
606 
607  call make_wsu(fstrmat, n, ndof, iw2, wsu)
608  end subroutine choose_slaves
609 
610  subroutine make_wsu(fstrMAT, n, ndof, iw2, wSU)
611  implicit none
612  type(fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
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
617 
618  if (fstrmat%num_lagrange == 0) return
619 
620  wsu=0.d0
621  do i=1,n
622  do idof=1,ndof
623  idx=(i-1)*ndof+idof
624  if (iw2(idx) > 0) then
625  js=fstrmat%indexU_lagrange(i-1)+1
626  je=fstrmat%indexU_lagrange(i)
627  do j=js,je
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)
631  endif
632  enddo
633  endif
634  enddo
635  enddo
636  !write(0,*) wSU
637  end subroutine make_wsu
638 
639  subroutine make_btmat(hecMAT, fstrMAT, iw2, wSL, BTmat, fg_paracon)
640  implicit none
641  type (hecmwst_matrix ), intent(inout) :: hecmat
642  type (fstrst_matrix_contact_lagrange), intent(inout) :: fstrmat
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
650 
651  ndof=hecmat%NDOF
652  if (fg_paracon) then
653  n=hecmat%NP
654  else
655  n=hecmat%N
656  endif
657  tmat%nr=n*ndof
658  tmat%nc=tmat%nr
659  if (fg_paracon) then
660  tmat%nnz=fstrmat%numL_lagrange*ndof-fstrmat%num_lagrange
661  else
662  tmat%nnz=tmat%nr+fstrmat%numL_lagrange*ndof-2*fstrmat%num_lagrange
663  endif
664  tmat%ndof=1
665 
666  allocate(tmat%index(0:tmat%nr))
667  allocate(tmat%item(tmat%nnz), tmat%A(tmat%nnz))
668  ! index
669  tmat%index(0)=0
670  do i=1,tmat%nr
671  if (iw2(i) > 0) then
672  nnz=ndof*(fstrmat%indexL_lagrange(iw2(i))-fstrmat%indexL_lagrange(iw2(i)-1))-1
673  else
674  if (fg_paracon) then
675  nnz = 0
676  else
677  nnz=1
678  endif
679  endif
680  tmat%index(i)=tmat%index(i-1)+nnz
681  enddo
682  if (tmat%nnz /= tmat%index(tmat%nr)) then
683  write(0,*) tmat%nnz, tmat%index(tmat%nr)
684  stop 'ERROR: Tmat%nnz wrong'
685  endif
686  ! item and A
687  do i=1,tmat%nr
688  l=tmat%index(i-1)+1
689  if (iw2(i) > 0) then
690  js=fstrmat%indexL_lagrange(iw2(i)-1)+1
691  je=fstrmat%indexL_lagrange(iw2(i))
692  factor=wsl(iw2(i))
693  do j=js,je
694  k=fstrmat%itemL_lagrange(j)
695  do jdof=1,ndof
696  kk=(k-1)*ndof+jdof
697  jj=(j-1)*ndof+jdof
698  if (kk==i) cycle
699  tmat%item(l)=kk
700  tmat%A(l)=fstrmat%AL_lagrange(jj)*factor
701  l=l+1
702  enddo
703  enddo
704  else
705  if (.not. fg_paracon) then
706  tmat%item(l)=i
707  tmat%A(l)=1.0d0
708  l=l+1
709  endif
710  endif
711  if (l /= tmat%index(i)+1) then
712  write(0,*) l, tmat%index(i)+1
713  stop 'ERROR: Tmat%index wrong'
714  endif
715  enddo
716  !call hecmw_localmat_write(Tmat, 0)
717  ! make 3x3-block version of Tmat
718  call hecmw_localmat_blocking(tmat, ndof, btmat)
719  call hecmw_localmat_free(tmat)
720  end subroutine make_btmat
721 
722  subroutine make_bttmat(hecMAT, fstrMAT, iw2, iwS, wSU, BTtmat, fg_paracon)
723  implicit none
724  type (hecmwst_matrix ), intent(inout) :: hecmat
725  type (fstrst_matrix_contact_lagrange), intent(inout) :: fstrmat
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
732 
733  ndof=hecmat%NDOF
734  if (fg_paracon) then
735  n=hecmat%NP
736  else
737  n=hecmat%N
738  endif
739  ttmat%nr=n*ndof
740  ttmat%nc=ttmat%nr
741  if (fg_paracon) then
742  ttmat%nnz=fstrmat%numU_lagrange*ndof-fstrmat%num_lagrange
743  else
744  ttmat%nnz=ttmat%nr+fstrmat%numU_lagrange*ndof-2*fstrmat%num_lagrange
745  endif
746  ttmat%ndof=1
747 
748  allocate(ttmat%index(0:ttmat%nr))
749  allocate(ttmat%item(ttmat%nnz), ttmat%A(ttmat%nnz))
750  ! index
751  ttmat%index(0)=0
752  do i=1,n
753  do idof=1,ndof
754  idx=(i-1)*ndof+idof
755  if (iw2(idx) <= 0) then
756  if (fstrmat%num_lagrange == 0) then
757  if (fg_paracon) then
758  nnz=0
759  else
760  nnz=1
761  endif
762  else
763  if (fg_paracon) then
764  nnz=fstrmat%indexU_lagrange(i)-fstrmat%indexU_lagrange(i-1)
765  else
766  nnz=fstrmat%indexU_lagrange(i)-fstrmat%indexU_lagrange(i-1)+1
767  endif
768  endif
769  else
770  nnz=0
771  endif
772  ttmat%index(idx)=ttmat%index(idx-1)+nnz
773  enddo
774  enddo
775  if (ttmat%nnz /= ttmat%index(ttmat%nr)) then
776  write(0,*) ttmat%nnz, ttmat%index(ttmat%nr)
777  stop 'ERROR: Ttmat%nnz wrong'
778  endif
779  ! item and A
780  do i=1,n
781  do idof=1,ndof
782  idx=(i-1)*ndof+idof
783  l=ttmat%index(idx-1)+1
784  if (iw2(idx) <= 0) then
785  if (.not. fg_paracon) then
786  ! one on diagonal
787  ttmat%item(l)=idx
788  ttmat%A(l)=1.0d0
789  l=l+1
790  endif
791  if (fstrmat%num_lagrange > 0) then
792  ! offdiagonal
793  js=fstrmat%indexU_lagrange(i-1)+1
794  je=fstrmat%indexU_lagrange(i)
795  do j=js,je
796  k=fstrmat%itemU_lagrange(j)
797  ttmat%item(l)=iws(k)
798  ttmat%A(l)=fstrmat%AU_lagrange((j-1)*ndof+idof)*wsu(k)
799  l=l+1
800  enddo
801  endif
802  else
803  ! no element
804  endif
805  if (l /= ttmat%index(idx)+1) then
806  write(0,*) l, ttmat%index(idx)+1
807  stop 'ERROR: Ttmat%index wrong'
808  endif
809  enddo
810  enddo
811  !call hecmw_localmat_write(Ttmat, 0)
812  ! make 3x3-block version of Tmat
813  call hecmw_localmat_blocking(ttmat, ndof, bttmat)
814  call hecmw_localmat_free(ttmat)
815  end subroutine make_bttmat
816 
817  subroutine make_contact_dof_list(hecMAT, fstrMAT, n_contact_dof, contact_dofs)
818  implicit none
819  type (hecmwst_matrix), intent(in) :: hecmat
820  type (fstrst_matrix_contact_lagrange), intent(in) :: fstrmat
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
827  n_contact_dof = 0
828  return
829  endif
830  ! lower
831  ndof = hecmat%NDOF
832  allocate(iw(hecmat%NP * ndof))
833  icnt = 0
834  do i = 1, fstrmat%num_lagrange
835  ls = fstrmat%indexL_lagrange(i-1)+1
836  le = fstrmat%indexL_lagrange(i)
837  do l = ls, le
838  j = fstrmat%itemL_lagrange(l)
839  ljdof1: do jdof = 1, ndof
840  val = fstrmat%AL_lagrange((l-1)*ndof+jdof)
841  !write(0,*) 'j,jdof,val',j,jdof,val
842  if (val == 0.0d0) cycle
843  idx = (j-1)*ndof+jdof
844  do k = 1, icnt
845  if (iw(k) == idx) cycle ljdof1
846  enddo
847  icnt = icnt + 1
848  iw(icnt) = idx
849  !write(0,*) 'icnt,idx',icnt,idx
850  enddo ljdof1
851  enddo
852  enddo
853  ! upper
854  do i = 1, hecmat%NP
855  ls = fstrmat%indexU_lagrange(i-1)+1
856  le = fstrmat%indexU_lagrange(i)
857  do l = ls, le
858  j = fstrmat%itemU_lagrange(l)
859  lidof1: do idof = 1, ndof
860  val = fstrmat%AU_lagrange((l-1)*ndof+idof)
861  !write(0,*) 'i,idof,val',i,idof,val
862  if (val == 0.0d0) cycle
863  idx = (i-1)*ndof+idof
864  do k = 1, icnt
865  if (iw(k) == idx) cycle lidof1
866  enddo
867  icnt = icnt + 1
868  iw(icnt) = idx
869  !write(0,*) 'icnt,idx',icnt,idx
870  enddo lidof1
871  enddo
872  enddo
873  call quick_sort(iw, 1, icnt)
874  allocate(contact_dofs(icnt))
875  contact_dofs(1:icnt) = iw(1:icnt)
876  n_contact_dof = icnt
877  deallocate(iw)
878  if (debug >= 2) write(0,*) ' DEBUG2: contact_dofs',contact_dofs(:)
879  end subroutine make_contact_dof_list
880 
881  subroutine make_new_b(hecMESH, hecMAT, BTtmat, iwS, wSL, num_lagrange, Bnew)
882  implicit none
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
892 
893  npndof=hecmat%NP*hecmat%NDOF
894  nndof=hecmat%N*hecmat%NDOF
895 
896  allocate(btmp(npndof))
897 
898  !!! BTtmat*(B+K*(-Bs^-1)*Blag)
899 
900  !B2=-Bs^-1*Blag
901  bnew=0.d0
902  do i=1,num_lagrange
903  bnew(iws(i))=wsl(i)*hecmat%B(npndof+i)
904  enddo
905  !Btmp=B+K*B2
906  call hecmw_matvec(hecmesh, hecmat, bnew, btmp)
907  do i=1,nndof
908  btmp(i)=hecmat%B(i)+btmp(i)
909  enddo
910  !B2=BTtmat*Btmp
911  call hecmw_localmat_mulvec(bttmat, btmp, bnew)
912 
913  deallocate(btmp)
914  end subroutine make_new_b
915 
916  subroutine assemble_b_paracon(hecMESH, hecMAT, conMAT, num_lagrange, conCOMM, Btot)
917  implicit none
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
924  ndof = hecmat%NDOF
925  nndof = hecmat%N * ndof
926  npndof = hecmat%NP * ndof
927  do i=1,npndof+num_lagrange
928  btot(i) = conmat%B(i)
929  enddo
930  call fstr_contact_comm_reduce_r(concomm, btot, hecmw_sum)
931  end subroutine assemble_b_paracon
932 
933  subroutine make_new_b_paracon(hecMESH, hecMAT, conMAT, Btot, hecMESHtmp, hecTKT, BTtmat, BKmat, &
934  iwS, wSL, num_lagrange, conCOMM, Bnew)
935  implicit none
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
947  ! SIZE:
948  ! Btot <=> hecMAT, hecMESH
949  ! Btmp <=> BKmat, hecMESHtmp
950  ! Bnew <=> hecTKT, hecMESHtmp
951  npndof = hecmat%NP*hecmat%NDOF
952  nndof = hecmat%N *hecmat%NDOF
953  npndof_new = hectkt%NP*hectkt%NDOF
954  allocate(btmp(npndof_new))
955  !
956  !! BTtmat*(B+K*(-Bs^-1)*Blag)
957  !
958  ! B2=-Bs^-1*Blag
959  bnew=0.d0
960  do i=1,num_lagrange
961  !Bnew(iwS(i))=wSL(i)*conMAT%B(npndof+i)
962  bnew(iws(i))=wsl(i)*btot(npndof+i)
963  enddo
964  ! send external contact dof => recv internal contact dof
965  call fstr_contact_comm_reduce_r(concomm, bnew, hecmw_sum)
966  ! Btmp=B+K*B2 (including update of Bnew)
967  call hecmw_update_3_r(hecmeshtmp, bnew, hectkt%NP)
968  call hecmw_localmat_mulvec(bkmat, bnew, btmp)
969  do i=1,nndof
970  btmp(i)=btot(i)+btmp(i)
971  enddo
972  ! B2=BTtmat*Btmp
973  call hecmw_update_3_r(hecmeshtmp, btmp, hectkt%NP)
974  call hecmw_localmat_mulvec(bttmat, btmp, bnew)
975  deallocate(btmp)
976  end subroutine make_new_b_paracon
977 
978  subroutine subst_blag(hecMAT, iwS, wSL, num_lagrange)
979  implicit none
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
985 
986  npndof=hecmat%NP*hecmat%NDOF
987  do i=1,num_lagrange
988  ilag=iws(i)
989  hecmat%X(ilag)=hecmat%X(ilag)-wsl(i)*hecmat%B(npndof+i)
990  enddo
991  end subroutine subst_blag
992 
993  subroutine comp_x_slave_paracon(hecMESH, hecMAT, Btot, hecMESHtmp, hecTKT, BTmat, &
994  num_lagrange, iwS, wSL, conCOMM, n_slave, slaves)
995  implicit none
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(:)
1009  ndof = hecmat%NDOF
1010  ndof2 = ndof*ndof
1011  npndof = hecmat%NP * ndof
1012  nndof = hecmat%N * ndof
1013  !!
1014  !! {X} = [BT] {Xp} - [-Bs^-1] {c}
1015  !!
1016  ! compute {X} = [BT] {Xp}
1017  call hecmw_update_3_r(hecmeshtmp, hectkt%X, hectkt%NP)
1018  call hecmw_localmat_mulvec(btmat, hectkt%X, hecmat%X)
1019  !
1020  ! compute {Xtmp} = [-Bs^-1] {c}
1021  allocate(xtmp(npndof))
1022  xtmp(:) = 0.0d0
1023  do ilag = 1, num_lagrange
1024  islave = iws(ilag)
1025  !Xtmp(islave) = wSL(ilag) * conMAT%B(npndof + ilag)
1026  xtmp(islave) = wsl(ilag) * btot(npndof + ilag)
1027  enddo
1028  !
1029  ! send external contact dof => recv internal contact dof
1030  call fstr_contact_comm_reduce_r(concomm, xtmp, hecmw_sum)
1031  !
1032  ! {X} = {X} - {Xtmp}
1033  do i = 1, n_slave
1034  islave = slaves(i)
1035  hecmat%X(islave) = hecmat%X(islave) - xtmp(islave)
1036  enddo
1037  deallocate(xtmp)
1038  end subroutine comp_x_slave_paracon
1039 
1040  subroutine comp_lag(hecMAT, iwS, wSU, num_lagrange)
1041  implicit none
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(:)
1049 
1050  ndof=hecmat%ndof
1051  ndof2=ndof*ndof
1052 
1053  xlag=>hecmat%X(hecmat%NP*ndof+1:hecmat%NP*ndof+num_lagrange)
1054 
1055  do ilag=1,num_lagrange
1056  is=iws(ilag)
1057  i=(is-1)/ndof+1
1058  idof=mod(is-1, ndof)+1
1059  xlag(ilag)=hecmat%B(is)
1060  !lower
1061  js=hecmat%indexL(i-1)+1
1062  je=hecmat%indexL(i)
1063  do j=js,je
1064  jj=hecmat%itemL(j)
1065  ij0=(j-1)*ndof2+(idof-1)*ndof
1066  j0=(jj-1)*ndof
1067  do jdof=1,ndof
1068  xlag(ilag)=xlag(ilag)-hecmat%AL(ij0+jdof)*hecmat%X(j0+jdof)
1069  enddo
1070  enddo
1071  !diag
1072  ij0=(i-1)*ndof2+(idof-1)*ndof
1073  j0=(i-1)*ndof
1074  do jdof=1,ndof
1075  xlag(ilag)=xlag(ilag)-hecmat%D(ij0+jdof)*hecmat%X(j0+jdof)
1076  enddo
1077  !upper
1078  js=hecmat%indexU(i-1)+1
1079  je=hecmat%indexU(i)
1080  do j=js,je
1081  jj=hecmat%itemU(j)
1082  ij0=(j-1)*ndof2+(idof-1)*ndof
1083  j0=(jj-1)*ndof
1084  do jdof=1,ndof
1085  xlag(ilag)=xlag(ilag)-hecmat%AU(ij0+jdof)*hecmat%X(j0+jdof)
1086  enddo
1087  enddo
1088  xlag(ilag)=-wsu(ilag)*xlag(ilag)
1089  enddo
1090  end subroutine comp_lag
1091 
1092  subroutine comp_lag_paracon(hecMESH, hecMAT, Btot, hecMESHtmp, hecTKT, BKmat, &
1093  n_slave, slaves, num_lagrange, iwS, wSU, conCOMM)
1094  implicit none
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(:)
1107  ndof=hecmat%ndof
1108  npndof = hecmat%NP * ndof
1109  nndof = hecmat%N * ndof
1110  npndof_new = hectkt%NP * ndof
1111  !! <SUMMARY>
1112  !! {lag} = [Bs^-T] ( {fs} - [Ksp Kss] {u} )
1113  !!
1114  ! 1. {Btmp} = [BKmat] {X}
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)
1119  !
1120  ! 2. {Btmp_s} = {fs} - {Btmp_s}
1121  !do i = 1, nndof
1122  ! Btmp(i) = Btot(i) - Btmp(i)
1123  !enddo
1124  do i = 1, n_slave
1125  islave = slaves(i)
1126  btmp(islave) = btot(islave) - btmp(islave)
1127  enddo
1128  !
1129  ! 3. send internal contact dof => recv external contact dof
1130  !call hecmw_update_3_R(hecMESH, Btmp, hecMAT%NP)
1131  ! Btmp(nndof+1:npndof) = 0.0d0
1132  call fstr_contact_comm_bcast_r(concomm, btmp)
1133  !
1134  ! 4. {lag} = - [-Bs^-T] {Btmp_s}
1135  xlag => hecmat%X(npndof+1:npndof+num_lagrange)
1136  do ilag = 1, num_lagrange
1137  islave = iws(ilag)
1138  !xlag(ilag)=-wSU(ilag)*(Btot(islave) - Btmp(islave))
1139  xlag(ilag)=-wsu(ilag)*btmp(islave)
1140  enddo
1141  deallocate(btmp)
1142  end subroutine comp_lag_paracon
1143 
1144  subroutine comp_lag_paracon2(hecMESH, hecMAT, conMAT, num_lagrange, iwS, wSU, conCOMM)
1145  implicit none
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(:)
1156  ndof=hecmat%ndof
1157  ndof2=ndof*ndof
1158  npndof = hecmat%NP * ndof
1159  nndof = hecmat%N * ndof
1160  !! <SUMMARY>
1161  !! {lag} = [Bs^-T] ( {fs} - [Ksp Kss] {u} )
1162  !!
1163  !
1164  ! {fs} = {hecMAT%B} + {conMAT%B}
1165  ! [K] {u} = [hecMAT] {u} + [conMAT] {u}
1166  !
1167  ! {Btmp} = {hecMAT%B} - [hecMAT] {u}
1168  allocate(btmp(npndof))
1169  call hecmw_matresid(hecmesh, hecmat, hecmat%X, hecmat%B, btmp)
1170  call fstr_contact_comm_bcast_r(concomm, btmp)
1171  !
1172  ! {Btmp} = {Btmp} + {conMAT%B} - [conMAT] {u}
1173  do ilag = 1, num_lagrange
1174  i = iws(ilag)
1175  btmp(i) = btmp(i) + conmat%B(i)
1176  irow = (i + ndof - 1) / ndof
1177  idof = i - ndof*(irow-1)
1178  ! lower
1179  js = conmat%indexL(irow-1)+1
1180  je = conmat%indexL(irow)
1181  do j = js, je
1182  jcol = conmat%itemL(j)
1183  do jdof = 1, ndof
1184  btmp(i) = btmp(i) - conmat%AL(ndof2*(j-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(jcol-1)+jdof)
1185  enddo
1186  enddo
1187  ! diag
1188  do jdof = 1, ndof
1189  btmp(i) = btmp(i) - conmat%D(ndof2*(irow-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(irow-1)+jdof)
1190  enddo
1191  ! upper
1192  js = conmat%indexU(irow-1)+1
1193  je = conmat%indexU(irow)
1194  do j = js, je
1195  jcol = conmat%itemU(j)
1196  do jdof = 1, ndof
1197  btmp(i) = btmp(i) - conmat%AU(ndof2*(j-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(jcol-1)+jdof)
1198  enddo
1199  enddo
1200  enddo
1201  !
1202  ! {lag} = - [-Bs^-T] {Btmp_s}
1203  xlag => hecmat%X(npndof+1:npndof+num_lagrange)
1204  do ilag = 1, num_lagrange
1205  islave = iws(ilag)
1206  !xlag(ilag)=-wSU(ilag)*(conMAT%B(islave) - Btmp(islave))
1207  xlag(ilag)=-wsu(ilag)*btmp(islave)
1208  enddo
1209  deallocate(btmp)
1210  end subroutine comp_lag_paracon2
1211 
1212  subroutine check_solution(hecMESH, hecMAT, fstrMAT, Btot, hecMESHtmp, hecTKT, BKmat, &
1213  conCOMM, n_slave, slaves)
1214  implicit none
1215  type (hecmwst_local_mesh), intent(in) :: hecmesh, hecmeshtmp
1216  type (hecmwst_matrix), intent(inout) :: hecmat, hectkt
1217  type (fstrst_matrix_contact_lagrange) , intent(in) :: fstrmat
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
1229  ndof = hecmat%NDOF
1230  nndof = hecmat%N * ndof
1231  npndof = hecmat%NP * ndof
1232  num_lagrange = fstrmat%num_lagrange
1233  !
1234  allocate(r(npndof + num_lagrange))
1235  r(:) = 0.0d0
1236  allocate(btmp(npndof))
1237  !
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)
1241  !
1242  !! {r} = {b} - [K] {x} - [Bt] {lag}
1243  !! {rlag} = {c} - [B] {x}
1244  !
1245  ! {r} = {b} - [K] {x}
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)
1250  !
1251  ! {r} = {r} - [Bt] {lag}
1252  btmp(:) = 0.0d0
1253  if (fstrmat%num_lagrange > 0) then
1254  do i = 1, hecmat%NP
1255  ls = fstrmat%indexU_lagrange(i-1)+1
1256  le = fstrmat%indexU_lagrange(i)
1257  do l = ls, le
1258  j = fstrmat%itemU_lagrange(l)
1259  do idof = 1, ndof
1260  btmp(ndof*(i-1)+idof) = btmp(ndof*(i-1)+idof) + fstrmat%AU_lagrange(ndof*(l-1)+idof) * xlag(j)
1261  enddo
1262  enddo
1263  enddo
1264  endif
1265  call fstr_contact_comm_reduce_r(concomm, btmp, hecmw_sum)
1266  r(1:nndof) = r(1:nndof) - btmp(1:nndof)
1267  !
1268  ! {rlag} = {c} - [B] {x}
1269  call hecmw_update_3_r(hecmesh, hecmat%X, hecmat%NP)
1270  do i = 1, num_lagrange
1271  rlag(i) = blag(i)
1272  ls = fstrmat%indexL_lagrange(i-1)+1
1273  le = fstrmat%indexL_lagrange(i)
1274  do l = ls, le
1275  j = fstrmat%itemL_lagrange(l)
1276  do jdof = 1, ndof
1277  rlag(i) = rlag(i) - fstrmat%AL_lagrange(ndof*(l-1)+jdof) * hecmat%X(ndof*(j-1)+jdof)
1278  enddo
1279  enddo
1280  enddo
1281  !
1282  ! residual in original system
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)
1291  endif
1292  if (n_slave > 0) then
1293  write(1000+myrank,*) 'R(slave):',slaves(:)
1294  write(1000+myrank,*) r(slaves(:))
1295  endif
1296  endif
1297  !
1298  call hecmw_innerproduct_r(hecmesh, ndof, r, r, rnrm2)
1299  call hecmw_innerproduct_r(hecmesh, ndof, btot, btot, bnrm2)
1300  rlagnrm2 = 0.0d0
1301  do i = 1, num_lagrange
1302  rlagnrm2 = rlagnrm2 + rlag(i)*rlag(i)
1303  enddo
1304  call hecmw_allreduce_r1(hecmesh, rlagnrm2, hecmw_sum)
1305  blagnrm2 = 0.0d0
1306  do i = 1, num_lagrange
1307  blagnrm2 = blagnrm2 + blag(i)*blag(i)
1308  enddo
1309  call hecmw_allreduce_r1(hecmesh, blagnrm2, hecmw_sum)
1310  !
1311  if (myrank == 0) then
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)
1314  endif
1315  end subroutine check_solution
1316 
1317  subroutine check_solution2(hecMESH, hecMAT, conMAT, fstrMAT, n_contact_dof, contact_dofs, &
1318  conCOMM, n_slave, slaves)
1319  implicit none
1320  type (hecmwst_local_mesh), intent(in) :: hecmesh
1321  type (hecmwst_matrix), intent(inout) :: hecmat
1322  type (hecmwst_matrix), intent(in) :: conmat
1323  type (fstrst_matrix_contact_lagrange) , intent(in) :: fstrmat
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
1335  ndof = hecmat%NDOF
1336  ndof2 = ndof*ndof
1337  nndof = hecmat%N * ndof
1338  npndof = hecmat%NP * ndof
1339  num_lagrange = fstrmat%num_lagrange
1340  !
1341  allocate(r(npndof + num_lagrange))
1342  r(:) = 0.0d0
1343  allocate(r_con(npndof))
1344  r_con(:) = 0.0d0
1345  !
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)
1349  !
1350  !! <SUMMARY>
1351  !! {r} = {b} - [K] {x} - [Bt] {lag}
1352  !! {rlag} = {c} - [B] {x}
1353  !
1354  ! 1. {r} = {r_org} + {r_con}
1355  ! {r_org} = {b_org} - [hecMAT] {x}
1356  ! {r_con} = {b_con} - [conMAT] {x} - [Bt] {lag}
1357  !
1358  ! 1.1 {r_org}
1359  ! {r} = {b} - [K] {x}
1360  call hecmw_matresid(hecmesh, hecmat, hecmat%X, hecmat%B, r)
1361  !
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(:))
1370  endif
1371  endif
1372  !
1373  ! 1.2 {r_con}
1374  ! 1.2.1 {r_con} = {b_con} - [conMAT]{x}
1375  !call hecmw_update_3_R(hecMESH, hecMAT%X, hecMAT%NP) ! X is already updated
1376  do i = 1, npndof
1377  r_con(i) = conmat%B(i)
1378  enddo
1379  do icon = 1, n_contact_dof
1380  i = contact_dofs(icon)
1381  irow = (i + ndof - 1) / ndof
1382  idof = i - ndof*(irow-1)
1383  ! lower
1384  js = conmat%indexL(irow-1)+1
1385  je = conmat%indexL(irow)
1386  do j = js, je
1387  jcol = conmat%itemL(j)
1388  do jdof = 1, ndof
1389  r_con(i) = r_con(i) - conmat%AL(ndof2*(j-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(jcol-1)+jdof)
1390  enddo
1391  enddo
1392  ! diag
1393  do jdof = 1, ndof
1394  r_con(i) = r_con(i) - conmat%D(ndof2*(irow-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(irow-1)+jdof)
1395  enddo
1396  ! upper
1397  js = conmat%indexU(irow-1)+1
1398  je = conmat%indexU(irow)
1399  do j = js, je
1400  jcol = conmat%itemU(j)
1401  do jdof = 1, ndof
1402  r_con(i) = r_con(i) - conmat%AU(ndof2*(j-1)+ndof*(idof-1)+jdof) * hecmat%X(ndof*(jcol-1)+jdof)
1403  enddo
1404  enddo
1405  enddo
1406  !
1407  ! 1.2.2 {r_con} = {r_con} - [Bt] {lag}
1408  if (fstrmat%num_lagrange > 0) then
1409  do i = 1, hecmat%NP
1410  ls = fstrmat%indexU_lagrange(i-1)+1
1411  le = fstrmat%indexU_lagrange(i)
1412  do l = ls, le
1413  j = fstrmat%itemU_lagrange(l)
1414  do idof = 1, ndof
1415  r_con(ndof*(i-1)+idof) = r_con(ndof*(i-1)+idof) - fstrmat%AU_lagrange(ndof*(l-1)+idof) * xlag(j)
1416  enddo
1417  enddo
1418  enddo
1419  endif
1420  !
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(:))
1431  endif
1432  endif
1433  !
1434  ! 1.2.3 send external part of {r_con}
1435  call fstr_contact_comm_reduce_r(concomm, r_con, hecmw_sum)
1436  !
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(:))
1445  endif
1446  endif
1447  !
1448  ! 1.3 {r} = {r_org} + {r_con}
1449  do i = 1,nndof
1450  r(i) = r(i) + r_con(i)
1451  enddo
1452  !
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(:))
1461  endif
1462  endif
1463  !
1464  ! 2. {rlag} = {c} - [B] {x}
1465  !call hecmw_update_3_R(hecMESH, hecMAT%X, hecMAT%NP) ! X is already updated
1466  do i = 1, num_lagrange
1467  rlag(i) = blag(i)
1468  ls = fstrmat%indexL_lagrange(i-1)+1
1469  le = fstrmat%indexL_lagrange(i)
1470  do l = ls, le
1471  j = fstrmat%itemL_lagrange(l)
1472  do jdof = 1, ndof
1473  rlag(i) = rlag(i) - fstrmat%AL_lagrange(ndof*(l-1)+jdof) * hecmat%X(ndof*(j-1)+jdof)
1474  enddo
1475  enddo
1476  enddo
1477  !
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)
1483  endif
1484  endif
1485  !
1486  call hecmw_innerproduct_r(hecmesh, ndof, r, r, rnrm2)
1487  rlagnrm2 = 0.0d0
1488  do i = 1, num_lagrange
1489  rlagnrm2 = rlagnrm2 + rlag(i)*rlag(i)
1490  enddo
1491  call hecmw_allreduce_r1(hecmesh, rlagnrm2, hecmw_sum)
1492  !
1493  if (myrank == 0) write(0,*) 'INFO: resid(x,lag,tot)',sqrt(rnrm2),sqrt(rlagnrm2),sqrt(rnrm2+rlagnrm2)
1494  end subroutine check_solution2
1495 
1496  subroutine mark_slave_dof(BTmat, mark, n_slave, slaves)
1497  implicit none
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
1503  ndof = btmat%ndof
1504  ndof2 = ndof*ndof
1505  mark(:) = 0
1506  do irow = 1, btmat%nr
1507  js = btmat%index(irow-1)+1
1508  je = btmat%index(irow)
1509  do j = js, je
1510  jcol = btmat%item(j)
1511  do idof = 1, ndof
1512  if (mark(ndof*(irow-1)+idof) == 1) cycle
1513  do jdof = 1, ndof
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
1517  exit
1518  endif
1519  enddo
1520  enddo
1521  enddo
1522  enddo
1523  n_slave = 0
1524  do i = 1, btmat%nr * ndof
1525  if (mark(i) /= 0) n_slave = n_slave + 1
1526  enddo
1527  if (debug >= 2) write(0,*) ' DEBUG2: n_slave',n_slave
1528  allocate(slaves(n_slave))
1529  n_slave = 0
1530  do i = 1, btmat%nr * ndof
1531  if (mark(i) /= 0) then
1532  n_slave = n_slave + 1
1533  slaves(n_slave) = i
1534  endif
1535  enddo
1536  if (debug >= 2) write(0,*) ' DEBUG2: slaves',slaves(:)
1537  end subroutine mark_slave_dof
1538 
1539  subroutine place_one_on_diag_of_unmarked_dof(BTmat, mark)
1540  implicit none
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
1545  ndof = btmat%ndof
1546  ndof2 = ndof*ndof
1547  ! Imat: unit matrix except for slave dofs
1548  imat%nr = btmat%nr
1549  imat%nc = btmat%nc
1550  imat%nnz = imat%nr
1551  imat%ndof = ndof
1552  allocate(imat%index(0:imat%nr))
1553  allocate(imat%item(imat%nnz))
1554  imat%index(0) = 0
1555  do i = 1, imat%nr
1556  imat%index(i) = i
1557  imat%item(i) = i
1558  enddo
1559  allocate(imat%A(ndof2 * imat%nnz))
1560  imat%A(:) = 0.0d0
1561  n_slave = 0
1562  n_other = 0
1563  do irow = 1, imat%nr
1564  do idof = 1, ndof
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
1568  else
1569  n_slave = n_slave + 1
1570  endif
1571  enddo
1572  enddo
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)
1577  btmat%nr = wmat%nr
1578  btmat%nc = wmat%nc
1579  btmat%nnz = wmat%nnz
1580  btmat%ndof = wmat%ndof
1581  btmat%index => wmat%index
1582  btmat%item => wmat%item
1583  btmat%A => wmat%A
1584  end subroutine place_one_on_diag_of_unmarked_dof
1585 
1586  subroutine place_num_on_diag_of_marked_dof(BTtKTmat, num, mark)
1587  implicit none
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
1594  ndof2 = ndof*ndof
1595  do irow = 1, bttktmat%nr
1596  js = bttktmat%index(irow-1)+1
1597  je = bttktmat%index(irow)
1598  do j = js, je
1599  jcol = bttktmat%item(j)
1600  if (irow /= jcol) cycle
1601  do idof = 1, ndof
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
1606  else
1607  if (bttktmat%A(ndof2*(j-1)+ndof*(idof-1)+idof) == 0.0d0) then
1608  !write(0,*) 'irow,idof',irow,idof
1609  n_error = n_error+1
1610  endif
1611  endif
1612  enddo
1613  enddo
1614  enddo
1615  if (n_error > 0) then
1616  write(0,*) 'n_error',n_error
1617  stop 'ERROR: zero diag on non-slave dof of BTtKTmat'
1618  endif
1619  end subroutine place_num_on_diag_of_marked_dof
1620 
1621  subroutine update_comm_table(hecMESH, BTmat, hecMESHtmp, BT_all)
1622  implicit none
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
1648  !!! PREPARATION FOR COMM_TABLE UPDATE
1649  call copy_mesh(hecmesh, hecmeshtmp)
1650  allocate(bt_exp(hecmesh%n_neighbor_pe))
1651  call extract_bt_exp(btmat, hecmesh, bt_exp)
1652 
1653  !!! UPDATE COMMUNICATION TABLE for Parallel Computation
1654  allocate(statuses(hecmw_status_size,2*hecmesh%n_neighbor_pe))
1655  allocate(requests(2*hecmesh%n_neighbor_pe))
1656 
1657  allocate(old_export_item(hecmesh%export_index(hecmesh%n_neighbor_pe)))
1658 
1659  n_send = 0
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'
1665  n_oldexp = 0
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)
1670  ! check location of old export nodes in original export list
1671  old_export => old_export_item(idx_0+1:idx_n)
1672  do k = 1,n_orgexp
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
1677  end if
1678  end do
1679  if (debug >= 1) write(0,*) 'DEBUG: making old_export done'
1680  ! gather new export nodes at the end of current export list
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'
1683  ! check consistency
1684  if (n_curexp /= n_orgexp - n_oldexp + n_newexp) &
1685  stop 'ERROR: unknown error(num of export nodes)' !!! ASSERTION
1686  ! make item_exp from item of BT_exp by converting column id to place in cur_export
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'
1689  ! add current export list to commtable
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))
1695  ! send current export info to neighbor pe
1696  sendbuf(1) = n_oldexp
1697  sendbuf(2) = n_newexp
1698  tag = 1001
1699  call hecmw_isend_int(sendbuf, 2, irank, tag, &
1700  hecmesh%MPI_COMM, requests(idom))
1701  if (n_oldexp > 0) then
1702  n_send = n_send + 1
1703  tag = 1002
1704  call hecmw_isend_int(old_export, n_oldexp, irank, tag, &
1705  hecmesh%MPI_COMM, requests(hecmesh%n_neighbor_pe+n_send))
1706  end if
1707  end do
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)
1711  ! receive current import info from neighbor pe
1712  tag = 1001
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))
1719  tag = 1002
1720  call hecmw_recv_int(old_import, n_oldimp, irank, tag, &
1721  hecmesh%MPI_COMM, statuses(:,1))
1722  end if
1723  !
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,*) 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)
1740 import => hecmeshtmp%import_item(hecmeshtmp%import_index(idom-1)+1:hecmeshtmp%import_index(idom)) !cur_
1741  end do
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)
1745 
1746  !!! Send BT_exp & Recv BT_imp; nnz and index
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
1751  tag = 1003
1752  call HECMW_ISEND_INT(sendbuf, 2, irank, tag, &
1753  hecMESH%MPI_COMM, requests(2*idom-1))
1754  tag = 1004
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))
1757  end do
1758  if (DEBUG >= 1) write(0,*) 'debug: isend bt_exp (nnz and index) done'
1759  BT_imp%nr = 0
1760  BT_imp%nc = hecMESHtmp%n_node - hecMESHtmp%nn_internal
1761  BT_imp%nnz = 0
1762  allocate(BT_imp%index(0:hecMESH%import_index(hecMESH%n_neighbor_pe)))
1763  BT_imp%index(0) = 0
1764  allocate(nnz_imp(hecMESH%n_neighbor_pe))
1765  do idom = 1,hecMESH%n_neighbor_pe
1766  irank = hecMESH%neighbor_pe(idom)
1767  tag = 1003
1768  call HECMW_RECV_INT(recvbuf, 2, irank, tag, &
1769  hecMESH%MPI_COMM, statuses(:,1))
1770  nr_imp = recvbuf(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' !!! ASSERTION
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))
1779  tag = 1004
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 !!! ASSERTION
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)
1786  stop
1787  endif
1788  do j = 1, nr_imp
1789  jj = hecMESH%import_item(idx_0+j) - hecMESH%nn_internal
1790  BT_imp%index(jj) = index_imp(j) - index_imp(j-1)
1791  end do
1792  deallocate(index_imp)
1793  end do
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)
1797  end do
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' !!! ASSERTION
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)
1803 
1804  !!! Send BT_exp & Recv BT_imp; item and val
1805  do idom = 1,hecMESH%n_neighbor_pe
1806  irank = hecMESH%neighbor_pe(idom)
1807  tag = 1005
1808  call HECMW_Isend_INT(BT_exp(idom)%item, BT_exp(idom)%nnz, &
1809  irank, tag, hecMESH%MPI_COMM, requests(2*idom-1))
1810  tag = 1006
1811  call HECMW_Isend_R(BT_exp(idom)%A, BT_exp(idom)%nnz * ndof2, &
1812  irank, tag, hecMESH%MPI_COMM, requests(2*idom))
1813  end do
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)))
1820  tag = 1005
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))
1824  tag = 1006
1825  call HECMW_Recv_R(val_imp, nnz_imp(idom) * ndof2, &
1826  irank, tag, hecMESH%MPI_COMM, statuses(:,1))
1827 
1828  ! convert column id of item_imp() to local id refering cur_import(:)
1829  idx_0_tmp = hecMESHtmp%import_index(idom-1)
1830  idx_n_tmp = hecMESHtmp%import_index(idom)
1831 import => hecmeshtmp%import_item(idx_0_tmp+1:idx_n_tmp) cur_
1832  n_curimp = idx_n_tmp - idx_0_tmp
1833  n_orgimp = idx_n - idx_0
1834  cnt = 0
1835  do j = 1, n_orgimp
1836  jj = hecMESH%import_item(idx_0+j) - hecMESH%nn_internal
1837  ks = BT_imp%index(jj-1)
1838  ke = BT_imp%index(jj)
1839  do k = ks+1, ke
1840  cnt = cnt + 1
1841  iimp = item_imp(cnt)
1842 .or. if (iimp <= 0 n_curimp < iimp) &
1843  stop 'error: received column id out of range' !!! ASSERTION
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)
1846  end do
1847  end do
1848  deallocate(item_imp, val_imp)
1849  end do
1850  deallocate(nnz_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)
1853 
1854  deallocate(statuses)
1855  deallocate(requests)
1856 
1857  ! make BT_all by combining BTmat and BT_exp
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))
1865  BT_all%index(0) = 0
1866  do i = 1, BTmat%nr
1867  BT_all%index(i) = BTmat%index(i)
1868  end do
1869  do i = 1, BT_imp%nr
1870  BT_all%index(BTmat%nr+i) = BT_all%index(BTmat%nr+i-1) + &
1871  BT_imp%index(i) - BT_imp%index(i-1)
1872  end do
1873  do i = 1, BTmat%nnz
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)
1876  end do
1877  do i = 1, BT_imp%nnz
1878  ii = BTmat%nnz + i
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)
1881  end do
1882  if (DEBUG >= 1) write(0,*) 'debug: making bt_all done'
1883 
1884  ! free BT_exp(:)
1885  do idom=1,hecMESH%n_neighbor_pe
1886  call hecmw_localmat_free(BT_exp(idom))
1887  end do
1888  deallocate(BT_exp)
1889  end subroutine update_comm_table
1890 
1891  subroutine copy_mesh(src, dst, fg_paracon)
1892  implicit none
1893  type (hecmwST_local_mesh), intent(in) :: src
1894  type (hecmwST_local_mesh), intent(out) :: dst
1895  logical, intent(in), optional :: fg_paracon
1896  dst%zero = src%zero
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 .and. if (present(fg_paracon) 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)
1922  else
1923  dst%import_index(:)= 0
1924  dst%export_index(:)= 0
1925  dst%import_item => null()
1926  dst%export_item => null()
1927  endif
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(:)
1932  !dst%mpc = src%mpc
1933  ! MPC is already set outside of here
1934  dst%mpc%n_mpc = 0
1935  dst%node => src%node
1936  end subroutine copy_mesh
1937 
1938  subroutine free_mesh(hecMESH, fg_paracon)
1939  implicit none
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 .and. if (present(fg_paracon) fg_paracon) then
1946  deallocate(hecMESH%import_item)
1947  deallocate(hecMESH%export_item)
1948  deallocate(hecMESH%global_node_ID)
1949  else
1950  if (associated(hecMESH%import_item)) deallocate(hecMESH%import_item)
1951  if (associated(hecMESH%export_item)) deallocate(hecMESH%export_item)
1952  endif
1953  deallocate(hecMESH%node_ID)
1954  deallocate(hecMESH%elem_type_item)
1955  !hecMESH%node => null()
1956  end subroutine free_mesh
1957 
1958  subroutine extract_BT_exp(BTmat, hecMESH, BT_exp)
1959  implicit none
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
1970  BT_exp(i)%nnz = 0
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
1979  end do
1980  allocate(BT_exp(i)%item(BT_exp(i)%nnz))
1981  allocate(BT_exp(i)%A(BT_exp(i)%nnz * ndof2))
1982  n = 0
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)
1986  n = n + 1
1987  !write(0,*) j, jrow, k, n
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)
1990  end do
1991  end do
1992  end do
1993  end subroutine extract_BT_exp
1994 
1995  subroutine extract_cols(BT_exp, cur_export, n_curexp)
1996  implicit none
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
2000  ! write(0,*) 'bt_exp%item(1:',BT_exp%nnz,')'
2001  ! write(0,*) BT_exp%item(1:BT_exp%nnz)
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)
2005  ! write(0,*) 'cur_export(1:',n_curexp,')'
2006  ! write(0,*) cur_export(1:n_curexp)
2007  end subroutine extract_cols
2008 
2009  subroutine reorder_current_export(cur_export, n_curexp, org_export, n_orgexp, n_newexp, nn_internal)
2010  implicit none
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
2019  n_newexp = 0
2020  allocate(new_export(n_curexp))
2021  do j = 1,n_curexp
2022  jnod = cur_export(j)
2023  if (jnod > nn_internal) &
2024  stop 'error: unknown error (jnod)' !!! ASSERTION
2025 .not. if ( is_included(org_export, n_orgexp, jnod)) then
2026  n_newexp = n_newexp + 1
2027  new_export(n_newexp) = jnod
2028  !write(0,*) 'found new export', jnod
2029  else if (n_newexp > 0) then
2030  cur_export(j - n_newexp) = jnod
2031  end if
2032  end do
2033  do j = 1,n_newexp
2034  cur_export(n_curexp - n_newexp + j) = new_export(j)
2035  end do
2036  deallocate(new_export)
2037  ! write(0,*) 'reordered cur_export(1:',n_curexp,')'
2038  ! write(0,*) cur_export(1:n_curexp)
2039  end subroutine reorder_current_export
2040 
2041  subroutine convert_BT_exp_col_id(BT_exp, cur_export, n_curexp)
2042  implicit none
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
2047  logical :: found
2048  ! make item_exp from item of BT_exp by converting column id to place in cur_export
2049  do i = 1, BT_exp%nnz
2050  icol = BT_exp%item(i)
2051  found = .false.
2052  do j = 1, n_curexp
2053  if (icol == cur_export(j)) then
2054  BT_exp%item(i) = j
2055  found = .true.
2056  exit
2057  end if
2058  end do
2059 .not. if ( found) then
2060  write(0,*) icol
2061  stop 'error: unknown error (item not found in cur_export)' !!! ASSERTION
2062  end if
2063  end do
2064  end subroutine convert_BT_exp_col_id
2065 
2066  subroutine append_commtable(n, index, item, idom, cur, ncur)
2067  implicit none
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(:)
2075  norg = index(n)
2076  allocate(tmp_item(norg))
2077  if (norg > 0) then
2078  tmp_item(:) = item(:)
2079  if (associated(item)) deallocate(item)
2080  end if
2081  allocate(item(norg + ncur))
2082  do j = idom,n
2083  index(j) = index(j) + ncur
2084  end do
2085  do j = 1,tmp_index(idom)
2086  item(j) = tmp_item(j)
2087  end do
2088  do j = 1,ncur
2089  item(tmp_index(idom)+j) = cur(j)
2090  end do
2091  do j = tmp_index(idom)+1,tmp_index(n)
2092  item(j+ncur) = tmp_item(j)
2093  end do
2094  deallocate(tmp_index, tmp_item)
2095  end subroutine append_commtable
2096 
2097  subroutine append_nodes(hecMESHtmp, n_newimp, i0)
2098  implicit none
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
2105 
2106  subroutine make_cur_import(org_import, n_orgimp, old_import, n_oldimp, &
2107  n_newimp, i0, cur_import)
2108  implicit none
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(:)
2112  ! integer(kind=kint), intent(out) :: n_curimp
2113  integer(kind=kint) :: ndel, i, j
2114  ndel = 0
2115  i = 1
2116 .and. do while (i <= n_orgimp ndel < n_oldimp)
2117  if (org_import(i) == old_import(ndel+1)) then
2118  ndel = ndel + 1
2119  else
2120  cur_import(i-ndel) = org_import(i)
2121  endif
2122  i = i + 1
2123  enddo
2124  if (ndel /= n_oldimp) stop 'error: unknown error (ndel)' !!! ASSERTION
2125  do j = i, n_orgimp
2126  cur_import(j-ndel) = org_import(j)
2127  enddo
2128  i = n_orgimp - ndel
2129  do j = 1, n_newimp
2130  cur_import(i + j) = i0+j
2131  end do
2132 import end subroutine make_cur_
2133 
2134  recursive subroutine quick_sort(array, id1, id2)
2135  implicit none
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)
2142  left = id1
2143  right = id2
2144  do
2145  do while (array(left) < pivot)
2146  left = left + 1
2147  end do
2148  do while (pivot < array(right))
2149  right = right - 1
2150  end do
2151  if (left >= right) exit
2152  tmp = array(left)
2153  array(left) = array(right)
2154  array(right) = tmp
2155  left = left + 1
2156  right = right - 1
2157  end do
2158  if (id1 < left-1) call quick_sort(array, id1, left-1)
2159  if (right+1 < id2) call quick_sort(array, right+1, id2)
2160  return
2161  end subroutine quick_sort
2162 
2163  subroutine unique(array, len, newlen)
2164  implicit none
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
2169  ndup = 0
2170  do i=2,len
2171  if (array(i) == array(i - 1 - ndup)) then
2172  ndup = ndup + 1
2173  else if (ndup > 0) then
2174  array(i - ndup) = array(i)
2175  endif
2176  end do
2177  newlen = len - ndup
2178  end subroutine unique
2179 
2180  function is_included(array, len, ival)
2181  implicit none
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.
2188  do i=1,len
2189  if (array(i) == ival) then
2190  is_included = .true.
2191  exit
2192  end if
2193  end do
2194  end function is_included
2195 
2196  !!
2197  !! Solve without elimination of Lagrange-multipliers
2198  !!
2199 
2200  subroutine solve_no_eliminate(hecMESH,hecMAT,fstrMAT)
2201  implicit none
2202  type (hecmwST_local_mesh), intent(in) :: hecMESH
2203  type (hecmwST_matrix ), intent(inout) :: hecMAT
2204  type (fstrST_matrix_contact_lagrange), intent(inout) :: fstrMAT !< type fstrST_matrix_contact_lagrange
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
2211 
2212  t1 = hecmw_wtime()
2213  if (DEBUG >= 1) write(0,*) 'debug: solve_no_eliminate, start', hecmw_wtime()-t1
2214 
2215  call hecmw_mat_init(hecMATLag)
2216 
2217  ndof = hecMAT%NDOF
2218  ndof2 = ndof*ndof
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
2223  !write(0,*) 'debug: hecmat: ndof,n,np=',hecMAT%NDOF,hecMAT%N,hecMAT%NP
2224  !write(0,*) 'debug: hecmatlag: ndof,n,np=',hecMATLag%NDOF,hecMATLag%N,hecMATLag%NP
2225 
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
2228 
2229  ! Upper: count num of blocks
2230  allocate(iwUr(hecMAT%N))
2231  allocate(iwUc(nb_lag))
2232  iwUr = 0
2233  do i = 1, hecMAT%N
2234  iwUc = 0
2235  ls=fstrMAT%indexU_lagrange(i-1)+1
2236  le=fstrMAT%indexU_lagrange(i)
2237  do l=ls,le
2238  j=fstrMAT%itemU_lagrange(l)
2239  jb_lag = (j+2)/3
2240  iwUc(jb_lag) = 1
2241  enddo
2242  do j=1,nb_lag
2243  if (iwUc(j) > 0) iwUr(i) = iwUr(i) + 1
2244  enddo
2245  !if (iwUr(i) > 0) write(0,*) 'debug: iwur(',i,')=',iwUr(i)
2246  enddo
2247 
2248  ! Lower: count num of blocks
2249  allocate(iwLr(nb_lag))
2250  allocate(iwLc(hecMAT%N))
2251  iwLr = 0
2252  do ib_lag = 1, nb_lag
2253  iwLc = 0
2254  i = hecMAT%N + ib_lag
2255  do idof = 1, ndof
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)
2260  do l=ls,le
2261  j=fstrMAT%itemL_lagrange(l)
2262  iwLc(j) = 1
2263  enddo
2264  enddo
2265  do j=1,hecMAT%N
2266  if (iwLc(j) > 0) iwLr(ib_lag) = iwLr(ib_lag) + 1
2267  enddo
2268  !if (iwLr(ib_lag) > 0) write(0,*) 'debug: iwlr(',ib_lag,')=',iwLr(ib_lag)
2269  enddo
2270 
2271  ! Upper: indexU
2272  allocate(hecMATLag%indexU(0:hecMATLag%NP))
2273  hecMATLag%indexU(0) = 0
2274  do i = 1, hecMAT%N
2275  hecMATLag%indexU(i) = hecMATLag%indexU(i-1) + &
2276  (hecMAT%indexU(i) - hecMAT%indexU(i-1)) + iwUr(i)
2277  enddo
2278  do i = hecMAT%N+1, hecMATLag%N
2279  hecMATLag%indexU(i) = hecMATLag%indexU(i-1)
2280  enddo
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))
2284  enddo
2285  hecMATLag%NPU = hecMATLag%indexU(hecMATLag%NP)
2286  !write(0,*) 'debug: hecmatlag%npu=',hecMATLag%NPU
2287 
2288  ! Lower: indexL
2289  allocate(hecMATLag%indexL(0:hecMATLag%NP))
2290  do i = 0, hecMAT%N
2291  hecMATLag%indexL(i) = hecMAT%indexL(i)
2292  enddo
2293  do i = hecMAT%N+1, hecMATLag%N
2294  hecMATLag%indexL(i) = hecMATLag%indexL(i-1) + iwLr(i-hecMAT%N)
2295  enddo
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))
2299  enddo
2300  hecMATLag%NPL = hecMATLag%indexL(hecMATLag%NP)
2301  !write(0,*) 'debug: hecmatlag%npl=',hecMATLag%NPL
2302 
2303  ! Upper: itemU and AU
2304  allocate(hecMATLag%itemU(hecMATLag%NPU))
2305  allocate(hecMATLag%AU(hecMATLag%NPU*ndof2))
2306  hecMATLag%AU = 0.d0
2307  do i = 1, hecMAT%N
2308  idx = hecMATLag%indexU(i-1)+1
2309  ! copy from hecMAT; internal
2310  ls = hecMAT%indexU(i-1)+1
2311  le = hecMAT%indexU(i)
2312  do l=ls,le
2313  ll = hecMAT%itemU(l)
2314  if (ll > hecMAT%N) cycle ! skip external
2315  hecMATLag%itemU(idx) = ll
2316  hecMATLag%AU((idx-1)*ndof2+1:idx*ndof2)=hecMAT%AU((l-1)*ndof2+1:l*ndof2)
2317  idx = idx + 1
2318  enddo
2319  ! Lag. itemU
2320  iwUc = 0
2321  ls=fstrMAT%indexU_lagrange(i-1)+1
2322  le=fstrMAT%indexU_lagrange(i)
2323  do l=ls,le
2324  j=fstrMAT%itemU_lagrange(l)
2325  jb_lag = (j+2)/3
2326  iwUc(jb_lag) = 1
2327  enddo
2328  idx_lag_s = idx
2329  do j=1,nb_lag
2330  if (iwUc(j) > 0) then
2331  hecMATLag%itemU(idx) = hecMAT%N + j
2332  idx = idx + 1
2333  endif
2334  enddo
2335  idx_lag_e = idx - 1
2336  ! Lag. AU
2337  ls=fstrMAT%indexU_lagrange(i-1)+1
2338  le=fstrMAT%indexU_lagrange(i)
2339  do l=ls,le
2340  j=fstrMAT%itemU_lagrange(l)
2341  jb_lag = (j+2)/3
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
2346  !if (hecMATLag%itemU(k) /= hecMAT%N + jb_lag) stop 'error itemu jb_lag'
2347  do idof = 1, ndof
2348  hecMATLag%AU((k-1)*ndof2+(idof-1)*ndof+jdof) = &
2349  fstrMAT%AU_lagrange((l-1)*ndof+idof)
2350  enddo
2351  enddo
2352  enddo
2353  ! copy from hecMAT; externl
2354  ls = hecMAT%indexU(i-1)+1
2355  le = hecMAT%indexU(i)
2356  do l=ls,le
2357  ll = hecMAT%itemU(l)
2358  if (ll <= hecMAT%N) cycle ! skip internal
2359  hecMATLag%itemU(idx) = ll + nb_lag
2360  hecMATLag%AU((idx-1)*ndof2+1:idx*ndof2)=hecMAT%AU((l-1)*ndof2+1:l*ndof2)
2361  idx = idx + 1
2362  enddo
2363  if (idx /= hecMATLag%indexU(i)+1) stop 'error idx indexu'
2364  enddo
2365  do i = hecMAT%N, hecMAT%NP
2366  idx = hecMATLag%indexU(i+nb_lag-1)+1
2367  ls=hecMAT%indexU(i-1)+1
2368  le=hecMAT%indexU(i)
2369  do l=ls,le
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)
2373  idx = idx + 1
2374  enddo
2375  enddo
2376 
2377  ! Lower: itemL and AL
2378  allocate(hecMATLag%itemL(hecMATLag%NPL))
2379  allocate(hecMATLag%AL(hecMATLag%NPL*ndof2))
2380  hecMATLag%AL = 0.d0
2381  do i = 1, hecMAT%indexL(hecMAT%N)
2382  hecMATLag%itemL(i) = hecMAT%itemL(i)
2383  enddo
2384  do i = 1, hecMAT%indexL(hecMAT%N)*ndof2
2385  hecMATLag%AL(i) = hecMAT%AL(i)
2386  enddo
2387  ! Lag. itemL
2388  idx = hecMAT%indexL(hecMAT%N) + 1
2389  do ib_lag = 1, nb_lag
2390  iwLc = 0
2391  i = hecMAT%N + ib_lag
2392  do idof = 1, ndof
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)
2397  do l=ls,le
2398  j=fstrMAT%itemL_lagrange(l)
2399  iwLc(j) = 1
2400  enddo
2401  enddo
2402  idx_lag_s = idx
2403  do j=1,hecMAT%N
2404  if (iwLc(j) > 0) then
2405  hecMATLag%itemL(idx) = j
2406  idx = idx + 1
2407  endif
2408  enddo
2409  idx_lag_e = idx - 1
2410  if (idx /= hecMATLag%indexL(hecMAT%N + ib_lag)+1) then
2411  stop 'error idx indexl'
2412  endif
2413  ! Lag. AL
2414  do idof = 1, ndof
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)
2419  do l=ls,le
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
2424  !if (hecMATLag%itemL(k) /= j) stop 'error iteml j'
2425  do jdof = 1, ndof
2426  hecMATLag%AL((k-1)*ndof2+(idof-1)*ndof+jdof) = &
2427  fstrMAT%AL_lagrange((l-1)*ndof+jdof)
2428  enddo
2429  enddo
2430  enddo
2431  enddo
2432  enddo
2433  do i = hecMAT%N+1, hecMAT%NP
2434  idx = hecMATLag%indexL(i+nb_lag-1)+1
2435  ls=hecMAT%indexL(i-1)+1
2436  le=hecMAT%indexL(i)
2437  do l=ls,le
2438  ll = hecMAT%itemL(l)
2439  if (ll <= hecMAT%N) then
2440  hecMATLag%itemL(idx) = ll
2441  else
2442  hecMATLag%itemL(idx) = ll + nb_lag
2443  endif
2444  hecMATLag%AL((idx-1)*ndof2+1:idx*ndof2)=hecMAT%AL((l-1)*ndof2+1:l*ndof2)
2445  idx = idx + 1
2446  enddo
2447  enddo
2448 
2449  deallocate(iwUr, iwUc, iwLr, iwLc)
2450 
2451  allocate(hecMATLag%D(hecMATLag%NP*ndof2))
2452  hecMATLag%D = 0.d0
2453  ! internal
2454  do i = 1, hecMAT%N*ndof2
2455  hecMATLag%D(i) = hecMAT%D(i)
2456  enddo
2457  ! Lag.
2458  do idof = ndof - ndofextra + 1, ndof
2459  hecMATLag%D((hecMATLag%N-1)*ndof2 + (idof-1)*ndof + idof) = 1.d0
2460  enddo
2461  ! external
2462  do i = hecMAT%N*ndof2+1, hecMAT%NP*ndof2
2463  hecMATLag%D(i + nb_lag*ndof2) = hecMAT%D(i)
2464  enddo
2465 
2466  allocate(hecMATLag%B(hecMATLag%NP*ndof))
2467  allocate(hecMATLag%X(hecMATLag%NP*ndof))
2468  hecMATLag%B = 0.d0
2469  hecMATLag%X = 0.d0
2470 
2471  ! internal
2472  do i = 1, hecMAT%N*ndof
2473  hecMATLag%B(i) = hecMAT%B(i)
2474  enddo
2475  ! Lag.
2476  do i = 1, fstrMAT%num_lagrange
2477  hecMATLag%B(hecMAT%N*ndof + i) = hecMAT%B(hecMAT%NP*ndof + i)
2478  enddo
2479  ! external
2480  do i = hecMAT%N*ndof+1, hecMAT%NP*ndof
2481  hecMATlag%B(i + nb_lag*ndof) = hecMAT%B(i)
2482  enddo
2483 
2484  hecMATLag%Iarray=hecMAT%Iarray
2485  hecMATLag%Rarray=hecMAT%Rarray
2486 
2487  if (DEBUG >= 1) write(0,*) 'debug: made hecmatlag', hecmw_wtime()-t1
2488 
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
2492  enddo
2493  endif
2494 
2495  call hecmw_solve(hecMESH, hecMATLag)
2496 
2497  hecMAT%Iarray=hecMATLag%Iarray
2498  hecMAT%Rarray=hecMATLag%Rarray
2499 
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
2503  enddo
2504  endif
2505 
2506  hecMAT%X = 0.d0
2507  ! internal
2508  do i = 1, hecMAT%N*ndof
2509  hecMAT%X(i) = hecMATLag%X(i)
2510  enddo
2511  ! external
2512  do i = hecMAT%N*ndof+1, hecMAT%NP*ndof
2513  hecMAT%X(i) = hecMATLag%X(i + nb_lag*ndof)
2514  enddo
2515  ! Lag.
2516  do i = 1, fstrMAT%num_lagrange
2517  hecMAT%X(hecMAT%NP*ndof + i) = hecMATLag%X(hecMAT%N*ndof + i)
2518  enddo
2519 
2520  call hecmw_mat_finalize(hecMATLag)
2521 
2522  if (DEBUG >= 1) write(0,*) 'debug: solve_no_eliminate end', hecmw_wtime()-t1
2523  end subroutine solve_no_eliminate
2524 
2525 end module m_solve_LINEQ_iter_contact
This module provides functions of reconstructing.
subroutine hecmw_solve(hecMESH, hecMAT)
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.f90:80
logical paracontactflag
PARALLEL CONTACT FLAG.
Definition: m_fstr.f90:84
This module provides interface of iteratie linear equation solver for contact problems using Lagrange...
subroutine, public solve_lineq_iter_contact_init(hecMESH, hecMAT, fstrMAT, is_sym)
subroutine, public solve_lineq_iter_contact(hecMESH, hecMAT, fstrMAT, istat, conMAT)
Structure for Lagrange multiplier-related part of stiffness matrix (Lagrange multiplier-related matri...