-
Notifications
You must be signed in to change notification settings - Fork 1
/
test_fortran.f90
105 lines (77 loc) · 2.34 KB
/
test_fortran.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
!=====================
!
!=====================
program main_test
use,intrinsic :: iso_c_binding
implicit none
type(C_PTR) :: plan
integer(C_INT) :: ndim
integer(C_INT) :: nhist
real(C_DOUBLE) :: tolerance
real(C_DOUBLE),allocatable :: xx(:)
real(C_DOUBLE),allocatable :: gradf(:)
real(C_DOUBLE),allocatable :: diag(:)
real(C_DOUBLE) :: ff
integer(C_INT) :: info
integer :: idim
integer :: iter
integer,parameter :: niter=20
#include "libmin.f03"
!=====================
! ndim = 1
ndim = 3
nhist = 5
tolerance = 1.0e-7
allocate(xx(ndim),gradf(ndim))
! plan = libmin_init(ndim,nhist)
! An estimate about the inverse Hessian diagonal
allocate(diag(ndim))
diag(1) = 1.0 / 1.00
diag(2) = 1.0 / 0.25
diag(3) = 1.0 / 6.00
plan = libmin_init_diag(ndim,nhist,diag)
deallocate(diag)
do idim=1,ndim
xx(idim) = 1.11258d0 + 1.5415 * idim
enddo
do iter=1,niter
ff = eval_function(ndim,xx)
gradf = eval_gradient(ndim,xx)
write(*,*) ' ===== iteration:',iter,' Function value: ',ff
write(*,*) 'Input x: ',xx(:)
if( NORM2(gradf) < tolerance ) then
write(*,*) 'Convergence reached'
write(*,*) 'Minarg f:',xx(:)
exit
endif
info = libmin_execute(plan,xx,ff,gradf)
write(*,*) 'Output x: ',xx(:)
if( info < 0 ) then
write(*,*) 'Problem happened',info
stop 'STOP'
endif
enddo
call libmin_destroy(plan)
deallocate(xx,gradf)
contains
function eval_function(ndim,xx) result(ff)
integer(C_INT),intent(in) :: ndim
real(C_DOUBLE),intent(in) :: xx(ndim)
real(C_DOUBLE) :: ff
!====
ff = 0.500_C_DOUBLE * ( xx(1) - 1.000_C_DOUBLE )**2 &
+ 0.125_C_DOUBLE * ( xx(2) - 3.000_C_DOUBLE )**2 &
+ 3.000_C_DOUBLE * ( xx(3) + 2.000_C_DOUBLE )**2
! ff = 0.50 * ( xx(1) - 0.500 )**2 + 0.0263 * ( xx(1) - 0.500)**4
end function eval_function
function eval_gradient(ndim,xx) result(gradf)
integer(C_INT),intent(in) :: ndim
real(C_DOUBLE),intent(in) :: xx(ndim)
real(C_DOUBLE) :: gradf(ndim)
!====
gradf(1) = 0.5000_C_DOUBLE * ( xx(1) - 1.000_C_DOUBLE ) * 2.0_C_DOUBLE
gradf(2) = 0.1250_C_DOUBLE * ( xx(2) - 3.000_C_DOUBLE ) * 2.0_C_DOUBLE
gradf(3) = 3.0000_C_DOUBLE * ( xx(3) + 2.000_C_DOUBLE ) * 2.0_C_DOUBLE
! gradf(1) = ( xx(1) - 0.500 ) + 4.0 * 0.0263 * ( xx(1) - 0.500)**3
end function eval_gradient
end program main_test