FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
m_matrix_partition_info.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!-------------------------------------------------------------------------------
5
7
8 use hecmw_util
11
12 private !default access control
13
16 public reovec
17
18 ! Handle matrix partitioning informations.
20 integer(kind=kint) :: ndm ! number of dm
21 type(child_matrix), pointer :: dm(:) ! partitioned matrix. 2^n number
22 integer(kind=kint) :: neqns_d ! number of eqns in D matrix
23 real(kind=kreal), pointer :: dsln(:,:) ! non-diagonal elements of dens D matrix which kept in proc0
24 real(kind=kreal), pointer :: diag(:,:) ! diagonal elements of dens D matrix
25 integer(kind=kint), pointer :: part_all(:) ! index of corresponding dm of a0 row
26 integer(kind=kint), pointer :: iperm_all(:) ! index in partitioned matrix of a0 row
27 end type
28
29 ! node of bisection partition tree
30 type matrix_partition_node
31 integer(kind=kint) :: id ! position in array
32 integer(kind=kint) :: depth ! depth in tree. root is 0.
33 integer(kind=kint) :: idm ! corresponding divided matrix no. root and tree node have 0. leaf node have 1..2^n
34
35 integer(kind=kint) :: neqns_a, nttbr_a ! size and number of non-zero element in a0.
36 integer(kind=kint) :: neqns_d ! size of D.
37 integer(kind=kint), pointer :: irow(:), jcol(:)
38 integer(kind=kint), pointer :: idx_g_a(:), idx_g_d(:) ! global indices of each point
39 end type
40
41 integer(kind=kint) :: maxdepth ! root is 0.
42 integer(kind=kint) :: iend_trunk ! position of last trunk node.
43 integer(kind=kint) :: ista_leaf ! position of first leaf node.
44 integer(kind=kint) :: iend_leaf ! position of last leaf node.
45
46 type(matrix_partition_node), pointer :: tree_array(:)
47
48 integer(kind=kint), parameter :: ilog = 16
49
50contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
51
52 subroutine matrix_partition_recursive_bisection(a0, ndiv, pmi)
53 implicit none
54 type(irjc_square_matrix), intent(inout) :: a0
55 integer(kind=kint), intent(in) :: ndiv
56 type(matrix_partition_info), intent(out) :: pmi
57
58 type(matrix_partition_node), pointer :: mdroot ! root of partition tree
59
60 integer(kind=kint) :: i, ndeg, nndeg, ndegt
61
62
63 ! make matrix partition tree using METIS library.
64 ! leaf of tree is according to divided child matrix.
65 ! trunk of tree is according to D region (parent process).
66 write(6,*) 'METIS 4.0 library is used for matrix partitioning.'
67 write(6,*) 'Copyright 1997, Regents of the University of Minnesota'
68 write(6,*) 'http://glaros.dtc.umn.edu/gkhome/metis/metis/overview'
69 pmi%ndm=2**ndiv
70 call matrix_partition_tree_initialize(ndiv)
71 call matrix_partition_tree_get_root(mdroot)
72 mdroot%neqns_a=a0%neqns
73 mdroot%nttbr_a=a0%nttbr
74 allocate(mdroot%idx_g_a(a0%neqns))
75 do i=1,a0%neqns
76 mdroot%idx_g_a(i)=i
77 end do
78 allocate(mdroot%irow(size(a0%irow)))
79 allocate(mdroot%jcol(size(a0%jcol)))
80 mdroot%irow=a0%irow
81 mdroot%jcol=a0%jcol
82 call make_matrix_partition_tree(mdroot)
83
84 ! get global partition array from matrix partition tree.
85 allocate(pmi%part_all(a0%neqns), pmi%iperm_all(a0%neqns))
86 call mkpart(pmi%part_all, pmi%iperm_all, pmi%neqns_d)
87
88 ndeg=a0%ndeg
89 nndeg = ndeg*ndeg
90 ndegt = (ndeg+1)*ndeg/2 ! triangle number. 6 for ndeg=3, 3 for ndeg=2
91 allocate(pmi%dsln(nndeg, pmi%neqns_d*(pmi%neqns_d - 1) / 2))
92 allocate(pmi%diag(ndegt, pmi%neqns_d))
93 allocate(pmi%dm(2**ndiv))
94 pmi%dsln=0
95 pmi%diag=0
96 call divmat(a0, pmi%part_all, pmi%iperm_all, pmi%dm, pmi%dsln, pmi%diag, pmi%neqns_d)
97
99
100 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
101
102 subroutine matrix_partition_tree_initialize(ndiv)
103 ! allocate 2^ndiv matrix_partition_node and set belongings in processes for each node.
104
105 implicit none
106 integer(kind=kint), intent(in) :: ndiv ! number of bisection
107
108 type(matrix_partition_node), pointer :: this
109
110 integer(kind=kint) :: tsize, i
111
112
113 ! allocate matrix divide tree
114 tsize=1
115 do i=1,ndiv
116 tsize=tsize+2**i
117 end do
118 allocate(tree_array(tsize))
119 iend_trunk = tsize - 2**ndiv
120 ista_leaf = iend_trunk+1
121 iend_leaf = tsize
122 maxdepth = ndiv
123 do i=1,tsize
124 if (i .le. iend_trunk) then
125 this=>tree_array(i)
126 this%id=i
127 this%idm=0
128 this%depth=depthinbinarytree(this%id)
129 else if (i .ge. ista_leaf) then
130 this=>tree_array(i)
131 this%id=i
132 this%idm=i - iend_trunk
133 this%depth=maxdepth
134 end if
135 end do
136 end subroutine matrix_partition_tree_initialize
137
138 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
139
140 subroutine matrix_partition_tree_get_root(root)
141 implicit none
142 type(matrix_partition_node),pointer :: root
143 root=>tree_array(1)
144 end subroutine matrix_partition_tree_get_root
145
146 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
147
148 subroutine matrix_partition_tree_get_left_child(this, left)
149 implicit none
150 type(matrix_partition_node), pointer :: this
151 type(matrix_partition_node), pointer :: left
152 integer(kind=kint) :: id_this, id_left
153 id_this = this%id
154 id_left = id_this*2
155 left => tree_array(id_left)
156 end subroutine matrix_partition_tree_get_left_child
157
158 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
159
160 subroutine matrix_partition_tree_get_right_child(this, right)
161 implicit none
162 type(matrix_partition_node), pointer :: this
163 type(matrix_partition_node), pointer :: right
164 integer(kind=kint) :: id_this, id_right
165 id_this = this%id
166 id_right = id_this*2+1
167 right => tree_array(id_right)
168 end subroutine matrix_partition_tree_get_right_child
169
170 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
171
172 recursive subroutine make_matrix_partition_tree(this)
173 ! make matrix partition tree according to graph structure of given matrix via matrix_partition_node.
174 implicit none
175 type(matrix_partition_node), pointer :: this, left, right
176 integer(kind=kint), allocatable :: ip_left(:), ip_right(:), ip_d(:), part(:), iperm(:)
177 integer(kind=kint) :: loc1, loc2, loc3, ipass
178 integer(kind=kint) :: i,j,k,l
179
180 allocate(ip_left(this%neqns_a), ip_right(this%neqns_a), ip_d(this%neqns_a)) ! local permtation index used by metis.
181 allocate(part(this%neqns_a), iperm(this%neqns_a))
182
183
184 if (this%depth .eq. maxdepth) then ! leaf node
185 return
186 end if
187
188 call matrix_partition_tree_get_left_child(this, left)
189 call matrix_partition_tree_get_right_child(this, right)
190
191 call bi_part_directive(this%neqns_a, this%nttbr_a, this%irow, this%jcol, left%neqns_a, right%neqns_a, this%neqns_d)
192
193 if ((left%neqns_a + right%neqns_a + this%neqns_d) .ne. this%neqns_a) then
194 write(6,*) 'Error: Fail to matrix partitioning in parallel direct solver.'
195 write(6,*) 'Matrix is too small.'
196 write(ilog,*) 'Error: Fail to matrix partitioning in parallel direct solver.'
197 write(ilog,*) 'Matrix is too small.'
199 end if
200
201 call get_part_result(left%neqns_a, ip_left, right%neqns_a, ip_right, this%neqns_d, ip_d)
202
203 ! set belongings
204 do i=1,left%neqns_a ! left child
205 part(ip_left(i))=1
206 end do
207 do i=1,right%neqns_a ! right child
208 part(ip_right(i))=2
209 end do
210 do i=1,this%neqns_d ! remain in this node
211 part(ip_d(i))=3
212 end do
213
214 ! set global indices of each point in terms of local array indices.
215 allocate(left%idx_g_a(left%neqns_a))
216 allocate(right%idx_g_a(right%neqns_a))
217 allocate(this%idx_g_d(this%neqns_d))
218
219 loc1=0
220 loc2=0
221 loc3=0
222 do i=1,this%neqns_a
223 if(part(i) .eq. 1) then
224 loc1=loc1+1
225 iperm(i)=loc1
226 left%idx_g_a(loc1)=this%idx_g_a(i)
227 end if
228 if(part(i) .eq. 2) then
229 loc2=loc2+1
230 iperm(i)=loc2
231 right%idx_g_a(loc2)=this%idx_g_a(i)
232 end if
233 if(part(i) .eq. 3) then
234 loc3=loc3+1
235 iperm(i)=loc3
236 this%idx_g_d(loc3)=this%idx_g_a(i)
237 end if
238 end do
239
240 ! devide a0 to a_left and a_right
241 do ipass=1,2
242 left%nttbr_a = 0
243 right%nttbr_a = 0
244
245 do l=1,this%nttbr_a
246 i=this%irow(l)
247 j=this%jcol(l)
248
249 if ((part(i) .eq. 1) .and. (part(j) .eq. 1)) then !left
250 left%nttbr_a=left%nttbr_a + 1
251 if(ipass .eq. 2) then
252 left%irow(left%nttbr_a)=iperm(i)
253 left%jcol(left%nttbr_a)=iperm(j)
254 end if
255 cycle
256 end if
257
258 if ((part(i) .eq. 2) .and. (part(j) .eq. 2)) then !right
259 right%nttbr_a=right%nttbr_a + 1
260 if(ipass .eq. 2) then
261 right%irow(right%nttbr_a)=iperm(i)
262 right%jcol(right%nttbr_a)=iperm(j)
263 end if
264 cycle
265 end if
266 end do
267
268 if (ipass .eq. 1) then
269 allocate(left%irow(left%nttbr_a), left%jcol(left%nttbr_a))
270 allocate(right%irow(right%nttbr_a), right%jcol(right%nttbr_a))
271 end if
272 end do
273
274 deallocate(ip_left, ip_right, ip_d, iperm, part, this%irow, this%jcol)
275
276 call make_matrix_partition_tree(left)
277 call make_matrix_partition_tree(right)
278 return
279 end subroutine make_matrix_partition_tree
280
281 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
282
283 subroutine mkpart(part_g, iperm_g, neqns_d)
284 ! make global partition array from matrix partition tree.
285 implicit none
286 integer(kind=kint), intent(out) :: part_g(:)
287 integer(kind=kint), intent(out) :: iperm_g(:)
288 integer(kind=kint), intent(out) :: neqns_d
289
290 type(matrix_partition_node), pointer :: node
291 integer(kind=kint) :: ncount
292 integer(kind=kint) :: i,j
293
294
295
296 ! for leaf
297 do i=iend_leaf, ista_leaf, -1
298 node=>tree_array(i)
299 ncount=0 ! node local numbering
300 do j=1,node%neqns_a
301 ncount=ncount+1
302 part_g(node%idx_g_a(j))=node%idm ! index of divided matrix. 1..2**ndiv
303 iperm_g(node%idx_g_a(j))=ncount ! divided matrix local numbering
304 end do
305 deallocate(node%idx_g_a)
306 end do
307
308 ! for trunk
309 ncount=0 !global numbering
310 do i=iend_trunk, 1, -1
311 node => tree_array(i)
312 do j=1,node%neqns_d
313 ncount=ncount+1
314 part_g(node%idx_g_d(j))=node%idm ! parent D. 0
315 iperm_g(node%idx_g_d(j))=ncount ! global numbering
316 end do
317 deallocate(node%idx_g_d)
318 end do
319
320 neqns_d = ncount ! size of D
321 deallocate(tree_array) ! tree_array is not used after here.
322 return
323 end subroutine mkpart
324
325 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
326
327 subroutine divmat(a0, part_g, iperm_g, dm, dsln, diag, neqns_d)
328 implicit none
329 ! divide a0 to n children and D region on parent, according to the global permtation information.
330
331 type(irjc_square_matrix), intent(in) :: a0
332 integer(kind=kint), intent(in) :: part_g(:), iperm_g(:)
333
334 type(child_matrix), intent(out), target :: dm(:)
335 real(kind=kreal), intent(out) :: dsln(:,:), diag(:,:) ! for dens D
336
337 integer(kind=kint), intent(in) :: neqns_d
338
339 ! internal
340 type(child_matrix), pointer :: dmc
341 integer(kind=kint) :: ndeg, ipass, ipos, itmp
342 integer(kind=kint) :: i,j,k,l,m,n, ii, jj
343
344 ! set neqns for each divided matrix
345 ndeg = a0%ndeg
346 do i=1,size(dm)
347 dm(i)%a%neqns = 0
348 end do
349
350 do i=1,a0%neqns
351 k=part_g(i)
352 if (k .ne. 0 ) then
353 dm(k)%a%neqns = dm(k)%a%neqns + 1
354 end if
355 end do
356
357 do i=1,size(dm)
358 dm(i)%ndeg = ndeg
359 dm(i)%a%ndeg = ndeg
360 dm(i)%c%ndeg = ndeg
361 dm(i)%c%nrows = neqns_d
362 dm(i)%c%ncols = dm(i)%a%neqns
363 dm(i)%ista_c = dm(i)%a%neqns + 1
364 end do
365
366 ! divide matrix element to child matrix
367 do ipass=1,2
368 do i=1,size(dm)
369 dm(i)%a%nttbr = 0
370 dm(i)%c%nttbr = 0
371 end do
372
373 do l=1,a0%nttbr
374 i=a0%irow(l)
375 j=a0%jcol(l)
376
377 if (part_g(i) .eq. part_g(j)) then ! same matrix
378 if ((part_g(i) .eq. 0) .or. (part_g(j) .eq. 0)) then ! D
379 if (iperm_g(i) .eq. iperm_g(j)) then ! diag in D
380 ipos=1
381 do m=1,ndeg
382 do n=1,m
383 diag(ipos, iperm_g(i))=a0%val(ndeg*(m-1)+n,l)
384 ipos=ipos+1
385 end do
386 end do
387
388 else if (iperm_g(i) .gt. iperm_g(j)) then ! dsln in D
389
390 !count index in dsln for specified i,j.
391 ii = iperm_g(i)
392 jj = iperm_g(j)
393
394 k = (ii-1)*(ii-2)/2 + jj
395
396 dsln(:,k)=a0%val(:,l)
397
398 else ! dsln in D with inverse
399 !count index in dsln for specified i,j.
400 jj = iperm_g(i)
401 ii = iperm_g(j)
402
403 k = (ii-1)*(ii-2)/2 + jj
404
405 do m=1,ndeg
406 do n=1,ndeg
407 dsln((m-1)*ndeg+n,k)=a0%val(m+(n-1)*ndeg,l)
408 end do
409 end do
410
411 end if
412 else !An
413 dmc=>dm(part_g(i))
414 dmc%a%nttbr = dmc%a%nttbr + 1
415 if(ipass.eq.2) then
416 dmc%a%irow(dmc%a%nttbr) = iperm_g(i)
417 dmc%a%jcol(dmc%a%nttbr) = iperm_g(j)
418 dmc%a%val(:,dmc%a%nttbr) = a0%val(:,l)
419 end if
420 end if
421 cycle
422 end if
423
424 if (part_g(i) .eq. 0) then !C run right in ndm part_g(j)
425 dmc=>dm(part_g(j))
426 dmc%c%nttbr = dmc%c%nttbr + 1
427 if(ipass.eq.2) then
428 dmc%c%irow(dmc%c%nttbr)=iperm_g(i)
429 dmc%c%jcol(dmc%c%nttbr)=iperm_g(j)
430 dmc%c%val(:,dmc%c%nttbr)=a0%val(:,l)
431 end if
432 cycle
433 end if
434
435 if (part_g(j) .eq. 0) then !C run down in ndm part_g(i). Invert val
436 dmc=>dm(part_g(i))
437 dmc%c%nttbr = dmc%c%nttbr + 1
438 if(ipass.eq.2) then
439 dmc%c%irow(dmc%c%nttbr)=iperm_g(j)
440 dmc%c%jcol(dmc%c%nttbr)=iperm_g(i)
441 do m=1,ndeg
442 do n=1,ndeg
443 dmc%c%val(n+ndeg*(m-1),dmc%c%nttbr)=a0%val(ndeg*(n-1)+m,l)
444 end do
445 end do
446 end if
447 cycle
448 end if
449
450 stop !never come here
451 end do
452
453 if (ipass .eq. 1) then
454 do i=1,size(dm)
455 allocate(dm(i)%a%irow(dm(i)%a%nttbr))
456 allocate(dm(i)%a%jcol(dm(i)%a%nttbr))
457 allocate(dm(i)%a%val(ndeg*ndeg,dm(i)%a%nttbr))
458
459 allocate(dm(i)%c%irow(dm(i)%c%nttbr))
460 allocate(dm(i)%c%jcol(dm(i)%c%nttbr))
461 allocate(dm(i)%c%val(ndeg*ndeg,dm(i)%c%nttbr))
462 end do
463 end if
464 end do
465
466 return
467 end subroutine divmat
468
469 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
470
471 subroutine reovec(r, iperm)
472 ! reorder vector as given iperm
473
474 implicit none
475
476 integer(kind=kint), dimension(:), intent(in) :: iperm ! permtation vector
477 real(kind=kreal), dimension(:,:), intent(inout) :: r ! reordered result will be return
478
479 real(kind=kreal), allocatable, dimension(:,:) :: tmp
480 integer(kind=kint) :: ndeg, neqns
481 integer(kind=kint), dimension(2) :: ishape
482 integer(kind=kint) :: i,j,k,l
483
484 ishape = shape(r)
485 ndeg = ishape(1)
486 neqns = ishape(2)
487 allocate(tmp(ndeg, neqns))
488
489
490 do i=1, neqns
491 do j=1, ndeg
492 tmp(j,iperm(i))=r(j,i)
493 end do
494 end do
495
496 r=tmp
497
498 deallocate(tmp)
499
500 return
501 end subroutine reovec
502
503 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
504
505 function depthinbinarytree(n)
506 ! return depth of given integer "n" in 1-ofset binary tree. root have 0 depth.
507 implicit none
508 integer(kind=kint) :: depthinbinarytree
509 integer(kind=kint), intent(in) :: n
510 integer(kind=kint) :: i
511 i=0
512 do while (n .ge. 2**(i+1))
513 i=i+1
514 end do
515 depthinbinarytree = i ! root is 0
516 end function
517
518 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
519
void bi_part_directive(int *neqns, int *nttbr, int *irow, int *jcol, int *num_graph1, int *num_graph2, int *num_separator)
Definition: matrix_repart.c:18
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)
subroutine, public reovec(r, iperm)
subroutine, public matrix_partition_recursive_bisection(a0, ndiv, pmi)
void get_part_result(int *num_graph1, int *igraph1, int *num_graph2, int *igraph2, int *num_separator, int *iseparator)