root / src / Hinvup_BFGS.f90 @ 3
Historique | Voir | Annoter | Télécharger (3,94 ko)
| 1 |
! This program reads an old hessian, old forces, new |
|---|---|
| 2 |
! forces and updates the hessian matrix inverse |
| 3 |
! using the BFGS scheme. |
| 4 |
! |
| 5 |
! This routines has been modified so that its input is now: |
| 6 |
SUBROUTINE Hinvup_BFGS(Nfree,ds,dgrad,Hess) |
| 7 |
|
| 8 |
IMPLICIT NONE |
| 9 |
|
| 10 |
INTEGER, PARAMETER :: KINT=KIND(1) |
| 11 |
INTEGER, PARAMETER :: KREAL=KIND(1.D0) |
| 12 |
|
| 13 |
INTEGER(KINT), INTENT(IN) :: NFREE |
| 14 |
REAL(KREAL), INTENT(INOUT) :: Hess(Nfree,Nfree) |
| 15 |
REAL(KREAL), INTENT(IN) :: ds(Nfree) |
| 16 |
REAL(KREAL), INTENT(IN) :: dGrad(Nfree) |
| 17 |
|
| 18 |
REAL(KREAL), ALLOCATABLE :: HessOld(:,:) , Wy(:) |
| 19 |
REAL(KREAL), ALLOCATABLE :: Hup(:,:),ssT(:,:) |
| 20 |
REAL(KREAL), ALLOCATABLE :: HupInt(:) |
| 21 |
REAL(KREAL), ALLOCATABLE :: Num(:,:), Num2(:,:) |
| 22 |
REAL(KREAL) :: Phi, Num3, Den, Den2,NormDS |
| 23 |
CHARACTER(120) :: fmt,fmt2,Line |
| 24 |
INTEGER(KINT) :: i,j,k, Idx |
| 25 |
|
| 26 |
LOGICAL debug,valid |
| 27 |
EXTERNAL valid |
| 28 |
|
| 29 |
debug=valid("Hinvup_BFGS")
|
| 30 |
|
| 31 |
if (debug) WRITE(*,*) "=========================== Entering Hinvup_BFGS ====================" |
| 32 |
|
| 33 |
ALLOCATE(HessOld(nfree,nfree), Hup(nfree,nfree), & |
| 34 |
HupInt(nfree), & |
| 35 |
Num(nfree,nfree), Num2(nfree,nfree),Wy(Nfree),ssT(Nfree,nFree) ) |
| 36 |
|
| 37 |
! Notations |
| 38 |
! s=Geom-GeomOld |
| 39 |
! y=Grad-GradOld |
| 40 |
! W= (H)-1 ie the inverse matrix that we are updating |
| 41 |
! W+=W - [s(yT)W+Wy(sT)]/(y,s)+ {1+[(y,Wy)/(y,s)]}(s(sT))/(y,s)
|
| 42 |
|
| 43 |
|
| 44 |
HessOld=Hess |
| 45 |
|
| 46 |
if (debug) THEN |
| 47 |
WRITE(*,*) 'DBG Hinvup_BFGS HessOld' |
| 48 |
Idx=min(12,Nfree) |
| 49 |
DO I=1,NFree |
| 50 |
WRITE(*,'(12(1X,F8.4))') HessOld(I,1:Idx) |
| 51 |
END DO |
| 52 |
WRITE(*,*) 'DBG Hinvup_BFGS ds' |
| 53 |
WRITE(*,'(12(1X,F8.4))') ds(1:Idx) |
| 54 |
WRITE(*,*) 'DBG Hinvup_BFGS dGrad' |
| 55 |
WRITE(*,'(12(1X,F8.4))') dGrad(1:Idx) |
| 56 |
END IF |
| 57 |
|
| 58 |
|
| 59 |
|
| 60 |
NormDS=sqrt(DOT_PRODUCT(DS,DS)) |
| 61 |
if (NormDS>=1e-6) THEN |
| 62 |
IF (debug) WRITE(*,*) 'NormDS=',NormDS |
| 63 |
IF (debug) WRITE(*,*) 'DS=',DS |
| 64 |
IF (debug) WRITE(*,*) 'DGrad=',DGrad |
| 65 |
! We calculate the denominator (DGrad,DS) |
| 66 |
den=DOT_PRODUCT(DGrad,DS) |
| 67 |
|
| 68 |
if (debug) WRITE(*,*) "(y,s)=(yT)s",den |
| 69 |
|
| 70 |
! We calculate Wy |
| 71 |
Wy=0 |
| 72 |
DO I=1,nfree |
| 73 |
DO J=1,NFree |
| 74 |
Wy(I)=Wy(I)+HessOld(i,j)*DGrad(j) |
| 75 |
END DO |
| 76 |
END DO |
| 77 |
|
| 78 |
if (debug) WRITE(*,*) "Wy=",Wy |
| 79 |
|
| 80 |
! We calculate Num=Wy(sT) |
| 81 |
Num=0. |
| 82 |
DO I=1,Nfree |
| 83 |
DO J=1,Nfree |
| 84 |
Num(I,J)=Wy(i)*Ds(j) |
| 85 |
END DO |
| 86 |
END DO |
| 87 |
|
| 88 |
if (debug) THEN |
| 89 |
WRITE(*,*) "Wy(sT)" |
| 90 |
DO I=1,NFree |
| 91 |
WRITE(*,*) Num(I,:) |
| 92 |
END DO |
| 93 |
END if |
| 94 |
|
| 95 |
! We calculate Num2=s(yT) |
| 96 |
Num2=0 |
| 97 |
Do I=1,NFree |
| 98 |
DO J=1,NFree |
| 99 |
Num2(I,J)=Ds(i)*DGrad(j) |
| 100 |
END DO |
| 101 |
END DO |
| 102 |
|
| 103 |
if (debug) THEN |
| 104 |
WRITE(*,*) "s(yT)" |
| 105 |
DO I=1,NFree |
| 106 |
WRITE(*,*) Num2(I,:) |
| 107 |
END DO |
| 108 |
END if |
| 109 |
|
| 110 |
|
| 111 |
! we add Num2*W to Num |
| 112 |
Do I=1,NFree |
| 113 |
DO J=1,NFree |
| 114 |
DO K=1,NFree |
| 115 |
Num(I,J)=Num(I,J)+Num2(I,K)*HessOld(K,J) |
| 116 |
END DO |
| 117 |
END DO |
| 118 |
END DO |
| 119 |
|
| 120 |
if (debug) THEN |
| 121 |
WRITE(*,*) "Num=s(yT)W+Wy(sT)" |
| 122 |
DO I=1,NFree |
| 123 |
WRITE(*,*) Num(I,:) |
| 124 |
END DO |
| 125 |
END if |
| 126 |
|
| 127 |
|
| 128 |
! we calculate ssT |
| 129 |
ssT=0. |
| 130 |
DO I=1,Nfree |
| 131 |
DO J=1,NFree |
| 132 |
ssT(I,J)=Ds(I)*Ds(J) |
| 133 |
END DO |
| 134 |
END DO |
| 135 |
|
| 136 |
if (debug) THEN |
| 137 |
WRITE(*,*) "s(sT)" |
| 138 |
DO I=1,NFree |
| 139 |
WRITE(*,*) ssT(I,:) |
| 140 |
END DO |
| 141 |
END if |
| 142 |
|
| 143 |
|
| 144 |
!We add everything |
| 145 |
Hup=(1+DOT_PRODUCT(DGRAD,Wy)/Den)*ssT/Den-Num/Den |
| 146 |
Hess=HessOld+Hup |
| 147 |
|
| 148 |
if (debug) THEN |
| 149 |
WRITE(*,*) 'DBG Hinvup_BFGS Hup' |
| 150 |
Idx=min(12,Nfree) |
| 151 |
DO I=1,NFree |
| 152 |
WRITE(*,'(12(1X,F8.4))') Hup(I,1:Idx) |
| 153 |
END DO |
| 154 |
WRITE(*,*) 'DBG Hinvup_BFGS Hess new' |
| 155 |
Idx=min(12,Nfree) |
| 156 |
DO I=1,NFree |
| 157 |
WRITE(*,'(12(1X,F8.4))') Hess(I,1:Idx) |
| 158 |
END DO |
| 159 |
END IF |
| 160 |
|
| 161 |
ELSE |
| 162 |
Hess=HessOld |
| 163 |
END IF |
| 164 |
! WRITE(*,*) 'New hess' |
| 165 |
! DO I=1,nfree |
| 166 |
! WRITE(*,fmt) Hess(i,:) |
| 167 |
! END DO |
| 168 |
|
| 169 |
DEALLOCATE(HessOld, Wy,ssT, Hup,HupInt, Num, Num2 ) |
| 170 |
|
| 171 |
if (debug) WRITE(*,*) "=========================== Hinvup_BFGS Over ====================" |
| 172 |
|
| 173 |
END SUBROUTINE Hinvup_BFGS |