FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
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 !-------------------------------------------------------------------------------
5 module m_crs_matrix_lag
6  use my_hecmw_util_lag
7 
8  ! compact row storage matrix
9 
10  private
11 
12  public crs_matrix
13  public irjctocrs
14  public symbolicirjctocrs
15  public crs_matrix_getvec
16 
17  type crs_matrix
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 
27 contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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 
324 end module m_crs_matrix_lag
m_crs_matrix_lag::irjctocrs
subroutine, public irjctocrs(ndeg, nttbr, irow, jcol, val, ncols, nrows, c)
Definition: m_crs_matrix_lag.f90:81
my_hecmw_util_lag
Definition: my_hecmw_util_f_lag.f90:6
m_fstr::idbg
integer(kind=kint), parameter idbg
Definition: m_fstr.f90:111
m_crs_matrix_lag::symbolicirjctocrs
subroutine, public symbolicirjctocrs(ndeg, nttbr, irow, jcol, ncols, nrows, c)
Definition: m_crs_matrix_lag.f90:29
hecmw_util::kreal
integer(kind=4), parameter kreal
Definition: hecmw_util_f.F90:16
m_crs_matrix_lag::crs_matrix_getvec
subroutine, public crs_matrix_getvec(c, k, v)
Definition: m_crs_matrix_lag.f90:293
m_crs_matrix_lag
Definition: m_crs_matrix_lag.f90:5