FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
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 
6 module m_crs_matrix
7  use hecmw_util
8 
9  ! compact row storage matrix
10 
11  private
12 
13  public crs_matrix
14  public irjctocrs
15  public symbolicirjctocrs
16  public crs_matrix_getvec
17 
18  type crs_matrix
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 
28 contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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 
325 end module m_crs_matrix
m_crs_matrix::irjctocrs
subroutine, public irjctocrs(ndeg, nttbr, irow, jcol, val, ncols, nrows, c)
Definition: m_crs_matrix.f90:82
m_crs_matrix::symbolicirjctocrs
subroutine, public symbolicirjctocrs(ndeg, nttbr, irow, jcol, ncols, nrows, c)
Definition: m_crs_matrix.f90:30
m_fstr::idbg
integer(kind=kint), parameter idbg
Definition: m_fstr.f90:111
hecmw_util
I/O and Utility.
Definition: hecmw_util_f.F90:7
m_crs_matrix
Definition: m_crs_matrix.f90:6
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
m_crs_matrix::crs_matrix_getvec
subroutine, public crs_matrix_getvec(c, k, v)
Definition: m_crs_matrix.f90:294