MODFLOW 6  version 6.6.0.dev0
USGS Modular Hydrologic Model
sparsekit.f90 File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine amux (n, x, y, a, ja, ia)
 
subroutine dvperm (n, x, perm)
 
subroutine cperm (nrow, a, ja, ia, ao, jao, iao, perm, job)
 
subroutine rperm (nrow, a, ja, ia, ao, jao, iao, perm, job)
 
subroutine dperm (nrow, a, ja, ia, ao, jao, iao, perm, qperm, job)
 

Function/Subroutine Documentation

◆ amux()

subroutine amux ( integer ( kind = 4 )  n,
real ( kind = 8 ), dimension(*)  x,
real ( kind = 8 ), dimension(n)  y,
real ( kind = 8 ), dimension(*)  a,
integer ( kind = 4 ), dimension(*)  ja,
integer ( kind = 4 ), dimension(*)  ia 
)

Definition at line 1 of file sparsekit.f90.

2 
3  !*****************************************************************************80
4  !
5  !! AMUX multiplies a CSR matrix A times a vector.
6  !
7  ! Discussion:
8  !
9  ! This routine multiplies a matrix by a vector using the dot product form.
10  ! Matrix A is stored in compressed sparse row storage.
11  !
12  ! Modified:
13  !
14  ! 07 January 2004
15  !
16  ! Author:
17  !
18  ! Youcef Saad
19  !
20  ! Parameters:
21  !
22  ! Input, integer ( kind = 4 ) N, the row dimension of the matrix.
23  !
24  ! Input, real X(*), and array of length equal to the column dimension
25  ! of A.
26  !
27  ! Input, real A(*), integer ( kind = 4 ) JA(*), IA(NROW+1), the matrix in CSR
28  ! Compressed Sparse Row format.
29  !
30  ! Output, real Y(N), the product A * X.
31  !
32  implicit none
33 
34  integer ( kind = 4 ) n
35 
36  real ( kind = 8 ) a(*)
37  integer ( kind = 4 ) i
38  integer ( kind = 4 ) ia(*)
39  integer ( kind = 4 ) ja(*)
40  integer ( kind = 4 ) k
41  real ( kind = 8 ) t
42  real ( kind = 8 ) x(*)
43  real ( kind = 8 ) y(n)
44 
45  do i = 1, n
46  !
47  ! Compute the inner product of row I with vector X.
48  !
49  t = 0.0d+00
50  do k = ia(i), ia(i+1)-1
51  t = t + a(k) * x(ja(k))
52  end do
53 
54  y(i) = t
55 
56  end do
57 
58  return
Here is the caller graph for this function:

◆ cperm()

subroutine cperm ( integer ( kind = 4 )  nrow,
real ( kind = 8 ), dimension(*)  a,
integer ( kind = 4 ), dimension(*)  ja,
integer ( kind = 4 ), dimension(nrow+1)  ia,
real ( kind = 8 ), dimension(*)  ao,
integer ( kind = 4 ), dimension(*)  jao,
integer ( kind = 4 ), dimension(nrow+1)  iao,
integer ( kind = 4 ), dimension(*)  perm,
integer ( kind = 4 )  job 
)

Definition at line 163 of file sparsekit.f90.

164 
165  !*****************************************************************************80
166  !
167  !! CPERM permutes the columns of a matrix.
168  !
169  ! Discussion:
170  !
171  ! This routine permutes the columns of a matrix a, ja, ia.
172  ! The result is written in the output matrix ao, jao, iao.
173  ! cperm computes B = A P, where P is a permutation matrix
174  ! that maps column j into column perm(j), i.e., on return
175  ! a(i,j) becomes a(i,perm(j)) in new matrix
176  !
177  ! if job=1 then ao, iao are not used.
178  ! This routine is in place: ja, jao can be the same.
179  ! If the matrix is initially sorted (by increasing column number)
180  ! then ao,jao,iao may not be on return.
181  !
182  ! Modified:
183  !
184  ! 07 January 2004
185  !
186  ! Author:
187  !
188  ! Youcef Saad
189  !
190  ! Parameters:
191  !
192  ! Input, integer ( kind = 4 ) NROW, the row dimension of the matrix.
193  !
194  ! Input, real A(*), integer ( kind = 4 ) JA(*), IA(NROW+1), the matrix in CSR
195  ! Compressed Sparse Row format.
196  !
197  ! perm = integer ( kind = 4 ) array of length ncol (number of columns of A
198  ! containing the permutation array the columns:
199  ! a(i,j) in the original matrix becomes a(i,perm(j))
200  ! in the output matrix.
201  !
202  ! job = integer ( kind = 4 ) indicating the work to be done:
203  ! job = 1 permute a, ja, ia into ao, jao, iao
204  ! (including the copying of real values ao and
205  ! the array iao).
206  ! job /= 1 : ignore real values ao and ignore iao.
207  !
208  !
209  ! on return:
210  !
211  ! ao, jao, iao = input matrix in a, ja, ia format (array ao not needed)
212  !
213  implicit none
214 
215  integer ( kind = 4 ) nrow
216 
217  real ( kind = 8 ) a(*)
218  real ( kind = 8 ) ao(*)
219  integer ( kind = 4 ) ia(nrow+1)
220  integer ( kind = 4 ) iao(nrow+1)
221  integer ( kind = 4 ) ja(*)
222  integer ( kind = 4 ) jao(*)
223  integer ( kind = 4 ) job
224  integer ( kind = 4 ) k
225  integer ( kind = 4 ) nnz
226  integer ( kind = 4 ) perm(*)
227 
228  nnz = ia(nrow+1)-1
229 
230  do k = 1, nnz
231  jao(k) = perm(ja(k))
232  end do
233  !
234  ! Done with the JA array. Return if no need to touch values.
235  !
236  if ( job /= 1 ) then
237  return
238  end if
239  !
240  ! Else get new pointers, and copy values too.
241  !
242  iao(1:nrow+1) = ia(1:nrow+1)
243 
244  ao(1:nnz) = a(1:nnz)
245 
246  return
Here is the caller graph for this function:

◆ dperm()

subroutine dperm ( integer ( kind = 4 )  nrow,
real ( kind = 8 ), dimension(*)  a,
integer ( kind = 4 ), dimension(*)  ja,
integer ( kind = 4 ), dimension(nrow+1)  ia,
real ( kind = 8 ), dimension(*)  ao,
integer ( kind = 4 ), dimension(*)  jao,
integer ( kind = 4 ), dimension(nrow+1)  iao,
integer ( kind = 4 ), dimension(nrow)  perm,
integer ( kind = 4 ), dimension(nrow)  qperm,
integer ( kind = 4 )  job 
)

Definition at line 353 of file sparsekit.f90.

354 
355  !*****************************************************************************80
356  !
357  !! DPERM permutes the rows and columns of a matrix stored in CSR format.
358  !
359  ! Discussion:
360  !
361  ! This routine computes P*A*Q, where P and Q are permutation matrices.
362  ! P maps row i into row perm(i) and Q maps column j into column qperm(j).
363  ! A(I,J) becomes A(perm(i),qperm(j)) in the new matrix.
364  !
365  ! In the particular case where Q is the transpose of P (symmetric
366  ! permutation of A) then qperm is not needed.
367  ! note that qperm should be of length ncol (number of columns) but this
368  ! is not checked.
369  !
370  ! The algorithm is "in place".
371  !
372  ! The column indices may not be sorted on return even if they are
373  ! sorted on entry.
374  !
375  ! In case job == 2 or job == 4, a and ao are never referred to
376  ! and can be dummy arguments.
377  !
378  ! Modified:
379  !
380  ! 07 January 2004
381  !
382  ! Author:
383  !
384  ! Youcef Saad
385  !
386  ! Parameters:
387  !
388  ! Input, integer ( kind = 4 ) NROW, the order of the matrix.
389  !
390  ! Input, real A(*), integer ( kind = 4 ) JA(*), IA(NROW+1), the matrix in CSR
391  ! Compressed Sparse Row format.
392  !
393  ! Input, integer ( kind = 4 ) PERM(NROW), the permutation array for the rows: PERM(I)
394  ! is the destination of row I in the permuted matrix; also the destination
395  ! of column I in case the permutation is symmetric (JOB <= 2).
396  !
397  ! Input, integer ( kind = 4 ) QPERM(NROW), the permutation array for the columns.
398  ! This should be provided only if JOB=3 or JOB=4, that is, only in
399  ! the case of a nonsymmetric permutation of rows and columns.
400  ! Otherwise QPERM is a dummy argument.
401  !
402  ! job = integer ( kind = 4 ) indicating the work to be done:
403  ! * job = 1,2 permutation is symmetric Ao :== P * A * transp(P)
404  ! job = 1 permute a, ja, ia into ao, jao, iao
405  ! job = 2 permute matrix ignoring real values.
406  ! * job = 3,4 permutation is non-symmetric Ao :== P * A * Q
407  ! job = 3 permute a, ja, ia into ao, jao, iao
408  ! job = 4 permute matrix ignoring real values.
409  !
410  ! Output, real AO(*), JAO(*), IAO(NROW+1), the permuted matrix in CSR
411  ! Compressed Sparse Row format.
412  !
413  implicit none
414 
415  integer ( kind = 4 ) nrow
416 
417  real ( kind = 8 ) a(*)
418  real ( kind = 8 ) ao(*)
419  integer ( kind = 4 ) ia(nrow+1)
420  integer ( kind = 4 ) iao(nrow+1)
421  integer ( kind = 4 ) ja(*)
422  integer ( kind = 4 ) jao(*)
423  integer ( kind = 4 ) job
424  integer ( kind = 4 ) locjob
425  integer ( kind = 4 ) perm(nrow)
426  integer ( kind = 4 ) qperm(nrow)
427  !
428  ! LOCJOB indicates whether or not real values must be copied.
429  !
430  locjob = mod( job, 2 )
431  !
432  ! Permute the rows first.
433  !
434  call rperm ( nrow, a, ja, ia, ao, jao, iao, perm, locjob )
435  !
436  ! Permute the columns.
437  !
438  locjob = 0
439 
440  if ( job <= 2 ) then
441  call cperm ( nrow, ao, jao, iao, ao, jao, iao, perm, locjob )
442  else
443  call cperm ( nrow, ao, jao, iao, ao, jao, iao, qperm, locjob )
444  end if
445 
446  return
subroutine rperm(nrow, a, ja, ia, ao, jao, iao, perm, job)
Definition: sparsekit.f90:250
subroutine cperm(nrow, a, ja, ia, ao, jao, iao, perm, job)
Definition: sparsekit.f90:164
Here is the call graph for this function:
Here is the caller graph for this function:

◆ dvperm()

subroutine dvperm ( integer ( kind = 4 )  n,
real ( kind = 8 ), dimension(n)  x,
integer ( kind = 4 ), dimension(n)  perm 
)

Definition at line 61 of file sparsekit.f90.

62 
63  !*****************************************************************************80
64  !
65  !! DVPERM performs an in-place permutation of a real vector.
66  !
67  ! Discussion:
68  !
69  ! This routine permutes a real vector X using a permutation PERM.
70  !
71  ! On return, the vector X satisfies,
72  !
73  ! x(perm(j)) :== x(j), j = 1,2,.., n
74  !
75  ! Modified:
76  !
77  ! 07 January 2004
78  !
79  ! Author:
80  !
81  ! Youcef Saad
82  !
83  ! Parameters:
84  !
85  ! Input, integer ( kind = 4 ) N, the length of X.
86  !
87  ! Input/output, real X(N), the vector to be permuted.
88  !
89  ! Input, integer ( kind = 4 ) PERM(N), the permutation.
90  !
91  implicit none
92 
93  integer ( kind = 4 ) n
94 
95  integer ( kind = 4 ) ii
96  integer ( kind = 4 ) init
97  integer ( kind = 4 ) k
98  integer ( kind = 4 ) next
99  integer ( kind = 4 ) perm(n)
100  real ( kind = 8 ) tmp
101  real ( kind = 8 ) tmp1
102  real ( kind = 8 ) x(n)
103 
104  init = 1
105  tmp = x(init)
106  ii = perm(init)
107  perm(init)= -perm(init)
108  k = 0
109  !
110  ! The main loop.
111  !
112  6 continue
113 
114  k = k + 1
115  !
116  ! Save the chased element.
117  !
118  tmp1 = x(ii)
119  x(ii) = tmp
120  next = perm(ii)
121 
122  if ( next < 0 ) then
123  go to 65
124  end if
125  !
126  ! Test for end.
127  !
128  if ( n < k ) then
129  perm(1:n) = -perm(1:n)
130  return
131  end if
132 
133  tmp = tmp1
134  perm(ii) = -perm(ii)
135  ii = next
136  !
137  ! End of the loop.
138  !
139  go to 6
140  !
141  ! Reinitialize cycle.
142  !
143  65 continue
144 
145  init = init + 1
146 
147  if ( n < init ) then
148  perm(1:n) = -perm(1:n)
149  return
150  end if
151 
152  if ( perm(init) < 0 ) then
153  go to 65
154  end if
155 
156  tmp = x(init)
157  ii = perm(init)
158  perm(init) = -perm(init)
159  go to 6
160 
subroutine init()
Definition: GridSorting.f90:24
Here is the caller graph for this function:

◆ rperm()

subroutine rperm ( integer ( kind = 4 )  nrow,
real ( kind = 8 ), dimension(*)  a,
integer ( kind = 4 ), dimension(*)  ja,
integer ( kind = 4 ), dimension(nrow+1)  ia,
real ( kind = 8 ), dimension(*)  ao,
integer ( kind = 4 ), dimension(*)  jao,
integer ( kind = 4 ), dimension(nrow+1)  iao,
integer ( kind = 4 ), dimension(nrow)  perm,
integer ( kind = 4 )  job 
)

Definition at line 249 of file sparsekit.f90.

250 
251  !*****************************************************************************80
252  !
253  !! RPERM permutes the rows of a matrix in CSR format.
254  !
255  ! Discussion:
256  !
257  ! This routine computes B = P*A where P is a permutation matrix.
258  ! the permutation P is defined through the array perm: for each j,
259  ! perm(j) represents the destination row number of row number j.
260  !
261  ! Modified:
262  !
263  ! 07 January 2004
264  !
265  ! Author:
266  !
267  ! Youcef Saad
268  !
269  ! Parameters:
270  !
271  ! Input, integer ( kind = 4 ) NROW, the row dimension of the matrix.
272  !
273  ! Input, real A(*), integer ( kind = 4 ) JA(*), IA(NROW+1), the matrix in CSR
274  ! Compressed Sparse Row format.
275  !
276  ! perm = integer ( kind = 4 ) array of length nrow containing the
277  ! permutation arrays
278  ! for the rows: perm(i) is the destination of row i in the
279  ! permuted matrix.
280  ! ---> a(i,j) in the original matrix becomes a(perm(i),j)
281  ! in the output matrix.
282  !
283  ! job = integer ( kind = 4 ) indicating the work to be done:
284  ! job = 1 permute a, ja, ia into ao, jao, iao
285  ! (including the copying of real values ao and
286  ! the array iao).
287  ! job /= 1 : ignore real values.
288  ! (in which case arrays a and ao are not needed nor
289  ! used).
290  !
291  ! Output, real AO(*), integer ( kind = 4 ) JAO(*), IAO(NROW+1), the permuted
292  ! matrix in CSR Compressed Sparse Row format.
293  !
294  ! note :
295  ! if (job/=1) then the arrays a and ao are not used.
296  !
297  implicit none
298 
299  integer ( kind = 4 ) nrow
300 
301  real ( kind = 8 ) a(*)
302  real ( kind = 8 ) ao(*)
303  integer ( kind = 4 ) i
304  integer ( kind = 4 ) ia(nrow+1)
305  integer ( kind = 4 ) iao(nrow+1)
306  integer ( kind = 4 ) ii
307  integer ( kind = 4 ) j
308  integer ( kind = 4 ) ja(*)
309  integer ( kind = 4 ) jao(*)
310  integer ( kind = 4 ) job
311  integer ( kind = 4 ) k
312  integer ( kind = 4 ) ko
313  integer ( kind = 4 ) perm(nrow)
314  logical values
315 
316  values = ( job == 1 )
317  !
318  ! Determine pointers for output matrix.
319  !
320  do j = 1, nrow
321  i = perm(j)
322  iao(i+1) = ia(j+1) - ia(j)
323  end do
324  !
325  ! Get pointers from lengths.
326  !
327  iao(1) = 1
328  do j = 1, nrow
329  iao(j+1) = iao(j+1) + iao(j)
330  end do
331  !
332  ! Copying.
333  !
334  do ii = 1, nrow
335  !
336  ! Old row = II is new row IPERM(II), and KO is the new pointer.
337  !
338  ko = iao(perm(ii))
339 
340  do k = ia(ii), ia(ii+1)-1
341  jao(ko) = ja(k)
342  if ( values ) then
343  ao(ko) = a(k)
344  end if
345  ko = ko + 1
346  end do
347 
348  end do
349 
350  return
Here is the caller graph for this function: