hdiff output

r29765/Dijinit.f90 2016-07-06 15:39:16.770215208 +0100 r29764/Dijinit.f90 2016-07-06 15:39:17.158220453 +0100
 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 17: !   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18: ! 18: !
 19:  19: 
 20: !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 20: !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 21: ! 21: !
 22: !  Dijkstra connection algorithm for pathsample. 22: !  Dijkstra connection algorithm for pathsample.
 23: ! 23: !
 24: SUBROUTINE DIJINIT(NWORST) 24: SUBROUTINE DIJINIT(NWORST)
 25: USE PORFUNCS 25: USE PORFUNCS
 26: USE COMMONS 26: USE COMMONS
 27: USE UTILS,ONLY : GETUNIT 
 28: IMPLICIT NONE 27: IMPLICIT NONE
 29:  28: 
 30: INTEGER J1, J2, J4, PARENT(NMIN), JMINW, NPERM, J5, LJ1, LJ2, NWORST, NSTEPS, NMINSTART, NMINEND, J6, MUNIT 29: INTEGER J1, J2, J4, PARENT(NMIN), JMINW, NPERM, J5, LJ1, LJ2, NWORST, NSTEPS, NMINSTART, NMINEND, J6
 31: INTEGER, ALLOCATABLE :: LOCATIONSTART(:), LOCATIONEND(:) 30: INTEGER, ALLOCATABLE :: LOCATIONSTART(:), LOCATIONEND(:)
 32: LOGICAL PERMANENT(NMIN), ISA(NMIN), ISB(NMIN), ISSTART(NMIN), NOTDONE 31: LOGICAL PERMANENT(NMIN), ISA(NMIN), ISB(NMIN), ISSTART(NMIN), NOTDONE
 33: DOUBLE PRECISION MINWEIGHT, DUMMY, TNEW, ELAPSED, PFTOTALSTART, HUGESAVE, THRESH 32: DOUBLE PRECISION MINWEIGHT, DUMMY, TNEW, ELAPSED, PFTOTALSTART, HUGESAVE, THRESH
 34: DOUBLE PRECISION MAXWEIGHT, SCALEFAC, PDMAX, PD 33: DOUBLE PRECISION MAXWEIGHT, SCALEFAC, PDMAX, PD
 35: INTEGER, DIMENSION(:), ALLOCATABLE :: PATHMINS 34: INTEGER, DIMENSION(:), ALLOCATABLE :: PATHMINS
 36: ! 35: !
 37: ! KIND=16 is not supported by Portland. If you want extra precision, uncomment the following line 36: ! KIND=16 is not supported by Portland. If you want extra precision, uncomment the following line
 38: ! and use NAG. 37: ! and use NAG.
 39: ! 38: !
 40: ! REAL(KIND=16) :: TMPWEIGHT, WEIGHT(NMIN) 39: ! REAL(KIND=16) :: TMPWEIGHT, WEIGHT(NMIN)
 41: REAL(KIND=8) :: TMPWEIGHT, WEIGHT(NMIN) 40: REAL(KIND=8) :: TMPWEIGHT, WEIGHT(NMIN)
 42:  41: 
 43:  42: 
 44: CALL CPU_TIME(ELAPSED) 43: CALL CPU_TIME(ELAPSED)
 45: DEALLOCATE(DMIN1,DMIN2) 44: DEALLOCATE(DMIN1,DMIN2)
 46: ALLOCATE(DMIN1(10000),DMIN2(10000)) ! surely this will be enough room! 45: ALLOCATE(DMIN1(10000),DMIN2(10000)) ! surely this will be enough room!
 47: IF (DIJPRUNET) THEN 46: IF (DIJPRUNET) ALLOCATE(PATHMINS(10000))
 48:    MUNIT=GETUNIT() 
 49:    ALLOCATE(PATHMINS(10000)) !change if more needed, but 10000 should be enough for most uses 
 50:    PRINT '(A)',' dijprune> Maximum number of minima on path is 10000.' 
 51: ENDIF 
 52: ISA(1:NMIN)=.FALSE. 47: ISA(1:NMIN)=.FALSE.
 53: ISB(1:NMIN)=.FALSE. 48: ISB(1:NMIN)=.FALSE.
 54: DO J1=1,NMINA 49: DO J1=1,NMINA
 55:    ISA(LOCATIONA(J1))=.TRUE. 50:    ISA(LOCATIONA(J1))=.TRUE.
 56: ENDDO 51: ENDDO
 57: DO J1=1,NMINB 52: DO J1=1,NMINB
 58:    ISB(LOCATIONB(J1))=.TRUE. 53:    ISB(LOCATIONB(J1))=.TRUE.
 59: ENDDO 54: ENDDO
 60:  55: 
 61: !!!!!!!!!!!!!!!!!!!   Dijkstra calculation    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 56: !!!!!!!!!!!!!!!!!!!   Dijkstra calculation    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 79:    LOCATIONEND(1:NMINB)=LOCATIONB(1:NMINB) 74:    LOCATIONEND(1:NMINB)=LOCATIONB(1:NMINB)
 80:    ISSTART(1:NMIN)=ISA(1:NMIN) 75:    ISSTART(1:NMIN)=ISA(1:NMIN)
 81:    PFTOTALSTART=PFTOTALA 76:    PFTOTALSTART=PFTOTALA
 82: ENDIF 77: ENDIF
 83:  78: 
 84: PDMAX=-1.0D0 79: PDMAX=-1.0D0
 85: DO J2=1,NMIN 80: DO J2=1,NMIN
 86:    DO J5=1,PAIRDISTMAX 81:    DO J5=1,PAIRDISTMAX
 87:       IF (PAIRDIST(J2,J5).GT.PDMAX) THEN 82:       IF (PAIRDIST(J2,J5).GT.PDMAX) THEN
 88:          PDMAX=PAIRDIST(J2,J5) 83:          PDMAX=PAIRDIST(J2,J5)
 89: !        IF (DEBUG) PRINT '(A,G20.10)','Dijinit> maximum neighbour metric value increased to',PDMAX  84: !        IF (DEBUG) PRINT '(A,G20.10)','Dijinit> maximum neighbour metric value increased to',PDMAX
 90: !        IF (DEBUG) PRINT '(A,2I8,G20.10)','Dijinit> J2,J5,PAIRDIST=',J2,J5,PAIRDIST(J2,J5) 85: !        IF (DEBUG) PRINT '(A,2I8,G20.10)','Dijinit> J2,J5,PAIRDIST=',J2,J5,PAIRDIST(J2,J5)
 91:       ENDIF 86:       ENDIF
 92:    ENDDO 87:    ENDDO
 93: ENDDO 88: ENDDO
 94: PRINT '(A,G20.10)','Dijinit> maximum neighbour metric value=',PDMAX 89: PRINT '(A,G20.10)','Dijinit> maximum neighbour metric value=',PDMAX
 95: ! 90: !
 96: !  Find largest weight for each B(A) minimum to all A(B) minima. 91: !  Find largest weight for each B(A) minimum to all A(B) minima.
 97: ! 92: !
 98: !  Added maximum weight condition via a scale factor. 93: !  Added maximum weight condition via a scale factor.
 99: !  Otherwise loss of precision can cause connections to be missed completely. DJW 29/7/08 94: !  Otherwise loss of precision can cause connections to be missed completely. DJW 29/7/08
104:    SCALEFAC=1.0D0 99:    SCALEFAC=1.0D0
105: 222   LJ1=LOCATIONSTART(J1)100: 222   LJ1=LOCATIONSTART(J1)
106:    WEIGHT(1:NMIN)=HUGE(1.0D0)101:    WEIGHT(1:NMIN)=HUGE(1.0D0)
107:    HUGESAVE=WEIGHT(1)102:    HUGESAVE=WEIGHT(1)
108:    WEIGHT(LJ1)=0.0D0103:    WEIGHT(LJ1)=0.0D0
109:    PERMANENT(1:NMIN)=.FALSE.104:    PERMANENT(1:NMIN)=.FALSE.
110:    PERMANENT(LJ1)=.TRUE.105:    PERMANENT(LJ1)=.TRUE.
111:    NPERM=1106:    NPERM=1
112:    PARENT(1:NMIN)=0 ! parent is initially undefined107:    PARENT(1:NMIN)=0 ! parent is initially undefined
113:    J4=LJ1108:    J4=LJ1
 109: 
114:    dijkstraloop: DO110:    dijkstraloop: DO
 111: 
115:       DO J2=1,NMIN112:       DO J2=1,NMIN
116:          IF (J2.EQ.J4) CYCLE113:          IF (J2.EQ.J4) CYCLE
117:          IF (PERMANENT(J2)) CYCLE114:          IF (PERMANENT(J2)) CYCLE
118:          PD=1.0D4*PDMAX115:          PD=1.0D4*PDMAX
119:          IF (.NOT.DIJPRUNET) THEN !for pruning the database all minima count not just the ones not searched yet116:          IF (.NOT.DIJPRUNET) THEN !for pruning the database all minima count not just the ones not searched yet
120:             DO J5=1,NPAIRDONE ! skip117:             DO J5=1,NPAIRDONE ! skip
121:                IF ((PAIR1(J5).EQ.J4).AND.(PAIR2(J5).EQ.J2)) GOTO 973118:                IF ((PAIR1(J5).EQ.J4).AND.(PAIR2(J5).EQ.J2)) GOTO 973
122:                IF ((PAIR1(J5).EQ.J2).AND.(PAIR2(J5).EQ.J4)) GOTO 973119:                IF ((PAIR1(J5).EQ.J2).AND.(PAIR2(J5).EQ.J4)) GOTO 973
123:             ENDDO120:             ENDDO
124:          ENDIF 121:          ENDIF 
200:                   NOTDONE=.FALSE.197:                   NOTDONE=.FALSE.
201:                ENDIF198:                ENDIF
202:             ENDIF199:             ENDIF
203:          ENDDO200:          ENDDO
204:          STOP !!! DJW201:          STOP !!! DJW
205:       ENDIF202:       ENDIF
206: 203: 
207:       J4=JMINW204:       J4=JMINW
208:       PERMANENT(J4)=.TRUE.205:       PERMANENT(J4)=.TRUE.
209:       NPERM=NPERM+1206:       NPERM=NPERM+1
210: !     PRINT '(A,2I8,G20.10)','permanent minimum J4,NPERM,WEIGHT=',J4,NPERM,WEIGHT(J4) 207: !     PRINT '(A,2I8,G20.10)','permanent minimum J4,NPERM,WEIGHT=',J4,NPERM,WEIGHT(J4)
211:       IF (WEIGHT(J4).GT.MAXWEIGHT) THEN208:       IF (WEIGHT(J4).GT.MAXWEIGHT) THEN
212:          SCALEFAC=SCALEFAC/10.0D0209:          SCALEFAC=SCALEFAC/10.0D0
213:          PRINT '(A,G20.10)','dijinit> Maximum weight is too large - scaling by ',SCALEFAC210:          PRINT '(A,G20.10)','dijinit> Maximum weight is too large - scaling by ',SCALEFAC
214:          GOTO 222211:          GOTO 222
215:       ENDIF212:       ENDIF
216: 213: 
217:       IF (NPERM.EQ.NMIN) EXIT dijkstraloop214:       IF (NPERM.EQ.NMIN) EXIT dijkstraloop
218: 215: 
219:    ENDDO dijkstraloop216:    ENDDO dijkstraloop
220: 217: 
222: ! 219: ! 
223: !  Summarise the best path for any A(B) and any B(A)220: !  Summarise the best path for any A(B) and any B(A)
224: !221: !
225: LJ2=LOCATIONEND(1)222: LJ2=LOCATIONEND(1)
226: LJ1=LOCATIONSTART(1)223: LJ1=LOCATIONSTART(1)
227: J5=LJ2224: J5=LJ2
228: NWORST=0225: NWORST=0
229: NSTEPS=0226: NSTEPS=0
230: PRINT '(A)','Dijinit> Summary of best path based on missing connection metric - note distance scaling is removed'227: PRINT '(A)','Dijinit> Summary of best path based on missing connection metric - note distance scaling is removed'
231: PRINT '(A)','    min1          energy        min2          energy             metric          edge weight            weight'228: PRINT '(A)','    min1          energy        min2          energy             metric          edge weight            weight'
232: IF (DIJPRUNET) OPEN(MUNIT,FILE='min.retain',POSITION='APPEND',ACTION='WRITE',STATUS='NEW')229: IF (DIJPRUNET) OPEN(121,FILE='min.retain',POSITION='APPEND',ACTION='WRITE',STATUS='NEW')
233: DO 230: DO 
234:    IF (PARENT(J5).EQ.0) THEN231:    IF (PARENT(J5).EQ.0) THEN
235:       PRINT '(A,I6,A)','Dijinit> ERROR - parent for J5=',J5,' is zero'232:       PRINT '(A,I6,A)','Dijinit> ERROR - parent for J5=',J5,' is zero'
236:       PRINT '(A)',     'Dijinit> Suggests all possible pairs have been tried!'233:       PRINT '(A)',     'Dijinit> Suggests all possible pairs have been tried!'
237:       STOP234:       STOP
238:    ENDIF235:    ENDIF
239: !  DUMMY=PAIRDIST(MAX(J5,PARENT(J5))*(MAX(J5,PARENT(J5))-1)/2+MIN(J5,PARENT(J5)))*SCALEFAC236: !  DUMMY=PAIRDIST(MAX(J5,PARENT(J5))*(MAX(J5,PARENT(J5))-1)/2+MIN(J5,PARENT(J5)))*SCALEFAC
240:    DUMMY=1.0D4*PDMAX*SCALEFAC237:    DUMMY=1.0D4*PDMAX*SCALEFAC
241:    DO J2=1,NPAIRDONE ! skip238:    DO J2=1,NPAIRDONE ! skip
242:       IF ((PAIR1(J2).EQ.J5).AND.(PAIR2(J2).EQ.PARENT(J5))) GOTO 864239:       IF ((PAIR1(J2).EQ.J5).AND.(PAIR2(J2).EQ.PARENT(J5))) GOTO 864
267:             ELSE264:             ELSE
268:                IF ((PARENT(J5).LE.NMINA+NMINB).AND.(PARENT(J5).GT.NMINA)) TMPWEIGHT=NMIN+1-J5265:                IF ((PARENT(J5).LE.NMINA+NMINB).AND.(PARENT(J5).GT.NMINA)) TMPWEIGHT=NMIN+1-J5
269:                IF ((J5.LE.NMINA+NMINB).AND.(J5.GT.NMINA)) TMPWEIGHT=NMIN+1-PARENT(J5)266:                IF ((J5.LE.NMINA+NMINB).AND.(J5.GT.NMINA)) TMPWEIGHT=NMIN+1-PARENT(J5)
270:             ENDIF 267:             ENDIF 
271:           ENDIF268:           ENDIF
272:       ELSEIF (EXPCOSTFUNCTION) THEN ! saves memory and CPU when endpoint separation is very large SAT269:       ELSEIF (EXPCOSTFUNCTION) THEN ! saves memory and CPU when endpoint separation is very large SAT
273:          IF (DUMMY.EQ.0.0D0) THEN270:          IF (DUMMY.EQ.0.0D0) THEN
274:             TMPWEIGHT=0.0D0271:             TMPWEIGHT=0.0D0
275:          ELSEIF (DUMMY.GT.300.0D0) THEN272:          ELSEIF (DUMMY.GT.300.0D0) THEN
276:              TMPWEIGHT=EXP(300.0D0)273:              TMPWEIGHT=EXP(300.0D0)
277:          ELSE 274:          ELSE
278:             TMPWEIGHT=EXP(DUMMY)275:             TMPWEIGHT=EXP(DUMMY)
279:          ENDIF276:          ENDIF
280:       ELSE ! compare higher powers to favour more small jumps over big ones DJW277:       ELSE ! compare higher powers to favour more small jumps over big ones DJW
281:          IF (DUMMY.EQ.0.0D0) THEN278:          IF (DUMMY.EQ.0.0D0) THEN
282:             TMPWEIGHT=0.0D0279:             TMPWEIGHT=0.0D0
283:          ELSEIF (INTERPCOSTFUNCTION) THEN280:          ELSEIF (INTERPCOSTFUNCTION) THEN
284:             TMPWEIGHT=DUMMY**COSTFUNCTIONPOWER281:             TMPWEIGHT=DUMMY**COSTFUNCTIONPOWER
285:          ELSEIF (COSTFUNCTIONPOWER.EQ.0) THEN282:          ELSEIF (COSTFUNCTIONPOWER.EQ.0) THEN
286:             TMPWEIGHT=1.0D0283:             TMPWEIGHT=1.0D0
287:          ELSEIF (COSTFUNCTIONPOWER.EQ.-1) THEN284:          ELSEIF (COSTFUNCTIONPOWER.EQ.-1) THEN
330: !751      CONTINUE327: !751      CONTINUE
331: !      ENDIF328: !      ENDIF
332:    ENDIF329:    ENDIF
333:    J5=PARENT(J5)330:    J5=PARENT(J5)
334:    IF (J5.EQ.LJ1) EXIT331:    IF (J5.EQ.LJ1) EXIT
335:    IF (J5.EQ.0) EXIT332:    IF (J5.EQ.0) EXIT
336: ENDDO333: ENDDO
337: PRINT '(2(A,I8))','Dijinit> Number of steps=',NSTEPS,' number of missing connections=',NWORST334: PRINT '(2(A,I8))','Dijinit> Number of steps=',NSTEPS,' number of missing connections=',NWORST
338: IF (DIJPRUNET) THEN !write the best path out to min.retain and then terminate335: IF (DIJPRUNET) THEN !write the best path out to min.retain and then terminate
339:    PATHMINS(NSTEPS+1)=J5336:    PATHMINS(NSTEPS+1)=J5
340:    WRITE(MUNIT,'(I8)') NSTEPS+1337:    WRITE(121,'(I8)') NSTEPS+1
341:    DO J6=1,NSTEPS+1338:    DO J6=1,NSTEPS+1
342:       WRITE(MUNIT,'(I8)') PATHMINS(J6)339:       WRITE(121,'(I8)') PATHMINS(J6)
343:    ENDDO340:    ENDDO
344:    CLOSE(MUNIT)341:    CLOSE(121)
345:    PRINT '(A)','Dijprune> Best path written to min.retain'342:    PRINT '(A)','Dijprune> Best path written to min.retain'
346:    DEALLOCATE(PATHMINS) 
347:    STOP343:    STOP
348: ENDIF344: ENDIF
349: IF (NWORST.EQ.0) THEN345: IF (NWORST.EQ.0) THEN
350:    PRINT '(A)','Dijinit> Connected path found'346:    PRINT '(A)','Dijinit> Connected path found'
351: !  IF (DIJCONT) THEN347: !  IF (DIJCONT) THEN
352: !     DIJINITT=.FALSE.348: !     DIJINITT=.FALSE.
353: !  ELSE349: !  ELSE
354:       STOP350:       STOP
355: !  ENDIF351: !  ENDIF
356: ENDIF352: ENDIF


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0