hdiff output

r31239/lbfgs.f90 2016-10-05 13:30:09.253543952 +0100 r31238/lbfgs.f90 2016-10-05 13:30:09.509547359 +0100
 45:  45: 
 46:      INTEGER,INTENT(IN) :: D,U  ! DIMENSIONALITY OF THE PROBLEM AND NUMBER OF UPDATES 46:      INTEGER,INTENT(IN) :: D,U  ! DIMENSIONALITY OF THE PROBLEM AND NUMBER OF UPDATES
 47:      INTEGER NPERSIST           ! number of persistent minima 47:      INTEGER NPERSIST           ! number of persistent minima
 48:      INTEGER PERSISTTHRESH      ! persistence threshold 48:      INTEGER PERSISTTHRESH      ! persistence threshold
 49:      LOGICAL PERSISTENT(NIMAGE+2) ! logical array to identify persistent minima 49:      LOGICAL PERSISTENT(NIMAGE+2) ! logical array to identify persistent minima
 50:      INTEGER NTIMESMIN(NIMAGE+2)  ! number of consecutive steps this image is a local minimum 50:      INTEGER NTIMESMIN(NIMAGE+2)  ! number of consecutive steps this image is a local minimum
 51:      LOGICAL NOTNEW  51:      LOGICAL NOTNEW 
 52:      DOUBLE PRECISION DIJ, DMAX, DS, DF, EWORST, EMAX, BESTE, ENOLD 52:      DOUBLE PRECISION DIJ, DMAX, DS, DF, EWORST, EMAX, BESTE, ENOLD
 53:      DOUBLE PRECISION, ALLOCATABLE :: REPTEMP(:) 53:      DOUBLE PRECISION, ALLOCATABLE :: REPTEMP(:)
 54:      INTEGER, ALLOCATABLE :: IREPTEMP(:) 54:      INTEGER, ALLOCATABLE :: IREPTEMP(:)
 55:      INTEGER MAXIM, NDECREASE, NFAIL, NRED 55:      INTEGER MAXIM, NDECREASE, NFAIL
 56:      LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM 56:      LOGICAL KNOWE, KNOWG, KNOWH, ADDATOM
 57:      COMMON /KNOWN/ KNOWE, KNOWG, KNOWH 57:      COMMON /KNOWN/ KNOWE, KNOWG, KNOWH
 58:  58: 
 59:      INTEGER :: J2,POINT,BOUND,NPT,CP,I,ISTAT,J1,JMINUS,JPLUS,J3,JDO,J4,NPEPFAIL,NBADTOTAL 59:      INTEGER :: J2,POINT,BOUND,NPT,CP,I,ISTAT,J1,JMINUS,JPLUS,J3,JDO,J4,NPEPFAIL,NBADTOTAL
 60:      INTEGER AT1(NATOMS),AT2(NATOMS),AT3(NATOMS),AT4(NATOMS) 60:      INTEGER AT1(NATOMS),AT2(NATOMS),AT3(NATOMS),AT4(NATOMS)
 61:      DOUBLE PRECISION :: YS,YY,SQ,YR,BETA,GNORM,DDOT,DUMMY,STPMIN,PREVGRAD,MEANSEP,EMINUS,EPLUS, & 61:      DOUBLE PRECISION :: YS,YY,SQ,YR,BETA,GNORM,DDOT,DUMMY,STPMIN,PREVGRAD,MEANSEP,EMINUS,EPLUS, &
 62:   &                      LCOORDS(NOPT), ENERGY, INVDTOACTIVE(NATOMS), STIME 62:   &                      LCOORDS(NOPT), ENERGY, INVDTOACTIVE(NATOMS), STIME
 63:      DOUBLE PRECISION,DIMENSION(D)     :: GTMP,DIAG,STP 63:      DOUBLE PRECISION,DIMENSION(D)     :: GTMP,DIAG,STP
 64:      DOUBLE PRECISION,DIMENSION(U)     :: RHO1,ALPHA 64:      DOUBLE PRECISION,DIMENSION(U)     :: RHO1,ALPHA
 65:      DOUBLE PRECISION,DIMENSION(0:U,D) :: SEARCHSTEP,GDIF 65:      DOUBLE PRECISION,DIMENSION(0:U,D) :: SEARCHSTEP,GDIF
117:           POINT = 0117:           POINT = 0
118:           IF (.NOT.DIAGCO) DIAG(1:D)=NEBDGUESS118:           IF (.NOT.DIAGCO) DIAG(1:D)=NEBDGUESS
119:           SEARCHSTEP(0,:)= -G*DIAG                           ! NR STEP FOR DIAGONAL INVERSE HESSIAN119:           SEARCHSTEP(0,:)= -G*DIAG                           ! NR STEP FOR DIAGONAL INVERSE HESSIAN
120:           GTMP           = SEARCHSTEP(0,:)120:           GTMP           = SEARCHSTEP(0,:)
121:           GNORM          = MAX(SQRT(DOT_PRODUCT(G,G)),1.0D-100)121:           GNORM          = MAX(SQRT(DOT_PRODUCT(G,G)),1.0D-100)
122:           STP            = MIN(1.0D0/GNORM, GNORM)           ! MAKE THE FIRST GUESS FOR THE STEP LENGTH CAUTIOUS122:           STP            = MIN(1.0D0/GNORM, GNORM)           ! MAKE THE FIRST GUESS FOR THE STEP LENGTH CAUTIOUS
123:      ELSE MAIN123:      ELSE MAIN
124:           BOUND=NITERDONE-1124:           BOUND=NITERDONE-1
125:           IF (NITERDONE.GT.NEBMUPDATE) BOUND=NEBMUPDATE125:           IF (NITERDONE.GT.NEBMUPDATE) BOUND=NEBMUPDATE
126:           YS=DOT_PRODUCT( GDIF(NPT/D,:), SEARCHSTEP(NPT/D,:)  )126:           YS=DOT_PRODUCT( GDIF(NPT/D,:), SEARCHSTEP(NPT/D,:)  )
127:           IF (ABS(YS).LT.1.0D-20) THEN127:           IF (YS==0.0D0) YS=1.0D0
128:               YS=SIGN(1.0D-20,YS) 
129:           ENDIF  
130:  128:  
131:           ! Update estimate of diagonal inverse Hessian elements.129:           ! Update estimate of diagonal inverse Hessian elements.
132:           ! We divide by both YS and YY at different points, so they had better not be zero!130:           ! We divide by both YS and YY at different points, so they had better not be zero!
133:           IF (.NOT.DIAGCO) THEN131:           IF (.NOT.DIAGCO) THEN
134:                YY=DOT_PRODUCT( GDIF(NPT/D,:) , GDIF(NPT/D,:) )132:                YY=DOT_PRODUCT( GDIF(NPT/D,:) , GDIF(NPT/D,:) )
135:                 IF (ABS(YY).LT.1.0D-20) THEN133:                IF (YY==0.0D0) YY=1.0D0
136:                     YY=SIGN(1.0D-20,YY)134: !              DIAG = ABS(YS/YY)
137:                 ENDIF 
138:                DIAG = YS/YY135:                DIAG = YS/YY
139:           ELSE136:           ELSE
140:                CALL CHECKINPUT137:                CALL CHECKINPUT
141:           ENDIF138:           ENDIF
142:  139:  
143:           ! COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980, "Updating quasi-Newton matrices with limited storage",140:           ! COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980, "Updating quasi-Newton matrices with limited storage",
144:           ! Mathematics of Computation, Vol.35, No.151, pp. 773-782141:           ! Mathematics of Computation, Vol.35, No.151, pp. 773-782
145:           CP= POINT; IF (POINT==0) CP = NEBMUPDATE142:           CP= POINT; IF (POINT==0) CP = NEBMUPDATE
146:           RHO1(CP)=1.0D0/YS143:           RHO1(CP)=1.0D0/YS
147:           GTMP = -G144:           GTMP = -G
148:           CP= POINT 145:           CP= POINT 
149:  146:  
150:           DO I= 1,BOUND 147:           DO I= 1,BOUND 
151: !              CP = CP - 1; IF (CP == -1) CP = M - 1148: !              CP = CP - 1; IF (CP == -1) CP = M - 1
152:                CP = CP - 1; IF (CP == -1) CP = NEBMUPDATE - 1149:                CP = CP - 1; IF (CP == -1) CP = NEBMUPDATE - 1
153:                SQ= DOT_PRODUCT( SEARCHSTEP(CP,:),GTMP )150:                SQ= DOT_PRODUCT( SEARCHSTEP(CP,:),GTMP )
154:                ALPHA(CP+1) = RHO1(CP+1) * SQ151:                ALPHA(CP+1) = RHO1(CP+1) * SQ
155:                GTMP        = -ALPHA(CP+1)*GDIF(CP,:) + GTMP152:                GTMP        = -ALPHA(CP+1)*GDIF(CP,:) + GTMP
156:           ENDDO153:           ENDDO
157: 154:               
158:           GTMP=DIAG*GTMP155:           GTMP=DIAG*GTMP
159: 156: 
160:           DO I=1,BOUND157:           DO I=1,BOUND
161:                YR= DOT_PRODUCT( GDIF(CP,:) , GTMP )158:                YR= DOT_PRODUCT( GDIF(CP,:) , GTMP )
162:                BETA= RHO1(CP+1)*YR159:                BETA= RHO1(CP+1)*YR
163:                BETA= ALPHA(CP+1)-BETA160:                BETA= ALPHA(CP+1)-BETA
164:                GTMP = BETA*SEARCHSTEP(CP,:) + GTMP161:                GTMP = BETA*SEARCHSTEP(CP,:) + GTMP
165:                CP=CP+1162:                CP=CP+1
166: !              IF (CP==M) CP=0163: !              IF (CP==M) CP=0
167:                IF (CP==NEBMUPDATE) CP=0164:                IF (CP==NEBMUPDATE) CP=0
168:           ENDDO165:           ENDDO
169: 166:               
170:           STP(1:D) = 1.0D0167:           STP(1:D) = 1.0D0
171:      ENDIF MAIN168:      ENDIF MAIN
172: 169: 
173:      !  Store the new search direction170:      !  Store the new search direction
174:      IF (NITERDONE.GT.1) SEARCHSTEP(POINT,:)=GTMP171:      IF (NITERDONE.GT.1) SEARCHSTEP(POINT,:)=GTMP
175:       172:       
176: !173: !
177: ! If the number of images has changed since G was declared then G is not the same174: ! If the number of images has changed since G was declared then G is not the same
178: ! size as Gtmp and Dot_Product cannot be used.175: ! size as Gtmp and Dot_Product cannot be used.
179: !176: !
266:      ENDIF263:      ENDIF
267:      NDECREASE=0264:      NDECREASE=0
268: 20   X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)265: 20   X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)
269: 266: 
270:      IF (PREVGRAD.LT.DNEBSWITCH) THEN267:      IF (PREVGRAD.LT.DNEBSWITCH) THEN
271:         CALL OLDNEBGRADIENT268:         CALL OLDNEBGRADIENT
272:      ELSE269:      ELSE
273:         CALL NEBGRADIENT270:         CALL NEBGRADIENT
274:      ENDIF271:      ENDIF
275: 272: 
276:      NRED = 0 
277:      DO WHILE ((ETOTAL-ENOLD).GT.NEBMAXERISE)273:      DO WHILE ((ETOTAL-ENOLD).GT.NEBMAXERISE)
278:         X(1:D) = X(1:D) - STP(1:D)*SEARCHSTEP(POINT,1:D)274:         X(1:D) = X(1:D) - STP(1:D)*SEARCHSTEP(POINT,1:D)
279:         STP(1:D) = STP(1:D)/10275:         STP(1:D) = STP(1:D)/10
280:         X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)276:         X(1:D) = X(1:D) + STP(1:D)*SEARCHSTEP(POINT,1:D)
281:         NRED = NRED + 1 
282:         IF (PREVGRAD.LT.DNEBSWITCH) THEN277:         IF (PREVGRAD.LT.DNEBSWITCH) THEN
283:            CALL OLDNEBGRADIENT278:            CALL OLDNEBGRADIENT
284:         ELSE279:         ELSE
285:            CALL NEBGRADIENT280:            CALL NEBGRADIENT
286:         END IF281:         END IF
287:         IF (NRED .EQ. 10) THEN 
288:             EXIT 
289:         END IF 
290:      END DO282:      END DO
291: 283: 
292:      ENOLD=ETOTAL284:      ENOLD=ETOTAL
293: 285: 
294:      IF (ETOTAL/NIMAGE.LT.COLDFUSIONLIMIT) THEN286:      IF (ETOTAL/NIMAGE.LT.COLDFUSIONLIMIT) THEN
295:         WRITE(*,'(A,2G20.10)') ' lbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/NIMAGE,COLDFUSIONLIMIT287:         WRITE(*,'(A,2G20.10)') ' lbfgs> Cold fusion diagnosed - step discarded, energy, limit=',ETOTAL/NIMAGE,COLDFUSIONLIMIT
296:         IF (DEBUG) CALL DUMPFILES("e")288:         IF (DEBUG) CALL DUMPFILES("e")
297:         IF (KADJUSTFRQ.GT.0) PRINT '(A,G20.10)',' lbfgs> Final DNEB force constant ',NEWNEBK(1)289:         IF (KADJUSTFRQ.GT.0) PRINT '(A,G20.10)',' lbfgs> Final DNEB force constant ',NEWNEBK(1)
298:         RETURN290:         RETURN
299:      ENDIF291:      ENDIF


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0