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