FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
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
25contains
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
381end module hecmw_matrix_contact
integer function, public hecmw_cmat_n_newpair(hecmat, symm)
subroutine, public hecmw_cmat_pack(cmat)
subroutine, public hecmw_cmat_substitute(dest, src)
subroutine, public hecmw_cmat_init(cmat)
subroutine, public hecmw_cmat_ass_bc(hecmat, inode, idof, rhs)
subroutine, public hecmw_cmat_lu(hecmat)
subroutine, public hecmw_cmat_finalize(cmat)
subroutine, public hecmw_cmat_clear(cmat)
subroutine, public hecmw_cmat_multvec_add(cmat, x, y, len)
subroutine, public hecmw_cmat_add(cmat, i, j, val)
subroutine, public hecmw_cmat_lu_free(hecmat)
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)