FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_matrix_ordering_CM.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 use hecmw_util
8 implicit none
9
10 private
13
14 integer(kind=kint), parameter :: DEBUG = 0
15
16contains
17
18 subroutine hecmw_matrix_ordering_cm(N, indexL, itemL, indexU, itemU, &
19 perm, iperm)
20 implicit none
21 integer(kind=kint), intent(in) :: n
22 integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
23 integer(kind=kint), intent(in) :: iteml(:), itemu(:)
24 integer(kind=kint), intent(out) :: perm(:), iperm(:)
25 integer(kind=kint), parameter :: nminmax = 5
26 integer(kind=kint) :: i, nmin
27 integer(kind=kint) :: mins(nminmax)
28 integer(kind=kint), allocatable :: degs(:)
29 integer(kind=kint), allocatable :: nlevel(:)
30 integer(kind=kint), allocatable :: lv_index(:,:), lv_item(:,:)
31 integer(kind=kint) :: nlevel_max, max_id
32 allocate(degs(n))
33 call count_degrees(n, indexl, indexu, itemu, degs)
34 call find_minimum_degrees(n, degs, nminmax, nmin, mins)
35 allocate(nlevel(nmin), lv_index(0:n,nmin), lv_item(n,nmin))
36 ! perform CM ordering starting from each minimum degree node
37 !$omp parallel default(none),private(i), &
38 !$omp& shared(nmin,N,indexL,itemL,indexU,itemU,degs,mins,nlevel,lv_index,lv_item)
39 !$omp do
40 do i=1,nmin
41 call ordering_cm_inner(n, indexl, iteml, indexu, itemu, degs, mins(i), &
42 nlevel(i), lv_index(:,i), lv_item(:,i))
43 if (debug > 0) write(*,*) 'DEBUG:: hecmw_matrix_ordering_CM: i, nstart, nlevel = ', i, mins(i), nlevel(i)
44 end do
45 !$omp end do
46 !$omp end parallel
47 deallocate(degs)
48 ! choose CM ordering with maximum nlevel
49 nlevel_max = nlevel(1)
50 max_id = 1
51 do i=2,nmin
52 if (nlevel(i) > nlevel_max) then
53 nlevel_max = nlevel(i)
54 max_id = i
55 end if
56 end do
57 if (debug > 0) write(*,*) 'DEBUG:: hecmw_matrix_ordering_CM: chose ordering',max_id
58 do i=1,n
59 perm(i) = lv_item(i,max_id)
60 iperm(perm(i)) = i
61 end do
62 deallocate(nlevel, lv_index, lv_item)
63 end subroutine hecmw_matrix_ordering_cm
64
65 subroutine hecmw_matrix_ordering_rcm(N, indexL, itemL, indexU, itemU, &
66 perm, iperm)
67 implicit none
68 integer(kind=kint), intent(in) :: n
69 integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
70 integer(kind=kint), intent(in) :: iteml(:), itemu(:)
71 integer(kind=kint), intent(out) :: perm(:), iperm(:)
72 call hecmw_matrix_ordering_cm(n, indexl, iteml, indexu, itemu, perm, iperm)
73 call reverse_ordering(n, perm, iperm)
74 if (debug > 0) then
75 call write_nonzero_profile(n, indexl, iteml, indexu, itemu, perm, iperm)
76 call write_perm(n, perm, iperm)
77 endif
78 end subroutine hecmw_matrix_ordering_rcm
79
80 subroutine ordering_cm_inner(N, indexL, itemL, indexU, itemU, degs, nstart, &
81 nlevel, lv_index, lv_item)
82 implicit none
83 integer(kind=kint), intent(in) :: n
84 integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
85 integer(kind=kint), intent(in) :: iteml(:), itemu(:)
86 integer(kind=kint), intent(in) :: degs(:)
87 integer(kind=kint), intent(in) :: nstart
88 integer(kind=kint), intent(out) :: nlevel
89 integer(kind=kint), intent(out) :: lv_index(0:)
90 integer(kind=kint), intent(out) :: lv_item(:)
91 integer(kind=kint), allocatable :: iwk(:)
92 integer(kind=kint) :: level, cntall, cnt, j, jnode, k, knode
93 allocate(iwk(n))
94 iwk(:) = 0
95 lv_index(0) = 0
96 ! first level
97 iwk(nstart) = 1
98 cntall = 1
99 lv_item(1) = nstart
100 lv_index(1) = 1
101 ! other levels
102 do level=2,n
103 cnt = 0
104 ! all nodes in previous level
105 prlv: do j = lv_index(level-2)+1, lv_index(level-1)
106 jnode = lv_item(j)
107 ! all connected nodes
108 do k = indexl(jnode-1)+1, indexl(jnode)
109 knode = iteml(k)
110 if (iwk(knode) == 0) then
111 iwk(knode) = level
112 cnt = cnt + 1
113 cntall = cntall + 1
114 lv_item(cntall) = knode
115 !if (cntall == N) exit PRLV
116 end if
117 end do
118 do k = indexu(jnode-1)+1, indexu(jnode)
119 knode = itemu(k)
120 if (knode > n) cycle
121 if (iwk(knode) == 0) then
122 iwk(knode) = level
123 cnt = cnt + 1
124 cntall = cntall + 1
125 lv_item(cntall) = knode
126 !if (cntall == N) exit PRLV
127 end if
128 end do
129 end do prlv
130 if (cnt == 0) then
131 if (debug > 0) write(*,*) 'DEBUG: choose any uncolored node..'
132 do knode = 1, n
133 if (iwk(knode) == 0) then
134 iwk(knode) = level
135 cnt = cnt + 1
136 cntall = cntall + 1
137 lv_item(cntall) = knode
138 exit
139 end if
140 end do
141 endif
142 lv_index(level) = cntall
143 call sort_nodes_by_degree(lv_item, lv_index(level-1)+1, lv_index(level), degs)
144 if (cntall == n) then
145 nlevel = level
146 exit
147 end if
148 end do
149 if (debug > 0) then
150 level = 0
151 do j = 1, n
152 jnode = lv_item(j)
153 if (jnode <= 0 .or. jnode > n) stop 'ERROR: ordering_CM_inner: out of range'
154 if (iwk(jnode) < 0) stop 'ERROR: ordering_CM_inner: duplicated node found'
155 if (iwk(jnode) < level) stop 'ERROR: ordering_CM_inner: not in level order'
156 level = iwk(jnode)
157 iwk(jnode) = -1
158 enddo
159 do j = 1, n
160 if (iwk(j) /= -1) stop 'ERROR: ordering_CM_inner: non-numbered node found'
161 enddo
162 endif
163 deallocate(iwk)
164 end subroutine ordering_cm_inner
165
166 subroutine count_degrees(N, indexL, indexU, itemU, degs)
167 implicit none
168 integer(kind=kint), intent(in) :: n
169 integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
170 integer(kind=kint), intent(in) :: itemu(:)
171 integer(kind=kint), intent(out) :: degs(:)
172 integer(kind=kint) :: i, j
173 do i=1,n
174 degs(i) = indexl(i) - indexl(i-1)
175 do j = indexu(i-1)+1, indexu(i)
176 if (itemu(j) <= n) degs(i) = degs(i) + 1
177 end do
178 end do
179 end subroutine count_degrees
180
181 subroutine find_minimum_degrees(N, degs, nminmax, nmin, mins)
182 implicit none
183 integer(kind=kint), intent(in) :: n
184 integer(kind=kint), intent(in) :: degs(:)
185 integer(kind=kint), intent(in) :: nminmax
186 integer(kind=kint), intent(out) :: nmin
187 integer(kind=kint), intent(out) :: mins(nminmax)
188 integer(kind=kint) :: degmin, i, j
189 degmin = n
190 nmin = 0
191 do i=1,n
192 ! skip unconnected nodes
193 if (degs(i) == 0) cycle
194 if (degs(i) < degmin) then
195 degmin = degs(i)
196 nmin = 1
197 mins(1) = i
198 else if (degs(i) == degmin) then
199 nmin = nmin + 1
200 if (nmin <= nminmax) mins(nmin) = i
201 end if
202 end do
203 if (debug > 0) write(*,*) 'DEBUG:: find_minimum_degrees: nmin, deg = ', nmin, degmin
204 if (nmin > nminmax) nmin = nminmax
205 end subroutine find_minimum_degrees
206
207 recursive subroutine sort_nodes_by_degree(lv_item, istart, iend, degs)
208 implicit none
209 integer(kind=kint), intent(inout) :: lv_item(:)
210 integer(kind=kint), intent(in) :: istart, iend
211 integer(kind=kint), intent(in) :: degs(:)
212 integer(kind=kint) :: pivot, left, right, itmp
213 if (istart >= iend) return
214 pivot = degs(lv_item((istart + iend) / 2))
215 left = istart
216 right = iend
217 do
218 do while (degs(lv_item(left)) < pivot)
219 left = left + 1
220 end do
221 do while (pivot < degs(lv_item(right)))
222 right = right - 1
223 end do
224 if (left >= right) exit
225 itmp = lv_item(left)
226 lv_item(left) = lv_item(right)
227 lv_item(right) = itmp
228 left = left + 1
229 right = right - 1
230 end do
231 call sort_nodes_by_degree(lv_item, istart, left-1, degs)
232 call sort_nodes_by_degree(lv_item, right+1, iend, degs)
233 end subroutine sort_nodes_by_degree
234
235 subroutine reverse_ordering(N, perm, iperm)
236 implicit none
237 integer(kind=kint), intent(in) :: n
238 integer(kind=kint), intent(inout) :: perm(:), iperm(:)
239 integer(kind=kint) :: i, n1
240 n1 = n + 1
241 do i=1,n
242 iperm(i) = n1 - iperm(i)
243 perm(iperm(i)) = i
244 end do
245 end subroutine reverse_ordering
246
247 subroutine write_nonzero_profile(N, indexL, itemL, indexU, itemU, perm, iperm)
248 implicit none
249 integer(kind=kint), intent(in) :: n
250 integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
251 integer(kind=kint), intent(in) :: iteml(:), itemu(:)
252 integer(kind=kint), intent(in) :: perm(:), iperm(:)
253 integer(kind=kint), parameter :: f_org = 901
254 integer(kind=kint), parameter :: f_new = 902
255 integer(kind=kint) :: i, j, irow, jcol
256 open(f_org, file='nzprof_org.txt', status='replace')
257 do irow = 1, n
258 i = irow
259 do j = indexl(i-1)+1, indexl(i)
260 jcol = iteml(j)
261 write(f_org,*) irow, jcol
262 end do
263 do j = indexu(i-1)+1, indexu(i)
264 jcol = itemu(j)
265 if (jcol > n) cycle
266 write(f_org,*) irow, jcol
267 end do
268 end do
269 close(f_org)
270 open(f_new, file='nzprof_new.txt', status='replace')
271 do irow = 1, n
272 i = perm(irow)
273 do j = indexl(i-1)+1, indexl(i)
274 jcol = iteml(j)
275 write(f_new,*) irow, iperm(jcol)
276 end do
277 do j = indexu(i-1)+1, indexu(i)
278 jcol = itemu(j)
279 if (jcol > n) cycle
280 write(f_new,*) irow, iperm(jcol)
281 end do
282 end do
283 close(f_new)
284 end subroutine write_nonzero_profile
285
286 subroutine write_perm(N, perm, iperm)
287 implicit none
288 integer(kind=kint), intent(in) :: n
289 integer(kind=kint), intent(in) :: perm(:), iperm(:)
290 integer(kind=kint), parameter :: f_perm = 903
291 integer(kind=kint) :: i
292 open(f_perm, file='perm_iperm.txt', status='replace')
293 do i = 1, n
294 write(f_perm,*) i, perm(i), iperm(i)
295 end do
296 close(f_perm)
297 end subroutine write_perm
298
I/O and Utility.
Definition: hecmw_util_f.F90:7
subroutine, public hecmw_matrix_ordering_cm(n, indexl, iteml, indexu, itemu, perm, iperm)
subroutine, public hecmw_matrix_ordering_rcm(n, indexl, iteml, indexu, itemu, perm, iperm)