Choral
newton_mod.f90
Go to the documentation of this file.
1 
5 
6 module newton_mod
7 
9  use basic_tools
10  use real_type
11  use abstract_interfaces, only: rntorn
12  use r1d_mod, only: norm_2
13  use krylov_mod
14 
15  implicit none
16  private
17 
18 
19  public :: solve, newton
20 
21  ! %----------------------------------------%
22  ! | |
23  ! | DERIVED TYPE |
24  ! | |
25  ! %----------------------------------------%
26  type newton
27 
28  real(RP) :: tol = 1e-5_rp, eps = 1e-3_rp
29  integer :: verb=0, itmax=100
30  type(krylov):: kry
31 
32  !! Data returned when solving
33  !!
34  real(RP) :: res
35  integer :: iter, feval=0
36  logical :: ierr = .false.
37 
38  end type newton
39 
40  interface newton
41  module procedure newton_create
42  end interface newton
43 
44  interface solve
45  module procedure newton_solve
46  end interface solve
47 
48 contains
49 
50 
53  function newton_create(tol, itmax, eps, verb) result(n)
54  type(newton) :: n
55  real(RP), intent(in), optional :: tol, eps
56  integer , intent(in), optional :: itmax, verb
57 
58  if (present(verb )) n%verb = verb
59  if (present(tol )) n%tol = tol
60  if (present(itmax)) n%itMax= itmax
61  if (present(eps )) n%eps = eps
62 
63  n%kry = krylov(type=kry_gmres, tol=1e-4_rp, &
64  & itmax=100, restart=3)
65 
66  end function newton_create
67 
68 
91  subroutine newton_solve(x, n, F)
92  real(RP), dimension(:), intent(inout) :: x
93  type(newton) , intent(inout) :: n
94  procedure(rntorn) :: F
95 
96  real(RP), dimension(size(x, 1)) :: Fx, z, x0
97  real(DP) :: cpu
98  integer :: ii
99 
100  n%ierr = .false.
101 
102  if (n%verb>0) then
103  write(*,*)'newton_mod: newton_solve: '
104  cpu = clock()
105  end if
106 
107  ! Initiation
108  !
109  n%iter = 0
110  x0 = x
111  call f(fx, x)
112  n%Feval = 1
113  n%res = residual(fx, x)
114  if (n%verb>1) write(*,*)' Iter', n%iter,', residual', n%res
115  if (n%res<n%tol) goto 10
116 
117  ! Newton loop
118  !
119  do ii=1, n%itmax
120 
121  z = x
122 
123  call solve(z, n%kry, fx, df)
124 
125  ! Updates
126  !
127  x = x - z
128 
129  call f(fx, x)
130  n%Feval = n%Feval + 1 + n%kry%AEval
131 
132  n%res = residual(fx, x) ! exit loop test
133  if (n%verb>1) &
134  & write(*,*)' Iter', ii,', residual', n%res
135 
136  ! Test loop end
137  !
138  if (n%res<n%tol) goto 10
139 
140  x0 = x
141 
142  end do
143 
144  ! Convergence failure flag
145  call warning("newton_mod: newton_solve: not converged", n%verb-1)
146  n%ierr = .true.
147 
148 10 continue
149  n%iter = ii
150  if (n%verb>0) then
151  write(*,*)" Iter ",n%iter
152  write(*,*)" Performed x-->F(x) ",n%Feval
153  write(*,*)" Residual ",n%res
154 
155  cpu = clock() - cpu
156  write(*,*)' CPU ', real(cpu, sp)
157  end if
158 
159  contains
160 
161 
162  ! z = Df(x).y \simeq (F(x + eps*y) - F(x))/eps
163  !
164  subroutine df(zz, yy)
165 
166  real(RP), dimension(:), intent(out) :: zz
167  real(RP), dimension(:), intent(in) :: yy
168 
169  call f(zz,x + n%eps*yy)
170  zz = (zz - fx)/n%eps
171 
172  end subroutine df
173 
174  end subroutine newton_solve
175 
176 
177 
180  function residual(Fx, x) result(res)
182  real(RP) :: res
183  real(RP), dimension(:), intent(in) :: Fx, x
184 
185  res = norm_2(x)
186  if (res < real_tol) res = 1._rp
187  res = norm_2(fx)/res
188 
189  end function residual
190 
191 end module newton_mod
192 
DEFINITION OF ABSTRACT_INTERFACES FOR THE LIBRARY CHORAL
BASIC TOOLS
Definition: basic_tools.F90:8
double precision function, public clock()
Returns the internal clock time
Definition: basic_tools.F90:94
integer, parameter kry_gmres
GmRes linear solver.
The type krylov defines the settings of a linear solver.
Definition: krylov_mod.f90:63
L2 vector norm.
Definition: R1d_mod.F90:82
integer, parameter, public sp
simple precision for real numbers
Definition: real_type.F90:60
NEWTON NON LINEAR SOLVER.
Definition: newton_mod.f90:6
REAL NUMBERS PRECISION IN CHORAL: selects simple/double/quad
Definition: real_type.F90:33
OPENMP OPERATIONS ON 1-DIMENSIONAL REAL ARRAYS
Definition: R1d_mod.F90:10
real(rp) function e(x, v1, v2)
real(rp) function residual(Fx, x)
Residual.
Definition: newton_mod.f90:181
type(newton) function newton_create(tol, itmax, eps, verb)
constructor
Definition: newton_mod.f90:54
CHORAL CONSTANTS
subroutine df(zz, yy)
Definition: newton_mod.f90:165
real(rp), parameter, public real_tol
Definition: real_type.F90:91
subroutine newton_solve(x, n, F)
Newton algorithm to solve F(x) = 0.
Definition: newton_mod.f90:92
DERIVED TYPE krylov: for the resolution of linear systems
Definition: krylov_mod.f90:25
subroutine, public warning(message, verb)
Warning message.