-
Notifications
You must be signed in to change notification settings - Fork 3
/
tao_minimizer.f90
264 lines (208 loc) · 7.91 KB
/
tao_minimizer.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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
!---------------------------------------------------------------------------
! !
! Copyright 2018 Anna Teruzzi, OGS, Trieste !
! !
! This file is part of 3DVarBio.
! 3DVarBio is based on OceanVar (Dobricic, 2006) !
! !
! 3DVarBio is free software: you can redistribute it and/or modify. !
! it under the terms of the GNU General Public License as published by !
! the Free Software Foundation, either version 3 of the License, or !
! (at your option) any later version. !
! !
! 3DVarBio is distributed in the hope that it will be useful, !
! but WITHOUT ANY WARRANTY; without even the implied warranty of !
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !
! GNU General Public License for more details. !
! !
! You should have received a copy of the GNU General Public License !
! along with OceanVar. If not, see <http://www.gnu.org/licenses/>. !
! !
!---------------------------------------------------------------------------
subroutine tao_minimizer
#include "petsc/finclude/petscvecdef.h"
use drv_str
use ctl_str
use mpi_str
use petscvec
implicit none
#include "petsc/finclude/petsctao.h"
PetscErrorCode :: ierr
Tao :: tao
Vec :: MyState ! array that stores the (temporary) state
PetscInt :: n, M, GlobalStart, MyEnd, iter
PetscReal :: fval, gnorm, cnorm, xdiff
PetscScalar :: MyTolerance
TaoConvergedReason :: reason
integer(i4) :: j
real(8) :: MaxGrad
! Working arrays
PetscInt, allocatable, dimension(:) :: loc
PetscScalar, allocatable, dimension(:) :: MyValues
PetscScalar, pointer :: xtmp(:)
external MyFuncAndGradient
if(MyId .eq. 0) then
print*,'PETSc-TAO lmvm minimizer configuration'
print*, ''
endif
if(MyId .eq. 0) then
write(drv%dia,*) ''
write(drv%dia,*) "Within tao_minimizer subroutine!"
endif
! Allocate working arrays
! TAO needs to know both n (the local size of the array)
! and M (global size of the array)
n = ctl%n
M = ctl%n_global
ALLOCATE(loc(n), MyValues(n))
! Create MyState array and fill it
! VecGetOwnergshipRange returns the indices related to the local
! section of the PETSc Vector (required in order to fill the initial state vector,
! further informations can be found in the PETSc user manual, Section 2.1 and 2.2)
call VecCreateMPI(Var3DCommunicator, n, M, MyState, ierr)
call VecGetOwnershipRange(MyState, GlobalStart, MyEnd, ierr)
if(drv%Verbose .eq. 1) &
print*, "MyState initialization by MyId ", MyId, "with indices: ", GlobalStart, MyEnd
if( ctl%n .ne. MyEnd - GlobalStart ) then
print*, ""
print*, "WARNING!!"
print*, "ctl%n .ne. GlobalStart - MyEnd"
print*, "ctl%n = ", ctl%n
print*, "GlobalStart = ", GlobalStart
print*, "MyEnd = ", MyEnd
print*, ""
endif
! Setting initial state vector to zero
do j = 1, ctl%n
loc(j) = GlobalStart + j - 1
MyValues(j) = 0.
end do
! Setting only local values (since each process can access at all entries of MyState)
call VecSetValues(MyState, ctl%n, loc, MyValues, INSERT_VALUES, ierr)
call VecAssemblyBegin(MyState, ierr)
call VecAssemblyEnd(MyState, ierr)
! Iteration counter initialization
drv%MyCounter = 0
! Create Tao object and set type LMVM (ones that use BFGS minimization algorithm)
call TaoCreate(Var3DCommunicator, tao, ierr)
CHKERRQ(ierr)
call TaoSetType(tao,"lmvm",ierr)
CHKERRQ(ierr)
! Set initial solution array, MyBounds and MyFuncAndGradient routines
call TaoSetInitialVector(tao, MyState, ierr)
CHKERRQ(ierr)
call TaoSetObjectiveAndGradientRoutine(tao, MyFuncAndGradient, PETSC_NULL_OBJECT, ierr)
CHKERRQ(ierr)
! Calling costf in order to compute
! the initial gradient. This will be used to
! set MyTolerance
call costf
MaxGrad = 0
do j=1,ctl%n
MaxGrad = max(MaxGrad, abs(ctl%g_c(j)))
end do
! since MaxGrad is the local maximum value, the MPI_Allreduce call
! is required to find the global maximum value
call MPI_Allreduce(MPI_IN_PLACE, MaxGrad, 1, MPI_REAL8, MPI_MAX, Var3DCommunicator, ierr)
MyTolerance = ctl%pgper * MaxGrad
if(MyId .eq. 0) then
print*, "Setting MyTolerance", MyTolerance
print*, ""
write(drv%dia,*) "Setting MyTolerance", MyTolerance
endif
! setting tolerances
call TaoSetTolerances(tao, MyTolerance, 1.d-4, ctl%pgper, ierr)
CHKERRQ(ierr)
! calling the solver to minimize the problem
call TaoSolve(tao, ierr)
CHKERRQ(ierr)
! printing solver info
if(MyId .eq. 0) then
print*, ''
print*, 'Tao Solver Info:'
print*, ''
endif
call TaoView(tao, PETSC_VIEWER_STDOUT_WORLD, ierr)
call TaoGetSolutionStatus(tao, iter, fval, gnorm, cnorm, xdiff, reason, ierr)
if(reason .lt. 0) then
if(MyId .eq. 0) then
print*, "TAO failed to find a solution"
print*, "Aborting.."
endif
call MPI_Barrier(Var3DCommunicator, ierr)
call MPI_Abort(Var3DCommunicator, -1, ierr)
endif
! Get the solution and copy into ctl%x_c array
call TaoGetSolutionVector(tao, MyState, ierr)
CHKERRQ(ierr)
call VecGetArrayReadF90(MyState, xtmp, ierr)
CHKERRQ(ierr)
do j = 1, ctl%n
ctl%x_c(j) = xtmp(j)
end do
call VecRestoreArrayReadF90(MyState, xtmp, ierr)
CHKERRQ(ierr)
! Deallocating variables
DEALLOCATE(loc, MyValues)
! free memory
call TaoDestroy(tao, ierr)
CHKERRQ(ierr)
call VecDestroy(MyState, ierr)
CHKERRQ(ierr)
if(MyId .eq. 0) then
write(drv%dia,*) 'Minimization done with ', drv%MyCounter
write(drv%dia,*) 'iterations'
write(drv%dia,*) ''
print*, ""
print*, "Minimization done with ", drv%MyCounter, "iterations"
print*, ""
endif
end subroutine tao_minimizer
!-------------------------------------------------!
! subroutine that, given a state !
! MyState (provided by Tao solver), !
! computes the value of cost function !
! in MyState and its gradient !
!-------------------------------------------------!
subroutine MyFuncAndGradient(tao, MyState, CostFunc, Grad, dummy, ierr)
#include "petsc/finclude/petscvecdef.h"
use set_knd
use drv_str
use ctl_str
use petscvec
use mpi_str
implicit none
#include "petsc/finclude/petsctao.h"
Tao :: tao
Vec :: MyState, Grad
PetscScalar :: CostFunc
PetscErrorCode :: ierr
integer(i4) :: dummy, j
! Working arrays
PetscScalar, pointer, dimension(:) :: my_grad
PetscScalar, pointer, dimension(:) :: xtmp
! read temporary state provided by Tao Solver
! and set it in ctl%x_c array in order to compute
! the actual value of Cost Function and the gradient
! with costf subroutine
call VecGetArrayReadF90(MyState, xtmp, ierr)
CHKERRQ(ierr)
do j=1,ctl%n
ctl%x_c(j) = xtmp(j)
end do
call VecRestoreArrayReadF90(MyState, xtmp, ierr)
CHKERRQ(ierr)
! compute function and gradient
call costf
! assign the Cost Function value computed by costf to CostFunc
CostFunc = ctl%f_c
call VecGetArrayF90(Grad, my_grad, ierr)
do j = 1, ctl%n
my_grad(j) = ctl%g_c(j)
end do
call VecRestoreArrayF90(Grad, my_grad, ierr)
! Update counter
drv%MyCounter = drv%MyCounter + 1
! Exit without errors
ierr = 0
end subroutine MyFuncAndGradient