Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2019 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_bfgs_b2
8 :
9 : CONTAINS
10 :
11 : ! The subroutine bfgs_b2 is an implementation of a BFGS algorithm with
12 : ! box constraints. It is motivated by the paper:
13 : ! Byrd et al., "A Limited Memory Algorithm for Bound Constrained Optimization",
14 : ! SIAM Journal on Scientific Computing 16, 1190 (1995)
15 : ! However, it is substantially different from that paper.
16 : !
17 : ! GM'2019
18 :
19 : ! Note: lastGradient and lastParameters have to be 0.0 if this is the first call to the subroutine and these
20 : ! vectors don't exist yet. Also historyLength has to be 0 in that case.
21 :
22 0 : SUBROUTINE bfgs_b2(vecLength,gradient,lastGradient,minConstraints,maxConstraints,enabledConstraints,parameters,&
23 0 : lastParameters,equalityLinCombi,equalityCriterion,maxHistoryLength,paramCorrections,&
24 0 : gradientCorrections,iStep,mixParam,l_converged,convCrit)
25 :
26 : USE m_constants
27 :
28 : IMPLICIT NONE
29 :
30 : INTEGER, INTENT(IN) :: vecLength
31 : INTEGER, INTENT(IN) :: maxHistoryLength
32 : INTEGER, INTENT(INOUT) :: iStep ! The index of the bfgs_b iteration
33 : REAL, INTENT(IN) :: mixParam
34 : REAL, INTENT(IN) :: equalityCriterion
35 : REAL, INTENT(INOUT) :: gradient(vecLength)
36 : REAL, INTENT(INOUT) :: lastGradient(vecLength)
37 : REAL, INTENT(IN) :: minConstraints(vecLength)
38 : REAL, INTENT(IN) :: maxConstraints(vecLength)
39 : LOGICAL, INTENT(IN) :: enabledConstraints(vecLength)
40 : REAL, INTENT(INOUT) :: parameters(vecLength)
41 : REAL, INTENT(INOUT) :: lastParameters(vecLength)
42 : REAL, INTENT(IN) :: equalityLinCombi(vecLength)
43 : REAL, INTENT(INOUT) :: paramCorrections(vecLength,maxHistoryLength)
44 : REAL, INTENT(INOUT) :: gradientCorrections(vecLength,maxHistoryLength)
45 : LOGICAL, INTENT(OUT) :: l_converged
46 : REAL, INTENT(IN) :: convCrit
47 :
48 : INTEGER :: i, j, stepIndex, numFreeParams, numFixedParams
49 : INTEGER :: historyLength
50 : REAL, PARAMETER :: eps = 1.0e-10
51 : REAL, PARAMETER :: largeEps = 1.0e-6
52 : REAL :: temp, norm, subspaceEqualityCrit, maxDisp, dampingFactor
53 : REAL :: ddot
54 :
55 : LOGICAL :: converged
56 :
57 0 : REAL :: lastGradCorrection(vecLength)
58 0 : REAL :: lastParamCorrection(vecLength)
59 0 : REAL :: prelimNextPoint(vecLength)
60 0 : REAL :: prelimGradient(vecLength)
61 0 : REAL :: tempVecA(vecLength)
62 0 : REAL :: subspaceGrad(vecLength)
63 0 : REAL :: subspaceNextPoint(vecLength)
64 0 : REAL :: subspaceLinCombi(vecLength)
65 0 : REAL :: diffPoint(vecLength)
66 0 : REAL :: diffGrad(vecLength)
67 :
68 0 : LOGICAL :: freeParams(vecLength)
69 :
70 0 : REAL, ALLOCATABLE :: hessMat(:,:), invHessMat(:,:)
71 0 : REAL, ALLOCATABLE :: tempMatB(:,:), tempMatC(:,:)
72 0 : REAL, ALLOCATABLE :: subspaceBasis(:,:)
73 0 : REAL, ALLOCATABLE :: subspaceInvHessMat(:,:)
74 :
75 0 : WRITE(*,*) 'bfgs_b2 (input): parameters(:vecLength): ', parameters(:vecLength)
76 0 : WRITE(*,*) 'bfgs_b2 (input): gradient(:vecLength): ', gradient(:vecLength)
77 :
78 0 : WRITE(5000,*) '========================================================================='
79 0 : WRITE(5000,*) 'ENTERED BFGS_B2: '
80 :
81 : !!! Test start
82 : ! DO i = 1, vecLength - 1
83 : ! IF(ABS(gradient(i)).LT.eps) CYCLE
84 : ! gradient(i) = sqrt(ABS(gradient(i)))*2.0*ATAN(gradient(i))/pi_const
85 : ! END DO
86 : !!! Test end
87 :
88 0 : l_converged = .FALSE.
89 :
90 : ! 0. Add last data pair to history iff it is not obviously destructive for the positive definiteness of the Hessian matrix
91 :
92 0 : IF(ANY(lastGradient(:).NE.0.0).OR.ANY(lastParameters(:).NE.0.0)) THEN
93 :
94 0 : WRITE(2501,'(a,12f10.6)') ' gradient: ', gradient(:)
95 0 : WRITE(2501,'(a,12f10.6)') 'lastGradient: ', lastGradient(:)
96 0 : WRITE(2502,'(a,12f10.6)') ' parameters: ', parameters(:)
97 0 : WRITE(2502,'(a,12f10.6)') 'lastParameters: ', lastParameters(:)
98 :
99 0 : lastGradCorrection(:) = gradient(:) - lastGradient(:)
100 0 : lastParamCorrection(:) = parameters(:) - lastParameters(:)
101 0 : temp = ddot(vecLength,lastParamCorrection(:),1,lastGradCorrection(:),1)
102 0 : IF (temp.GT.eps) THEN
103 : ! Add data pair to history
104 0 : istep = istep + 1
105 0 : stepIndex = MOD(iStep - 1,maxHistoryLength) + 1
106 0 : paramCorrections(:,stepIndex) = lastParamCorrection(:)
107 0 : gradientCorrections(:,stepIndex) = lastGradCorrection(:)
108 : END IF
109 :
110 : ! Check for convergence and return if converged
111 0 : WRITE(*,*) 'MAXVAL(ABS(lastGradCorrection(:))): ', MAXVAL(ABS(lastGradCorrection(:)))
112 0 : WRITE(*,*) 'MAXVAL(ABS(lastParamCorrection(:))): ', MAXVAL(ABS(lastParamCorrection(:)))
113 0 : WRITE(*,*) 'MAXVAL(ABS(lastGradCorrection(:vecLength-1))): ', MAXVAL(ABS(lastGradCorrection(:vecLength-1)))
114 0 : WRITE(*,*) 'MAXVAL(ABS(lastParamCorrection(:vecLength-1))): ', MAXVAL(ABS(lastParamCorrection(:vecLength-1)))
115 : ! IF(ALL(ABS(lastGradCorrection(:)).LT.convCrit).AND.ALL(ABS(lastParamCorrection(:)).LT.convCrit)) THEN
116 0 : IF(ALL(ABS(lastParamCorrection(:vecLength-1)).LT.convCrit)) THEN
117 0 : l_converged = .TRUE.
118 : RETURN
119 : END IF
120 : END IF
121 :
122 : ! 1. Construct Hessian matrix and its inverse.
123 :
124 0 : historyLength = MIN(istep,maxHistoryLength)
125 :
126 0 : ALLOCATE(hessMat(vecLength,vecLength))
127 0 : ALLOCATE(invHessMat(vecLength,vecLength))
128 0 : ALLOCATE(tempMatB(vecLength,vecLength))
129 :
130 0 : hessMat(:,:) = 0.0
131 0 : invHessMat(:,:) = 0.0
132 0 : tempMatB(:,:) = 0.0
133 :
134 0 : DO i = 1, vecLength
135 0 : hessMat(i,i) = 1.0
136 0 : invHessMat(i,i) = 1.0
137 : END DO
138 :
139 0 : IF(historyLength.EQ.0) THEN
140 0 : DO i = 1, vecLength
141 0 : hessMat(i,i) = 1.0 / mixParam
142 0 : invHessMat(i,i) = mixParam
143 : END DO
144 : END IF
145 :
146 0 : DO i = 1, historyLength
147 0 : stepIndex = MOD(iStep - historyLength + i - 1,maxHistoryLength) + 1
148 :
149 : ! First the Hessian matrix
150 : ! second correction term
151 0 : tempMatB(:,:) = 0.0
152 0 : tempVecA(:) = 0.0
153 0 : CALL dgemv('N',vecLength,vecLength,1.0,hessMat(:,:),vecLength,paramCorrections(:,stepIndex),1,0.0,tempVecA(:),1)
154 0 : norm = ddot(vecLength,paramCorrections(:,stepIndex),1,tempVecA(:),1)
155 0 : CALL dgemm('N','N',vecLength,vecLength,1,-1.0/norm,tempVecA(:),vecLength,tempVecA(:),1,0.0,tempMatB,vecLength)
156 :
157 : ! first correction term
158 0 : norm = ddot(vecLength,gradientCorrections(:,stepIndex),1,paramCorrections(:,stepIndex),1)
159 0 : CALL dgemm('N','N',vecLength,vecLength,1,1.0/norm,gradientCorrections(:,stepIndex),vecLength,gradientCorrections(:,stepIndex),1,1.0,tempMatB,vecLength)
160 0 : hessMat(:,:) = hessMat(:,:) + tempMatB(:,:)
161 :
162 : ! Now the inverse of the Hessian matrix
163 0 : tempMatB(:,:) = 0.0
164 0 : tempVecA(:) = 0.0
165 0 : CALL dgemv('N',vecLength,vecLength,1.0,invHessMat(:,:),vecLength,gradientCorrections(:,stepIndex),1,0.0,tempVecA(:),1)
166 0 : norm = ddot(vecLength,paramCorrections(:,stepIndex),1,gradientCorrections(:,stepIndex),1)
167 :
168 : ! second correction term
169 0 : CALL dgemm('N','N',vecLength,vecLength,1,1.0,tempVecA(:),vecLength,paramCorrections(:,stepIndex),1,0.0,tempMatB,vecLength)
170 :
171 0 : tempMatB(:,:) = tempMatB(:,:) + TRANSPOSE(tempMatB(:,:))
172 0 : tempMatB(:,:) = -tempMatB(:,:) / norm
173 :
174 : ! first correction term
175 0 : temp = norm + ddot(vecLength,gradientCorrections(:,stepIndex),1,tempVecA(:),1)
176 0 : CALL dgemm('N','N',vecLength,vecLength,1,temp/(norm*norm),paramCorrections(:,stepIndex),vecLength,paramCorrections(:,stepIndex),1,1.0,tempMatB,vecLength)
177 0 : invHessMat(:,:) = invHessMat(:,:) + tempMatB(:,:)
178 :
179 : !Check whether invHessMat is the inverse of hessMat (for debugging only)
180 0 : tempMatB(:,:) = 0.0
181 0 : CALL dgemm('N','N',vecLength,vecLength,vecLength,1.0,hessMat(:,:),vecLength,invHessMat(:,:),vecLength,0.0,tempMatB,vecLength)
182 :
183 : END DO
184 :
185 : ! 2. Determine next point based on the gradient, the inverse of the Hessian, and the box constraints
186 :
187 0 : freeParams(:) = .TRUE.
188 0 : prelimNextPoint(:) = parameters(:)
189 0 : prelimGradient(:) = gradient(:)
190 :
191 0 : WRITE(5300,*) '=========================================================='
192 0 : WRITE(5300,*) 'bfgs_b2: Starting next point search: '
193 0 : WRITE(5300,*) '=========================================================='
194 :
195 0 : converged = .FALSE.
196 0 : DO WHILE (.NOT.converged)
197 :
198 0 : WRITE(5300,*) 'Point A: '
199 0 : WRITE(5300,*) 'prelimNextPoint: '
200 0 : WRITE(5300,'(5f20.10)') prelimNextPoint(:)
201 0 : WRITE(5300,*) 'prelimGradient: '
202 0 : WRITE(5300,'(5f20.10)') prelimGradient(:)
203 :
204 :
205 : ! Determine free parameters (parameters not on boundaries with a gradient pointing outwards)
206 0 : numFixedParams = 0
207 0 : DO i = 1, vecLength
208 0 : IF(enabledConstraints(i)) THEN
209 0 : IF((prelimNextPoint(i).LE.minConstraints(i)+eps).AND.(prelimGradient(i).GT.0.0)) THEN
210 0 : freeParams(i) = .FALSE.
211 0 : numFixedParams = numFixedParams + 1
212 : END IF
213 0 : IF((prelimNextPoint(i).GE.maxConstraints(i)-eps).AND.(prelimGradient(i).LT.0.0)) THEN
214 0 : freeParams(i) = .FALSE.
215 0 : numFixedParams = numFixedParams + 1
216 : END IF
217 : END IF
218 : END DO
219 0 : numFreeParams = vecLength - numFixedParams
220 :
221 : ! Perform basis transformation to obtain inverse Hessian matrix for the free parameter subspace
222 :
223 0 : ALLOCATE (subspaceBasis(vecLength,numFreeParams))
224 0 : subspaceBasis = 0.0
225 0 : subspaceGrad = 0.0
226 : j = 1
227 0 : DO i = 1, vecLength
228 0 : IF(freeParams(i)) THEN
229 0 : subspaceBasis(i,j) = 1.0
230 0 : subspaceGrad(j) = prelimGradient(i)
231 0 : j = j + 1
232 : END IF
233 : END DO
234 :
235 0 : ALLOCATE (tempMatC(vecLength,numFreeParams))
236 0 : ALLOCATE (subspaceInvHessMat(numFreeParams,numFreeParams))
237 0 : tempMatC(:,:) = 0.0
238 0 : subspaceInvHessMat(:,:) = 0.0
239 :
240 0 : WRITE(5000,*) '========================================================================='
241 0 : WRITE(5000,*) 'prelimNextPoint: '
242 0 : WRITE(5000,'(5f20.10)') prelimNextPoint(:)
243 0 : WRITE(5000,*) 'prelimGradient: '
244 0 : WRITE(5000,'(5f20.10)') prelimGradient(:)
245 0 : WRITE(5000,*) 'subspaceBasis: '
246 0 : DO i = 1, SIZE(subspaceBasis,2)
247 0 : WRITE(5000,'(5f20.10)') subspaceBasis(:,i)
248 : END DO
249 :
250 0 : CALL dgemm('N','N',vecLength,numFreeParams,vecLength,1.0,invHessMat(:,:),vecLength,subspaceBasis(:,:),vecLength,0.0,tempMatC(:,:),vecLength)
251 0 : CALL dgemm('N','N',numFreeParams,numFreeParams,vecLength,1.0,TRANSPOSE(subspaceBasis(:,:)),numFreeParams,tempMatC(:,:),vecLength,0.0,subspaceInvHessMat(:,:),numFreeParams)
252 :
253 0 : DEALLOCATE(tempMatC)
254 0 : DEALLOCATE(subspaceBasis)
255 :
256 : ! Calculate next point (ignoring constraints)
257 :
258 0 : tempVecA(:) = 0.0
259 :
260 0 : WRITE(5000,*) 'numFreeParams: ', numFreeParams
261 0 : WRITE(5000,*) 'subspaceInvHessMat: '
262 0 : DO i = 1, SIZE(subspaceInvHessMat,2)
263 0 : WRITE(5000,'(5f20.10)') subspaceInvHessMat(:,i)
264 : END DO
265 0 : WRITE(5000,*) '-subspaceGrad: '
266 0 : WRITE(5000,'(5f20.10)') -subspaceGrad(:numFreeParams)
267 0 : WRITE(5000,*) 'tempVecA: '
268 0 : WRITE(5000,'(5f20.10)') tempVecA(:numFreeParams)
269 :
270 0 : CALL dgemv('N',numFreeParams,numFreeParams,1.0,subspaceInvHessMat(:,:),numFreeParams,-subspaceGrad(:numFreeParams),1,0.0,tempVecA(:numFreeParams),1)
271 :
272 0 : WRITE(5000,*) '----------------------------'
273 0 : WRITE(5000,*) 'tempVecA: '
274 0 : WRITE(5000,'(5f20.10)') tempVecA(:numFreeParams)
275 0 : WRITE(5000,*) 'prelimNextPoint: '
276 0 : WRITE(5000,'(5f20.10)') prelimNextPoint(:)
277 :
278 : !!! TEST START
279 : ! temp = SQRT(ddot(numFreeParams,tempVecA(:numFreeParams),1,tempVecA(:numFreeParams),1))
280 0 : temp = MAXVAL(ABS(tempVecA(:numFreeParams-1)))
281 0 : maxDisp = 0.05
282 0 : IF (temp.GT.maxDisp) tempVecA(:numFreeParams) = tempVecA(:numFreeParams) / (temp/maxDisp)
283 : dampingFactor = 0.5
284 : !!! TEST END
285 :
286 : j = 1
287 0 : DO i = 1, vecLength
288 0 : IF(freeParams(i)) THEN
289 0 : prelimNextPoint(i) = prelimNextPoint(i) + dampingFactor * tempVecA(j)
290 : ! prelimNextPoint(i) = prelimNextPoint(i) + tempVecA(j)
291 0 : j = j + 1
292 : END IF
293 : END DO
294 :
295 0 : WRITE(5000,*) 'prelimNextPoint: '
296 0 : WRITE(5000,'(5f20.10)') prelimNextPoint(:)
297 :
298 0 : DEALLOCATE(subspaceInvHessMat)
299 :
300 : ! Constrain next point to box and check convergence
301 0 : converged = .TRUE.
302 0 : DO i = 1, vecLength
303 0 : IF (enabledConstraints(i).AND.freeParams(i)) THEN
304 : IF(prelimNextPoint(i).LT.minConstraints(i)-eps) THEN
305 : ! converged = .FALSE.
306 : END IF
307 : IF(prelimNextPoint(i).GT.maxConstraints(i)+eps) THEN
308 : ! converged = .FALSE.
309 : END IF
310 0 : IF(prelimNextPoint(i).LT.minConstraints(i)) THEN
311 0 : prelimNextPoint(i) = minConstraints(i)
312 : END IF
313 0 : IF(prelimNextPoint(i).GT.maxConstraints(i)) THEN
314 0 : prelimNextPoint(i) = maxConstraints(i)
315 : END IF
316 : END IF
317 : END DO
318 :
319 0 : WRITE(5000,*) 'prelimNextPoint: '
320 0 : WRITE(5000,'(5f20.10)') prelimNextPoint(:)
321 0 : WRITE(5000,*) 'prelimGradient: '
322 0 : WRITE(5000,'(5f20.10)') prelimGradient(:)
323 :
324 : ! If not converged calculate prelimGradient for next point according to quadratic model
325 : IF(.NOT.converged) THEN
326 : diffPoint(:) = prelimNextPoint(:) - parameters(:)
327 : diffGrad(:) = 0.0
328 : CALL dgemv('N',vecLength,vecLength,1.0,hessMat(:,:),vecLength,diffPoint(:),1,0.0,diffGrad(:),1)
329 : prelimGradient(:) = prelimGradient(:) + diffGrad(:)
330 : END IF
331 :
332 : ! Test: Set the gradient entry for the Lagrange multiplier explicitly
333 : ! (I hope this makes the approach more stable. The here introduced data is not automatically obtained)
334 0 : temp = ddot(vecLength,equalityLinCombi(:vecLength),1,prelimNextPoint(:vecLength),1)
335 0 : prelimGradient(vecLength) = temp - equalityCriterion ! derivative with respect to the Lagrange parameter
336 :
337 0 : WRITE(5000,*) 'prelimGradient: '
338 0 : WRITE(5000,'(5f20.10)') prelimGradient(:)
339 :
340 : END DO
341 :
342 : ! 3. Project prelimNextPoint onto the allowed subspace given by the equality constraints
343 : ! Note: I only consider a single equality constraint so this part becomes simpler
344 :
345 : ! Determine free parameters (parameters not on boundaries with a gradient pointing outwards)
346 :
347 : ! WRITE(5200,*) '=========================================================='
348 : ! WRITE(5200,*) 'bfgs_b2: Starting charge constraint ensurance: '
349 : ! WRITE(5200,*) '=========================================================='
350 :
351 : converged = .FALSE.
352 0 : DO WHILE (.NOT.converged)
353 :
354 : ! WRITE(5200,*) 'Point A: '
355 : ! WRITE(5200,*) 'prelimNextPoint: '
356 : ! WRITE(5200,'(5f20.10)') prelimNextPoint(:)
357 :
358 0 : freeParams(:) = .TRUE.
359 0 : numFixedParams = 0
360 0 : DO i = 1, vecLength
361 0 : IF(enabledConstraints(i)) THEN
362 0 : IF(prelimNextPoint(i).LT.minConstraints(i) - eps) THEN
363 0 : freeParams(i) = .FALSE.
364 0 : numFixedParams = numFixedParams + 1
365 : END IF
366 0 : IF(prelimNextPoint(i).GT.maxConstraints(i) + eps) THEN
367 0 : freeParams(i) = .FALSE.
368 0 : numFixedParams = numFixedParams + 1
369 : END IF
370 0 : IF(prelimNextPoint(i).LT.minConstraints(i)) THEN
371 0 : prelimNextPoint(i) = minConstraints(i)
372 : END IF
373 0 : IF(prelimNextPoint(i).GT.maxConstraints(i)) THEN
374 0 : prelimNextPoint(i) = maxConstraints(i)
375 : END IF
376 : END IF
377 : END DO
378 0 : numFreeParams = vecLength - numFixedParams
379 :
380 0 : subspaceNextPoint(:) = 0.0
381 0 : subspaceLinCombi(:) = 0.0
382 0 : subspaceEqualityCrit = equalityCriterion
383 0 : j = 1
384 0 : DO i = 1, vecLength
385 0 : IF(freeParams(i)) THEN
386 0 : subspaceNextPoint(j) = prelimNextPoint(i)
387 0 : subspaceLinCombi(j) = equalityLinCombi(i)
388 0 : j = j + 1
389 : ELSE
390 0 : subspaceEqualityCrit = subspaceEqualityCrit - prelimNextPoint(i) * equalityLinCombi(i)
391 : END IF
392 : END DO
393 :
394 0 : temp = ddot(numFreeParams,subspaceLinCombi(:),1,subspaceNextPoint(:),1)
395 0 : converged = .TRUE.
396 0 : IF(ABS(temp - subspaceEqualityCrit).GT.largeEps) THEN
397 0 : converged = .FALSE.
398 :
399 : ! WRITE(5100,*) 'subspaceNextPoint(:): '
400 : ! WRITE(5100,*) subspaceNextPoint(:)
401 : ! WRITE(5100,*) 'subspaceEqualityCrit: ', subspaceEqualityCrit
402 : ! WRITE(5100,*) 'temp: ', temp
403 :
404 0 : subspaceNextPoint(:) = subspaceNextPoint(:) * subspaceEqualityCrit / temp
405 : END IF
406 0 : j = 1
407 0 : DO i = 1, vecLength
408 0 : IF(freeParams(i)) THEN
409 0 : prelimNextPoint(i) = subspaceNextPoint(j)
410 0 : j = j + 1
411 : END IF
412 : END DO
413 :
414 : ! Also ensure that the equality criterion is fulfilled for the whole space.
415 : ! (This is only for pathological cases)
416 :
417 0 : temp = ddot(vecLength,equalityLinCombi(:vecLength),1,prelimNextPoint(:vecLength),1)
418 0 : IF(ABS(temp - equalityCriterion).GT.largeEps) THEN
419 0 : converged = .FALSE.
420 0 : prelimNextPoint(:vecLength) = prelimNextPoint(:vecLength) * equalityCriterion / temp
421 : END IF
422 :
423 : END DO
424 :
425 : ! 4. Update lastParameters, lastGradient, parameters
426 :
427 0 : lastParameters(:) = parameters(:)
428 0 : lastGradient(:) = gradient(:)
429 0 : parameters(:) = prelimNextPoint(:)
430 : ! parameters(:) = 0.25*(prelimNextPoint(:) - lastParameters(:)) + lastParameters(:)
431 :
432 0 : DEALLOCATE(tempMatB)
433 0 : DEALLOCATE(invHessMat)
434 0 : DEALLOCATE(hessMat)
435 :
436 0 : WRITE(*,*) 'bfgs_b2 (output): parameters(:vecLength): ', parameters(:vecLength)
437 :
438 : END SUBROUTINE bfgs_b2
439 :
440 : END MODULE m_bfgs_b2
|
