FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
m_crs_matrix.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
9 ! compact row storage matrix
10
11 private
12
13 public crs_matrix
14 public irjctocrs
17
19 integer(kind=kint) :: nrows ! number of rows in whole matrix.
20 integer(kind=kint) :: ncols ! number of columns in whole matrix.
21 integer(kind=kint) :: nttbr ! number of non zero elements in whole matrix.
22 integer(kind=kint) :: ndeg ! degree of freedom of each element
23 integer(kind=kint), pointer :: ia(:) ! size is ia(nrows+1); ia(k)..ia(k+1)-1 elements in ja(:) show columns of k'th row in whole matrix.
24 integer(kind=kint), pointer :: ja(:) ! size is ja(nttbr); store column indices of non zero elements.
25 real(kind=kreal), pointer :: val(:,:) ! size is val(ndeg*ndeg, nttbr). store matrix elements.
26 end type crs_matrix
27
28contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
29 subroutine symbolicirjctocrs(ndeg, nttbr, irow, jcol, ncols, nrows, c)
30
31 implicit none
32
33 integer(kind=kint), intent(in) :: ndeg ! degree of freedom of each element
34 integer(kind=kint), intent(in) :: nttbr ! number of non zero elements in whole matrix.
35 integer(kind=kint), intent(in) :: irow(:) ! irjc matrix style row number of element
36 integer(kind=kint), intent(in) :: jcol(:) ! irjc matrix style column number of element
37 integer(kind=kint), intent(in) :: ncols ! C column size
38 integer(kind=kint), intent(in) :: nrows ! C row size
39
40 type (crs_matrix), intent(out) :: c
41
42 ! internal !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
43 integer(kind=kint), allocatable :: istat(:), jrapt(:), jcolno(:)
44 integer(kind=kint) :: ntt
45
46 integer(kind=kint) :: ipass, i,j,k,l,m,n
47
48 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
49
50 allocate(istat(nrows))
51 allocate(jrapt(nttbr))
52 allocate(jcolno(nttbr))
53
54 istat=0
55 jrapt=0
56 jcolno=0
57 ntt=0
58
59 do l=1,nttbr
60 i=irow(l)
61 j=jcol(l)
62 call reserv(i,j,istat,jrapt,jcolno,ntt)
63 end do
64
65 ! allocation and set c%ia, c%ja
66 c%nrows=nrows
67 c%ncols=ncols
68 c%nttbr=ntt
69 c%ndeg=ndeg
70
71 call crsallocation(c)
72
73 call stiajac(c,istat,jrapt,jcolno)
74
75 deallocate(istat,jrapt,jcolno)
76
77 end subroutine symbolicirjctocrs
78
79 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
80
81 subroutine irjctocrs(ndeg, nttbr, irow, jcol, val, ncols, nrows, c)
82
83 implicit none
84
85 integer(kind=kint), intent(in) :: ndeg ! degree of freedom of each element
86 integer(kind=kint), intent(in) :: nttbr ! number of non zero elements in whole matrix.
87 integer(kind=kint), intent(in) :: irow(:) ! irjc matrix style row number of element
88 integer(kind=kint), intent(in) :: jcol(:) ! irjc matrix style column number of element
89 real(kind=kreal), intent(in) :: val(:,:) ! store matrix element
90 integer(kind=kint), intent(in) :: ncols ! C column size
91 integer(kind=kint), intent(in) :: nrows ! C row size
92
93 type (crs_matrix), intent(out) :: c
94
95 ! internal !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
96 integer(kind=kint), allocatable :: istat(:), jrapt(:), jcolno(:)
97 integer(kind=kint) :: ntt
98
99 integer(kind=kint) :: ipass, i,j,k,l,m,n
100
101 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
102
103 allocate(istat(nrows))
104 allocate(jrapt(nttbr))
105 allocate(jcolno(nttbr))
106
107 istat=0
108 jrapt=0
109 jcolno=0
110 ntt=0
111
112 do ipass=1,2
113 do l=1,nttbr
114 i=irow(l)
115 j=jcol(l)
116
117 if(ipass.eq.1)then
118 call reserv(i,j,istat,jrapt,jcolno,ntt)
119 endif
120 if(ipass.eq.2)then
121 call stval(c,i,j,val(:,l),0)
122 endif
123 end do
124
125 ! allocation and set c%ia, c%ja
126 if(ipass.eq.1)then
127 c%nrows=nrows
128 c%ncols=ncols
129 c%nttbr=ntt
130 c%ndeg=ndeg
131
132 call crsallocation(c)
133
134 call stiajac(c,istat,jrapt,jcolno)
135
136 deallocate(istat,jrapt,jcolno)
137 endif
138 end do
139
140 end subroutine irjctocrs
141
142 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
143
144 subroutine reserv(i,j,istat,jrapt,jcolno,k)
145 ! count number of non zero elements
146 implicit none
147 integer(kind=kint) :: jrapt(:)
148 integer(kind=kint) :: jcolno(:)
149 integer(kind=kint) :: istat(:)
150 integer(kind=kint) :: i,j,k,l, locr, loc
151 locr=0
152 loc=istat(i)
153 110 continue
154 if(loc.eq.0) goto 120
155 if(jcolno(loc).eq.j) then
156 goto 100
157 elseif(jcolno(loc).gt.j) then
158 goto 130
159 endif
160 locr=loc
161 loc=jrapt(loc)
162 goto 110
163 120 continue
164 k=k+1
165 if(locr.eq.0) then
166 istat(i)=k
167 else
168 jrapt(locr)=k
169 endif
170 jcolno(k)=j
171 goto 150
172 130 continue
173 k=k+1
174 if(locr.eq.0) then
175 istat(i)=k
176 else
177 jrapt(locr)=k
178 endif
179 jrapt(k)=loc
180 jcolno(k)=j
181 150 continue
182 100 continue
183 return
184 end subroutine reserv
185
186 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
187
188 subroutine crsallocation(c)
189 implicit none
190 type (crs_matrix) :: c
191
192 if (0 >= c%nrows) then
193 call errtrp('set positive nrows')
194 else if (0 >= c%ndeg) then
195 call errtrp('set positive ndeg')
196 else if (0 >= c%nttbr) then
197 call errtrp('set positive nttbr')
198 end if
199
200 allocate(c%ia(c%nrows+1))
201 allocate(c%ja(c%nttbr))
202 allocate(c%val(c%ndeg*c%ndeg,c%nttbr))
203
204 return
205 end subroutine crsallocation
206
207 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
208
209 subroutine stiajac(c,istat,jrapt,jcolno)
210 implicit none
211 ! set ia, ja
212 type (crs_matrix) :: c
213 integer(kind=kint) :: istat(:)
214 integer(kind=kint) :: jrapt(:)
215 integer(kind=kint) :: jcolno(:)
216 integer(kind=kint) :: i,j,k,l, ii, loc, idbg
217
218 if (0 >= c%ncols) then
219 call errtrp('set positive ncols')
220 else if (0 >= c%nrows) then
221 call errtrp('set positive nrows')
222 else if (0 >= c%nttbr) then
223 call errtrp('set positive nttbr')
224 else if (.not. associated(c%ia)) then
225 call errtrp('ia is not allocated')
226 else if (.not. associated(c%ja)) then
227 call errtrp('ja is not allocated')
228 else if (c%nrows+1 /= size(c%ia)) then
229 call errtrp('ia size unmatched with nrows')
230 else if (c%nttbr /= size(c%ja)) then
231 call errtrp('ja size unmatched with nttbr')
232 else if ((c%ndeg*c%ndeg /= size(c%val, dim=1)) .or. (c%nttbr /= size(c%val, dim=2))) then
233 call errtrp('ja size unmatched with ndeg, nttbr')
234 end if
235
236 idbg=0
237 c%ia(1)=1
238 l=0
239 do 100 k=1,c%nrows
240 loc=istat(k)
241 110 continue
242 if(loc.eq.0) goto 120
243 ii=jcolno(loc)
244 l=l+1
245 c%ja(l)=ii
246 loc=jrapt(loc)
247 goto 110
248 120 c%ia(k+1)=l+1
249 100 continue
250 if(idbg.ne.0) then
251 write(20,*) 'c%ia '
252 write(20,60) (c%ia(i),i=1,c%nrows+1) ! ; call flush(20)
253 write(20,*) 'c%ja '
254 write(20,60) (c%ja(i),i=1,c%ja(c%nrows+1))
255 end if
256 60 format(10i7)
257 return
258 end subroutine stiajac
259
260 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
261
262 subroutine stval(c,i,j,val,itrans)
263 ! store matrix element
264 implicit none
265
266 type (crs_matrix) :: c
267 real(kind=kreal), dimension(:) :: val
268 integer(kind=kint) :: ndeg,itrans
269 integer(kind=kint) :: i,j,k,m,n
270
271 ndeg=c%ndeg
272
273 do k=c%ia(i),c%ia(i+1)-1
274 if(j.eq.c%ja(k)) then
275 if(itrans.eq.0) then
276 c%val(:,k)=val(:)
277 else
278 do m=1,ndeg
279 do n=1,ndeg
280 c%val(m + (n-1)*ndeg,k)=val((m-1)*ndeg + n)
281 end do
282 end do
283 end if
284 return
285 end if
286 end do
287 write(6,*) "something wrong"
288 stop
289 end subroutine stval
290
291 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
292
293 subroutine crs_matrix_getvec(c, k, v)
294 ! return k'th row vector in whole matrix
295 implicit none
296 type (crs_matrix), intent(in) :: c
297 integer(kind=kint), intent(in) :: k ! given in "ndeg opened" numbering.
298 real(kind=kreal), dimension(:), intent(out) :: v ! output as "ndeg opened" dens vector
299
300 integer(kind=kint) :: ndeg, nndeg, i, idx, jcol, iofset
301
302 ndeg=c%ndeg
303 nndeg=ndeg*ndeg
304 v=0
305 jcol = (k+ndeg-1) / ndeg ! row number in sparse matrix
306 iofset = mod(k+ndeg-1, ndeg) ! offset in val. 0offset
307
308 do i=c%ia(jcol),c%ia(jcol+1)-1
309 idx=c%ja(i) ! row number in sparse matrix
310 v(ndeg*(idx-1)+1:ndeg*(idx-1)+ndeg) &
311 & =c%val(ndeg*iofset + 1 : ndeg*iofset + ndeg ,i)
312 end do
313
314 return
315 end subroutine crs_matrix_getvec
316
317 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
318
319 subroutine errtrp(mes)
320 character(*) mes
321 write(*,*) 'Error in m_crs_matrix: ', mes
322 stop
323 end subroutine errtrp
324
325end module m_crs_matrix
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=4), parameter kreal
subroutine, public crs_matrix_getvec(c, k, v)
subroutine, public symbolicirjctocrs(ndeg, nttbr, irow, jcol, ncols, nrows, c)
subroutine, public irjctocrs(ndeg, nttbr, irow, jcol, val, ncols, nrows, c)