FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
m_cclsmatrix.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! compact column storage matrix
8
10 integer :: neqns ! number of rows in whole matrix.
11 integer :: ncol ! number of columns in whole matrix.
12 integer :: nttbr ! number of non zero elements in whole matrix.
13 integer :: ndeg ! degree of freedom of each element
14 integer, pointer :: ia(:) ! size is ia(ncol+1). ia(k)..ia(k+1)-1 elements in ja(:) show rows of k'th column in whole matrix.
15 integer, pointer :: ja(:) ! size is ja(nttbr). store row indices of non zero elements.
16 real(8), pointer :: val(:,:) ! size is val(ndeg*ndeg, nttbr). store matrix elements.
17end type ccls_matrix
18
19contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20 subroutine reserv(i,j,jstat,irpt,irowno,k)
21 ! count number of non zero elements
22 implicit none
23! integer irpt(*),irowno(*),jstat(*)
24 integer, dimension(:) :: irpt
25 integer, dimension(:) :: irowno
26 integer, dimension(:) :: jstat
27 integer i,j,k,l, locr, loc
28 locr=0
29 loc=jstat(j)
30 110 continue
31 if(loc.eq.0) goto 120
32 if(irowno(loc).eq.i) then
33 goto 100
34 elseif(irowno(loc).gt.i) then
35 goto 130
36 endif
37 locr=loc
38 loc=irpt(loc)
39 goto 110
40 120 continue
41 k=k+1
42 if(locr.eq.0) then
43 jstat(j)=k
44 else
45 irpt(locr)=k
46 endif
47 irowno(k)=i
48 goto 150
49 130 continue
50 k=k+1
51 if(locr.eq.0) then
52 jstat(j)=k
53 else
54 irpt(locr)=k
55 endif
56 irpt(k)=loc
57 irowno(k)=i
58 150 continue
59 100 continue
60 return
61 end subroutine reserv
62
63 subroutine cclsallocation(c)
64 implicit none
65 type (ccls_matrix) :: c
66
67 if (0 >= c%ncol) then
68 call m_cclsmatrix_errtrp('set positive ncol')
69 else if (0 >= c%ndeg) then
70 call m_cclsmatrix_errtrp('set positive ndeg')
71 else if (0 >= c%nttbr) then
72 call m_cclsmatrix_errtrp('set positive nttbr')
73 end if
74
75 allocate(c%ia(c%ncol+1))
76 allocate(c%ja(c%nttbr))
77 allocate(c%val(c%ndeg*c%ndeg,c%nttbr))
78
79 return
80 end subroutine cclsallocation
81
82 subroutine stiajac(c,jstat,irpt,irowno)
83 implicit none
84! set ia, ja
85 type (ccls_matrix) :: c
86! integer jstat(*),irpt(*),irowno(*)
87 integer, dimension(:) :: jstat
88 integer, dimension(:) :: irpt
89 integer, dimension(:) :: irowno
90 integer i,j,k,l, ii, loc, idbg
91
92 if (0 >= c%neqns) then
93 call m_cclsmatrix_errtrp('set positive neqns')
94 else if (0 >= c%ncol) then
95 call m_cclsmatrix_errtrp('set positive ncol')
96 else if (0 >= c%nttbr) then
97 call m_cclsmatrix_errtrp('set positive nttbr')
98 else if (.false. == associated(c%ia)) then
99 call m_cclsmatrix_errtrp('ia is not allocated')
100 else if (.false. == associated(c%ja)) then
101 call m_cclsmatrix_errtrp('ja is not allocated')
102 else if (c%ncol+1 /= size(c%ia)) then
103 call m_cclsmatrix_errtrp('ia size unmatched with ncol')
104 else if (c%nttbr /= size(c%ja)) then
105 call m_cclsmatrix_errtrp('ja size unmatched with nttbr')
106 else if ((c%ndeg*c%ndeg /= size(c%val, dim=1)) .or. (c%nttbr /= size(c%val, dim=2))) then
107 call m_cclsmatrix_errtrp('ja size unmatched with ndeg, nttbr')
108 end if
109
110 idbg=0
111 c%ia(1)=1
112 l=0
113 do 100 k=1,c%ncol
114 loc=jstat(k)
115 110 continue
116 if(loc.eq.0) goto 120
117 ii=irowno(loc)
118 l=l+1
119 c%ja(l)=ii
120 loc=irpt(loc)
121 goto 110
122 120 c%ia(k+1)=l+1
123 100 continue
124 if(idbg.ne.0) then
125 write(20,*) 'c%ia '
126 write(20,60) (c%ia(i),i=1,c%ncol+1) ! ; call flush(20)
127 write(20,*) 'c%ja '
128 write(20,60) (c%ja(i),i=1,c%ja(c%ncol+1))
129 end if
130 60 format(10i7)
131 return
132 end subroutine stiajac
133
134subroutine stval(c,i,j,val,itrans)
135! store matrix element
136implicit none
137
138type (ccls_matrix) :: c
139real(8), dimension(c%ndeg*c%ndeg) :: val
140integer ndeg,itrans
141integer i,j,k,m,n
142
143ndeg=c%ndeg
144
145do k=c%ia(j),c%ia(j+1)-1
146 if(i.eq.c%ja(k)) then
147 if(itrans.eq.0) then
148 c%val(:,k)=val(:)
149 else
150 do m=1,ndeg
151 do n=1,ndeg
152 c%val(m + (n-1)*ndeg,k)=val((m-1)*ndeg + n)
153 end do
154 end do
155 end if
156 return
157 end if
158end do
159write(6,*) "something wrong"
160stop
161end subroutine stval
162
163subroutine m_cclsmatrix_getvec(c, k, v)
164! return k'th column vector in whole matrix
165implicit none
166type (ccls_matrix), intent(in) :: c
167integer, intent(in) :: k ! given in "ndeg opened" numbering.
168!real(8), dimension(c%ndeg*c%neqns), intent(out) :: v ! output as "ndeg opened" dens vector
169real(8), dimension(:), intent(out) :: v ! output as "ndeg opened" dens vector
170
171integer ndeg, nndeg, i, idx, jcol, iofset
172
173ndeg=c%ndeg
174nndeg=ndeg*ndeg
175v=0
176jcol = (k+ndeg-1) / ndeg ! column number in sparse matrix
177iofset = mod(k+ndeg-1, ndeg) ! offset in val. 0offset
178
179do i=c%ia(jcol),c%ia(jcol+1)-1
180 idx=c%ja(i) ! row number in sparse matrix
181 v(ndeg*(idx-1)+1:ndeg*(idx-1)+ndeg) &
182& =c%val(ndeg*iofset + 1 : ndeg*iofset + ndeg ,i)
183end do
184
185return
186end subroutine m_cclsmatrix_getvec
187
188subroutine m_cclsmatrix_errtrp(mes)
189character(*) mes
190write(*,*) 'Error in m_cclsmatrix: ', mes
191stop
192end subroutine m_cclsmatrix_errtrp
193
194end module m_cclsmatrix
subroutine cclsallocation(c)
subroutine m_cclsmatrix_getvec(c, k, v)
subroutine stval(c, i, j, val, itrans)
subroutine m_cclsmatrix_errtrp(mes)
subroutine reserv(i, j, jstat, irpt, irowno, k)
subroutine stiajac(c, jstat, irpt, irowno)