1!------------------------------------------------------------------------------- 2! Copyright (c) 2019 FrontISTR Commons 3! This software is released under the MIT License, see LICENSE.txt 4!------------------------------------------------------------------------------- 5!> This module provides interface of iteratie linear equation solver for 6!! contact problems using Lagrange multiplier. 7module m_solve_LINEQ_iter_contact 8 use m_fstr 9 use fstr_matrix_con_contact 10 use hecmw_solver 11 12 private 13 public :: solve_LINEQ_iter_contact_init 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 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 !< type fstrST_matrix_contact_lagrange 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 !< type fstrST_matrix_contact_lagrange 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 !< type fstrST_matrix_contact_lagrange 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 !< type fstrST_matrix_contact_lagrange 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 !< type fstrST_matrix_contact_lagrange 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 !< type fstrST_matrix_contact_lagrange 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 !< 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 2525end module m_solve_LINEQ_iter_contact 2526