LCOV - code coverage report
Current view: top level - rdmft - bfgs_b2.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 163 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          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             :    IMPLICIT NONE
      27             : 
      28             :    INTEGER,        INTENT(IN)    :: vecLength
      29             :    INTEGER,        INTENT(IN)    :: maxHistoryLength
      30             :    INTEGER,        INTENT(INOUT) :: iStep ! The index of the bfgs_b iteration
      31             :    REAL,           INTENT(IN)    :: mixParam
      32             :    REAL,           INTENT(IN)    :: equalityCriterion
      33             :    REAL,           INTENT(IN)    :: gradient(vecLength)
      34             :    REAL,           INTENT(INOUT) :: lastGradient(vecLength)
      35             :    REAL,           INTENT(IN)    :: minConstraints(vecLength)
      36             :    REAL,           INTENT(IN)    :: maxConstraints(vecLength)
      37             :    LOGICAL,        INTENT(IN)    :: enabledConstraints(vecLength)
      38             :    REAL,           INTENT(INOUT) :: parameters(vecLength)
      39             :    REAL,           INTENT(INOUT) :: lastParameters(vecLength)
      40             :    REAL,           INTENT(IN)    :: equalityLinCombi(vecLength)
      41             :    REAL,           INTENT(INOUT) :: paramCorrections(vecLength,maxHistoryLength)
      42             :    REAL,           INTENT(INOUT) :: gradientCorrections(vecLength,maxHistoryLength)
      43             :    LOGICAL,        INTENT(OUT)   :: l_converged
      44             :    REAL,           INTENT(IN)    :: convCrit
      45             : 
      46             :    INTEGER                :: i, j, stepIndex, numFreeParams, numFixedParams
      47             :    INTEGER                :: historyLength
      48             :    REAL, PARAMETER        :: eps = 1.0e-10
      49             :    REAL, PARAMETER        :: largeEps = 1.0e-6
      50             :    REAL                   :: temp, norm, subspaceEqualityCrit
      51             :    REAL                   :: ddot
      52             : 
      53             :    LOGICAL                :: converged
      54             : 
      55           0 :    REAL                   :: lastGradCorrection(vecLength)
      56           0 :    REAL                   :: lastParamCorrection(vecLength)
      57           0 :    REAL                   :: prelimNextPoint(vecLength)
      58           0 :    REAL                   :: prelimGradient(vecLength)
      59           0 :    REAL                   :: tempVecA(vecLength)
      60           0 :    REAL                   :: subspaceGrad(vecLength)
      61           0 :    REAL                   :: subspaceNextPoint(vecLength)
      62           0 :    REAL                   :: subspaceLinCombi(vecLength)
      63           0 :    REAL                   :: diffPoint(vecLength)
      64           0 :    REAL                   :: diffGrad(vecLength)
      65             : 
      66           0 :    LOGICAL                :: freeParams(vecLength)
      67             : 
      68           0 :    REAL, ALLOCATABLE      :: hessMat(:,:), invHessMat(:,:)
      69           0 :    REAL, ALLOCATABLE      :: tempMatB(:,:), tempMatC(:,:)
      70           0 :    REAL, ALLOCATABLE      :: subspaceBasis(:,:)
      71           0 :    REAL, ALLOCATABLE      :: subspaceInvHessMat(:,:)
      72             : 
      73           0 :    WRITE(*,*) 'bfgs_b2 (input): parameters(:vecLength): ', parameters(:vecLength)
      74           0 :    WRITE(*,*) 'bfgs_b2 (input): gradient(:vecLength): ', gradient(:vecLength)
      75             : 
      76           0 :    l_converged = .FALSE.
      77             : 
      78             :    ! 0. Add last data pair to history iff it is not obviously destructive for the positive definiteness of the Hessian matrix
      79             : 
      80           0 :    IF(ANY(lastGradient(:).NE.0.0).OR.ANY(lastParameters(:).NE.0.0)) THEN
      81           0 :       lastGradCorrection(:) = gradient(:) - lastGradient(:)
      82           0 :       lastParamCorrection(:) = parameters(:) - lastParameters(:)
      83           0 :       temp = ddot(vecLength,lastParamCorrection(:),1,lastGradCorrection(:),1)
      84           0 :       IF (temp.GT.eps) THEN
      85             :          ! Add data pair to history
      86           0 :          istep = istep + 1
      87           0 :          stepIndex = MOD(iStep - 1,maxHistoryLength) + 1
      88           0 :          paramCorrections(:,stepIndex) = lastParamCorrection(:)
      89           0 :          gradientCorrections(:,stepIndex) = lastGradCorrection(:)
      90             :       END IF
      91             : 
      92             :       ! Check for convergence and return if converged
      93           0 :       IF(ALL(ABS(lastGradCorrection(:)).LT.convCrit).AND.ALL(ABS(lastParamCorrection(:)).LT.convCrit)) THEN
      94           0 :          l_converged = .TRUE.
      95           0 :          RETURN
      96             :       END IF
      97             :    END IF
      98             : 
      99             :    ! 1. Construct Hessian matrix and its inverse.
     100             : 
     101           0 :    historyLength = MIN(istep,maxHistoryLength)
     102             : 
     103           0 :    ALLOCATE(hessMat(vecLength,vecLength))
     104           0 :    ALLOCATE(invHessMat(vecLength,vecLength))
     105           0 :    ALLOCATE(tempMatB(vecLength,vecLength))
     106             : 
     107           0 :    hessMat(:,:) = 0.0
     108           0 :    invHessMat(:,:) = 0.0
     109           0 :    tempMatB(:,:) = 0.0
     110             : 
     111           0 :    DO i = 1, vecLength
     112           0 :       hessMat(i,i) = 1.0
     113           0 :       invHessMat(i,i) = 1.0
     114             :    END DO
     115             : 
     116           0 :    IF(historyLength.EQ.0) THEN
     117           0 :       DO i = 1, vecLength
     118           0 :          hessMat(i,i) = 1.0 / mixParam
     119           0 :          invHessMat(i,i) = mixParam
     120             :       END DO
     121             :    END IF
     122             : 
     123           0 :    DO i = 1, historyLength
     124           0 :       stepIndex = MOD(iStep - historyLength + i - 1,maxHistoryLength) + 1
     125             : 
     126             :       ! First the Hessian matrix
     127             :       ! second correction term
     128           0 :       tempMatB(:,:) = 0.0
     129           0 :       tempVecA(:) = 0.0
     130           0 :       CALL dgemv('N',vecLength,vecLength,1.0,hessMat(:,:),vecLength,paramCorrections(:,stepIndex),1,0.0,tempVecA(:),1)
     131           0 :       norm = ddot(vecLength,paramCorrections(:,stepIndex),1,tempVecA(:),1)
     132           0 :       CALL dgemm('N','N',vecLength,vecLength,1,-1.0/norm,tempVecA(:),vecLength,tempVecA(:),1,0.0,tempMatB,vecLength)
     133             : 
     134             :       ! first correction term
     135           0 :       norm = ddot(vecLength,gradientCorrections(:,stepIndex),1,paramCorrections(:,stepIndex),1)
     136           0 :       CALL dgemm('N','N',vecLength,vecLength,1,1.0/norm,gradientCorrections(:,stepIndex),vecLength,gradientCorrections(:,stepIndex),1,1.0,tempMatB,vecLength)
     137           0 :       hessMat(:,:) = hessMat(:,:) + tempMatB(:,:)
     138             : 
     139             :       ! Now the inverse of the Hessian matrix
     140           0 :       tempMatB(:,:) = 0.0
     141           0 :       tempVecA(:) = 0.0
     142           0 :       CALL dgemv('N',vecLength,vecLength,1.0,invHessMat(:,:),vecLength,gradientCorrections(:,stepIndex),1,0.0,tempVecA(:),1)
     143           0 :       norm = ddot(vecLength,paramCorrections(:,stepIndex),1,gradientCorrections(:,stepIndex),1)
     144             : 
     145             :       ! second correction term
     146           0 :       CALL dgemm('N','N',vecLength,vecLength,1,1.0,tempVecA(:),vecLength,paramCorrections(:,stepIndex),1,0.0,tempMatB,vecLength)
     147             : 
     148           0 :       tempMatB(:,:) = tempMatB(:,:) + TRANSPOSE(tempMatB(:,:))
     149           0 :       tempMatB(:,:) = -tempMatB(:,:) / norm
     150             : 
     151             :       ! first correction term
     152           0 :       temp = norm + ddot(vecLength,gradientCorrections(:,stepIndex),1,tempVecA(:),1)
     153           0 :       CALL dgemm('N','N',vecLength,vecLength,1,temp/(norm*norm),paramCorrections(:,stepIndex),vecLength,paramCorrections(:,stepIndex),1,1.0,tempMatB,vecLength)
     154           0 :       invHessMat(:,:) = invHessMat(:,:) + tempMatB(:,:)
     155             : 
     156             :       !Check whether invHessMat is the inverse of hessMat (for debugging only)
     157           0 :       tempMatB(:,:) = 0.0
     158           0 :       CALL dgemm('N','N',vecLength,vecLength,vecLength,1.0,hessMat(:,:),vecLength,invHessMat(:,:),vecLength,0.0,tempMatB,vecLength)      
     159             : 
     160             :    END DO
     161             : 
     162             :    ! 2. Determine next point based on the gradient, the inverse of the Hessian, and the box constraints
     163             : 
     164           0 :    freeParams(:) = .TRUE.
     165           0 :    prelimNextPoint(:) = parameters(:)
     166           0 :    prelimGradient(:) = gradient(:)
     167             : 
     168             :    converged = .FALSE.
     169           0 :    DO WHILE (.NOT.converged)
     170             : 
     171             :       ! Determine free parameters (parameters not on boundaries with a gradient pointing outwards)
     172             :       numFixedParams = 0
     173           0 :       DO i = 1, vecLength
     174           0 :          IF(enabledConstraints(i)) THEN
     175           0 :             IF((prelimNextPoint(i).LE.minConstraints(i)+eps).AND.(prelimGradient(i).GT.0.0)) THEN
     176           0 :                freeParams(i) = .FALSE.
     177           0 :                numFixedParams = numFixedParams + 1
     178             :             END IF
     179           0 :             IF((prelimNextPoint(i).GE.maxConstraints(i)-eps).AND.(prelimGradient(i).LT.0.0)) THEN
     180           0 :                freeParams(i) = .FALSE.
     181           0 :                numFixedParams = numFixedParams + 1
     182             :             END IF
     183             :          END IF
     184             :       END DO
     185           0 :       numFreeParams = vecLength - numFixedParams
     186             : 
     187             :       ! Perform basis transformation to obtain inverse Hessian matrix for the free parameter subspace
     188             : 
     189           0 :       ALLOCATE (subspaceBasis(vecLength,numFreeParams))
     190           0 :       subspaceBasis = 0.0
     191           0 :       subspaceGrad = 0.0
     192             :       j = 1
     193           0 :       DO i = 1, vecLength
     194           0 :          IF(freeParams(i)) THEN
     195           0 :             subspaceBasis(i,j) = 1.0
     196           0 :             subspaceGrad(j) = prelimGradient(i)
     197           0 :             j = j + 1
     198             :          END IF
     199             :       END DO
     200             : 
     201           0 :       ALLOCATE (tempMatC(vecLength,numFreeParams))
     202           0 :       ALLOCATE (subspaceInvHessMat(numFreeParams,numFreeParams))
     203           0 :       tempMatC(:,:) = 0.0
     204           0 :       subspaceInvHessMat(:,:) = 0.0
     205             : 
     206           0 :       CALL dgemm('N','N',vecLength,numFreeParams,vecLength,1.0,invHessMat(:,:),vecLength,subspaceBasis(:,:),vecLength,0.0,tempMatC(:,:),vecLength)
     207           0 :       CALL dgemm('N','N',numFreeParams,numFreeParams,vecLength,1.0,TRANSPOSE(subspaceBasis(:,:)),numFreeParams,tempMatC(:,:),vecLength,0.0,subspaceInvHessMat(:,:),numFreeParams)
     208             :       
     209           0 :       DEALLOCATE(tempMatC)
     210           0 :       DEALLOCATE(subspaceBasis)
     211             : 
     212             :       ! Calculate next point (ignoring constraints)
     213             : 
     214           0 :       tempVecA(:) = 0.0
     215           0 :       CALL dgemv('N',numFreeParams,numFreeParams,1.0,subspaceInvHessMat(:,:),numFreeParams,-subspaceGrad(:numFreeParams),1,0.0,tempVecA(:numFreeParams),1)
     216             : 
     217           0 :       j = 1
     218           0 :       DO i = 1, vecLength
     219           0 :          IF(freeParams(i)) THEN
     220           0 :             prelimNextPoint(i) = prelimNextPoint(i) + tempVecA(j)
     221           0 :             j = j + 1
     222             :          END IF
     223             :       END DO
     224             : 
     225           0 :       DEALLOCATE(subspaceInvHessMat)
     226             : 
     227             :       ! Constrain next point to box and check convergence
     228           0 :       converged = .TRUE.
     229           0 :       DO i = 1, vecLength
     230           0 :          IF (enabledConstraints(i).AND.freeParams(i)) THEN
     231           0 :             IF(prelimNextPoint(i).LT.minConstraints(i)-eps) THEN
     232           0 :                converged = .FALSE.
     233             :             END IF
     234           0 :             IF(prelimNextPoint(i).GT.maxConstraints(i)+eps) THEN
     235           0 :                converged = .FALSE.
     236             :             END IF
     237           0 :             IF(prelimNextPoint(i).LT.minConstraints(i)) THEN
     238           0 :                prelimNextPoint(i) = minConstraints(i)
     239             :             END IF
     240           0 :             IF(prelimNextPoint(i).GT.maxConstraints(i)) THEN
     241           0 :                prelimNextPoint(i) = maxConstraints(i)
     242             :             END IF
     243             :          END IF
     244             :       END DO
     245             : 
     246             : 
     247             :       ! If not converged calculate prelimGradient for next point according to quadratic model
     248           0 :       IF(.NOT.converged) THEN
     249           0 :          diffPoint(:) = prelimNextPoint(:) - parameters(:)
     250           0 :          diffGrad(:) = 0.0
     251           0 :          CALL dgemv('N',vecLength,vecLength,1.0,hessMat(:,:),vecLength,diffPoint(:),1,0.0,diffGrad(:),1)
     252           0 :          prelimGradient(:) = prelimGradient(:) + diffGrad(:)
     253             :       END IF
     254             : 
     255             :    END DO
     256             : 
     257             :    ! 3. Project prelimNextPoint onto the allowed subspace given by the equality constraints
     258             :    ! Note: I only consider a single equality constraint so this part becomes simpler
     259             : 
     260             :    ! Determine free parameters (parameters not on boundaries with a gradient pointing outwards)
     261             : 
     262             :    converged = .FALSE.
     263           0 :    DO WHILE (.NOT.converged)
     264             : 
     265           0 :       freeParams(:) = .TRUE.
     266             :       numFixedParams = 0
     267           0 :       DO i = 1, vecLength
     268           0 :          IF(enabledConstraints(i)) THEN
     269           0 :             IF(prelimNextPoint(i).LT.minConstraints(i) - eps) THEN
     270           0 :                freeParams(i) = .FALSE.
     271           0 :                numFixedParams = numFixedParams + 1
     272             :             END IF
     273           0 :             IF(prelimNextPoint(i).GT.maxConstraints(i) + eps) THEN
     274           0 :                freeParams(i) = .FALSE.
     275           0 :                numFixedParams = numFixedParams + 1
     276             :             END IF
     277           0 :             IF(prelimNextPoint(i).LT.minConstraints(i)) THEN
     278           0 :                prelimNextPoint(i) = minConstraints(i)
     279             :             END IF
     280           0 :             IF(prelimNextPoint(i).GT.maxConstraints(i)) THEN
     281           0 :                prelimNextPoint(i) = maxConstraints(i)
     282             :             END IF
     283             :          END IF
     284             :       END DO
     285           0 :       numFreeParams = vecLength - numFixedParams
     286             : 
     287           0 :       subspaceNextPoint(:) = 0.0
     288           0 :       subspaceLinCombi(:) = 0.0
     289           0 :       subspaceEqualityCrit = equalityCriterion
     290           0 :       j = 1
     291           0 :       DO i = 1, vecLength
     292           0 :          IF(freeParams(i)) THEN
     293           0 :             subspaceNextPoint(j) = prelimNextPoint(i)
     294           0 :             subspaceLinCombi(j) = equalityLinCombi(i)
     295           0 :             j = j + 1
     296             :          ELSE
     297           0 :             subspaceEqualityCrit = subspaceEqualityCrit - prelimNextPoint(i) * equalityLinCombi(i)
     298             :          END IF
     299             :       END DO
     300           0 :       temp = ddot(numFreeParams,subspaceLinCombi(:),1,subspaceNextPoint(:),1)
     301           0 :       converged = .TRUE.
     302           0 :       IF(ABS(temp - subspaceEqualityCrit).GT.largeEps) THEN
     303             :          converged = .FALSE.
     304           0 :          subspaceNextPoint(:) = subspaceNextPoint(:) * subspaceEqualityCrit / temp
     305             :       END IF
     306           0 :       j = 1
     307           0 :       DO i = 1, vecLength
     308           0 :          IF(freeParams(i)) THEN
     309           0 :             prelimNextPoint(i) = subspaceNextPoint(j)
     310           0 :             j = j + 1
     311             :          END IF
     312             :       END DO
     313             :    END DO
     314             : 
     315             :    ! 4. Update lastParameters, lastGradient, parameters
     316             : 
     317           0 :    lastParameters(:) = parameters(:)
     318           0 :    lastGradient(:) = gradient(:)
     319           0 :    parameters(:) = prelimNextPoint(:)
     320             : 
     321           0 :    DEALLOCATE(tempMatB)
     322           0 :    DEALLOCATE(invHessMat)
     323           0 :    DEALLOCATE(hessMat)
     324             : 
     325           0 :    WRITE(*,*) 'bfgs_b2 (output): parameters(:vecLength): ', parameters(:vecLength)
     326             : 
     327             : END SUBROUTINE bfgs_b2
     328             : 
     329             : END MODULE m_bfgs_b2

Generated by: LCOV version 1.13