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