SICOPOLIS V5-dev  Revision 1264
sico_maths_m_grad.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : s i c o _ s l e _ s o l v e r s
4 !
5 !> @file
6 !!
7 !! Solvers for systems of linear equations used by SICOPOLIS.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2018 Ralf Greve, Tatsuru Sato
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Solvers for systems of linear equations used by SICOPOLIS.
34 !<------------------------------------------------------------------------------
36 
37 use sico_types_m
38 
39  private:: transpose_csr
40  public :: sor_spr_grad, tri_sle_grad
41 #ifdef BUILD_LIS
42  public :: sico_lis_solver_grad
43 #endif
44 
45  interface sor_sprs_grad
46  module procedure sor_sprs_grad_local
47  end interface
48 
49  interface tri_sle_grad
50  module procedure tri_sle_grad_local
51  end interface
52 
53 #ifdef BUILD_LIS
54  interface sico_lis_solver_grad
55  module procedure sico_lis_solver_grad_local
56  end interface
57 #endif
58 contains
59 
60 subroutine transpose_csr(a_value, a_index, a_diag_index, a_ptr, &
61  nnz, nmax, &
62  b_value, b_index, b_diag_index, b_ptr)
63 
64 implicit none
65 
66 integer(i4b), intent(in) :: nnz, nmax
67 integer(i4b), dimension(nmax+1), intent(in) :: a_ptr
68 integer(i4b), dimension(nnz), intent(in) :: a_index
69 integer(i4b), dimension(nmax), intent(in) :: a_diag_index
70 real(dp), dimension(nnz), intent(in) :: a_value
71 integer(i4b), dimension(nmax+1), intent(out) :: b_ptr
72 integer(i4b), dimension(nnz), intent(out) :: b_index
73 integer(i4b), dimension(nmax), intent(out) :: b_diag_index
74 real(dp), dimension(nnz), intent(out) :: b_value
75 
76 integer(i4b) :: nr, k
77 integer(i4b), dimension(nmax) :: b_ptr_plus
78 
79  b_ptr = 0
80  do nr=1, nmax
81  do k=a_ptr(nr), a_ptr(nr+1)-1
82  b_ptr(a_index(k)+1) = b_ptr(a_index(k)+1) + 1
83  end do
84  end do
85  b_ptr(1) = 1
86  do nr=1, nmax
87  b_ptr(nr+1) = b_ptr(nr) + b_ptr(nr+1)
88  end do
89  b_ptr_plus = 0
90  do nr=1, nmax
91  do k=a_ptr(nr), a_ptr(nr+1)-1
92  b_index(b_ptr(a_index(k)) + b_ptr_plus(a_index(k))) = nr
93  if(a_index(k).eq.nr)then
94  b_diag_index(nr) = b_ptr(a_index(k)) + b_ptr_plus(a_index(k))
95  end if
96  b_value(b_ptr(a_index(k)) + b_ptr_plus(a_index(k))) = a_value(k)
97  b_ptr_plus(a_index(k)) = b_ptr_plus(a_index(k)) + 1
98  end do
99  end do
100 end subroutine transpose_csr
101 
102 subroutine transpose_csr_no_diagonal(a_value, a_index, a_ptr, &
103  nnz, nmax, &
104  b_value, b_index, b_ptr)
106 implicit none
107 
108 integer(i4b), intent(in) :: nnz, nmax
109 integer(i4b), dimension(nmax+1), intent(in) :: a_ptr
110 integer(i4b), dimension(nnz), intent(in) :: a_index
111 real(dp), dimension(nnz), intent(in) :: a_value
112 integer(i4b), dimension(nmax+1), intent(out) :: b_ptr
113 integer(i4b), dimension(nnz), intent(out) :: b_index
114 real(dp), dimension(nnz), intent(out) :: b_value
115 
116 integer(i4b) :: nr, k
117 integer(i4b), dimension(nmax) :: b_ptr_plus
118 
119  b_ptr = 0
120  do nr=1, nmax
121  do k=a_ptr(nr), a_ptr(nr+1)-1
122  b_ptr(a_index(k)+1) = b_ptr(a_index(k)+1) + 1
123  end do
124  end do
125  b_ptr(1) = 1
126  do nr=1, nmax
127  b_ptr(nr+1) = b_ptr(nr) + b_ptr(nr+1)
128  end do
129  b_ptr_plus = 0
130  do nr=1, nmax
131  do k=a_ptr(nr), a_ptr(nr+1)-1
132  b_index(b_ptr(a_index(k)) + b_ptr_plus(a_index(k))) = nr
133  b_value(b_ptr(a_index(k)) + b_ptr_plus(a_index(k))) = a_value(k)
134  b_ptr_plus(a_index(k)) = b_ptr_plus(a_index(k)) + 1
135  end do
136  end do
137 end subroutine transpose_csr_no_diagonal
138 !-------------------------------------------------------------------------------
139 !> SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b
140 !! [matrix storage: compressed sparse row CSR,
141 !! represented by arrays lgs_a_value(values), lgs_a_index (indices)
142 !! and lgs_a_ptr (pointers)].
143 !<------------------------------------------------------------------------------
144 subroutine sor_sprs_grad_local(lgs_a_value, lgs_a_value_b, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, &
145  lgs_b_value, lgs_b_value_b, &
146  nnz, nmax, n_sprs, omega, eps_sor, lgs_x_value, lgs_x_value_b, ierr)
148 use sico_maths_m, only : sor_sprs
149 !use sico_maths_m, only : sor_sprs_local
150 implicit none
151 
152 integer(i4b), intent(in) :: nnz, nmax, n_sprs
153 real(dp), intent(in) :: omega, eps_sor
154 integer(i4b), dimension(nmax+1), intent(in) :: lgs_a_ptr
155 integer(i4b), dimension(nnz), intent(in) :: lgs_a_index
156 integer(i4b), dimension(nmax), intent(in) :: lgs_a_diag_index
157 real(dp), dimension(nnz), intent(in) :: lgs_a_value
158 real(dp), dimension(nmax), intent(in) :: lgs_b_value
159 
160 integer(i4b), intent(out) :: ierr
161 real(dp), dimension(nmax), intent(inout) :: lgs_x_value
162 
163 real(dp), dimension(nnz), intent(inout) :: lgs_a_value_b
164 real(dp), dimension(nmax), intent(inout) :: lgs_b_value_b
165 real(dp), dimension(nmax), intent(inout) :: lgs_x_value_b
166 
167 integer(i4b), dimension(nmax+1) :: lgs_aT_ptr
168 integer(i4b), dimension(nnz) :: lgs_aT_index
169 integer(i4b), dimension(nmax) :: lgs_aT_diag_index
170 real(dp), dimension(nnz) :: lgs_aT_value
171 real(dp), dimension(nmax) :: incrbb
172 integer(i4b) :: iter
173 integer(i4b) :: nr, k
174 real(dp), allocatable, dimension(:) :: lgs_x_value_prev
175 real(dp) :: b_nr
176 logical :: flag_convergence
177 logical :: inif, infor
178 
179  incrbb = 0.0
180  call transpose_csr(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, &
181  nnz, nmax, &
182  lgs_at_value, lgs_at_index, lgs_at_diag_index, lgs_at_ptr)
183 
184  call sor_sprs(lgs_at_value, lgs_at_index, lgs_at_diag_index, lgs_at_ptr, &
185  lgs_x_value_b, &
186  nnz, nmax, n_sprs, omega, eps_sor, incrbb, ierr)
187 
188  DO nr=1,nmax
189  lgs_b_value_b(nr) = lgs_b_value_b(nr) + incrbb(nr)
190  ENDDO
191 
192  call sor_sprs(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, &
193  lgs_b_value, &
194  nnz, nmax, n_sprs, omega, eps_sor, lgs_x_value, ierr)
195  do nr=1, nmax
196  inif = .false.
197  infor = .false.
198  if(nr .eq. lgs_a_diag_index(nr)) then
199  inif = .true.
200  lgs_a_value_b(lgs_a_diag_index(nr)) = lgs_a_value_b(lgs_a_diag_index(nr)) -&
201  lgs_x_value(lgs_a_index(nr)) * incrbb(nr)
202  end if
203  do k=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
204  infor = .true.
205  lgs_a_value_b(k) = lgs_a_value_b(k) - lgs_x_value(lgs_a_index(k)) * incrbb(nr)
206  end do
207  end do
208  lgs_x_value_b = 0.0
209 
210 end subroutine sor_sprs_grad_local
211 
212 !-------------------------------------------------------------------------------
213 !> Solution of a system of linear equations Ax=b with tridiagonal matrix A.
214 !! @param[in] a0 a0(j) is element A_(j,j-1) of Matrix A
215 !! @param[in] a1 a1(j) is element A_(j,j) of Matrix A
216 !! @param[in] a2 a2(j) is element A_(j,j+1) of Matrix A
217 !! @param[in] b inhomogeneity vector
218 !! @param[in] nrows size of matrix A (indices run from 0 (!!!) to nrows)
219 !! @param[out] x Solution vector.
220 !<------------------------------------------------------------------------------
221 subroutine tri_sle_grad_local(a0, a0_b, a1, a1_b, a2, a2_b, x, x_b, b, b_b, nrows)
223 use sico_maths_m, only : tri_sle
224 !use sico_maths_m, only : tri_sle_local
225 implicit none
226 
227 integer(i4b), intent(in) :: nrows
228 real(dp), dimension(0:*), intent(in) :: a0, a2
229 real(dp), dimension(0:*), intent(inout) :: a1, b
230 
231 real(dp), dimension(0:*), intent(out) :: x
232 
233 
234 real(dp), dimension(0:*), intent(inout) :: a0_b, a2_b
235 real(dp), dimension(0:*), intent(inout) :: a1_b, b_b
236 real(dp), dimension(0:*), intent(inout) :: x_b
237 real(dp), dimension(0:nrows) :: a0T, a1T, a2T
238 real(dp), dimension(0:nrows) :: incrbb
239 integer(i4b) :: i
240  a0t(0) = 0.0
241  a0t(1:nrows) = a2(0:nrows-1)
242  a1t = a1(0:nrows)
243  a2t(0:nrows-1) = a0(1:nrows)
244  a2t(nrows) = 0.0
245  call tri_sle(a0t, a1t, a2t, incrbb, x_b, nrows)
246  DO i=0,nrows
247  b_b(i) = b_b(i) + incrbb(i)
248  ENDDO
249  call tri_sle(a0, a1, a2, x, b, nrows)
250  DO i=0,nrows
251  a0_b(i) = a0_b(i) - x(i-1)*incrbb(i)
252  ENDDO
253  DO i=0,nrows
254  a1_b(i) = a1_b(i) - x(i)*incrbb(i)
255  ENDDO
256  DO i=0,nrows
257  a2_b(i) = a2_b(i) - x(i+1)*incrbb(i)
258  ENDDO
259  DO i=0,nrows
260  x_b(i) = 0
261  ENDDO
262 end subroutine tri_sle_grad_local
263 
264 !-------------------------------------------------------------------------------
265 #ifdef BUILD_LIS
266 subroutine sico_lis_solver_grad_local(nmax, nnz, &
267  lgs_a_ptr, lgs_a_index, &
268  lgs_a_value, lgs_a_value_b, lgs_b_value,&
269  lgs_b_value_b, lgs_x_value, lgs_x_value_b)
270 
271 use sico_maths_m, only :sico_lis_solver
272 implicit none
273 
274 integer(i4b) :: nc, nr
275 integer(i4b), intent(in) :: nmax
276 integer(i4b), intent(in) :: nnz
277 integer(i4b), dimension(nmax+1), intent(in) :: lgs_a_ptr
278 integer(i4b), dimension(nnz), intent(in) :: lgs_a_index
279 
280 real(dp), dimension(nnz), intent(in) :: lgs_a_value
281 real(dp), dimension(nmax), intent(in) :: lgs_b_value
282 real(dp), dimension(nmax), intent(inout) :: lgs_x_value
283 
284 real(dp), dimension(nnz), intent(inout) :: lgs_a_value_b
285 real(dp), dimension(nmax), intent(inout) :: lgs_b_value_b
286 real(dp), dimension(nmax), intent(inout) :: lgs_x_value_b
287 
288 integer(i4b), dimension(nmax+1) :: lgs_aT_ptr
289 integer(i4b), dimension(nnz) :: lgs_aT_index
290 real(dp), dimension(nnz) :: lgs_aT_value
291 real(dp), dimension(nmax) :: incrbb
292 
293  call transpose_csr_no_diagonal(lgs_a_value, lgs_a_index, lgs_a_ptr, &
294  nnz, nmax, &
295  lgs_at_value, lgs_at_index, lgs_at_ptr)
296  call sico_lis_solver(nmax, nnz, &
297  lgs_at_ptr, lgs_at_index, &
298  lgs_at_value, lgs_x_value_b, incrbb)
299  DO nr=1,nmax
300  lgs_b_value_b(nr) = lgs_b_value_b(nr) + incrbb(nr)
301  ENDDO
302 
303  call sico_lis_solver(nmax, nnz, &
304  lgs_a_ptr, lgs_a_index, &
305  lgs_a_value, lgs_b_value, lgs_x_value)
306 
307  do nr=1, nmax
308  do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
309  lgs_a_value_b(nc) = lgs_a_value_b(nc) - lgs_x_value(lgs_a_index(nc)) * incrbb(nr)
310  end do
311  end do
312  lgs_x_value_b = 0.0
313 end subroutine sico_lis_solver_grad_local
314 #endif
315 end module sico_maths_m_grad
316 !
Solvers for systems of linear equations used by SICOPOLIS.
subroutine transpose_csr_no_diagonal(a_value, a_index, a_ptr, nnz, nmax, b_value, b_index, b_ptr)
subroutine sor_sprs_grad_local(lgs_a_value, lgs_a_value_b, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, lgs_b_value, lgs_b_value_b, nnz, nmax, n_sprs, omega, eps_sor, lgs_x_value, lgs_x_value_b, ierr)
SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b [matrix storage: compressed sparse row ...
Several mathematical tools used by SICOPOLIS.
subroutine tri_sle(a0, a1, a2, x, b, nrows)
Solution of a system of linear equations Ax=b with tridiagonal matrix A.
subroutine sor_sprs(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr,
SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b [matrix storage: compressed sparse row ...
Declarations of kind types for SICOPOLIS.
subroutine tri_sle_grad_local(a0, a0_b, a1, a1_b, a2, a2_b, x, x_b, b, b_b, nrows)
Solution of a system of linear equations Ax=b with tridiagonal matrix A.