MODFLOW 6  version 6.6.0.dev0
USGS Modular Hydrologic Model
Sparse.f90
Go to the documentation of this file.
2 !cdl (8/12/2013) sparse matrix module for managing the structure
3 !of a matrix. Module uses FORTRAN 2003 extensions to manage
4 !the data structures in an object oriented fashion.
5 
6  use kindmodule, only: dp, i4b, lgp
7  implicit none
8 
9  type rowtype
10  integer(I4B) :: nnz ! number of nonzero entries in the row
11  integer(I4B), allocatable, dimension(:) :: icolarray ! array of column numbers
12  end type rowtype
13 
14  type, public :: sparsematrix
15  integer(I4B) :: offset !< global offset for first row in this matrix (default = 0)
16  integer(I4B) :: nrow !< number of rows in the matrix
17  integer(I4B) :: ncol !< number of columns in the matrix
18  integer(I4B) :: nnz !< number of nonzero matrix entries
19  integer(I4B) :: nnz_od !< number of off-diagonal nonzero matrix entries
20  type(rowtype), allocatable, dimension(:) :: row !< one rowtype for each matrix row
21  contains
23  procedure :: addconnection
24  procedure :: filliaja
25  procedure :: sort
26  procedure :: destroy
27 
28  procedure, private :: initializefixed
29  procedure, private :: initialize
30  end type sparsematrix
31 
32 contains
33 
34  subroutine destroy(this)
35  class(sparsematrix), intent(inout) :: this
36  deallocate (this%row)
37  end subroutine destroy
38 
39  subroutine initialize(this, nrow, ncol, rowmaxnnz)
40  !initial the sparse matrix. This subroutine
41  !acts a method for a sparse matrix by initializing
42  !the row data. It presently requires one maximum
43  !value for all rows, however, this can be changed
44  !to a vector of maximum values with one value for
45  !each row.
46  ! -- dummy
47  class(sparsematrix), intent(inout) :: this
48  integer(I4B), intent(in) :: nrow, ncol
49  integer(I4B), intent(in), dimension(nrow) :: rowmaxnnz
50  ! -- local
51  integer(I4B) :: i
52  ! -- code
53  this%offset = 0
54  this%nrow = nrow
55  this%ncol = ncol
56  this%nnz = 0
57  this%nnz_od = 0
58  allocate (this%row(nrow))
59  do i = 1, nrow
60  allocate (this%row(i)%icolarray(rowmaxnnz(i)))
61  this%row(i)%icolarray = 0
62  this%row(i)%nnz = 0
63  end do
64  end subroutine initialize
65 
66  ! overload
67  subroutine initializefixed(this, nrow, ncol, maxnnz)
68  implicit none
69  class(sparsematrix), intent(inout) :: this
70  integer(I4B), intent(in) :: nrow, ncol
71  integer(I4B), intent(in) :: maxnnz
72  ! local
73  integer(I4B), dimension(:), allocatable :: rowmaxnnz
74  integer(I4B) :: i
75 
76  allocate (rowmaxnnz(nrow))
77 
78  do i = 1, nrow
79  rowmaxnnz(i) = maxnnz
80  end do
81 
82  call this%initialize(nrow, ncol, rowmaxnnz)
83  deallocate (rowmaxnnz)
84 
85  end subroutine initializefixed
86 
87  subroutine filliaja(this, ia, ja, ierror, sort)
88  !allocate and fill the ia and ja arrays using information
89  !from the sparsematrix.
90  !ierror is returned as:
91  ! 0 if no error
92  ! 1 if ia is not the correct size
93  ! 2 if ja is not the correct size
94  ! 3 if both ia and ja are not correct size
95  ! -- dummy
96  class(sparsematrix), intent(inout) :: this
97  integer(I4B), dimension(:), intent(inout) :: ia, ja
98  integer(I4B), intent(inout) :: ierror
99  logical, intent(in), optional :: sort
100  ! -- local
101  logical :: sortja
102  integer(I4B) :: i, j, ipos
103  ! -- code
104  !
105  ! -- process optional dummy variables
106  if (present(sort)) then
107  sortja = sort
108  else
109  sortja = .false.
110  end if
111  !
112  ! -- initialize error variable
113  ierror = 0
114  !
115  ! -- check for error conditions
116  if (ubound(ia, dim=1) /= this%nrow + 1) then
117  ierror = 1
118  end if
119  if (ubound(ja, dim=1) /= this%nnz) then
120  ierror = ierror + 2
121  end if
122  if (ierror /= 0) then
123  return
124  end if
125  !
126  ! -- sort this
127  if (sortja) then
128  call this%sort()
129  end if
130  !
131  ! -- fill ia and ja
132  ipos = 1
133  ia(1) = ipos
134  do i = 1, this%nrow
135  do j = 1, this%row(i)%nnz
136  ja(ipos) = this%row(i)%icolarray(j)
137  ipos = ipos + 1
138  end do
139  ia(i + 1) = ipos
140  end do
141  end subroutine filliaja
142 
143  subroutine addconnection(this, i, j, inodup, iaddop)
144  !add a connection to the sparsematrix. if inodup
145  !(for no duplicates) is 1, then j is added only
146  !if it is unique.
147  ! -- dummy
148  class(sparsematrix), intent(inout) :: this
149  integer(I4B), intent(in) :: i, j, inodup
150  integer(I4B), optional, intent(inout) :: iaddop
151  ! -- local
152  integer(I4B) :: irow_local
153  integer(I4B) :: iadded
154  ! -- code
155  !
156  ! -- when distributed system, reduce row numbers to local range
157  irow_local = i - this%offset
158  !
159  call insert(j, this%row(irow_local), inodup, iadded)
160  this%nnz = this%nnz + iadded
161  if (j < this%offset + 1 .or. j > this%offset + this%nrow) then
162  ! count the off-diagonal entries separately
163  this%nnz_od = this%nnz_od + iadded
164  end if
165  if (present(iaddop)) iaddop = iadded
166  end subroutine addconnection
167 
168  subroutine insert(j, thisrow, inodup, iadded)
169  !insert j into thisrow (for row i)
170  !inodup=1 means do not include duplicate connections
171  !iadded is 1 if a new entry was added (meaning that nnz for the row was increased)
172  !iadded is 0 if duplicate and not added. Used to track total number of connections
173  ! -- dummy
174  integer(I4B), intent(in) :: j, inodup
175  type(rowtype), intent(inout) :: thisrow
176  integer(I4B), allocatable, dimension(:) :: iwk
177  integer(I4B), intent(inout) :: iadded
178  ! -- local
179  integer(I4B) :: jj, maxnnz
180  ! -- code
181  iadded = 0
182  maxnnz = ubound(thisrow%icolarray, dim=1)
183  if (thisrow%icolarray(1) == 0) then
184  thisrow%icolarray(1) = j
185  thisrow%nnz = thisrow%nnz + 1
186  iadded = 1
187  return
188  end if
189  if (thisrow%nnz == maxnnz) then
190  ! -- increase size of the row
191  allocate (iwk(thisrow%nnz))
192  iwk = thisrow%icolarray
193  deallocate (thisrow%icolarray)
194  ! -- Specify how to increase the size of the icolarray. Adding 1
195  ! will be most memory conservative, but may be a little slower
196  ! due to frequent allocate/deallocate. Another option would be
197  ! to double the size: maxnnz=maxnnz*2
198  maxnnz = maxnnz + 1
199  allocate (thisrow%icolarray(maxnnz))
200  thisrow%icolarray(1:thisrow%nnz) = iwk(1:thisrow%nnz)
201  thisrow%icolarray(thisrow%nnz + 1:maxnnz) = 0
202  end if
203  if (inodup == 1) then
204  do jj = 1, thisrow%nnz
205  if (thisrow%icolarray(jj) == j) then
206  return
207  end if
208  end do
209  end if
210  !
211  ! -- add the connection to end
212  thisrow%nnz = thisrow%nnz + 1
213  thisrow%icolarray(thisrow%nnz) = j
214  iadded = 1
215  end subroutine insert
216 
217  subroutine sort(this, with_csr)
218  !sort the icolarray for each row, but do not include
219  !the diagonal position in the sort so that it stays in front
220  ! -- dummy
221  class(sparsematrix), intent(inout) :: this
222  logical(LGP), optional :: with_csr
223  ! -- local
224  integer(I4B) :: i, nval, start_idx
225  ! -- code
226  start_idx = 2
227  if (present(with_csr)) then
228  if (with_csr) then
229  ! CSR: don't put diagonal up front
230  start_idx = 1
231  end if
232  end if
233 
234  do i = 1, this%nrow
235  nval = this%row(i)%nnz
236  call sortintarray(nval - start_idx + 1, &
237  this%row(i)%icolarray(start_idx:nval))
238  end do
239  end subroutine sort
240 
241  subroutine sortintarray(nval, iarray)
242  !simple subroutine for sorting an array
243  !in place. It is not the fastest sort function
244  !but should suffice for relatively short nodelists.
245  ! -- dummy
246  integer(I4B), intent(in) :: nval
247  integer(I4B), intent(inout), dimension(nval) :: iarray
248  ! -- local
249  integer(I4B) :: i, j, itemp
250  ! -- code
251  do i = 1, nval - 1
252  do j = i + 1, nval
253  if (iarray(i) > iarray(j)) then
254  itemp = iarray(j)
255  iarray(j) = iarray(i)
256  iarray(i) = itemp
257  end if
258  end do
259  end do
260  end subroutine sortintarray
261 
262  subroutine csr_diagsum(ia, flowja)
263  !Add up the off diagonal terms and put the sum in the
264  !diagonal position
265  ! -- dummy
266  integer(I4B), dimension(:), contiguous :: ia
267  real(DP), dimension(:), contiguous :: flowja
268  ! -- local
269  integer(I4B) :: nodes
270  integer(I4B) :: n
271  integer(I4B) :: iposdiag
272  integer(I4B) :: ipos
273  ! -- code
274  nodes = size(ia) - 1
275  do n = 1, nodes
276  iposdiag = ia(n)
277  do ipos = ia(n) + 1, ia(n + 1) - 1
278  flowja(iposdiag) = flowja(iposdiag) + flowja(ipos)
279  end do
280  end do
281  end subroutine csr_diagsum
282 
283 end module sparsemodule
subroutine init()
Definition: GridSorting.f90:24
This module defines variable data types.
Definition: kind.f90:8
subroutine sort(this, with_csr)
Definition: Sparse.f90:218
subroutine csr_diagsum(ia, flowja)
Definition: Sparse.f90:263
subroutine destroy(this)
Definition: Sparse.f90:35
subroutine filliaja(this, ia, ja, ierror, sort)
Definition: Sparse.f90:88
subroutine initialize(this, nrow, ncol, rowmaxnnz)
Definition: Sparse.f90:40
subroutine initializefixed(this, nrow, ncol, maxnnz)
Definition: Sparse.f90:68
subroutine sortintarray(nval, iarray)
Definition: Sparse.f90:242
subroutine addconnection(this, i, j, inodup, iaddop)
Definition: Sparse.f90:144
subroutine insert(j, thisrow, inodup, iadded)
Definition: Sparse.f90:169