FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_matrix_reorder.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
9 implicit none
10
11 private
17
18contains
19
20 subroutine hecmw_matrix_reorder_profile(N, perm, iperm, &
21 indexL, indexU, itemL, itemU, indexLp, indexUp, itemLp, itemUp)
22 implicit none
23 integer(kind=kint), intent(in) :: n
24 integer(kind=kint), intent(in) :: perm(:), iperm(:)
25 integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
26 integer(kind=kint), intent(in) :: iteml(:), itemu(:)
27 integer(kind=kint), intent(out) :: indexlp(0:), indexup(0:)
28 integer(kind=kint), intent(out) :: itemlp(:), itemup(:)
29 integer(kind=kint) :: cntl, cntu, inew, iold, j, jold, jnew
30 cntl = 0
31 cntu = 0
32 indexlp(0) = 0
33 indexup(0) = 0
34 do inew=1,n
35 iold = perm(inew)
36 ! original L
37 do j = indexl(iold-1)+1, indexl(iold)
38 jold = iteml(j)
39 jnew = iperm(jold)
40 if (jnew < inew) then
41 cntl = cntl + 1
42 itemlp(cntl) = jnew
43 else
44 cntu = cntu + 1
45 itemup(cntu) = jnew
46 end if
47 end do
48 ! original U
49 do j = indexu(iold-1)+1, indexu(iold)
50 jold = itemu(j)
51 if (jold > n) cycle
52 jnew = iperm(jold)
53 if (jnew < inew) then
54 cntl = cntl + 1
55 itemlp(cntl) = jnew
56 else
57 cntu = cntu + 1
58 itemup(cntu) = jnew
59 end if
60 end do
61 indexlp(inew) = cntl
62 indexup(inew) = cntu
63 call hecmw_qsort_int_array(itemlp, indexlp(inew-1)+1, indexlp(inew))
64 call hecmw_qsort_int_array(itemup, indexup(inew-1)+1, indexup(inew))
65 end do
66 end subroutine hecmw_matrix_reorder_profile
67
68 subroutine hecmw_matrix_reorder_values(N, NDOF, perm, iperm, &
69 indexL, indexU, itemL, itemU, AL, AU, D, &
70 indexLp, indexUp, itemLp, itemUp, ALp, AUp, Dp)
71 implicit none
72 integer(kind=kint), intent(in) :: n, ndof
73 integer(kind=kint), intent(in) :: perm(:), iperm(:)
74 integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
75 integer(kind=kint), intent(in) :: iteml(:), itemu(:)
76 real(kind=kreal), intent(in) :: al(:), au(:), d(:)
77 integer(kind=kint), intent(in) :: indexlp(0:), indexup(0:)
78 integer(kind=kint), intent(in) :: itemlp(:), itemup(:)
79 real(kind=kreal), intent(out) :: alp(:), aup(:), dp(:)
80 dp = 0.d0
81 alp = 0.d0
82 aup = 0.d0
83 ! call reorder_diag(N, NDOF, perm, D, Dp)
84 ! call reorder_off_diag(N, NDOF, perm, &
85 ! indexL, indexU, itemL, itemU, AL, AU, &
86 ! indexLp, itemLp, ALp)
87 ! call reorder_off_diag(N, NDOF, perm, &
88 ! indexL, indexU, itemL, itemU, AL, AU, &
89 ! indexUp, itemUp, AUp)
90 call reorder_diag2(n, ndof, iperm, d, dp)
91 call reorder_off_diag2(n, ndof, iperm, &
92 indexl, iteml, al, &
93 indexlp, indexup, itemlp, itemup, alp, aup)
94 call reorder_off_diag2(n, ndof, iperm, &
95 indexu, itemu, au, &
96 indexlp, indexup, itemlp, itemup, alp, aup)
97 end subroutine hecmw_matrix_reorder_values
98
99 subroutine hecmw_matrix_reorder_vector(N, NDOF, perm, X, Xp)
100 implicit none
101 integer(kind=kint), intent(in) :: n, ndof
102 integer(kind=kint), intent(in) :: perm(:)
103 real(kind=kreal), intent(in) :: x(:)
104 real(kind=kreal), intent(out) :: xp(:)
105 integer(kind=kint) :: inew, iold, j0new, j0old, j
106 !$omp parallel default(none),private(inew,iold,j0new,j0old,j), &
107 !$omp shared(N,perm,NDOF,Xp,X)
108 !$omp do
109 do inew=1,n
110 iold = perm(inew)
111 j0new = (inew-1)*ndof
112 j0old = (iold-1)*ndof
113 do j=1,ndof
114 xp(j0new + j) = x(j0old + j)
115 end do
116 end do
117 !$omp end do
118 !$omp end parallel
119 end subroutine hecmw_matrix_reorder_vector
120
121 subroutine hecmw_matrix_reorder_back_vector(N, NDOF, perm, Xp, X)
122 implicit none
123 integer(kind=kint), intent(in) :: n, ndof
124 integer(kind=kint), intent(in) :: perm(:)
125 real(kind=kreal), intent(in) :: xp(:)
126 real(kind=kreal), intent(out) :: x(:)
127 integer(kind=kint) :: inew, iold, j0new, j0old, j
128 !$omp parallel default(none),private(inew,iold,j0new,j0old,j), &
129 !$omp& shared(N,perm,NDOF,X,Xp)
130 !$omp do
131 do inew=1,n
132 iold = perm(inew)
133 j0new = (inew-1)*ndof
134 j0old = (iold-1)*ndof
135 do j=1,ndof
136 x(j0old + j) = xp(j0new + j)
137 end do
138 end do
139 !$omp end do
140 !$omp end parallel
142
143 subroutine hecmw_matrix_reorder_renum_item(N, perm, indexXp, itemXp)
144 implicit none
145 integer(kind=kint), intent(in) :: n
146 integer(kind=kint), intent(in) :: perm(:)
147 integer(kind=kint), intent(in) :: indexxp(0:)
148 integer(kind=kint), intent(inout) :: itemxp(:)
149 integer(kind=kint) :: npx, i
150 npx = indexxp(n)
151 !$omp parallel default(none),private(i),shared(NPX,itemXp,perm)
152 !$omp do
153 do i=1,npx
154 itemxp(i) = perm( itemxp(i) )
155 end do
156 !$omp end do
157 !$omp end parallel
159
160 subroutine reorder_diag(N, NDOF, perm, D, Dp)
161 implicit none
162 integer(kind=kint), intent(in) :: n, ndof
163 integer(kind=kint), intent(in) :: perm(:)
164 real(kind=kreal), intent(in) :: d(:)
165 real(kind=kreal), intent(out) :: dp(:)
166 integer(kind=kint) :: ndof2, inew, iold, j0new, j0old, j
167 ndof2 = ndof*ndof
168 ! diagonal
169 do inew=1,n
170 iold = perm(inew)
171 j0new = (inew-1)*ndof2
172 j0old = (iold-1)*ndof2
173 do j=1,ndof2
174 dp(j0new + j) = d(j0old + j)
175 end do
176 end do
177 end subroutine reorder_diag
178
179 subroutine reorder_diag2(N, NDOF, iperm, D, Dp)
180 implicit none
181 integer(kind=kint), intent(in) :: n, ndof
182 integer(kind=kint), intent(in) :: iperm(:)
183 real(kind=kreal), intent(in) :: d(:)
184 real(kind=kreal), intent(out) :: dp(:)
185 integer(kind=kint) :: ndof2, inew, iold, j0new, j0old, j
186 ndof2 = ndof*ndof
187 ! diagonal
188 !$omp parallel default(none),private(iold,inew,j0old,j0new,j), &
189 !$omp& shared(N,iperm,NDOF2,Dp,D)
190 !$omp do
191 do iold=1,n
192 inew = iperm(iold)
193 j0old = (iold-1)*ndof2
194 j0new = (inew-1)*ndof2
195 do j=1,ndof2
196 dp(j0new + j) = d(j0old + j)
197 end do
198 end do
199 !$omp end do
200 !$omp end parallel
201 end subroutine reorder_diag2
202
203 subroutine reorder_off_diag(N, NDOF, perm, &
204 indexL, indexU, itemL, itemU, AL, AU, &
205 indexXp, itemXp, AXp)
206 implicit none
207 integer(kind=kint), intent(in) :: n, ndof
208 integer(kind=kint), intent(in) :: perm(:)
209 integer(kind=kint), intent(in) :: indexl(0:), indexu(0:)
210 integer(kind=kint), intent(in) :: iteml(:), itemu(:)
211 real(kind=kreal), intent(in) :: al(:), au(:)
212 integer(kind=kint), intent(in) :: indexxp(0:)
213 integer(kind=kint), intent(in) :: itemxp(:)
214 real(kind=kreal), intent(out) :: axp(:)
215 integer(kind=kint) :: ndof2, inew, iold
216 integer(kind=kint) :: jsoldl, jeoldl, jsoldu, jeoldu
217 integer(kind=kint) :: jnew, knew, kold, jold, l0new, l0old, l
218 ndof2 = ndof*ndof
219 ! new L
220 do inew=1,n
221 iold = perm(inew)
222 jsoldl = indexl(iold-1)+1
223 jeoldl = indexl(iold)
224 jsoldu = indexu(iold-1)+1
225 jeoldu = indexu(iold)
226 do jnew = indexxp(inew-1)+1, indexxp(inew)
227 knew = itemxp(jnew)
228 kold = perm(knew)
229 if (kold < iold) then
230 call hecmw_bsearch_int_array(iteml, jsoldl, jeoldl, kold, jold)
231 if (jold < 0) then
232 write(0,*) 'DEBUG:: jold < 0 in reorder_off_diag'
233 cycle
234 end if
235 l0new = (jnew-1)*ndof2
236 l0old = (jold-1)*ndof2
237 do l=1,ndof2
238 axp(l0new + l) = al(l0old + l)
239 end do
240 else
241 call hecmw_bsearch_int_array(itemu, jsoldu, jeoldu, kold, jold)
242 if (jold < 0) then
243 write(0,*) 'DEBUG:: jold < 0 in reorder_off_diag'
244 cycle
245 end if
246 l0new = (jnew-1)*ndof2
247 l0old = (jold-1)*ndof2
248 do l=1,ndof2
249 axp(l0new + l) = au(l0old + l)
250 end do
251 end if
252 end do
253 end do
254 end subroutine reorder_off_diag
255
256 subroutine reorder_off_diag2(N, NDOF, iperm, &
257 indexX, itemX, AX, &
258 indexLp, indexUp, itemLp, itemUp, ALp, AUp)
259 implicit none
260 integer(kind=kint), intent(in) :: n, ndof
261 integer(kind=kint), intent(in) :: iperm(:)
262 integer(kind=kint), intent(in) :: indexx(0:)
263 integer(kind=kint), intent(in) :: itemx(:)
264 real(kind=kreal), intent(in) :: ax(:)
265 integer(kind=kint), intent(in) :: indexlp(0:), indexup(0:)
266 integer(kind=kint), intent(in) :: itemlp(:), itemup(:)
267 real(kind=kreal), intent(out) :: alp(:), aup(:)
268 integer(kind=kint) :: ndof2, iold, inew
269 integer(kind=kint) :: jsnewl, jenewl, jsnewu, jenewu
270 integer(kind=kint) :: jold, kold, knew, jnew, l0old, l0new, l
271 ndof2 = ndof*ndof
272 ! new L
273 !$omp parallel default(none), &
274 !$omp private(iold,inew,jsnewL,jenewL,jsnewU,jenewU,jold,kold,knew,jnew,l0old,l0new,l), &
275 !$omp shared(N,iperm,indexLp,indexUp,indexX,itemX,itemLp,NDOF2,ALp,AX,itemUp,AUp)
276 !$omp do
277 do iold=1,n
278 inew = iperm(iold)
279 jsnewl = indexlp(inew-1)+1
280 jenewl = indexlp(inew)
281 jsnewu = indexup(inew-1)+1
282 jenewu = indexup(inew)
283 do jold = indexx(iold-1)+1, indexx(iold)
284 kold = itemx(jold)
285 if (kold > n) cycle
286 knew = iperm(kold)
287 if (knew < inew) then
288 call hecmw_bsearch_int_array(itemlp, jsnewl, jenewl, knew, jnew)
289 if (jnew < 0) then
290 write(0,*) 'ERROR:: jnew < 0 in reorder_off_diag2'
292 end if
293 l0old = (jold-1)*ndof2
294 l0new = (jnew-1)*ndof2
295 do l=1,ndof2
296 alp(l0new + l) = ax(l0old + l)
297 end do
298 else
299 call hecmw_bsearch_int_array(itemup, jsnewu, jenewu, knew, jnew)
300 if (jnew < 0) then
301 write(0,*) 'ERROR:: jnew < 0 in reorder_off_diag2'
303 end if
304 l0old = (jold-1)*ndof2
305 l0new = (jnew-1)*ndof2
306 do l=1,ndof2
307 aup(l0new + l) = ax(l0old + l)
308 end do
309 end if
310 end do
311 end do
312 !$omp end do
313 !$omp end parallel
314 end subroutine reorder_off_diag2
315
316end module hecmw_matrix_reorder
subroutine, public hecmw_bsearch_int_array(array, istart, iend, val, idx)
recursive subroutine, public hecmw_qsort_int_array(array, istart, iend)
subroutine, public hecmw_matrix_reorder_renum_item(n, perm, indexxp, itemxp)
subroutine, public hecmw_matrix_reorder_values(n, ndof, perm, iperm, indexl, indexu, iteml, itemu, al, au, d, indexlp, indexup, itemlp, itemup, alp, aup, dp)
subroutine, public hecmw_matrix_reorder_back_vector(n, ndof, perm, xp, x)
subroutine, public hecmw_matrix_reorder_profile(n, perm, iperm, indexl, indexu, iteml, itemu, indexlp, indexup, itemlp, itemup)
subroutine, public hecmw_matrix_reorder_vector(n, ndof, perm, x, xp)
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)