FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
hecmw_matrix_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 !-------------------------------------------------------------------------------
6  use hecmw_util
7  implicit none
8 
9  private
10  public :: hecmw_cmat_init
11  public :: hecmw_cmat_finalize
12  public :: hecmw_cmat_clear
13  public :: hecmw_cmat_add
14  public :: hecmw_cmat_pack
15  public :: hecmw_cmat_multvec_add
16  public :: hecmw_cmat_ass_bc
17  public :: hecmw_cmat_n_newpair
18  public :: hecmw_cmat_lu
19  public :: hecmw_cmat_lu_free
20  public :: hecmw_cmat_substitute
21 
22  integer(kind=kint), parameter :: CMAT_MAX_VAL_INIT = 128
23  integer(kind=kint), parameter :: CMAT_MAX_VAL_GROW = 2
24 
25 contains
26 
27  subroutine hecmw_cmat_init( cmat )
28  type(hecmwst_matrix_contact) :: cmat
29 
30  nullify( cmat%pair )
31  cmat%max_val = 0
32  call hecmw_cmat_clear( cmat)
33  end subroutine hecmw_cmat_init
34 
35  subroutine hecmw_cmat_finalize( cmat )
36  type(hecmwst_matrix_contact) :: cmat
37 
38  if( cmat%max_val > 0 ) deallocate( cmat%pair )
39  call hecmw_cmat_init( cmat )
40  end subroutine hecmw_cmat_finalize
41 
42  subroutine hecmw_cmat_clear( cmat )
43  type(hecmwst_matrix_contact) :: cmat
44 
45  cmat%n_val = 0
46  cmat%checked = .true.
47  cmat%sorted = .true.
48  cmat%max_row = 0
49  cmat%max_col = 0
50  end subroutine hecmw_cmat_clear
51 
52  function compair_pair_by_index( p1, p2 )
53  integer(kind=kint) :: compair_pair_by_index
54  type(hecmwst_index_value_pair) :: p1, p2
55 
56  if( p1%i < p2%i ) then
57  compair_pair_by_index = -1
58  else if( p1%i > p2%i ) then
59  compair_pair_by_index = 1
60  else if( p1%j < p2%j ) then
61  compair_pair_by_index = -1
62  else if( p1%j > p2%j ) then
63  compair_pair_by_index = 1
64  else
65  compair_pair_by_index = 0
66  endif
67  end function compair_pair_by_index
68 
69  subroutine cmat_resize( cmat, newlen )
70  type(hecmwst_matrix_contact) :: cmat
71  integer(kind=kint) :: newlen
72  type(hecmwst_index_value_pair), allocatable :: temp(:)
73 
74  integer(kind=kint) :: i
75 
76  if( newlen <= cmat%max_val ) return
77 
78  if( cmat%max_val > 0 ) then
79  allocate( temp( cmat%n_val ) )
80  do i = 1, cmat%n_val
81  temp(i)%i = cmat%pair(i)%i
82  temp(i)%j = cmat%pair(i)%j
83  temp(i)%val = cmat%pair(i)%val
84  enddo
85  deallocate( cmat%pair )
86  endif
87 
88  allocate( cmat%pair( newlen ) )
89 
90  if( cmat%max_val > 0 ) then
91  do i = 1, cmat%n_val
92  cmat%pair(i)%i = temp(i)%i
93  cmat%pair(i)%j = temp(i)%j
94  cmat%pair(i)%val = temp(i)%val
95  enddo
96  deallocate( temp )
97  endif
98 
99  cmat%max_val = newlen
100  end subroutine cmat_resize
101 
102  subroutine cmat_grow( cmat )
103  type(hecmwst_matrix_contact) :: cmat
104  integer(kind=kint) :: newlen
105 
106  if( cmat%max_val == 0 ) then
107  newlen = cmat_max_val_init
108  else
109  newlen = cmat%max_val * cmat_max_val_grow
110  endif
111  call cmat_resize( cmat, newlen )
112  end subroutine cmat_grow
113 
114  subroutine hecmw_cmat_add( cmat, i, j, val )
115  type(hecmwst_matrix_contact) :: cmat
116  integer(kind=kint) :: i
117  integer(kind=kint) :: j
118  real(kind=kreal), dimension(:,:) :: val
119  integer(kind=kint) :: cmp
120 
121  if( cmat%n_val == cmat%max_val ) then
122  call cmat_grow( cmat )
123  endif
124 
125  cmat%pair( cmat%n_val+1 )%i = i
126  cmat%pair( cmat%n_val+1 )%j = j
127  cmat%pair( cmat%n_val+1 )%val = val
128 
129  if( cmat%n_val > 0 .and. cmat%sorted ) then
130  cmp = compair_pair_by_index( cmat%pair( cmat%n_val ), cmat%pair( cmat%n_val+1 ) )
131  if( cmp > 0 ) then
132  cmat%checked = .false.
133  cmat%sorted = .false.
134  endif
135 
136  if( cmat%checked .and. cmp == 0 ) then
137  cmat%checked = .false.
138  endif
139  endif
140 
141  cmat%n_val = cmat%n_val + 1
142 
143  if( cmat%max_row < i ) cmat%max_row = i
144  if( cmat%max_col < j ) cmat%max_col = j
145  end subroutine hecmw_cmat_add
146 
147  recursive subroutine sort_pair_by_index( pair, first, last )
148  type(hecmwst_index_value_pair), pointer :: pair(:)
149  integer(kind=kint) :: first
150  integer(kind=kint) :: last
151  integer(kind=kint) :: left, right
152  type(hecmwst_index_value_pair) :: pivot, temp
153 
154  pivot = pair( (first + last) / 2 )
155  left = first
156  right = last
157  do
158  do while( compair_pair_by_index( pair(left), pivot ) < 0 )
159  left = left + 1
160  enddo
161  do while( compair_pair_by_index( pivot, pair(right) ) < 0 )
162  right = right - 1
163  enddo
164  if ( left >= right ) exit
165 
166  temp = pair(left)
167  pair(left) = pair(right)
168  pair(right) = temp
169 
170  left = left + 1
171  right = right - 1
172  enddo
173  if( first < left - 1 ) call sort_pair_by_index( pair, first, left - 1 )
174  if( right + 1 < last ) call sort_pair_by_index( pair, right + 1, last )
175  end subroutine sort_pair_by_index
176 
177  subroutine hecmw_cmat_pack( cmat )
178  type(hecmwst_matrix_contact) :: cmat
179  integer(kind=kint) :: i
180  integer(kind=kint) :: n_dup, cmp
181 
182  if( cmat%checked .or. cmat%n_val <= 1 ) return
183 
184  if( .not. cmat%sorted ) then
185  call sort_pair_by_index( cmat%pair, 1, cmat%n_val )
186  cmat%sorted = .true.
187  endif
188 
189  n_dup = 0
190  do i = 2,cmat%n_val
191  cmp = compair_pair_by_index( cmat%pair( i-1-n_dup ), cmat%pair( i ) )
192  if( cmp == 0 ) then
193  n_dup = n_dup + 1
194  cmat%pair( i-n_dup )%val = cmat%pair( i-n_dup )%val + cmat%pair( i )%val
195  else
196  cmat%pair( i-n_dup )%i = cmat%pair( i )%i
197  cmat%pair( i-n_dup )%j = cmat%pair( i )%j
198  cmat%pair( i-n_dup )%val = cmat%pair( i )%val
199  endif
200  enddo
201  cmat%n_val = cmat%n_val - n_dup
202  cmat%checked = .true.
203  end subroutine hecmw_cmat_pack
204 
205  subroutine hecmw_cmat_multvec_add( cmat, X, Y, len )
206  type(hecmwst_matrix_contact) :: cmat
207  real(kind=kreal) :: x(:), y(:)
208  integer(kind=kint) :: len
209  integer(kind=kint) :: i, ii, jj, k, l
210  integer, parameter :: ndof = 3
211 
212  if( cmat%max_row > len .or. cmat%max_col > len ) then
213  write(*,*) 'ERROR: cmat_multvec_add: vector too short'
215  endif
216 
217  do i = 1, cmat%n_val
218  ii = ndof * (cmat%pair(i)%i - 1)
219  jj = ndof * (cmat%pair(i)%j - 1)
220  do l = 1, ndof
221  do k = 1, ndof
222  y(ii + k) = &
223  & y(ii + k) + cmat%pair(i)%val(k, l) * x(jj + l)
224  enddo
225  enddo
226  enddo
227  end subroutine hecmw_cmat_multvec_add
228 
229 
230  subroutine hecmw_cmat_ass_bc(hecMAT, inode, idof, RHS )
231  type (hecmwst_matrix) :: hecmat
232  integer(kind=kint) :: inode, idof
233  real(kind=kreal) :: rhs
234  integer(kind=kint) :: nval, i, cnode
235  integer(kind=kint), parameter :: ndof=3
236 
237  ! NDOF = hecMAT%NDOF
238  if( ndof < idof ) return
239 
240  !-- modify rhs
241  if( rhs /= 0.d0 ) then
242  do nval=1,hecmat%cmat%n_val
243  if( hecmat%cmat%pair(nval)%j /= inode ) cycle
244  cnode = hecmat%cmat%pair(nval)%i
245  if( cnode == inode ) then
246  do i=1,ndof
247  if( i==idof ) cycle
248  hecmat%B(ndof*(cnode-1)+i) = hecmat%B(ndof*(cnode-1)+i) &
249  - hecmat%cmat%pair(nval)%val(i,idof)*rhs
250  enddo
251  else
252  do i=1, ndof
253  hecmat%B(ndof*(cnode-1)+i) = hecmat%B(ndof*(cnode-1)+i) &
254  - hecmat%cmat%pair(nval)%val(i,idof)*rhs
255  enddo
256  endif
257  enddo
258  endif
259 
260  !-- clear cmat
261  do nval=1,hecmat%cmat%n_val
262  if( hecmat%cmat%pair(nval)%i /= inode ) cycle
263 
264  hecmat%cmat%pair(nval)%val(idof,:)=0.d0
265  enddo
266 
267  do nval=1,hecmat%cmat%n_val
268  if( hecmat%cmat%pair(nval)%j /= inode ) cycle
269 
270  hecmat%cmat%pair(nval)%val(:,idof)=0.d0
271  enddo
272 
273  end subroutine hecmw_cmat_ass_bc
274 
275 
276  integer function hecmw_cmat_n_newpair(hecMAT, symm )
277  type(hecmwst_matrix), intent(in) :: hecmat
278  integer, intent(in) :: symm
279  integer(kind=kint) :: nval, ind, jnd, k, m, mnd
280  logical :: isfind
281 
283  do nval=1,hecmat%cmat%n_val
284  ind = hecmat%cmat%pair(nval)%i
285  jnd = hecmat%cmat%pair(nval)%j
286 
287  if( (symm==1) .and. (ind<jnd) ) cycle
288  if( ind==jnd ) cycle
289  isfind = .false.
290  do k=1,hecmat%NP
291  if( ind/=k .and. jnd/=k ) cycle
292  do m= hecmat%indexL(k-1)+1, hecmat%indexL(k)
293  mnd= hecmat%itemL(m)
294  if( (ind==k .and. jnd==mnd) .or. (ind==mnd .and. jnd==k) ) then
295  isfind = .true.
296  endif
297  if ( isfind ) exit
298  enddo
299  if( isfind ) exit
300  enddo
301  if( .not. isfind ) hecmw_cmat_n_newpair = hecmw_cmat_n_newpair+1
302  enddo
303  end function hecmw_cmat_n_newpair
304 
305 
306  subroutine hecmw_cmat_lu( hecMAT )
307  type (hecmwst_matrix) :: hecmat
308  integer(kind=kint) :: i, j, k, l, countcal, countcau
309 
310  allocate (hecmat%indexCL(0:hecmat%NP), hecmat%indexCU(0:hecmat%NP))
311 
312  hecmat%indexCL = 0
313  hecmat%indexCU = 0
314 
315  do i= 1, hecmat%NP
316  do j= 1, hecmat%cmat%n_val
317  if (hecmat%cmat%pair(j)%i == i .and. hecmat%cmat%pair(j)%j < i) then
318  hecmat%indexCL(i) = hecmat%indexCL(i) + 1
319  endif
320  if (hecmat%cmat%pair(j)%i == i .and. hecmat%cmat%pair(j)%j > i) then
321  hecmat%indexCU(i) = hecmat%indexCU(i) + 1
322  endif
323  enddo
324  hecmat%indexCL(i) = hecmat%indexCL(i) + hecmat%indexCL(i-1)
325  hecmat%indexCU(i) = hecmat%indexCU(i) + hecmat%indexCU(i-1)
326  enddo
327 
328  hecmat%NPCL = hecmat%indexCL(hecmat%NP)
329  hecmat%NPCU = hecmat%indexCU(hecmat%NP)
330 
331  allocate (hecmat%itemCL(hecmat%NPCL), hecmat%itemCU(hecmat%NPCU))
332  allocate (hecmat%CAL(9*hecmat%NPCL), hecmat%CAU(9*hecmat%NPCU))
333 
334  countcal = 0
335  countcau = 0
336  do i= 1, hecmat%NP
337  do j= 1, hecmat%cmat%n_val
338  if (hecmat%cmat%pair(j)%i == i .and. hecmat%cmat%pair(j)%j < i) then
339  countcal = countcal + 1
340  hecmat%itemCL(countcal) = hecmat%cmat%pair(j)%j
341  do k= 1, 3
342  do l= 1, 3
343  hecmat%CAL(9*(countcal-1)+3*(k-1)+l) = hecmat%cmat%pair(j)%val(k,l)
344  enddo
345  enddo
346  endif
347  if (hecmat%cmat%pair(j)%i == i .and. hecmat%cmat%pair(j)%j > i) then
348  countcau = countcau + 1
349  hecmat%itemCU(countcau) = hecmat%cmat%pair(j)%j
350  do k= 1, 3
351  do l= 1, 3
352  hecmat%CAU(9*(countcau-1)+3*(k-1)+l) = hecmat%cmat%pair(j)%val(k,l)
353  enddo
354  enddo
355  endif
356  enddo
357  enddo
358  end subroutine hecmw_cmat_lu
359 
360 
361  subroutine hecmw_cmat_lu_free( hecMAT )
362  type (hecmwst_matrix) :: hecmat
363 
364  deallocate (hecmat%indexCL, hecmat%itemCL, hecmat%CAL)
365  deallocate (hecmat%indexCU, hecmat%itemCU, hecmat%CAU)
366  end subroutine hecmw_cmat_lu_free
367 
368  subroutine hecmw_cmat_substitute( dest, src )
369  implicit none
370  type(hecmwst_matrix_contact) :: dest
371  type(hecmwst_matrix_contact) :: src
372  dest%n_val = src%n_val
373  dest%max_val = src%max_val
374  if (associated(src%pair)) dest%pair => src%pair
375  dest%checked = src%checked
376  dest%sorted = src%sorted
377  dest%max_row = src%max_row
378  dest%max_col = src%max_col
379  end subroutine hecmw_cmat_substitute
380 
381 end module hecmw_matrix_contact
subroutine, public hecmw_cmat_pack(cmat)
subroutine, public hecmw_cmat_substitute(dest, src)
subroutine, public hecmw_cmat_init(cmat)
integer function, public hecmw_cmat_n_newpair(hecMAT, symm)
subroutine, public hecmw_cmat_lu(hecMAT)
subroutine, public hecmw_cmat_finalize(cmat)
subroutine, public hecmw_cmat_lu_free(hecMAT)
subroutine, public hecmw_cmat_clear(cmat)
subroutine, public hecmw_cmat_ass_bc(hecMAT, inode, idof, RHS)
subroutine, public hecmw_cmat_add(cmat, i, j, val)
subroutine, public hecmw_cmat_multvec_add(cmat, X, Y, len)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=kint) function hecmw_comm_get_comm()
integer(kind=4), parameter kreal
subroutine hecmw_abort(comm)