SICOPOLIS V5-dev  Revision 1264
sico_maths_m_stub.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 #if (CALCTHK!=3 && CALCTHK!=4 && CALCTHK!=6 )
40  public :: sor_sprs, tri_sle, bilinint, my_erfc
41 #else
42  public :: tri_sle, bilinint, my_erfc, sico_lis_solver
43 #endif
44 #if (CALCTHK!=3 && CALCTHK!=4 && CALCTHK!=6 )
45  interface sor_sprs
46  module procedure sor_sprs_stub
47  end interface
48 #endif
49  interface tri_sle
50  module procedure tri_sle_stub
51  end interface
52 
53  interface bilinint
54  module procedure bilinint_stub
55  end interface
56 
57  interface my_erfc
58  module procedure my_erfc_stub
59  end interface
60 
61 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
62 !#if (CALCTHK==3 || CALCTHK==4 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
63  interface sico_lis_solver
64  module procedure sico_lis_solver_stub
65  end interface
66 #endif
67 
68 contains
69 
70 !-------------------------------------------------------------------------------
71 !> SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b
72 !! [matrix storage: compressed sparse row CSR,
73 !! represented by arrays lgs_a_value(values), lgs_a_index (indices)
74 !! and lgs_a_ptr (pointers)].
75 !<------------------------------------------------------------------------------
76 #if (CALCTHK!=3 && CALCTHK!=4 && CALCTHK!=6 )
77 subroutine sor_sprs_stub(lgs_a_value, lgs_a_index, lgs_a_diag_index,&
78  lgs_a_ptr, &
79  lgs_b_value, &
80  nnz, nmax,&
81 #ifdef ALLOW_OPENAD
82  n_sprs,&
83 #endif
84  omega, eps_sor, lgs_x_value, ierr)
85  !$openad xxx template oad_template_sor_sprs.f90
86 implicit none
87 
88 integer(i4b), intent(in) :: nnz, nmax, n_sprs
89 real(dp), intent(in) :: omega, eps_sor
90 integer(i4b), dimension(nmax+1), intent(inout) :: lgs_a_ptr
91 integer(i4b), dimension(nnz), intent(inout) :: lgs_a_index
92 integer(i4b), dimension(nmax), intent(inout) :: lgs_a_diag_index
93 real(dp), dimension(nnz), intent(inout) :: lgs_a_value
94 real(dp), dimension(nmax), intent(inout) :: lgs_b_value
95 
96 integer(i4b), intent(out) :: ierr
97 real(dp), dimension(nmax), intent(inout) :: lgs_x_value
98 
99 integer(i4b) :: iter
100 integer(i4b) :: nr, k
101 real(dp), allocatable, dimension(:) :: lgs_x_value_prev
102 real(dp) :: b_nr
103 logical :: flag_convergence
104 
105 lgs_x_value_prev(0) = lgs_x_value(0) + omega + eps_sor + lgs_a_value(0) + lgs_b_value(0)
106 b_nr = lgs_x_value_prev(0)
107 lgs_x_value(0) = b_nr
108 
109 end subroutine sor_sprs_stub
110 #endif
111 
112 !-------------------------------------------------------------------------------
113 !> Solution of a system of linear equations Ax=b with tridiagonal matrix A.
114 !! @param[in] a0 a0(j) is element A_(j,j-1) of Matrix A
115 !! @param[in] a1 a1(j) is element A_(j,j) of Matrix A
116 !! @param[in] a2 a2(j) is element A_(j,j+1) of Matrix A
117 !! @param[in] b inhomogeneity vector
118 !! @param[in] nrows size of matrix A (indices run from 0 (!!!) to nrows)
119 !! @param[out] x Solution vector.
120 !<------------------------------------------------------------------------------
121 subroutine tri_sle_stub(a0, a1, a2, x, b, nrows)
123 implicit none
124 
125 integer(i4b), intent(in) :: nrows
126 
127 real(dp), dimension(0:KCMAX+KTMAX+KRMAX+IMAX+JMAX), intent(in) :: a0, a2
128 real(dp), dimension(0:KCMAX+KTMAX+KRMAX+IMAX+JMAX), intent(inout) :: a1, b
129 
130 real(dp), dimension(0:KCMAX+KTMAX+KRMAX+IMAX+JMAX), intent(out) :: x
131 real(dp), dimension(0:nrows) :: help_x
132 integer(i4b) :: n
133 
134 !-------- Generate an upper triangular matrix
135 ! ('obere Dreiecksmatrix') --------
136 !x(1) = a0(0) + a1(0) + a2(0) + b(0)
137 !x(0) = x(1)
138 !b(0) = x(1)
139 !a1(0) = x(1)
140 
141 do n=1, nrows
142  a1(n) = a1(n) - a0(n)/a1(n-1)*a2(n-1)
143 end do
144 
145 do n=1, nrows
146  b(n) = b(n) - a0(n)/a1(n-1)*b(n-1)
147 ! a0(n) = 0.0_dp , not needed in the following, therefore
148 ! not set
149 end do
150 
151 !-------- Iterative solution of the new system --------
152 
153 ! x(nrows) = b(nrows)/a1(nrows)
154 
155 ! do n=nrows-1, 0, -1
156 ! x(n) = (b(n)-a2(n)*x(n+1))/a1(n)
157 ! end do
158 
159 help_x(0) = b(nrows)/a1(nrows)
160 
161 do n=1, nrows
162  help_x(n) = b(nrows-n)/a1(nrows-n) &
163  -a2(nrows-n)/a1(nrows-n)*help_x(n-1)
164 end do
165 
166 do n=0, nrows
167  x(n) = help_x(nrows-n)
168 end do
169 
170 end subroutine tri_sle_stub
171 
172 !-------------------------------------------------------------------------------
173 !> Bilinear interpolation.
174 !<------------------------------------------------------------------------------
175  subroutine bilinint_stub(x1, x2, y1, y2, z11, z12, z21, z22, x, y, retval)
177  implicit none
178 
179  real(dp), intent(in) :: x1, x2, y1, y2, z11, z12, z21, z22, x, y
180 
181  real(dp) :: t, u
182  real(dp), intent(out) :: retval
183 
184  real(dp), parameter :: I = 1.0_dp
185 
186  t = (x-x1)/(x2-x1)
187  u = (y-y1)/(y2-y1)
188 
189  retval = (i-t)*(i-u)*z11 + (i-t)*u*z12 + t*(i-u)*z21 + t*u*z22
190 
191  end subroutine bilinint_stub
192 
193 !-------------------------------------------------------------------------------
194 !> Computation of the complementary error function erfc(x) = 1-erf(x)
195 !! with a fractional error everywhere less than 1.2 x 10^(-7)
196 !! (formula by Press et al., 'Numerical Recipes in Fortran 77').
197 !<------------------------------------------------------------------------------
198  subroutine my_erfc_stub(x, retval)
200  implicit none
201 
202  real(dp), intent(in) :: x
203  real(dp), intent(out) :: retval
204 
205  real(dp) :: t, z
206 
207  z = abs(x)
208  t = 1.0_dp/(1.0_dp+0.5_dp*z)
209 
210  retval = t * exp( -z*z -1.26551223_dp &
211  + t * ( 1.00002368_dp &
212  + t * ( 0.37409196_dp &
213  + t * ( 0.09678418_dp &
214  + t * ( -0.18628806_dp &
215  + t * ( 0.27886807_dp &
216  + t * ( -1.13520398_dp &
217  + t * ( 1.48851587_dp &
218  + t * ( -0.82215223_dp &
219  + t * 0.17087277_dp ) ) ) ) ) ) ) ) )
220 
221  if (x < 0.0_dp) retval = 2.0_dp-retval
222 
223  end subroutine my_erfc_stub
224 
225 !-------------------------------------------------------------------------------
226 #ifdef ALLOW_OPENAD /* OAD VERSION */
227  subroutine my_erfc(x, retval)
228 
229  implicit none
230 
231  real(dp), intent(in) :: x
232  real(dp), intent(out) :: retval
233 
234  real(dp) :: t, z
235 
236  z = abs(x)
237  t = 1.0_dp/(1.0_dp+0.5_dp*z)
238 
239  retval = t * exp( -z*z -1.26551223_dp &
240  + t * ( 1.00002368_dp &
241  + t * ( 0.37409196_dp &
242  + t * ( 0.09678418_dp &
243  + t * ( -0.18628806_dp &
244  + t * ( 0.27886807_dp &
245  + t * ( -1.13520398_dp &
246  + t * ( 1.48851587_dp &
247  + t * ( -0.82215223_dp &
248  + t * 0.17087277_dp ) ) ) ) ) ) ) ) )
249 
250  if (x < 0.0_dp) retval = 2.0_dp-retval
251 
252  end subroutine my_erfc
253 #endif /* OAD VERSION */
254 
255 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
256 !#if (CALCTHK==3 || CALCTHK==4 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
257 subroutine sico_lis_solver_stub(nmax, nnz, &
258  lgs_a_ptr, lgs_a_index, &
259  lgs_a_value, lgs_b_value, lgs_x_value)
260 
261  !$openad xxx template oad_template_sico_lis_solver.f90
262 implicit none
263 
264 integer(i4b), intent(in) :: nmax
265 integer(i4b), intent(in) :: nnz
266 integer(i4b), dimension(nmax+1),intent(inout) :: lgs_a_ptr
267 integer(i4b), dimension(nnz), intent(inout) :: lgs_a_index
268 real(dp), dimension(nnz), intent(inout) :: lgs_a_value
269 real(dp), dimension(nmax), intent(inout) :: lgs_b_value
270 real(dp), dimension(nmax), intent(inout) :: lgs_x_value
271 
272 ! ------ Settings for Lis
273 
274 lgs_x_value(1) = lgs_b_value(1) / lgs_a_value(1)
275 
276 end subroutine sico_lis_solver_stub
277 #endif
278 
279 end module sico_maths_m_stub
280 !
Declarations of kind types for SICOPOLIS.
Solvers for systems of linear equations used by SICOPOLIS.
subroutine my_erfc_stub(x, retval)
Computation of the complementary error function erfc(x) = 1-erf(x) with a fractional error everywhere...
subroutine sor_sprs_stub(lgs_a_value, lgs_a_index, lgs_a_diag_index, lgs_a_ptr, lgs_b_value, nnz, nmax,
SOR solver for a system of linear equations lgs_a*lgs_x=lgs_b [matrix storage: compressed sparse row ...
subroutine tri_sle_stub(a0, a1, a2, x, b, nrows)
Solution of a system of linear equations Ax=b with tridiagonal matrix A.
subroutine bilinint_stub(x1, x2, y1, y2, z11, z12, z21, z22, x, y, retval)
Bilinear interpolation.