FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
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
21contains
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.
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,*) '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)
1740 !cur_import => hecMESHtmp%import_item(hecMESHtmp%import_index(idom-1)+1:hecMESHtmp%import_index(idom))
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 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
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 if (iimp <= 0 .or. 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 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)
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 if (present(fg_paracon) .and. 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 if (.not. 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 if (.not. 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 do while (i <= n_orgimp .and. 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 end subroutine make_cur_import
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
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
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...