hdiff output

r33314/bnbalign.f90 2017-09-15 20:30:17.849825210 +0100 r33313/bnbalign.f90 2017-09-15 20:30:19.173842914 +0100
  3: MODULE GOPERMDIST  3: MODULE GOPERMDIST
  4:   4: 
  5: ! SAVECOORDSA(3*NATOMS,NSTRUCTS) stores the centred candidate structures  5: ! SAVECOORDSA(3*NATOMS,NSTRUCTS) stores the centred candidate structures
  6: ! SAVECOORDSB(3*NATOMS) stores the centred target structure  6: ! SAVECOORDSB(3*NATOMS) stores the centred target structure
  7:   7: 
  8: ! PERMCOORDSB(3,NATOMS,NPERMGROUP) stores the structures for the k-d tree  8: ! PERMCOORDSB(3,NATOMS,NPERMGROUP) stores the structures for the k-d tree
  9:   9: 
 10:  10: 
 11: USE COMMONS, ONLY : PERMGROUP, NPERMSIZE, NPERMGROUP, BESTPERM, MYUNIT, & 11: USE COMMONS, ONLY : PERMGROUP, NPERMSIZE, NPERMGROUP, BESTPERM, MYUNIT, &
 12:  & NSETS, SETS, OHCELLT, PERMINVOPT, PERMDIST, PERMOPT, BOXLX, BOXLY, BOXLZ 12:  & NSETS, SETS, OHCELLT, PERMINVOPT, PERMDIST, PERMOPT, BOXLX, BOXLY, BOXLZ
 13: USE PREC, ONLY: INT64 
 14: USE PRIORITYQUEUE, ONLY: QUEUE 13: USE PRIORITYQUEUE, ONLY: QUEUE
 15:  14: 
 16: IMPLICIT NONE 15: IMPLICIT NONE
 17:  16: 
 18: INTEGER, SAVE :: NATOMS, NCALC, NLAP, NQUENCH, NBAD 17: INTEGER, SAVE :: NATOMS, NCALC, NLAP, NQUENCH, NBAD
 19: INTEGER, SAVE :: PMAXNEI = 60 ! Number of nearest neighbours to store 18: INTEGER, SAVE :: PMAXNEI = 60 ! Number of nearest neighbours to store
 20: DOUBLE PRECISION, PARAMETER :: PSCALE = 1.D6 ! Scale for linear assignment problem 19: DOUBLE PRECISION, PARAMETER :: PSCALE = 1.D6 ! Scale for linear assignment problem
 21: DOUBLE PRECISION, PARAMETER :: PI = 3.141592653589793D0 20: DOUBLE PRECISION, PARAMETER :: PI = 3.141592653589793D0
 22: ! Absolute Tolerance, Relative Tolerance, Relative Tolerance for MINPERMDIST quench 21: ! Absolute Tolerance, Relative Tolerance, Relative Tolerance for MINPERMDIST quench
 23: DOUBLE PRECISION, SAVE :: ATOL=1D-8, RTOL=1D-1, MPRTOL=1.D-1 22: DOUBLE PRECISION, SAVE :: ATOL=1D-8, RTOL=1D-1, MPRTOL=1.D-1
502:                                                 & DUMMYLDISTS(:M*NPERM,J1), &501:                                                 & DUMMYLDISTS(:M*NPERM,J1), &
503:                  & MATMUL(FVECS(:,I),DUMMYDOTDISP(:,:M*NPERM,J1)).GT.0.D0)502:                  & MATMUL(FVECS(:,I),DUMMYDOTDISP(:,:M*NPERM,J1)).GT.0.D0)
504:             END DO503:             END DO
505: 504: 
506:             CALL PERMNEARESTNEIGHBOURDISTS(DUMMYLDISTS2,DUMMYIDX,NATOMS, &505:             CALL PERMNEARESTNEIGHBOURDISTS(DUMMYLDISTS2,DUMMYIDX,NATOMS, &
507:              & PMAXNEI,DUMMYNEARIDX,DUMMYNEARLDISTS,NPERMGROUP)506:              & PMAXNEI,DUMMYNEARIDX,DUMMYNEARLDISTS,NPERMGROUP)
508: 507: 
509:             D = SUM(DUMMYNEARLDISTS)508:             D = SUM(DUMMYNEARLDISTS)
510:             ESTLOWER = MIN(D, ESTLOWER)509:             ESTLOWER = MIN(D, ESTLOWER)
511: 510: 
512:             IF(DEBUG) WRITE(MYUNIT, "(A,I16,A,G20.5)") &511:             IF(DEBUG) WRITE(MYUNIT, "(A,I16,A,G10.5)") &
513:      & "gopermdist> estimating for face         = ", I, &512:      & "gopermdist> estimating for face         = ", I, &
514:      & "         lower bound = ", D**0.5513:      & "         lower bound = ", D**0.5
515:         END DO514:         END DO
516:         ESTLOWER = SQRT(ESTLOWER)515:         ESTLOWER = SQRT(ESTLOWER)
517: 516: 
518:     ELSE517:     ELSE
519:         CALL PERMNEARESTNEIGHBOURDISTS(DUMMYLDISTS,DUMMYIDX,NATOMS,PMAXNEI, &518:         CALL PERMNEARESTNEIGHBOURDISTS(DUMMYLDISTS,DUMMYIDX,NATOMS,PMAXNEI, &
520:          & DUMMYNEARIDX,DUMMYNEARLDISTS,NPERMGROUP)519:          & DUMMYNEARIDX,DUMMYNEARLDISTS,NPERMGROUP)
521: 520: 
522:         ESTLOWER = SUM(DUMMYNEARLDISTS)**0.5521:         ESTLOWER = SUM(DUMMYNEARLDISTS)**0.5
755: ! DISTANCE RETURN INACCURATE754: ! DISTANCE RETURN INACCURATE
756: IMPLICIT NONE755: IMPLICIT NONE
757: 756: 
758: INTEGER, INTENT(IN) :: NATOMS,NPERMGROUP,MAXNEI,NIDX(MAXNEI*NATOMS,NPERMGROUP)757: INTEGER, INTENT(IN) :: NATOMS,NPERMGROUP,MAXNEI,NIDX(MAXNEI*NATOMS,NPERMGROUP)
759: DOUBLE PRECISION, INTENT(IN) :: NDISTS(MAXNEI*NATOMS,NPERMGROUP)758: DOUBLE PRECISION, INTENT(IN) :: NDISTS(MAXNEI*NATOMS,NPERMGROUP)
760: 759: 
761: DOUBLE PRECISION, INTENT(OUT) :: DIST760: DOUBLE PRECISION, INTENT(OUT) :: DIST
762: INTEGER, INTENT(OUT) :: PERM(NATOMS), INFO761: INTEGER, INTENT(OUT) :: PERM(NATOMS), INFO
763: 762: 
764: ! COULD SET THESE AS MODULE VARIABLES763: ! COULD SET THESE AS MODULE VARIABLES
765: INTEGER(KIND=INT64) :: KK(NATOMS*MAXNEI), CC(NATOMS*MAXNEI)764: INTEGER*8 :: KK(NATOMS*MAXNEI), CC(NATOMS*MAXNEI)
766: INTEGER(KIND=INT64) :: FIRST(NATOMS+1), X(NATOMS), Y(NATOMS)765: INTEGER*8 :: FIRST(NATOMS+1), X(NATOMS), Y(NATOMS)
767: INTEGER(KIND=INT64) :: U(NATOMS), V(NATOMS), N8, SZ8, H766: INTEGER*8 :: U(NATOMS), V(NATOMS), N8, SZ8, H
768: INTEGER N,M,I,J,K,K1,I1,J1,NDUMMY767: INTEGER N,M,I,J,K,K1,I1,J1,NDUMMY
769: 768: 
770: DIST = 0.D0769: DIST = 0.D0
771: INFO=0770: INFO=0
772: 771: 
773: NDUMMY=0772: NDUMMY=0
774: 773: 
775: DO J1=1,NPERMGROUP774: DO J1=1,NPERMGROUP
776: 775: 
777:     N = NPERMSIZE(J1)776:     N = NPERMSIZE(J1)
1172:       double precision DUMMY1171:       double precision DUMMY
1173: 1172: 
1174: !     Parameters1173: !     Parameters
1175: !       scale : Precision1174: !       scale : Precision
1176: !       maxnei: Maximum number of closest neighbours1175: !       maxnei: Maximum number of closest neighbours
1177:       double precision scale, d, h1176:       double precision scale, d, h
1178: 1177: 
1179:       parameter (scale = 1.0d6   )1178:       parameter (scale = 1.0d6   )
1180: !      parameter (maxnei = 60     )1179: !      parameter (maxnei = 60     )
1181: 1180: 
1182:       integer(kind=INT64) first(n+1)!, x(n), y(n)1181:       integer*8 first(n+1)!, x(n), y(n)
1183: !      integer(kind=INT64) u(n), v(n)1182: !      integer*8 u(n), v(n)
1184:       integer   m, i, j, k, l, l2, t, a1183:       integer   m, i, j, k, l, l2, t, a
1185:       integer(kind=INT64) n8, sz81184:       integer*8 n8, sz8
1186:       integer J11185:       integer J1
1187: 1186: 
1188: !     Distance function1187: !     Distance function
1189:       double precision permdist1188:       double precision permdist
1190: 1189: 
1191:       s(1)=sx1190:       s(1)=sx
1192:       s(2)=sy1191:       s(2)=sy
1193:       s(3)=sz1192:       s(3)=sz
1194:       m = maxnei1193:       m = maxnei
1195:       if(n .le. maxnei) m = n1194:       if(n .le. maxnei) m = n


r33314/DSOFT.f90 2017-09-15 20:30:17.629822268 +0100 r33313/DSOFT.f90 2017-09-15 20:30:18.953839976 +0100
 23: !    Kostelec, P. J., & Rockmore, D. N. (2008). FFTs on the rotation group.  23: !    Kostelec, P. J., & Rockmore, D. N. (2008). FFTs on the rotation group. 
 24: !    Journal of Fourier Analysis and Applications, 14(2), 145–179.  24: !    Journal of Fourier Analysis and Applications, 14(2), 145–179. 
 25: !    http://doi.org/10.1007/s00041-008-9013-5 25: !    http://doi.org/10.1007/s00041-008-9013-5
 26:  26: 
 27: !    Dependencies: 27: !    Dependencies:
 28: !        1. FFTW 28: !        1. FFTW
 29:  29: 
 30: MODULE DSOFT 30: MODULE DSOFT
 31:  31: 
 32: USE FFTW3 32: USE FFTW3
 33: USE PREC, ONLY: INT64, REAL64 
 34:  33: 
 35: IMPLICIT NONE 34: IMPLICIT NONE
 36:  35: 
 37: INTEGER(KIND=INT64), SAVE :: BW 36: INTEGER*8, SAVE :: BW
 38: DOUBLE PRECISION, PARAMETER :: PI = 3.141592653589793D0 37: DOUBLE PRECISION, PARAMETER :: PI = 3.141592653589793D0
 39: DOUBLE PRECISION, SAVE, ALLOCATABLE :: WEIGHTS(:), WIGNERD(:,:,:,:) 38: DOUBLE PRECISION, SAVE, ALLOCATABLE :: WEIGHTS(:), WIGNERD(:,:,:,:)
 40:  39: 
 41: CONTAINS 40: CONTAINS
 42:  41: 
 43: SUBROUTINE SETBANDWIDTH(BANDWIDTH) 42: SUBROUTINE SETBANDWIDTH(BANDWIDTH)
 44:  43: 
 45: IMPLICIT NONE 44: IMPLICIT NONE
 46: INTEGER(KIND=INT64), INTENT(IN) :: BANDWIDTH 45: INTEGER*8, INTENT(IN) :: BANDWIDTH
 47:  46: 
 48: ! Check if bandwidth has already been calculated 47: ! Check if bandwidth has already been calculated
 49: IF (BW.NE.BANDWIDTH) THEN 48: IF (BW.NE.BANDWIDTH) THEN
 50:     IF (ALLOCATED(WEIGHTS)) THEN 49:     IF (ALLOCATED(WEIGHTS)) THEN
 51:         DEALLOCATE(WEIGHTS) 50:         DEALLOCATE(WEIGHTS)
 52:     ENDIF 51:     ENDIF
 53:     IF (ALLOCATED(WIGNERD)) THEN 52:     IF (ALLOCATED(WIGNERD)) THEN
 54:         DEALLOCATE(WIGNERD) 53:         DEALLOCATE(WIGNERD)
 55:     ENDIF 54:     ENDIF
 56:     ALLOCATE(WEIGHTS(2*BANDWIDTH)) 55:     ALLOCATE(WEIGHTS(2*BANDWIDTH))
 59:     CALL CALCWIGNERD(BANDWIDTH) 58:     CALL CALCWIGNERD(BANDWIDTH)
 60: ENDIF 59: ENDIF
 61:  60: 
 62: BW = BANDWIDTH 61: BW = BANDWIDTH
 63:  62: 
 64: END SUBROUTINE SETBANDWIDTH 63: END SUBROUTINE SETBANDWIDTH
 65:  64: 
 66: SUBROUTINE MAKEWEIGHTS(BANDWIDTH) 65: SUBROUTINE MAKEWEIGHTS(BANDWIDTH)
 67:  66: 
 68: IMPLICIT NONE 67: IMPLICIT NONE
 69: INTEGER(KIND=INT64), INTENT(IN) :: BANDWIDTH 68: INTEGER*8, INTENT(IN) :: BANDWIDTH
 70:  69: 
 71: DOUBLE PRECISION FUDGE, SINJ 70: DOUBLE PRECISION FUDGE, SINJ
 72: INTEGER(KIND=INT64) J, K 71: INTEGER*8 J, K
 73:  72: 
 74: FUDGE = PI / 4 / BANDWIDTH 73: FUDGE = PI / 4 / BANDWIDTH
 75:  74: 
 76: DO J=1, BANDWIDTH*2 75: DO J=1, BANDWIDTH*2
 77:     WEIGHTS(J) = 0 76:     WEIGHTS(J) = 0
 78:     SINJ = 2.D0 * SIN((2*J-1)*FUDGE) / BANDWIDTH   77:     SINJ = 2.D0 * SIN((2*J-1)*FUDGE) / BANDWIDTH  
 79:     DO K=1,BANDWIDTH   78:     DO K=1,BANDWIDTH  
 80:     WEIGHTS(J) = WEIGHTS(J) + SINJ * SIN((2*J-1)*(2*K-1)*FUDGE) / (2*K - 1)  79:     WEIGHTS(J) = WEIGHTS(J) + SINJ * SIN((2*J-1)*(2*K-1)*FUDGE) / (2*K - 1) 
 81:     ENDDO 80:     ENDDO
 82: ENDDO 81: ENDDO
 88: ! The Wigner little d elements are calculated with a recurrence relation 87: ! The Wigner little d elements are calculated with a recurrence relation
 89: ! This subroutine calculates the appropriate coefficients of the recurrent  88: ! This subroutine calculates the appropriate coefficients of the recurrent 
 90: ! relation. For more information see: 89: ! relation. For more information see:
 91: ! 90: !
 92: !    Kostelec, P. J., & Rockmore, D. N. (2008). FFTs on the rotation group.  91: !    Kostelec, P. J., & Rockmore, D. N. (2008). FFTs on the rotation group. 
 93: !    Journal of Fourier Analysis and Applications, 14(2), 145–179.  92: !    Journal of Fourier Analysis and Applications, 14(2), 145–179. 
 94: !    http://doi.org/10.1007/s00041-008-9013-5 93: !    http://doi.org/10.1007/s00041-008-9013-5
 95:  94: 
 96: IMPLICIT NONE 95: IMPLICIT NONE
 97:  96: 
 98: INTEGER(KIND=INT64), INTENT(IN) :: J,M1,M2 97: INTEGER*8, INTENT(IN) :: J,M1,M2
 99: DOUBLE PRECISION, INTENT(OUT) :: A, B, C 98: DOUBLE PRECISION, INTENT(OUT) :: A, B, C
100:  99: 
101: DOUBLE PRECISION T1,T2,T3,T4,T5,DJ,DM1,DM2100: DOUBLE PRECISION T1,T2,T3,T4,T5,DJ,DM1,DM2
102: 101: 
103: DJ = REAL(J,8)102: DJ = REAL(J,8)
104: DM1 = REAL(M1,8)103: DM1 = REAL(M1,8)
105: DM2 = REAL(M2,8)104: DM2 = REAL(M2,8)
106: 105: 
107: T1 = ((2.D0*DJ +3.D0)/(2.D0*DJ + 1.D0))**0.5D0106: T1 = ((2.D0*DJ +3.D0)/(2.D0*DJ + 1.D0))**0.5D0
108: T3 = (DJ+1.D0)*(2.D0*DJ+1.D0)107: T3 = (DJ+1.D0)*(2.D0*DJ+1.D0)
131: ! stores result in WIGNERD(k, l, m1, m2) in the SOFT module.130: ! stores result in WIGNERD(k, l, m1, m2) in the SOFT module.
132: !131: !
133: ! Follows method described in:132: ! Follows method described in:
134: !133: !
135: !    Kostelec, P. J., & Rockmore, D. N. (2008). FFTs on the rotation group. 134: !    Kostelec, P. J., & Rockmore, D. N. (2008). FFTs on the rotation group. 
136: !    Journal of Fourier Analysis and Applications, 14(2), 145–179. 135: !    Journal of Fourier Analysis and Applications, 14(2), 145–179. 
137: !    http://doi.org/10.1007/s00041-008-9013-5136: !    http://doi.org/10.1007/s00041-008-9013-5
138: 137: 
139: IMPLICIT NONE138: IMPLICIT NONE
140: 139: 
141: INTEGER(KIND=INT64), INTENT(IN) :: BANDWIDTH140: INTEGER*8, INTENT(IN) :: BANDWIDTH
142: 141: 
143: DOUBLE PRECISION COSB(2*BANDWIDTH+1), COSB2(2*BANDWIDTH+1), SINB2(2*BANDWIDTH+1)142: DOUBLE PRECISION COSB(2*BANDWIDTH+1), COSB2(2*BANDWIDTH+1), SINB2(2*BANDWIDTH+1)
144: DOUBLE PRECISION SINCOSB2(2*BANDWIDTH+1), SINDIVCOSB2(2*BANDWIDTH+1)143: DOUBLE PRECISION SINCOSB2(2*BANDWIDTH+1), SINDIVCOSB2(2*BANDWIDTH+1)
145: DOUBLE PRECISION, DIMENSION(0:3*BANDWIDTH-1) :: FACTORIALS144: DOUBLE PRECISION, DIMENSION(0:3*BANDWIDTH-1) :: FACTORIALS
146: DOUBLE PRECISION FACTOR, FUDGE, BETA, A, B, C, JM1(2*BANDWIDTH+1), T1,T2,T3,T4145: DOUBLE PRECISION FACTOR, FUDGE, BETA, A, B, C, JM1(2*BANDWIDTH+1), T1,T2,T3,T4
147: INTEGER(KIND=INT64) I,J,M1,M2,IND1,IND2,MAXM146: INTEGER*8 I,J,M1,M2,IND1,IND2,MAXM
148: 147: 
149: FUDGE = PI / 4 / BANDWIDTH148: FUDGE = PI / 4 / BANDWIDTH
150: DO I=1,2*BANDWIDTH149: DO I=1,2*BANDWIDTH
151:     BETA = FUDGE * (2*I-1)150:     BETA = FUDGE * (2*I-1)
152:     COSB(I) = COS(BETA)151:     COSB(I) = COS(BETA)
153:     COSB2(I) = COS(BETA/2)152:     COSB2(I) = COS(BETA/2)
154:     SINB2(I) = SIN(BETA/2)153:     SINB2(I) = SIN(BETA/2)
155:     SINCOSB2 = SINB2*COSB2154:     SINCOSB2 = SINB2*COSB2
156:     SINDIVCOSB2 = SINB2/COSB2155:     SINDIVCOSB2 = SINB2/COSB2
157: ENDDO 156: ENDDO 
200: END SUBROUTINE CALCWIGNERD199: END SUBROUTINE CALCWIGNERD
201: 200: 
202: 201: 
203: SUBROUTINE SOFT(INPUT, OUTPUT, BANDWIDTH)202: SUBROUTINE SOFT(INPUT, OUTPUT, BANDWIDTH)
204: 203: 
205: ! Performs discrete SO3 Fourier Analysis for a real input array for a function204: ! Performs discrete SO3 Fourier Analysis for a real input array for a function
206: ! defined on SO(3) returns a complex array of the Fourier Coefficients.205: ! defined on SO(3) returns a complex array of the Fourier Coefficients.
207: 206: 
208: IMPLICIT NONE207: IMPLICIT NONE
209: 208: 
210: INTEGER(KIND=INT64), INTENT(IN) :: BANDWIDTH209: INTEGER*8, INTENT(IN) :: BANDWIDTH
211: DOUBLE PRECISION, INTENT(IN) :: INPUT(2*BANDWIDTH,2*BANDWIDTH,2*BANDWIDTH)210: DOUBLE PRECISION, INTENT(IN) :: INPUT(2*BANDWIDTH,2*BANDWIDTH,2*BANDWIDTH)
212: COMPLEX(KIND=REAL64), INTENT(OUT) :: OUTPUT(BANDWIDTH, 2*BANDWIDTH-1, 2*BANDWIDTH-1)211: COMPLEX*16, INTENT(OUT) :: OUTPUT(BANDWIDTH, 2*BANDWIDTH-1, 2*BANDWIDTH-1)
213: 212: 
214: !INCLUDE "fftw3.f90"213: !INCLUDE "fftw3.f90"
215: COMPLEX(KIND=REAL64) IN1D(2*BANDWIDTH), OUT1D(2*BANDWIDTH), TEMP(2*BANDWIDTH, 2*BANDWIDTH, 2*BANDWIDTH)214: COMPLEX*16 IN1D(2*BANDWIDTH), OUT1D(2*BANDWIDTH), TEMP(2*BANDWIDTH, 2*BANDWIDTH, 2*BANDWIDTH)
216: INTEGER(KIND=INT64) PLAN, K1,K2,K3,M1,M2,I1,I2,IND1,IND2,J,MAXM215: INTEGER*8 PLAN, K1,K2,K3,M1,M2,I1,I2,IND1,IND2,J,MAXM
217: 216: 
218: 217: 
219: CALL SETBANDWIDTH(BANDWIDTH)218: CALL SETBANDWIDTH(BANDWIDTH)
220: 219: 
221: CALL DFFTW_PLAN_DFT_1D(PLAN, (2*BANDWIDTH), IN1D, OUT1D, FFTW_FORWARD, FFTW_ESTIMATE)220: CALL DFFTW_PLAN_DFT_1D(PLAN, (2*BANDWIDTH), IN1D, OUT1D, FFTW_FORWARD, FFTW_ESTIMATE)
222: 221: 
223: ! Do FFT on axis 1222: ! Do FFT on axis 1
224: DO K1=1,2*BANDWIDTH223: DO K1=1,2*BANDWIDTH
225:     DO K2=1,2*BANDWIDTH224:     DO K2=1,2*BANDWIDTH
226:         DO K3=1,2*BANDWIDTH225:         DO K3=1,2*BANDWIDTH
227:             IN1D(K3) = CMPLX(INPUT(K3,K2,K1), 0.D0, REAL64)226:             IN1D(K3) = CMPLX(INPUT(K3,K2,K1),0.D0, SELECTED_REAL_KIND(15,307))
228:         ENDDO227:         ENDDO
229:         CALL DFFTW_EXECUTE_(PLAN, IN1D, OUT1D)228:         CALL DFFTW_EXECUTE_(PLAN, IN1D, OUT1D)
230:         DO K3=1,2*BANDWIDTH229:         DO K3=1,2*BANDWIDTH
231:             TEMP(K3,K2,K1) = OUT1D(K3)230:             TEMP(K3,K2,K1) = OUT1D(K3)
232:         ENDDO231:         ENDDO
233:     ENDDO232:     ENDDO
234: ENDDO233: ENDDO
235: 234: 
236: ! Do FFT on axis 3235: ! Do FFT on axis 3
237: DO K1=1,2*BANDWIDTH236: DO K1=1,2*BANDWIDTH
240:             IN1D(K3) = TEMP(K2,K1,K3)239:             IN1D(K3) = TEMP(K2,K1,K3)
241:         ENDDO240:         ENDDO
242:         CALL DFFTW_EXECUTE_(PLAN, IN1D, OUT1D)241:         CALL DFFTW_EXECUTE_(PLAN, IN1D, OUT1D)
243:         DO K3=1,2*BANDWIDTH242:         DO K3=1,2*BANDWIDTH
244:             TEMP(K2,K1,K3) = OUT1D(K3)/(2*BANDWIDTH)**2243:             TEMP(K2,K1,K3) = OUT1D(K3)/(2*BANDWIDTH)**2
245:         ENDDO244:         ENDDO
246:     ENDDO245:     ENDDO
247: ENDDO246: ENDDO
248: 247: 
249: ! Perform Discrete Wigner Transform248: ! Perform Discrete Wigner Transform
250: OUTPUT = CMPLX(0.D0, 0.D0, REAL64)249: OUTPUT = CMPLX(0.D0,0.D0,8)
251: DO M2=-BANDWIDTH-1,BANDWIDTH-1250: DO M2=-BANDWIDTH-1,BANDWIDTH-1
252:     I2 = MODULO(M2, 2*BANDWIDTH) + 1251:     I2 = MODULO(M2, 2*BANDWIDTH) + 1
253:     IND2 = MODULO(M2, 2*BANDWIDTH-1) + 1252:     IND2 = MODULO(M2, 2*BANDWIDTH-1) + 1
254:     DO M1=-BANDWIDTH-1,BANDWIDTH-1253:     DO M1=-BANDWIDTH-1,BANDWIDTH-1
255:         I1 = MODULO(M1, 2*BANDWIDTH) + 1254:         I1 = MODULO(M1, 2*BANDWIDTH) + 1
256:         IND1 = MODULO(M1, 2*BANDWIDTH-1) + 1255:         IND1 = MODULO(M1, 2*BANDWIDTH-1) + 1
257:         MAXM = MAX(ABS(M1),ABS(M2))256:         MAXM = MAX(ABS(M1),ABS(M2))
258:         DO J=MAXM, BANDWIDTH-1257:         DO J=MAXM, BANDWIDTH-1
259:             DO K1=1,2*BANDWIDTH258:             DO K1=1,2*BANDWIDTH
260:                 OUTPUT(J+1,IND1,IND2) = OUTPUT(J+1,IND1,IND2) + WIGNERD(K1,J+1,IND1,IND2)*WEIGHTS(K1)*TEMP(I1,K1,I2)259:                 OUTPUT(J+1,IND1,IND2) = OUTPUT(J+1,IND1,IND2) + WIGNERD(K1,J+1,IND1,IND2)*WEIGHTS(K1)*TEMP(I1,K1,I2)
267: 266: 
268: END SUBROUTINE SOFT267: END SUBROUTINE SOFT
269: 268: 
270: SUBROUTINE ISOFT(INPUT, OUTPUT, BANDWIDTH)269: SUBROUTINE ISOFT(INPUT, OUTPUT, BANDWIDTH)
271: 270: 
272: ! Performs SO3 Fourier Synthesis for a complex input array of Fourier Coefficients271: ! Performs SO3 Fourier Synthesis for a complex input array of Fourier Coefficients
273: ! Generates a complex output array.272: ! Generates a complex output array.
274: 273: 
275: IMPLICIT NONE274: IMPLICIT NONE
276: 275: 
277: INTEGER(KIND=INT64), INTENT(IN) :: BANDWIDTH276: INTEGER*8, INTENT(IN) :: BANDWIDTH
278: COMPLEX(KIND=REAL64), INTENT(IN) :: INPUT(BANDWIDTH, 2*BANDWIDTH-1, 2*BANDWIDTH-1)277: COMPLEX*16, INTENT(IN) :: INPUT(BANDWIDTH, 2*BANDWIDTH-1, 2*BANDWIDTH-1)
279: COMPLEX(KIND=REAL64), INTENT(OUT) :: OUTPUT(2*BANDWIDTH,2*BANDWIDTH,2*BANDWIDTH)278: COMPLEX*16, INTENT(OUT) :: OUTPUT(2*BANDWIDTH,2*BANDWIDTH,2*BANDWIDTH)
280: 279: 
281: !INCLUDE "fftw3.f90"280: !INCLUDE "fftw3.f90"
282: COMPLEX(KIND=REAL64) IN1D(2*BANDWIDTH), OUT1D(2*BANDWIDTH), TEMP(2*BANDWIDTH, 2*BANDWIDTH, 2*BANDWIDTH)281: COMPLEX*16 IN1D(2*BANDWIDTH), OUT1D(2*BANDWIDTH), TEMP(2*BANDWIDTH, 2*BANDWIDTH, 2*BANDWIDTH)
283: INTEGER(KIND=INT64) PLAN, K1,K2,K3,M1,M2,I1,I2,IND1,IND2,J,MAXM282: INTEGER*8 PLAN, K1,K2,K3,M1,M2,I1,I2,IND1,IND2,J,MAXM
284: 283: 
285: CALL SETBANDWIDTH(BANDWIDTH)284: CALL SETBANDWIDTH(BANDWIDTH)
286: 285: 
287: CALL DFFTW_PLAN_DFT_1D(PLAN, (2*BANDWIDTH), IN1D, OUT1D, FFTW_BACKWARD, FFTW_ESTIMATE)286: CALL DFFTW_PLAN_DFT_1D(PLAN, (2*BANDWIDTH), IN1D, OUT1D, FFTW_BACKWARD, FFTW_ESTIMATE)
288: 287: 
289: ! Discrete inverse Wigner Transform288: ! Discrete inverse Wigner Transform
290: TEMP = CMPLX(0.D0, 0.D0, REAL64)289: TEMP = CMPLX(0.D0,0.D0,8)
291: DO M2=-BANDWIDTH-1,BANDWIDTH-1290: DO M2=-BANDWIDTH-1,BANDWIDTH-1
292:     I2 = MODULO(M2, 2*BANDWIDTH) + 1291:     I2 = MODULO(M2, 2*BANDWIDTH) + 1
293:     IND2 = MODULO(M2, 2*BANDWIDTH-1) + 1292:     IND2 = MODULO(M2, 2*BANDWIDTH-1) + 1
294:     DO M1=-BANDWIDTH-1,BANDWIDTH-1293:     DO M1=-BANDWIDTH-1,BANDWIDTH-1
295:         I1 = MODULO(M1, 2*BANDWIDTH) + 1294:         I1 = MODULO(M1, 2*BANDWIDTH) + 1
296:         IND1 = MODULO(M1, 2*BANDWIDTH-1) + 1295:         IND1 = MODULO(M1, 2*BANDWIDTH-1) + 1
297:         MAXM = MAX(ABS(M1),ABS(M2))296:         MAXM = MAX(ABS(M1),ABS(M2))
298:         DO K1=1,2*BANDWIDTH297:         DO K1=1,2*BANDWIDTH
299:             DO J=MAXM, BANDWIDTH-1298:             DO J=MAXM, BANDWIDTH-1
300:                 TEMP(I1,K1,I2) = TEMP(I1,K1,I2) + WIGNERD(K1,J+1,IND1,IND2)*INPUT(J+1,IND1,IND2)299:                 TEMP(I1,K1,I2) = TEMP(I1,K1,I2) + WIGNERD(K1,J+1,IND1,IND2)*INPUT(J+1,IND1,IND2)


r33314/fastbulk.f90 2017-09-15 20:30:18.073828204 +0100 r33313/fastbulk.f90 2017-09-15 20:30:19.401845963 +0100
 91: !    FASTOVERLAPUTILS (fastutils.f90) depends on (minperm.f90) 91: !    FASTOVERLAPUTILS (fastutils.f90) depends on (minperm.f90)
 92: !        Helper Module Needed for Peak Fitting and FFT routines 92: !        Helper Module Needed for Peak Fitting and FFT routines
 93:  93: 
 94: !*********************************************************************** 94: !***********************************************************************
 95:  95: 
 96: MODULE BULKFASTOVERLAP 96: MODULE BULKFASTOVERLAP
 97:  97: 
 98: USE COMMONS, ONLY : PERMGROUP, NPERMSIZE, NPERMGROUP, NATOMS, MYUNIT, NSETS, SETS, & 98: USE COMMONS, ONLY : PERMGROUP, NPERMSIZE, NPERMGROUP, NATOMS, MYUNIT, NSETS, SETS, &
 99:  & BOXLX, BOXLY, BOXLZ 99:  & BOXLX, BOXLY, BOXLZ
100: USE FASTOVERLAPUTILS, ONLY : DUMMYA, DUMMYB, XBESTA, XBESTASAVE100: USE FASTOVERLAPUTILS, ONLY : DUMMYA, DUMMYB, XBESTA, XBESTASAVE
101: USE PREC, ONLY: REAL64 
102: 101: 
103: IMPLICIT NONE102: IMPLICIT NONE
104: 103: 
105: ! If this is set to a value other than zero, algorithm will use this value104: ! If this is set to a value other than zero, algorithm will use this value
106: ! else it will set KWIDTH = 1/3 average interatomic separation.105: ! else it will set KWIDTH = 1/3 average interatomic separation.
107: DOUBLE PRECISION, SAVE :: KWIDTH=0.D0106: DOUBLE PRECISION, SAVE :: KWIDTH=0.D0
108: LOGICAL, SAVE :: OHCELLTSAVE107: LOGICAL, SAVE :: OHCELLTSAVE
109: DOUBLE PRECISION, SAVE :: OHOPSMAT(3,3,48)108: DOUBLE PRECISION, SAVE :: OHOPSMAT(3,3,48)
110: 109: 
111: DATA OHOPSMAT / &110: DATA OHOPSMAT / &
340: SUBROUTINE ALIGNGROUP(COORDS1LIST,N1LIST,COORDS2LIST,N2LIST,NATOMS,DEBUG, &339: SUBROUTINE ALIGNGROUP(COORDS1LIST,N1LIST,COORDS2LIST,N2LIST,NATOMS,DEBUG, &
341:     & BOXLX,BOXLY,BOXLZ,KWIDTH,NDISPLACEMENTS,NWAVE,NFSPACE,DISTMAT,ALIGNEDCOORDS2,SYM)340:     & BOXLX,BOXLY,BOXLZ,KWIDTH,NDISPLACEMENTS,NWAVE,NFSPACE,DISTMAT,ALIGNEDCOORDS2,SYM)
342: 341: 
343: IMPLICIT NONE342: IMPLICIT NONE
344: INTEGER, INTENT(IN) :: N1LIST, N2LIST, NATOMS, NDISPLACEMENTS, NFSPACE, NWAVE343: INTEGER, INTENT(IN) :: N1LIST, N2LIST, NATOMS, NDISPLACEMENTS, NFSPACE, NWAVE
345: LOGICAL, INTENT(IN) :: DEBUG,SYM344: LOGICAL, INTENT(IN) :: DEBUG,SYM
346: DOUBLE PRECISION, INTENT(IN) :: BOXLX, BOXLY, BOXLZ, KWIDTH345: DOUBLE PRECISION, INTENT(IN) :: BOXLX, BOXLY, BOXLZ, KWIDTH
347: DOUBLE PRECISION, INTENT(INOUT) :: COORDS1LIST(3*NATOMS,N1LIST), COORDS2LIST(3*NATOMS,N2LIST)346: DOUBLE PRECISION, INTENT(INOUT) :: COORDS1LIST(3*NATOMS,N1LIST), COORDS2LIST(3*NATOMS,N2LIST)
348: DOUBLE PRECISION, INTENT(OUT) :: DISTMAT(N1LIST,N2LIST), ALIGNEDCOORDS2(3*NATOMS,N1LIST,N2LIST)347: DOUBLE PRECISION, INTENT(OUT) :: DISTMAT(N1LIST,N2LIST), ALIGNEDCOORDS2(3*NATOMS,N1LIST,N2LIST)
349: 348: 
350: COMPLEX(KIND=REAL64) FCOEFF1(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP,N1LIST), &349: COMPLEX*16 FCOEFF1(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP,N1LIST), &
351:     & FCOEFF2(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP,N2LIST), FCOEFFS(NFSPACE,NFSPACE,NFSPACE)350:     & FCOEFF2(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP,N2LIST), FCOEFFS(NFSPACE,NFSPACE,NFSPACE)
352: DOUBLE PRECISION WAVEK(3, 2*NWAVE+1,2*NWAVE+1,2*NWAVE+1), K2(2*NWAVE+1,2*NWAVE+1,2*NWAVE+1), DIST2351: DOUBLE PRECISION WAVEK(3, 2*NWAVE+1,2*NWAVE+1,2*NWAVE+1), K2(2*NWAVE+1,2*NWAVE+1,2*NWAVE+1), DIST2
353: INTEGER I,J,K,JX,JY,JZ,NDISPS352: INTEGER I,J,K,JX,JY,JZ,NDISPS
354: 353: 
355: IF (DEBUG) WRITE(MYUNIT,'(A)') 'fastoverlap> starting group alignment'354: IF (DEBUG) WRITE(MYUNIT,'(A)') 'fastoverlap> starting group alignment'
356: IF (DEBUG) WRITE(MYUNIT,'(A,I5,A,I5)') 'fastoverlap> aligning ', N1LIST, ' structures with ', N2LIST355: IF (DEBUG) WRITE(MYUNIT,'(A,I5,A,I5)') 'fastoverlap> aligning ', N1LIST, ' structures with ', N2LIST
357: 356: 
358: CALL SETWAVEK(NWAVE,WAVEK,BOXLX,BOXLY,BOXLZ)357: CALL SETWAVEK(NWAVE,WAVEK,BOXLX,BOXLY,BOXLZ)
359: DO JZ=1,2*NWAVE+1358: DO JZ=1,2*NWAVE+1
360:     DO JY=1,2*NWAVE+1359:     DO JY=1,2*NWAVE+1
417: IMPLICIT NONE416: IMPLICIT NONE
418: 417: 
419: INTEGER, INTENT(IN) :: NATOMS, NDISPLACEMENTS, NFSPACE, NWAVE418: INTEGER, INTENT(IN) :: NATOMS, NDISPLACEMENTS, NFSPACE, NWAVE
420: LOGICAL, INTENT(IN) :: DEBUG419: LOGICAL, INTENT(IN) :: DEBUG
421: DOUBLE PRECISION, INTENT(IN) :: BOXLX, BOXLY, BOXLZ, KWIDTH420: DOUBLE PRECISION, INTENT(IN) :: BOXLX, BOXLY, BOXLZ, KWIDTH
422: DOUBLE PRECISION, INTENT(INOUT) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)421: DOUBLE PRECISION, INTENT(INOUT) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)
423: DOUBLE PRECISION, INTENT(OUT) :: DISTANCE, DIST2422: DOUBLE PRECISION, INTENT(OUT) :: DISTANCE, DIST2
424: 423: 
425: DOUBLE PRECISION WAVEK(3, 2*NWAVE+1,2*NWAVE+1,2*NWAVE+1), K2, DISTSAVE424: DOUBLE PRECISION WAVEK(3, 2*NWAVE+1,2*NWAVE+1,2*NWAVE+1), K2, DISTSAVE
426: DOUBLE PRECISION SAVEA(3*NATOMS), SAVEB(3*NATOMS)425: DOUBLE PRECISION SAVEA(3*NATOMS), SAVEB(3*NATOMS)
427: COMPLEX(KIND=REAL64) FCOEFFS(NFSPACE,NFSPACE,NFSPACE), FCOEFFA(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP), &426: COMPLEX*16 FCOEFFS(NFSPACE,NFSPACE,NFSPACE), FCOEFFA(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP), &
428:  & FCOEFFB(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP), FCOEFFDUMMYA(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP)427:  & FCOEFFB(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP), FCOEFFDUMMYA(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP)
429: INTEGER J, JX, JY, JZ, OPNUM, NDISPS, JXL, JYL, JZL, JXH, JYH, JZH, JXI, JYI, JZI428: INTEGER J, JX, JY, JZ, OPNUM, NDISPS, JXL, JYL, JZL, JXH, JYH, JZH, JXI, JYI, JZI
430: 429: 
431: CALL CHECKKEYWORDS()430: CALL CHECKKEYWORDS()
432: OHCELLTSAVE = OHCELLT431: OHCELLTSAVE = OHCELLT
433: OHCELLT = .FALSE.432: OHCELLT = .FALSE.
434: 433: 
435: ! Calculating Fourier Coefficients of COORDSA and COORDSB434: ! Calculating Fourier Coefficients of COORDSA and COORDSB
436: CALL SETWAVEK(NWAVE,WAVEK,BOXLX,BOXLY,BOXLZ)435: CALL SETWAVEK(NWAVE,WAVEK,BOXLX,BOXLY,BOXLZ)
437: CALL PERIODICFOURIERPERM(COORDSA,NATOMS,NWAVE,NFSPACE,WAVEK,FCOEFFA,NPERMGROUP)436: CALL PERIODICFOURIERPERM(COORDSA,NATOMS,NWAVE,NFSPACE,WAVEK,FCOEFFA,NPERMGROUP)
438: CALL PERIODICFOURIERPERM(COORDSB,NATOMS,NWAVE,NFSPACE,WAVEK,FCOEFFB,NPERMGROUP)437: CALL PERIODICFOURIERPERM(COORDSB,NATOMS,NWAVE,NFSPACE,WAVEK,FCOEFFB,NPERMGROUP)
439: 438: 
440: !FCOEFFS = CMPLX(0.D0, 0.D0, REAL64)439: !FCOEFFS = DCMPLX(0.D0, 0.D0)
441: FCOEFFA = CONJG(FCOEFFA)440: FCOEFFA = CONJG(FCOEFFA)
442: 441: 
443: ! Calculating Fourier Coefficients of overlap integral442: ! Calculating Fourier Coefficients of overlap integral
444: DO JZ=1,2*NWAVE+1443: DO JZ=1,2*NWAVE+1
445:     DO JY=1,2*NWAVE+1444:     DO JY=1,2*NWAVE+1
446:         DO JX=1,2*NWAVE+1445:         DO JX=1,2*NWAVE+1
447:             K2 = EXP(-0.5D0 * (WAVEK(1,JX,JY,JZ)**2 + WAVEK(2,JX,JY,JZ)**2 + WAVEK(3,JX,JY,JZ)**2)*KWIDTH**2)446:             K2 = EXP(-0.5D0 * (WAVEK(1,JX,JY,JZ)**2 + WAVEK(2,JX,JY,JZ)**2 + WAVEK(3,JX,JY,JZ)**2)*KWIDTH**2)
448:             FCOEFFA(JX,JY,JZ,:) = FCOEFFA(JX,JY,JZ,:) * K2447:             FCOEFFA(JX,JY,JZ,:) = FCOEFFA(JX,JY,JZ,:) * K2
449:             FCOEFFB(JX,JY,JZ,:) = FCOEFFB(JX,JY,JZ,:) * K2448:             FCOEFFB(JX,JY,JZ,:) = FCOEFFB(JX,JY,JZ,:) * K2
450:             !FCOEFFS(JX,JY,JZ) = SUM(FCOEFFA(JX,JY,JZ,:)*FCOEFFB(JX,JY,JZ,:))449:             !FCOEFFS(JX,JY,JZ) = SUM(FCOEFFA(JX,JY,JZ,:)*FCOEFFB(JX,JY,JZ,:))
461: 460: 
462: IF (OHCELLTSAVE) THEN461: IF (OHCELLTSAVE) THEN
463:     DISTSAVE = HUGE(DISTSAVE)462:     DISTSAVE = HUGE(DISTSAVE)
464:     DO OPNUM=1,48463:     DO OPNUM=1,48
465:         IF (DEBUG) WRITE(MYUNIT,'(A,I2)') 'fastoverlap> Trying Oh symmetry operation number ',OPNUM464:         IF (DEBUG) WRITE(MYUNIT,'(A,I2)') 'fastoverlap> Trying Oh symmetry operation number ',OPNUM
466:         CALL OHOPS(COORDSA,SAVEA,OPNUM,NATOMS)465:         CALL OHOPS(COORDSA,SAVEA,OPNUM,NATOMS)
467:         ! Applying octahedral symmetry operation to FCOEFFA466:         ! Applying octahedral symmetry operation to FCOEFFA
468:         CALL OHTRANSFORMCOEFFS(FCOEFFA, FCOEFFDUMMYA, NWAVE, NFSPACE-NWAVE-1, NPERMGROUP, OPNUM)467:         CALL OHTRANSFORMCOEFFS(FCOEFFA, FCOEFFDUMMYA, NWAVE, NFSPACE-NWAVE-1, NPERMGROUP, OPNUM)
469: 468: 
470:         ! Recalculating Fourier Coefficients469:         ! Recalculating Fourier Coefficients
471: !        FCOEFFS = CMPLX(0.D0, 0.D0, REAL64)470: !        FCOEFFS = DCMPLX(0.D0, 0.D0)
472: !        DO J=1,NPERMGROUP471: !        DO J=1,NPERMGROUP
473: !            DO JZ=1,2*NWAVE+1472: !            DO JZ=1,2*NWAVE+1
474: !                DO JY=1,2*NWAVE+1473: !                DO JY=1,2*NWAVE+1
475: !                    DO JX=1,2*NWAVE+1474: !                    DO JX=1,2*NWAVE+1
476: !                        FCOEFFS(JX,JY,JZ) = FCOEFFS(JX,JY,JZ) + &475: !                        FCOEFFS(JX,JY,JZ) = FCOEFFS(JX,JY,JZ) + &
477: !                        & FCOEFFDUMMYA(JX,JY,JZ,J)*FCOEFFB(JX,JY,JZ,J)476: !                        & FCOEFFDUMMYA(JX,JY,JZ,J)*FCOEFFB(JX,JY,JZ,J)
478: !                    ENDDO477: !                    ENDDO
479: !                ENDDO478: !                ENDDO
480: !            ENDDO479: !            ENDDO
481: !        ENDDO480: !        ENDDO
517: END SUBROUTINE ALIGN1516: END SUBROUTINE ALIGN1
518: 517: 
519: SUBROUTINE ALIGNCOEFFS(COORDSB,COORDSA,NATOMS,DEBUG,FCOEFFS,NFSPACE,LX,LY,LZ,DISTANCE,DIST2,NDISPS)518: SUBROUTINE ALIGNCOEFFS(COORDSB,COORDSA,NATOMS,DEBUG,FCOEFFS,NFSPACE,LX,LY,LZ,DISTANCE,DIST2,NDISPS)
520: 519: 
521: USE FASTOVERLAPUTILS, ONLY : FFT3D, FINDPEAKS520: USE FASTOVERLAPUTILS, ONLY : FFT3D, FINDPEAKS
522: IMPLICIT NONE521: IMPLICIT NONE
523: 522: 
524: INTEGER, INTENT(INOUT) :: NDISPS523: INTEGER, INTENT(INOUT) :: NDISPS
525: INTEGER, INTENT(IN) :: NATOMS, NFSPACE524: INTEGER, INTENT(IN) :: NATOMS, NFSPACE
526: LOGICAL, INTENT(IN) :: DEBUG525: LOGICAL, INTENT(IN) :: DEBUG
527: COMPLEX(KIND=REAL64), INTENT(IN) ::  FCOEFFS(NFSPACE,NFSPACE,NFSPACE)526: COMPLEX*16, INTENT(IN) ::  FCOEFFS(NFSPACE,NFSPACE,NFSPACE)
528: DOUBLE PRECISION, INTENT(IN) :: LX, LY, LZ527: DOUBLE PRECISION, INTENT(IN) :: LX, LY, LZ
529: DOUBLE PRECISION, INTENT(INOUT) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)528: DOUBLE PRECISION, INTENT(INOUT) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)
530: DOUBLE PRECISION, INTENT(OUT) :: DISTANCE, DIST2529: DOUBLE PRECISION, INTENT(OUT) :: DISTANCE, DIST2
531: 530: 
532: COMPLEX(KIND=REAL64) FSPACECMPLX(NFSPACE,NFSPACE,NFSPACE)531: COMPLEX*16 FSPACECMPLX(NFSPACE,NFSPACE,NFSPACE)
533: DOUBLE PRECISION FSPACE(NFSPACE,NFSPACE,NFSPACE), DISPS(NDISPS,3), R(3,3), BESTDIST532: DOUBLE PRECISION FSPACE(NFSPACE,NFSPACE,NFSPACE), DISPS(NDISPS,3), R(3,3), BESTDIST
534: DOUBLE PRECISION AMPLITUDES(NDISPS)533: DOUBLE PRECISION AMPLITUDES(NDISPS)
535: INTEGER J, J1534: INTEGER J, J1
536: 535: 
537: BOXLX = LX; BOXLY = LY; BOXLZ = LZ536: BOXLX = LX; BOXLY = LY; BOXLZ = LZ
538: 537: 
539: CALL FFT3D(NFSPACE,NFSPACE,NFSPACE,FCOEFFS,FSPACECMPLX)538: CALL FFT3D(NFSPACE,NFSPACE,NFSPACE,FCOEFFS,FSPACECMPLX)
540: FSPACE = ABS(FSPACECMPLX)539: FSPACE = ABS(FSPACECMPLX)
541: 540: 
542: CALL FINDPEAKS(FSPACE, DISPS, AMPLITUDES, NDISPS, DEBUG)541: CALL FINDPEAKS(FSPACE, DISPS, AMPLITUDES, NDISPS, DEBUG)
615: ! NATOMS: system size614: ! NATOMS: system size
616: ! NWAVE: number of wavevectors modes, FCOEFF will have (2*NWAVE+1)^3 elements615: ! NWAVE: number of wavevectors modes, FCOEFF will have (2*NWAVE+1)^3 elements
617: ! COORDS: coordinate vector616: ! COORDS: coordinate vector
618: ! WAVEK: wavevectors617: ! WAVEK: wavevectors
619: ! FCOEFF: fourier coefficients of coordinates618: ! FCOEFF: fourier coefficients of coordinates
620: 619: 
621: IMPLICIT NONE620: IMPLICIT NONE
622: 621: 
623: INTEGER, INTENT(IN) :: NATOMS, NWAVE, NCOEFF622: INTEGER, INTENT(IN) :: NATOMS, NWAVE, NCOEFF
624: DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS), WAVEK(3,2*NWAVE+1,2*NWAVE+1,2*NWAVE+1)623: DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS), WAVEK(3,2*NWAVE+1,2*NWAVE+1,2*NWAVE+1)
625: COMPLEX(KIND=REAL64), INTENT(OUT) :: FCOEFF(NCOEFF,NCOEFF,NCOEFF)624: COMPLEX*16, INTENT(OUT) :: FCOEFF(NCOEFF,NCOEFF,NCOEFF)
626: 625: 
627: INTEGER IX,IY,IZ, J, K626: INTEGER IX,IY,IZ, J, K
628: DOUBLE PRECISION :: KR627: DOUBLE PRECISION :: KR
629: 628: 
630: !FCOEFF = CMPLX(0.d0, 0.d0, REAL64)629: !FCOEFF = DCMPLX(0.d0,0.d0)
631: FCOEFF = CMPLX(0.0D0, 0.0D0, REAL64)630: FCOEFF = CMPLX(0.0D0, 0.0D0, kind=kind(1.0D0))
632: DO IX=1,2*NWAVE+1631: DO IX=1,2*NWAVE+1
633:     DO IY=1,2*NWAVE+1632:     DO IY=1,2*NWAVE+1
634:         DO IZ=1,2*NWAVE+1633:         DO IZ=1,2*NWAVE+1
635: !            FCOEFF(IX,IY,IZ) = CMPLX(0.d0, 0.d0, REAL64)634: !            FCOEFF(IX,IY,IZ) = DCMPLX(0.d0,0.d0)
636:             DO J=1, NATOMS635:             DO J=1, NATOMS
637:                 KR=0.d0636:                 KR=0.d0
638:                 DO K=1,3637:                 DO K=1,3
639:                     KR = KR + COORDS(3*J-3+K) * WAVEK(K,IX,IY,IZ)638:                     KR = KR + COORDS(3*J-3+K) * WAVEK(K,IX,IY,IZ)
640:                 ENDDO639:                 ENDDO
641:                 FCOEFF(IX,IY,IZ) = FCOEFF(IX,IY,IZ) + EXP(CMPLX(0.0D0, -KR, REAL64))640:                 FCOEFF(IX,IY,IZ) = FCOEFF(IX,IY,IZ) + EXP(CMPLX(0.0D0, -KR, kind=kind(1.0D0)))
642:             ENDDO641:             ENDDO
643:         ENDDO642:         ENDDO
644:     ENDDO643:     ENDDO
645: ENDDO644: ENDDO
646: 645: 
647: END SUBROUTINE PERIODICFOURIER646: END SUBROUTINE PERIODICFOURIER
648: 647: 
649: SUBROUTINE PERIODICFOURIERPERM(COORDS,NATOMS,NWAVE,NCOEFF,WAVEK,FCOEFF,NPERMGROUP)!,PERMGROUP,NPERMSIZE,NPERMGROUP)648: SUBROUTINE PERIODICFOURIERPERM(COORDS,NATOMS,NWAVE,NCOEFF,WAVEK,FCOEFF,NPERMGROUP)!,PERMGROUP,NPERMSIZE,NPERMGROUP)
650: ! Calculates Fourier coefficients of the different permutations of a structure.649: ! Calculates Fourier coefficients of the different permutations of a structure.
651: 650: 
652: IMPLICIT NONE651: IMPLICIT NONE
653: 652: 
654: INTEGER, INTENT(IN) :: NPERMGROUP653: INTEGER, INTENT(IN) :: NPERMGROUP
655: INTEGER, INTENT(IN) :: NATOMS, NWAVE, NCOEFF654: INTEGER, INTENT(IN) :: NATOMS, NWAVE, NCOEFF
656: !INTEGER, INTENT(IN) :: PERMGROUP(NATOMS), NPERMSIZE(NPERMGROUP)655: !INTEGER, INTENT(IN) :: PERMGROUP(NATOMS), NPERMSIZE(NPERMGROUP)
657: DOUBLE PRECISION, INTENT(IN) :: COORDS(NATOMS*3),  WAVEK(3,2*NWAVE+1,2*NWAVE+1,2*NWAVE+1)656: DOUBLE PRECISION, INTENT(IN) :: COORDS(NATOMS*3),  WAVEK(3,2*NWAVE+1,2*NWAVE+1,2*NWAVE+1)
658: !DOUBLE PRECISION, INTENT(IN) :: BOXLX,BOXLY,BOXLZ657: !DOUBLE PRECISION, INTENT(IN) :: BOXLX,BOXLY,BOXLZ
659: COMPLEX(KIND=REAL64), INTENT(OUT) :: FCOEFF(NCOEFF,NCOEFF,NCOEFF,NPERMGROUP)658: COMPLEX*16, INTENT(OUT) :: FCOEFF(NCOEFF,NCOEFF,NCOEFF,NPERMGROUP)
660: 659: 
661: COMPLEX(KIND=REAL64) FCOEFFDUMMY(NCOEFF,NCOEFF,NCOEFF)660: COMPLEX*16 FCOEFFDUMMY(NCOEFF,NCOEFF,NCOEFF)
662: DOUBLE PRECISION PDUMMY(3*NATOMS)661: DOUBLE PRECISION PDUMMY(3*NATOMS)
663: INTEGER NDUMMY, J1, J2, PATOMS662: INTEGER NDUMMY, J1, J2, PATOMS
664: 663: 
665: NDUMMY=1664: NDUMMY=1
666: 665: 
667: DO J1=1,NPERMGROUP666: DO J1=1,NPERMGROUP
668:     PATOMS=NPERMSIZE(J1)667:     PATOMS=NPERMSIZE(J1)
669:     DO J2=1,PATOMS668:     DO J2=1,PATOMS
670:         PDUMMY(3*(J2-1)+1)=COORDS(3*(PERMGROUP(NDUMMY+J2-1)-1)+1)669:         PDUMMY(3*(J2-1)+1)=COORDS(3*(PERMGROUP(NDUMMY+J2-1)-1)+1)
671:         PDUMMY(3*(J2-1)+2)=COORDS(3*(PERMGROUP(NDUMMY+J2-1)-1)+2)670:         PDUMMY(3*(J2-1)+2)=COORDS(3*(PERMGROUP(NDUMMY+J2-1)-1)+2)
676:     NDUMMY=NDUMMY+NPERMSIZE(J1)675:     NDUMMY=NDUMMY+NPERMSIZE(J1)
677: ENDDO676: ENDDO
678: 677: 
679: END SUBROUTINE PERIODICFOURIERPERM678: END SUBROUTINE PERIODICFOURIERPERM
680: 679: 
681: SUBROUTINE DOTFOURIERCOEFFS(FCOEFFB,FCOEFFA,NWAVE,NCOEFF,FCOEFFS,NPERMGROUP)680: SUBROUTINE DOTFOURIERCOEFFS(FCOEFFB,FCOEFFA,NWAVE,NCOEFF,FCOEFFS,NPERMGROUP)
682: 681: 
683: IMPLICIT NONE682: IMPLICIT NONE
684: 683: 
685: INTEGER, INTENT(IN) :: NPERMGROUP, NWAVE, NCOEFF684: INTEGER, INTENT(IN) :: NPERMGROUP, NWAVE, NCOEFF
686: COMPLEX(KIND=REAL64), INTENT(IN) :: FCOEFFA(NCOEFF,NCOEFF,NCOEFF,NPERMGROUP),FCOEFFB(NCOEFF,NCOEFF,NCOEFF,NPERMGROUP)685: COMPLEX*16, INTENT(IN) :: FCOEFFA(NCOEFF,NCOEFF,NCOEFF,NPERMGROUP),FCOEFFB(NCOEFF,NCOEFF,NCOEFF,NPERMGROUP)
687: COMPLEX(KIND=REAL64), INTENT(OUT) :: FCOEFFS(NCOEFF,NCOEFF,NCOEFF)686: COMPLEX*16, INTENT(OUT) :: FCOEFFS(NCOEFF,NCOEFF,NCOEFF)
688: 687: 
689: INTEGER J688: INTEGER J
690: 689: 
691: FCOEFFS = CMPLX(0.D0, 0.D0, REAL64)690: FCOEFFS = CMPLX(0.D0,0.D0)
692: 691: 
693: DO J=1,NPERMGROUP692: DO J=1,NPERMGROUP
694:     FCOEFFS = FCOEFFS + FCOEFFA(:,:,:,J)*FCOEFFB(:,:,:,J)693:     FCOEFFS = FCOEFFS + FCOEFFA(:,:,:,J)*FCOEFFB(:,:,:,J)
695: END DO694: END DO
696: 695: 
697: END SUBROUTINE DOTFOURIERCOEFFS696: END SUBROUTINE DOTFOURIERCOEFFS
698: 697: 
699: SUBROUTINE CALCFSPACE(NATOMS,COORDSA,COORDSB,NWAVE,WAVEK,KWIDTH,NFSPACE,FSPACE)!,NPERMGROUP)698: SUBROUTINE CALCFSPACE(NATOMS,COORDSA,COORDSB,NWAVE,WAVEK,KWIDTH,NFSPACE,FSPACE)!,NPERMGROUP)
700: !699: !
701: ! Calculate FASTOVERLAP real space array700: ! Calculate FASTOVERLAP real space array
707: 706: 
708: IMPLICIT NONE707: IMPLICIT NONE
709: 708: 
710: INTEGER, INTENT(IN) :: NATOMS, NWAVE, NFSPACE!, NPERMGROUP709: INTEGER, INTENT(IN) :: NATOMS, NWAVE, NFSPACE!, NPERMGROUP
711: DOUBLE PRECISION, INTENT(IN) :: KWIDTH710: DOUBLE PRECISION, INTENT(IN) :: KWIDTH
712: DOUBLE PRECISION, INTENT(IN) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)711: DOUBLE PRECISION, INTENT(IN) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)
713: DOUBLE PRECISION, INTENT(IN) :: WAVEK(3, 2*NWAVE+1,2*NWAVE+1,2*NWAVE+1)712: DOUBLE PRECISION, INTENT(IN) :: WAVEK(3, 2*NWAVE+1,2*NWAVE+1,2*NWAVE+1)
714: 713: 
715: DOUBLE PRECISION, INTENT(OUT) :: FSPACE(NFSPACE, NFSPACE, NFSPACE)714: DOUBLE PRECISION, INTENT(OUT) :: FSPACE(NFSPACE, NFSPACE, NFSPACE)
716: 715: 
717: COMPLEX(KIND=REAL64) FCOEFFA(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP), FCOEFFB(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP), COEFF716: COMPLEX*16 FCOEFFA(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP), FCOEFFB(NFSPACE,NFSPACE,NFSPACE,NPERMGROUP), COEFF
718: COMPLEX(KIND=REAL64) FCOEFF(NFSPACE,NFSPACE,NFSPACE)717: COMPLEX*16 FCOEFF(NFSPACE,NFSPACE,NFSPACE)
719: COMPLEX(KIND=REAL64) FSPACECMPLX(NFSPACE,NFSPACE,NFSPACE)718: COMPLEX*16 FSPACECMPLX(NFSPACE,NFSPACE,NFSPACE)
720: 719: 
721: INTEGER I, JX, JY, JZ720: INTEGER I, JX, JY, JZ
722: DOUBLE PRECISION K2721: DOUBLE PRECISION K2
723: 722: 
724: CALL PERIODICFOURIERPERM(COORDSA,NATOMS,NWAVE,NFSPACE,WAVEK,FCOEFFA,NPERMGROUP)!,PERMGROUP,NPERMSIZE,NPERMGROUP)723: CALL PERIODICFOURIERPERM(COORDSA,NATOMS,NWAVE,NFSPACE,WAVEK,FCOEFFA,NPERMGROUP)!,PERMGROUP,NPERMSIZE,NPERMGROUP)
725: CALL PERIODICFOURIERPERM(COORDSB,NATOMS,NWAVE,NFSPACE,WAVEK,FCOEFFB,NPERMGROUP)!,PERMGROUP,NPERMSIZE,NPERMGROUP)724: CALL PERIODICFOURIERPERM(COORDSB,NATOMS,NWAVE,NFSPACE,WAVEK,FCOEFFB,NPERMGROUP)!,PERMGROUP,NPERMSIZE,NPERMGROUP)
726: 725: 
727: FCOEFF = CMPLX(0.D0, 0.D0, REAL64)726: FCOEFF = CMPLX(0.D0, 0.D0, kind=kind(1.0D0))
728: FCOEFFB = CONJG(FCOEFFB)727: FCOEFFB = CONJG(FCOEFFB)
729: 728: 
730: DO JX=1,2*NWAVE+1729: DO JX=1,2*NWAVE+1
731:     DO JY=1,2*NWAVE+1730:     DO JY=1,2*NWAVE+1
732:         DO JZ=1,2*NWAVE+1731:         DO JZ=1,2*NWAVE+1
733:             COEFF = CMPLX(0.D0, 0.D0, REAL64)732:             COEFF = CMPLX(0.D0, 0.D0,kind=kind(1.0D0))
734:             K2 = -(WAVEK(1,JX,JY,JZ)**2 + WAVEK(2,JX,JY,JZ)**2 + WAVEK(3,JX,JY,JZ)**2)*KWIDTH**2733:             K2 = -(WAVEK(1,JX,JY,JZ)**2 + WAVEK(2,JX,JY,JZ)**2 + WAVEK(3,JX,JY,JZ)**2)*KWIDTH**2
735:             COEFF = SUM(FCOEFFA(JX,JY,JZ,:)*FCOEFFB(JX,JY,JZ,:))*EXP(K2)734:             COEFF = SUM(FCOEFFA(JX,JY,JZ,:)*FCOEFFB(JX,JY,JZ,:))*EXP(K2)
736:             FCOEFF(JX,JY,JZ) = COEFF735:             FCOEFF(JX,JY,JZ) = COEFF
737:         ENDDO736:         ENDDO
738:     ENDDO737:     ENDDO
739: ENDDO738: ENDDO
740: 739: 
741: !Set average overlap to 0740: !Set average overlap to 0
742: FCOEFF(NWAVE+1,NWAVE+1,NWAVE+1)=(0.d0,0.d0)741: FCOEFF(NWAVE+1,NWAVE+1,NWAVE+1)=(0.d0,0.d0)
743: 742: 
1168: !    J, Iind = ohopsmat[:,:,i].T.nonzero()1167: !    J, Iind = ohopsmat[:,:,i].T.nonzero()
1169: !    Jind = (Iind + 1) * ohopsmat[Iind,J,i].astype(int)1168: !    Jind = (Iind + 1) * ohopsmat[Iind,J,i].astype(int)
1170: !    print '    CASE ({})'.format(i+1)1169: !    print '    CASE ({})'.format(i+1)
1171: !    print prestring1170: !    print prestring
1172: !    print arraystr.format(Jstr, Js[Jind])1171: !    print arraystr.format(Jstr, Js[Jind])
1173: !    print poststring1172: !    print poststring
1174: !print 'END SELECT'1173: !print 'END SELECT'
1175: !1174: !
1176: IMPLICIT NONE1175: IMPLICIT NONE
1177: INTEGER, INTENT(IN) :: NF2, NWAVE, NPERMGROUP, OPNUM1176: INTEGER, INTENT(IN) :: NF2, NWAVE, NPERMGROUP, OPNUM
1178: COMPLEX(KIND=REAL64), INTENT(IN) :: FCOEFF(-NWAVE:NF2,-NWAVE:NF2,-NWAVE:NF2,NPERMGROUP)1177: COMPLEX*16, INTENT(IN) :: FCOEFF(-NWAVE:NF2,-NWAVE:NF2,-NWAVE:NF2,NPERMGROUP)
1179: COMPLEX(KIND=REAL64), INTENT(OUT) :: FCOEFFDUMMY(-NWAVE:NF2,-NWAVE:NF2,-NWAVE:NF2,NPERMGROUP)1178: COMPLEX*16, INTENT(OUT) :: FCOEFFDUMMY(-NWAVE:NF2,-NWAVE:NF2,-NWAVE:NF2,NPERMGROUP)
1180: 1179: 
1181: INTEGER JX, JY, JZ, J1180: INTEGER JX, JY, JZ, J
1182: 1181: 
1183: SELECT CASE (OPNUM)1182: SELECT CASE (OPNUM)
1184:     CASE (1)1183:     CASE (1)
1185:         DO J=1,NPERMGROUP1184:         DO J=1,NPERMGROUP
1186:             DO JZ=-NWAVE,NWAVE1185:             DO JZ=-NWAVE,NWAVE
1187:                 DO JY=-NWAVE,NWAVE1186:                 DO JY=-NWAVE,NWAVE
1188:                     DO JX=-NWAVE,NWAVE1187:                     DO JX=-NWAVE,NWAVE
1189:                         FCOEFFDUMMY(JX,JY,JZ,J) = FCOEFF(JX,JY,JZ,J)1188:                         FCOEFFDUMMY(JX,JY,JZ,J) = FCOEFF(JX,JY,JZ,J)


r33314/fastclusters.f90 2017-09-15 20:30:18.289831094 +0100 r33313/fastclusters.f90 2017-09-15 20:30:19.621848906 +0100
107: !        Helper Module Needed for Peak Fitting and FFT routines107: !        Helper Module Needed for Peak Fitting and FFT routines
108: !    DSOFT (DSOFT.f90) 108: !    DSOFT (DSOFT.f90) 
109: !        Module for performing discrete SO(3) transforms, depends on fftw.109: !        Module for performing discrete SO(3) transforms, depends on fftw.
110: 110: 
111: !***********************************************************************111: !***********************************************************************
112: 112: 
113: MODULE CLUSTERFASTOVERLAP113: MODULE CLUSTERFASTOVERLAP
114: 114: 
115: USE COMMONS, ONLY : PERMGROUP, NPERMSIZE, NPERMGROUP, NATOMS, BESTPERM, MYUNIT115: USE COMMONS, ONLY : PERMGROUP, NPERMSIZE, NPERMGROUP, NATOMS, BESTPERM, MYUNIT
116: USE FASTOVERLAPUTILS, ONLY : DUMMYA, DUMMYB, XBESTA, XBESTASAVE116: USE FASTOVERLAPUTILS, ONLY : DUMMYA, DUMMYB, XBESTA, XBESTASAVE
117: USE PREC, ONLY: INT64, REAL64 
118: 117: 
119: LOGICAL, SAVE :: PERMINVOPTSAVE, NOINVERSIONSAVE118: LOGICAL, SAVE :: PERMINVOPTSAVE, NOINVERSIONSAVE
120: 119: 
121: DOUBLE PRECISION, PARAMETER :: PI = 3.141592653589793D0120: DOUBLE PRECISION, PARAMETER :: PI = 3.141592653589793D0
122: 121: 
123: CONTAINS122: CONTAINS
124: 123: 
125: SUBROUTINE FOM_ALIGN_CLUSTERS(COORDSB, COORDSA, NATOMS, DEBUG, L, KWIDTH, DISTANCE, DIST2, RMATBEST, NROTATIONS)124: SUBROUTINE FOM_ALIGN_CLUSTERS(COORDSB, COORDSA, NATOMS, DEBUG, L, KWIDTH, DISTANCE, DIST2, RMATBEST, NROTATIONS)
126: 125: 
127: !  COORDSA becomes the optimal alignment of the optimal permutation(-inversion)126: !  COORDSA becomes the optimal alignment of the optimal permutation(-inversion)
144: USE FASTOVERLAPUTILS, ONLY: SETNATOMS143: USE FASTOVERLAPUTILS, ONLY: SETNATOMS
145: IMPLICIT NONE144: IMPLICIT NONE
146: 145: 
147: INTEGER, INTENT(IN) :: NATOMS, L146: INTEGER, INTENT(IN) :: NATOMS, L
148: INTEGER, INTENT(IN) :: NROTATIONS147: INTEGER, INTENT(IN) :: NROTATIONS
149: LOGICAL, INTENT(IN) :: DEBUG148: LOGICAL, INTENT(IN) :: DEBUG
150: DOUBLE PRECISION, INTENT(INOUT) :: KWIDTH ! Gaussian Kernel width149: DOUBLE PRECISION, INTENT(INOUT) :: KWIDTH ! Gaussian Kernel width
151: DOUBLE PRECISION, INTENT(INOUT) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)150: DOUBLE PRECISION, INTENT(INOUT) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)
152: DOUBLE PRECISION, INTENT(OUT) :: DISTANCE, DIST2, RMATBEST(3,3)151: DOUBLE PRECISION, INTENT(OUT) :: DISTANCE, DIST2, RMATBEST(3,3)
153: 152: 
154: COMPLEX(KIND=REAL64) PIMML(-L:L,-L:L,0:L)153: COMPLEX*16 PIMML(-L:L,-L:L,0:L)
155: COMPLEX(KIND=REAL64) IMML(-L:L,-L:L,0:L), YMLA(-L:L,0:L,NATOMS), YMLB(-L:L,0:L,NATOMS)154: COMPLEX*16 IMML(-L:L,-L:L,0:L), YMLA(-L:L,0:L,NATOMS), YMLB(-L:L,0:L,NATOMS)
156: 155: 
157: DOUBLE PRECISION SAVEA(3*NATOMS),SAVEB(3*NATOMS),COMA(3),COMB(3)156: DOUBLE PRECISION SAVEA(3*NATOMS),SAVEB(3*NATOMS),COMA(3),COMB(3)
158: DOUBLE PRECISION ANGLES(NROTATIONS,3), DISTSAVE, RMATSAVE(3,3), WORSTRAD, DIST2SAVE157: DOUBLE PRECISION ANGLES(NROTATIONS,3), DISTSAVE, RMATSAVE(3,3), WORSTRAD, DIST2SAVE
159: INTEGER J,J1,J2,M1,M2,IND2,NROT,NDUMMY,INVERT,PATOMS158: INTEGER J,J1,J2,M1,M2,IND2,NROT,NDUMMY,INVERT,PATOMS
160: INTEGER SAVEPERM(NATOMS), KEEPPERM(NATOMS)159: INTEGER SAVEPERM(NATOMS), KEEPPERM(NATOMS)
161: 160: 
162: ! Checking keywords are set properly161: ! Checking keywords are set properly
163: CALL CHECKKEYWORDS()162: CALL CHECKKEYWORDS()
164: 163: 
165: ! Allocate arrays164: ! Allocate arrays
184: ENDDO183: ENDDO
185: COMA = COMA/NATOMS184: COMA = COMA/NATOMS
186: COMB = COMB/NATOMS185: COMB = COMB/NATOMS
187: DO J=1,NATOMS186: DO J=1,NATOMS
188:     COORDSA(3*J-2:3*J) = COORDSA(3*J-2:3*J) - COMA187:     COORDSA(3*J-2:3*J) = COORDSA(3*J-2:3*J) - COMA
189:     COORDSB(3*J-2:3*J) = COORDSB(3*J-2:3*J) - COMB188:     COORDSB(3*J-2:3*J) = COORDSB(3*J-2:3*J) - COMB
190: ENDDO189: ENDDO
191: 190: 
192: 191: 
193: ! Calculating overlap integral separately for each permutation group192: ! Calculating overlap integral separately for each permutation group
194: IMML = CMPLX(0.D0, 0.D0, REAL64)193: IMML = CMPLX(0.D0,0.D0,8)
195: NDUMMY=1194: NDUMMY=1
196: DO J1=1,NPERMGROUP195: DO J1=1,NPERMGROUP
197:     PATOMS=INT(NPERMSIZE(J1),4)196:     PATOMS=INT(NPERMSIZE(J1),4)
198:     DO J2=1,PATOMS197:     DO J2=1,PATOMS
199:         IND2 = PERMGROUP(NDUMMY+J2-1)198:         IND2 = PERMGROUP(NDUMMY+J2-1)
200:         SAVEA(3*J2-2:3*J2)=COORDSA(3*IND2-2:3*IND2)199:         SAVEA(3*J2-2:3*J2)=COORDSA(3*IND2-2:3*IND2)
201:         SAVEB(3*J2-2:3*J2)=COORDSB(3*IND2-2:3*IND2)200:         SAVEB(3*J2-2:3*J2)=COORDSB(3*IND2-2:3*IND2)
202:     ENDDO201:     ENDDO
203:     CALL FOURIERCOEFFS(SAVEB,SAVEA,PATOMS,L,KWIDTH,PIMML,YMLB,YMLA)202:     CALL FOURIERCOEFFS(SAVEB,SAVEA,PATOMS,L,KWIDTH,PIMML,YMLB,YMLA)
204:     DO J=0,L203:     DO J=0,L
300: USE FASTOVERLAPUTILS, ONLY: SETNATOMS299: USE FASTOVERLAPUTILS, ONLY: SETNATOMS
301: IMPLICIT NONE300: IMPLICIT NONE
302: 301: 
303: INTEGER, INTENT(IN) :: NATOMS, N, L302: INTEGER, INTENT(IN) :: NATOMS, N, L
304: INTEGER, INTENT(IN) :: NROTATIONS303: INTEGER, INTENT(IN) :: NROTATIONS
305: LOGICAL, INTENT(IN) :: DEBUG304: LOGICAL, INTENT(IN) :: DEBUG
306: DOUBLE PRECISION, INTENT(IN) :: HWIDTH, KWIDTH305: DOUBLE PRECISION, INTENT(IN) :: HWIDTH, KWIDTH
307: DOUBLE PRECISION, INTENT(INOUT) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)306: DOUBLE PRECISION, INTENT(INOUT) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)
308: DOUBLE PRECISION, INTENT(OUT) :: DISTANCE, DIST2, RMATBEST(3,3)307: DOUBLE PRECISION, INTENT(OUT) :: DISTANCE, DIST2, RMATBEST(3,3)
309: 308: 
310: COMPLEX(KIND=REAL64) PIMML(-L:L,-L:L,0:L)309: COMPLEX*16 PIMML(-L:L,-L:L,0:L)
311: COMPLEX(KIND=REAL64) IMML(-L:L,-L:L,0:L), YMLA(-L:L,0:L,NATOMS), YMLB(-L:L,0:L,NATOMS)310: COMPLEX*16 IMML(-L:L,-L:L,0:L), YMLA(-L:L,0:L,NATOMS), YMLB(-L:L,0:L,NATOMS)
312: COMPLEX(KIND=REAL64) COEFFSA(0:N,-L:L,0:L,NPERMGROUP), COEFFSB(0:N,-L:L,0:L,NPERMGROUP)311: COMPLEX*16 COEFFSA(0:N,-L:L,0:L,NPERMGROUP), COEFFSB(0:N,-L:L,0:L,NPERMGROUP)
313: 312: 
314: DOUBLE PRECISION SAVEA(3*NATOMS),SAVEB(3*NATOMS)313: DOUBLE PRECISION SAVEA(3*NATOMS),SAVEB(3*NATOMS)
315: DOUBLE PRECISION ANGLES(NROTATIONS,3), DISTSAVE, RMATSAVE(3,3), WORSTRAD, DIST2SAVE314: DOUBLE PRECISION ANGLES(NROTATIONS,3), DISTSAVE, RMATSAVE(3,3), WORSTRAD, DIST2SAVE
316: INTEGER J,J1,J2,M1,M2,IND2,NROT,NDUMMY,INVERT,PATOMS315: INTEGER J,J1,J2,M1,M2,IND2,NROT,NDUMMY,INVERT,PATOMS
317: INTEGER SAVEPERM(NATOMS), KEEPPERM(NATOMS)316: INTEGER SAVEPERM(NATOMS), KEEPPERM(NATOMS)
318: 317: 
319: 318: 
320: ! Checking keywords are set properly319: ! Checking keywords are set properly
321: CALL CHECKKEYWORDS()320: CALL CHECKKEYWORDS()
322: 321: 
323: ! Allocate arrays322: ! Allocate arrays
324: CALL SETNATOMS(NATOMS)323: CALL SETNATOMS(NATOMS)
325: 324: 
326: ! Setting keywords for fastoverlap use of minpermdist, will be reset when exiting program325: ! Setting keywords for fastoverlap use of minpermdist, will be reset when exiting program
327: PERMINVOPTSAVE = PERMINVOPT326: PERMINVOPTSAVE = PERMINVOPT
328: NOINVERSIONSAVE = NOINVERSION327: NOINVERSIONSAVE = NOINVERSION
329: PERMINVOPT = .FALSE.328: PERMINVOPT = .FALSE.
330: NOINVERSION = .TRUE.329: NOINVERSION = .TRUE.
331: 330: 
332: ! Calculating overlap integral separately for each permutation group331: ! Calculating overlap integral separately for each permutation group
333: IMML = CMPLX(0.D0, 0.D0, REAL64)332: IMML = CMPLX(0.D0,0.D0,8)
334: NDUMMY=1333: NDUMMY=1
335: DO J1=1,NPERMGROUP334: DO J1=1,NPERMGROUP
336:     PATOMS=INT(NPERMSIZE(J1),4)335:     PATOMS=INT(NPERMSIZE(J1),4)
337:     DO J2=1,PATOMS336:     DO J2=1,PATOMS
338:         IND2 = PERMGROUP(NDUMMY+J2-1)337:         IND2 = PERMGROUP(NDUMMY+J2-1)
339:         SAVEA(3*J2-2:3*J2)=COORDSA(3*IND2-2:3*IND2)338:         SAVEA(3*J2-2:3*J2)=COORDSA(3*IND2-2:3*IND2)
340:         SAVEB(3*J2-2:3*J2)=COORDSB(3*IND2-2:3*IND2)339:         SAVEB(3*J2-2:3*J2)=COORDSB(3*IND2-2:3*IND2)
341:     ENDDO340:     ENDDO
342:     CALL HARMONICCOEFFS(SAVEA, PATOMS, COEFFSA(:,:,:,J1), N, L, HWIDTH, KWIDTH)341:     CALL HARMONICCOEFFS(SAVEA, PATOMS, COEFFSA(:,:,:,J1), N, L, HWIDTH, KWIDTH)
343:     CALL HARMONICCOEFFS(SAVEB, PATOMS, COEFFSB(:,:,:,J1), N, L, HWIDTH, KWIDTH)342:     CALL HARMONICCOEFFS(SAVEB, PATOMS, COEFFSB(:,:,:,J1), N, L, HWIDTH, KWIDTH)
357: 356: 
358: IF (PERMINVOPTSAVE.AND.(.NOT.(CHRMMT.OR.AMBERT.OR.AMBER12T))) THEN 357: IF (PERMINVOPTSAVE.AND.(.NOT.(CHRMMT.OR.AMBERT.OR.AMBER12T))) THEN 
359:     IF (DEBUG) WRITE(MYUNIT,'(A)') 'fastoverlap> inverting geometry for comparison with target'358:     IF (DEBUG) WRITE(MYUNIT,'(A)') 'fastoverlap> inverting geometry for comparison with target'
360:     ! Saving non inverted configuration359:     ! Saving non inverted configuration
361:     XBESTASAVE(1:3*NATOMS) = SAVEA(1:3*NATOMS)360:     XBESTASAVE(1:3*NATOMS) = SAVEA(1:3*NATOMS)
362:     KEEPPERM(1:NATOMS) = BESTPERM(1:NATOMS)361:     KEEPPERM(1:NATOMS) = BESTPERM(1:NATOMS)
363:     SAVEA = -COORDSA(1:3*NATOMS)362:     SAVEA = -COORDSA(1:3*NATOMS)
364:     NROT = NROTATIONS363:     NROT = NROTATIONS
365: 364: 
366:     ! Recalculating Fourier Coefficients for inverted COORDSA365:     ! Recalculating Fourier Coefficients for inverted COORDSA
367:     IMML = CMPLX(0.D0, 0.D0, REAL64)366:     IMML = CMPLX(0.D0,0.D0,8)
368:     NDUMMY=1367:     NDUMMY=1
369:     DO J1=1,NPERMGROUP368:     DO J1=1,NPERMGROUP
370:         DO J=0,L369:         DO J=0,L
371:             COEFFSA(:,:,J,J1) = COEFFSA(:,:,J,J1) * (-1)**(J)370:             COEFFSA(:,:,J,J1) = COEFFSA(:,:,J,J1) * (-1)**(J)
372:         ENDDO371:         ENDDO
373:         CALL DOTHARMONICCOEFFS(COEFFSB(:,:,:,J1), COEFFSA(:,:,:,J1), N, L, PIMML)372:         CALL DOTHARMONICCOEFFS(COEFFSB(:,:,:,J1), COEFFSA(:,:,:,J1), N, L, PIMML)
374:         DO J=0,L373:         DO J=0,L
375:             DO M2=-J,J374:             DO M2=-J,J
376:                 DO M1=-J,J375:                 DO M1=-J,J
377:                 IMML(M1,M2,J) = IMML(M1,M2,J) + PIMML(M1,M2,J)376:                 IMML(M1,M2,J) = IMML(M1,M2,J) + PIMML(M1,M2,J)
423: USE COMMONS, ONLY: PERMOPT, PERMINVOPT422: USE COMMONS, ONLY: PERMOPT, PERMINVOPT
424: 423: 
425: IMPLICIT NONE424: IMPLICIT NONE
426: 425: 
427: INTEGER, INTENT(IN) :: NATOMS, L426: INTEGER, INTENT(IN) :: NATOMS, L
428: INTEGER, INTENT(INOUT) :: NROTATIONS427: INTEGER, INTENT(INOUT) :: NROTATIONS
429: LOGICAL, INTENT(IN) :: DEBUG428: LOGICAL, INTENT(IN) :: DEBUG
430: DOUBLE PRECISION, INTENT(INOUT) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)429: DOUBLE PRECISION, INTENT(INOUT) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS)
431: DOUBLE PRECISION, INTENT(OUT) :: ANGLES(NROTATIONS,3)430: DOUBLE PRECISION, INTENT(OUT) :: ANGLES(NROTATIONS,3)
432: DOUBLE PRECISION, INTENT(OUT) :: DISTANCE, DIST2, RMATBEST(3,3)431: DOUBLE PRECISION, INTENT(OUT) :: DISTANCE, DIST2, RMATBEST(3,3)
433: COMPLEX(KIND=REAL64), INTENT(IN) :: IMML(-L:L,-L:L,0:L)432: COMPLEX*16, INTENT(IN) :: IMML(-L:L,-L:L,0:L)
434: 433: 
435: COMPLEX(KIND=REAL64) ILMM(0:L,0:2*L,0:2*L)434: COMPLEX*16 ILMM(0:L,0:2*L,0:2*L)
436: DOUBLE PRECISION OVERLAP(2*L+2,2*L+2,2*L+2)435: DOUBLE PRECISION OVERLAP(2*L+2,2*L+2,2*L+2)
437: DOUBLE PRECISION AMPLITUDES(NROTATIONS), BESTDIST, RMATSAVE(3,3), RMAT(3,3), WORSTRAD436: DOUBLE PRECISION AMPLITUDES(NROTATIONS), BESTDIST, RMATSAVE(3,3), RMAT(3,3), WORSTRAD
438: INTEGER J, J1437: INTEGER J, J1
439: 438: 
440: 439: 
441: CALL CALCOVERLAP(IMML, OVERLAP, L, ILMM)440: CALL CALCOVERLAP(IMML, OVERLAP, L, ILMM)
442: CALL FINDROTATIONS(OVERLAP, L, ANGLES, AMPLITUDES, NROTATIONS, DEBUG)441: CALL FINDROTATIONS(OVERLAP, L, ANGLES, AMPLITUDES, NROTATIONS, DEBUG)
443: IF (DEBUG) WRITE(MYUNIT,'(A,I3,A)') 'fastoverlap> found ', NROTATIONS, ' candidate rotations'442: IF (DEBUG) WRITE(MYUNIT,'(A,I3,A)') 'fastoverlap> found ', NROTATIONS, ' candidate rotations'
444: 443: 
445: 444: 
564: ! up to L, returns R, the distance COORD is from origin563: ! up to L, returns R, the distance COORD is from origin
565: ! Calculates value of Legendre Polynomial Recursively564: ! Calculates value of Legendre Polynomial Recursively
566: 565: 
567: ! UNSTABLE WHEN Z CLOSE TO 0 OR 1566: ! UNSTABLE WHEN Z CLOSE TO 0 OR 1
568: 567: 
569: IMPLICIT NONE568: IMPLICIT NONE
570: 569: 
571: DOUBLE PRECISION, INTENT(IN) :: COORD(3)570: DOUBLE PRECISION, INTENT(IN) :: COORD(3)
572: INTEGER, INTENT(IN) :: L571: INTEGER, INTENT(IN) :: L
573: DOUBLE PRECISION, INTENT(OUT) :: R572: DOUBLE PRECISION, INTENT(OUT) :: R
574: COMPLEX(KIND=REAL64), INTENT(OUT) :: YML(-L:L,0:L)573: COMPLEX*16, INTENT(OUT) :: YML(-L:L,0:L)
575: 574: 
576: INTEGER J, M, INDM1, INDM0, INDM2575: INTEGER J, M, INDM1, INDM0, INDM2
577: DOUBLE PRECISION THETA, PHI, Z, FACTORIALS(0:2*L), SQRTZ, SQRTMJ576: DOUBLE PRECISION THETA, PHI, Z, FACTORIALS(0:2*L), SQRTZ, SQRTMJ
578: COMPLEX(KIND=REAL64) EXPIM(-L:L)577: COMPLEX*16 EXPIM(-L:L)
579: 578: 
580: R = (COORD(1)**2+COORD(2)**2+COORD(3)**2)**0.5579: R = (COORD(1)**2+COORD(2)**2+COORD(3)**2)**0.5
581: PHI = ATAN2(COORD(2), COORD(1))580: PHI = ATAN2(COORD(2), COORD(1))
582: Z = COORD(3)/R581: Z = COORD(3)/R
583: SQRTZ = SQRT(1.D0-Z**2)582: SQRTZ = SQRT(1.D0-Z**2)
584: 583: 
585: !Calculating Associate Legendre Function584: !Calculating Associate Legendre Function
586: YML = CMPLX(0.D0, 0.D0, REAL64)585: YML = CMPLX(0.D0,0.D0, 8)
587: YML(0,0) = (4*PI)**(-0.5)586: YML(0,0) = (4*PI)**(-0.5)
588: 587: 
589: ! Initialising Recurrence for Associated Legendre Polynomials588: ! Initialising Recurrence for Associated Legendre Polynomials
590: ! Calculating normalised Legendre Polynomials for better numerical stability589: ! Calculating normalised Legendre Polynomials for better numerical stability
591: ! Pnorm^m_l = \sqrt{(l-m)!/(l+m)!} P^m_l590: ! Pnorm^m_l = \sqrt{(l-m)!/(l+m)!} P^m_l
592: DO J=0, L-1591: DO J=0, L-1
593:     YML(J+1,J+1) = - SQRT((2.D0*J+1.D0)/(2.D0*J+2.D0)) * SQRTZ* YML(J,J)592:     YML(J+1,J+1) = - SQRT((2.D0*J+1.D0)/(2.D0*J+2.D0)) * SQRTZ* YML(J,J)
594:     ! Calculating first recurrence term593:     ! Calculating first recurrence term
595:     YML(J, J+1) = -SQRT(2.D0*(J+1))*Z/SQRTZ * YML(J+1, J+1)594:     YML(J, J+1) = -SQRT(2.D0*(J+1))*Z/SQRTZ * YML(J+1, J+1)
596: ENDDO595: ENDDO
598: ! Recurrence for normalised Associated Legendre Polynomials597: ! Recurrence for normalised Associated Legendre Polynomials
599: DO J=1,L598: DO J=1,L
600:     DO M=J-1,-J+1,-1599:     DO M=J-1,-J+1,-1
601:         SQRTMJ = SQRT((J+M)*(J-M+1.D0))600:         SQRTMJ = SQRT((J+M)*(J-M+1.D0))
602:         YML(M-1, J) = -2*M*Z/SQRTMJ/SQRTZ * YML(M, J) - SQRT((J-M)*(J+M+1.D0))/SQRTMJ * YML(M+1,J)601:         YML(M-1, J) = -2*M*Z/SQRTMJ/SQRTZ * YML(M, J) - SQRT((J-M)*(J+M+1.D0))/SQRTMJ * YML(M+1,J)
603:     ENDDO602:     ENDDO
604: ENDDO603: ENDDO
605: 604: 
606: ! Calculating exp(imPHI) component605: ! Calculating exp(imPHI) component
607: DO M=-L,L606: DO M=-L,L
608:     EXPIM(M) = EXP(CMPLX(0.D0, M*PHI, REAL64))607:     EXPIM(M) = EXP(CMPLX(0.D0, M*PHI, 8))
609: ENDDO608: ENDDO
610: 609: 
611: ! Calculate Spherical Harmonics610: ! Calculate Spherical Harmonics
612: DO J=1,L611: DO J=1,L
613:     DO M=-J,J612:     DO M=-J,J
614:         INDM0 = MODULO(M, 2*L+1)613:         INDM0 = MODULO(M, 2*L+1)
615:         YML(M,J) = EXPIM(M)*YML(M,J) * SQRT((2.D0*J+1.D0))614:         YML(M,J) = EXPIM(M)*YML(M,J) * SQRT((2.D0*J+1.D0))
616:     ENDDO615:     ENDDO
617: ENDDO616: ENDDO
618: 617: 
622: 621: 
623: ! Calculates the Spherical Harmonics associated with coordinate COORD622: ! Calculates the Spherical Harmonics associated with coordinate COORD
624: ! up to L, returns R, the distance COORD is from origin623: ! up to L, returns R, the distance COORD is from origin
625: ! Calculates value of Legendre Polynomial Recursively624: ! Calculates value of Legendre Polynomial Recursively
626: 625: 
627: IMPLICIT NONE626: IMPLICIT NONE
628: 627: 
629: DOUBLE PRECISION, INTENT(IN) :: COORD(3)628: DOUBLE PRECISION, INTENT(IN) :: COORD(3)
630: INTEGER, INTENT(IN) :: L629: INTEGER, INTENT(IN) :: L
631: DOUBLE PRECISION, INTENT(OUT) :: R630: DOUBLE PRECISION, INTENT(OUT) :: R
632: COMPLEX(KIND=REAL64), INTENT(OUT) :: YML(-L:L,0:L)631: COMPLEX*16, INTENT(OUT) :: YML(-L:L,0:L)
633: 632: 
634: INTEGER J, M, INDM1, INDM0, INDM2, ISIG633: INTEGER J, M, INDM1, INDM0, INDM2, ISIG
635: DOUBLE PRECISION THETA, PHI, Z, FACTORIALS(0:2*L), SQRTZ, SQRTMJ, PLM(0:L), IPN(0:L), FACT634: DOUBLE PRECISION THETA, PHI, Z, FACTORIALS(0:2*L), SQRTZ, SQRTMJ, PLM(0:L), IPN(0:L), FACT
636: COMPLEX(KIND=REAL64) EXPIM(-L:L)635: COMPLEX*16 EXPIM(-L:L)
637: 636: 
638: R = (COORD(1)**2+COORD(2)**2+COORD(3)**2)**0.5637: R = (COORD(1)**2+COORD(2)**2+COORD(3)**2)**0.5
639: PHI = ATAN2(COORD(2), COORD(1))638: PHI = ATAN2(COORD(2), COORD(1))
640: Z = COORD(3)/R639: Z = COORD(3)/R
641: 640: 
642: !Calculating Associate Legendre Function641: !Calculating Associate Legendre Function
643: YML = CMPLX(0.D0, 0.D0, REAL64)642: YML = CMPLX(0.D0,0.D0, 8)
644: YML(0,0) = (4*PI)**(-0.5)643: YML(0,0) = (4*PI)**(-0.5)
645: 644: 
646: FACT = (2*PI)**(-0.5)645: FACT = (2*PI)**(-0.5)
647: 646: 
648: DO J=0, L647: DO J=0, L
649:     ! Calculate Normalised Legendre Polynomial648:     ! Calculate Normalised Legendre Polynomial
650:     CALL XDNRMP(J,0,J,Z,1,PLM(0:J),IPN(0:J),ISIG)649:     CALL XDNRMP(J,0,J,Z,1,PLM(0:J),IPN(0:J),ISIG)
651:     YML(0:J,J) = PLM(0:J) * FACT650:     YML(0:J,J) = PLM(0:J) * FACT
652:     DO M=1,J651:     DO M=1,J
653:         YML(-M,J) = YML(M,J)652:         YML(-M,J) = YML(M,J)
654:         YML(M,J) = YML(-M,J) * (-1)**M653:         YML(M,J) = YML(-M,J) * (-1)**M
655:     ENDDO654:     ENDDO
656: ENDDO655: ENDDO
657: 656: 
658: ! Calculating exp(imPHI) component657: ! Calculating exp(imPHI) component
659: DO M=-L,L658: DO M=-L,L
660:     EXPIM(M) = EXP(CMPLX(0.D0, M*PHI, REAL64))659:     EXPIM(M) = EXP(CMPLX(0.D0, M*PHI, 8))
661: ENDDO660: ENDDO
662: 661: 
663: ! Calculate Spherical Harmonics662: ! Calculate Spherical Harmonics
664: DO J=1,L663: DO J=1,L
665:     DO M=-J,J664:     DO M=-J,J
666:         INDM0 = MODULO(M, 2*L+1)665:         INDM0 = MODULO(M, 2*L+1)
667:         YML(M,J) = EXPIM(M)*YML(M,J) !* SQRT((2.D0*J+1.D0))666:         YML(M,J) = EXPIM(M)*YML(M,J) !* SQRT((2.D0*J+1.D0))
668:     ENDDO667:     ENDDO
669: ENDDO668: ENDDO
670: 669: 
675: !674: !
676: ! For a set of Gaussian Kernels of width KWIDTH at COORDS, 675: ! For a set of Gaussian Kernels of width KWIDTH at COORDS, 
677: ! this will calculate the coefficients of the isotropic quantum harmonic basis676: ! this will calculate the coefficients of the isotropic quantum harmonic basis
678: ! cnlm with length scale HWIDTH up to N and L.677: ! cnlm with length scale HWIDTH up to N and L.
679: !678: !
680: 679: 
681: IMPLICIT NONE680: IMPLICIT NONE
682: 681: 
683: INTEGER, INTENT(IN) :: NATOMS, N, L682: INTEGER, INTENT(IN) :: NATOMS, N, L
684: DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS), HWIDTH, KWIDTH683: DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS), HWIDTH, KWIDTH
685: COMPLEX(KIND=REAL64), INTENT(OUT) :: CNML(0:N,-L:L,0:L)684: COMPLEX*16, INTENT(OUT) :: CNML(0:N,-L:L,0:L)
686: 685: 
687: COMPLEX(KIND=REAL64) :: YML(-L:L,0:L)686: COMPLEX*16 :: YML(-L:L,0:L)
688: DOUBLE PRECISION HARMCOEFFS(0:2*N+L,0:N,0:L), DNL(0:N,0:L+2*N), RJ687: DOUBLE PRECISION HARMCOEFFS(0:2*N+L,0:N,0:L), DNL(0:N,0:L+2*N), RJ
689: INTEGER I,J,K,SI,M,INDM, S688: INTEGER I,J,K,SI,M,INDM, S
690: 689: 
691: CNML = CMPLX(0.D0, 0.D0, REAL64)690: CNML = CMPLX(0.D0,0.D0,8)
692: 691: 
693: DO K=1,NATOMS692: DO K=1,NATOMS
694:     CALL RYML(COORDS(3*K-2:3*K), RJ, YML, L)693:     CALL RYML(COORDS(3*K-2:3*K), RJ, YML, L)
695:     CALL HARMONICNL(N,L+2*N,RJ,KWIDTH,HWIDTH,DNL)694:     CALL HARMONICNL(N,L+2*N,RJ,KWIDTH,HWIDTH,DNL)
696:     DO J=0,L695:     DO J=0,L
697:         DO M=-J,J696:         DO M=-J,J
698:             INDM = MODULO(M,2*L+1)697:             INDM = MODULO(M,2*L+1)
699:             DO I=0,N698:             DO I=0,N
700:                 CNML(I,M,J) = CNML(I,M,J) + DNL(I,J) * CONJG(YML(M,J))699:                 CNML(I,M,J) = CNML(I,M,J) + DNL(I,J) * CONJG(YML(M,J))
701:             ENDDO700:             ENDDO
711: ! For a set of Gaussian Kernels of width KWIDTH at COORDS, 710: ! For a set of Gaussian Kernels of width KWIDTH at COORDS, 
712: ! this will calculate the coefficients of the isotropic quantum harmonic basis711: ! this will calculate the coefficients of the isotropic quantum harmonic basis
713: ! cnlm with length scale HWIDTH up to N and L.712: ! cnlm with length scale HWIDTH up to N and L.
714: ! Returns coefficients of the different permutations groups713: ! Returns coefficients of the different permutations groups
715: !714: !
716: 715: 
717: IMPLICIT NONE716: IMPLICIT NONE
718: 717: 
719: INTEGER, INTENT(IN) :: NATOMS, N, L, NPERMGROUP718: INTEGER, INTENT(IN) :: NATOMS, N, L, NPERMGROUP
720: DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS), HWIDTH, KWIDTH719: DOUBLE PRECISION, INTENT(IN) :: COORDS(3*NATOMS), HWIDTH, KWIDTH
721: COMPLEX(KIND=REAL64), INTENT(OUT) :: CNML(0:N,-L:L,0:L,1:NPERMGROUP)720: COMPLEX*16, INTENT(OUT) :: CNML(0:N,-L:L,0:L,1:NPERMGROUP)
722: 721: 
723: DOUBLE PRECISION DUMMY(3*NATOMS)722: DOUBLE PRECISION DUMMY(3*NATOMS)
724: INTEGER J1, J2, IND2, NDUMMY, PATOMS723: INTEGER J1, J2, IND2, NDUMMY, PATOMS
725: 724: 
726: ! Calculating overlap integral separately for each permutation group725: ! Calculating overlap integral separately for each permutation group
727: NDUMMY=1726: NDUMMY=1
728: DO J1=1,NPERMGROUP727: DO J1=1,NPERMGROUP
729:     PATOMS=NPERMSIZE(J1)728:     PATOMS=NPERMSIZE(J1)
730:     DO J2=1,PATOMS729:     DO J2=1,PATOMS
731:         IND2 = PERMGROUP(NDUMMY+J2-1)730:         IND2 = PERMGROUP(NDUMMY+J2-1)
736: ENDDO735: ENDDO
737: 736: 
738: END SUBROUTINE HARMONICCOEFFSPERM737: END SUBROUTINE HARMONICCOEFFSPERM
739: 738: 
740: SUBROUTINE HARMONICCOEFFSMULTI(COORDSLIST,NATOMS,NLIST,CNMLLIST,N,L,HWIDTH,KWIDTH,NPERMGROUP)739: SUBROUTINE HARMONICCOEFFSMULTI(COORDSLIST,NATOMS,NLIST,CNMLLIST,N,L,HWIDTH,KWIDTH,NPERMGROUP)
741: 740: 
742: IMPLICIT NONE741: IMPLICIT NONE
743: 742: 
744: INTEGER, INTENT(IN) :: NATOMS, NLIST, N, L, NPERMGROUP743: INTEGER, INTENT(IN) :: NATOMS, NLIST, N, L, NPERMGROUP
745: DOUBLE PRECISION, INTENT(IN) :: COORDSLIST(3*NATOMS, NLIST), HWIDTH, KWIDTH744: DOUBLE PRECISION, INTENT(IN) :: COORDSLIST(3*NATOMS, NLIST), HWIDTH, KWIDTH
746: COMPLEX(KIND=REAL64), INTENT(OUT) :: CNMLLIST(0:N,-L:L,0:L,1:NPERMGROUP, NLIST)745: COMPLEX*16, INTENT(OUT) :: CNMLLIST(0:N,-L:L,0:L,1:NPERMGROUP, NLIST)
747: 746: 
748: INTEGER I747: INTEGER I
749: 748: 
750: !write(*,*) NATOMS, NLIST, N, L, NPERMGROUP749: !write(*,*) NATOMS, NLIST, N, L, NPERMGROUP
751: !WRITE(*,*) SHAPE(CNMLLIST), SHAPE(COORDSLIST)750: !WRITE(*,*) SHAPE(CNMLLIST), SHAPE(COORDSLIST)
752: 751: 
753: DO I=1,NLIST752: DO I=1,NLIST
754:     CALL HARMONICCOEFFSPERM(COORDSLIST(:,I),NATOMS,CNMLLIST(:,:,:,:,I),N,L,HWIDTH,KWIDTH,NPERMGROUP)753:     CALL HARMONICCOEFFSPERM(COORDSLIST(:,I),NATOMS,CNMLLIST(:,:,:,:,I),N,L,HWIDTH,KWIDTH,NPERMGROUP)
755: ENDDO754: ENDDO
756: 755: 
757: END SUBROUTINE HARMONICCOEFFSMULTI756: END SUBROUTINE HARMONICCOEFFSMULTI
758: 757: 
759: SUBROUTINE DOTHARMONICCOEFFS(C1NML, C2NML, N, L, IMML)758: SUBROUTINE DOTHARMONICCOEFFS(C1NML, C2NML, N, L, IMML)
760: 759: 
761: IMPLICIT NONE760: IMPLICIT NONE
762: 761: 
763: INTEGER, INTENT(IN) :: N, L762: INTEGER, INTENT(IN) :: N, L
764: COMPLEX(KIND=REAL64), INTENT(IN) :: C1NML(0:N,-L:L,0:L), C2NML(0:N,-L:L,0:L)763: COMPLEX*16, INTENT(IN) :: C1NML(0:N,-L:L,0:L), C2NML(0:N,-L:L,0:L)
765: COMPLEX(KIND=REAL64), INTENT(OUT) :: IMML(-L:L,-L:L,0:L)764: COMPLEX*16, INTENT(OUT) :: IMML(-L:L,-L:L,0:L)
766: 765: 
767: INTEGER I, J, M1, M2, INDM1, INDM2766: INTEGER I, J, M1, M2, INDM1, INDM2
768: 767: 
769: IMML = CMPLX(0.D0, 0.D0, REAL64)768: IMML = CMPLX(0.D0,0.D0,8)
770: 769: 
771: DO J=0,L770: DO J=0,L
772:     DO M2=-J,J771:     DO M2=-J,J
773:         DO M1=-J,J772:         DO M1=-J,J
774:             DO I=0,N773:             DO I=0,N
775:                 IMML(M1,M2,J) = IMML(M1,M2,J) + CONJG(C1NML(I,M1,J))*C2NML(I,M2,J)774:                 IMML(M1,M2,J) = IMML(M1,M2,J) + CONJG(C1NML(I,M1,J))*C2NML(I,M2,J)
776:             ENDDO775:             ENDDO
777:         ENDDO776:         ENDDO
778:     ENDDO777:     ENDDO
779: ENDDO778: ENDDO
780: 779: 
781: END SUBROUTINE DOTHARMONICCOEFFS780: END SUBROUTINE DOTHARMONICCOEFFS
782: 781: 
783: SUBROUTINE DOTHARMONICCOEFFSPERM(C1NML, C2NML, N, L, IMML, NPERMGROUP)782: SUBROUTINE DOTHARMONICCOEFFSPERM(C1NML, C2NML, N, L, IMML, NPERMGROUP)
784: 783: 
785: IMPLICIT NONE784: IMPLICIT NONE
786: 785: 
787: INTEGER, INTENT(IN) :: N, L, NPERMGROUP786: INTEGER, INTENT(IN) :: N, L, NPERMGROUP
788: COMPLEX(KIND=REAL64), INTENT(IN) :: C1NML(0:N,-L:L,0:L,NPERMGROUP), C2NML(0:N,-L:L,0:L,NPERMGROUP)787: COMPLEX*16, INTENT(IN) :: C1NML(0:N,-L:L,0:L,NPERMGROUP), C2NML(0:N,-L:L,0:L,NPERMGROUP)
789: COMPLEX(KIND=REAL64), INTENT(OUT) :: IMML(-L:L,-L:L,0:L)788: COMPLEX*16, INTENT(OUT) :: IMML(-L:L,-L:L,0:L)
790: 789: 
791: INTEGER I, J, M1, M2, K, INDM1, INDM2790: INTEGER I, J, M1, M2, K, INDM1, INDM2
792: 791: 
793: IMML = CMPLX(0.D0, 0.D0, REAL64)792: IMML = CMPLX(0.D0,0.D0,8)
794: 793: 
795: DO K=1,NPERMGROUP794: DO K=1,NPERMGROUP
796:     DO J=0,L795:     DO J=0,L
797:         DO M2=-J,J796:         DO M2=-J,J
798:             DO M1=-J,J797:             DO M1=-J,J
799:                 DO I=0,N798:                 DO I=0,N
800:                     IMML(M1,M2,J) = IMML(M1,M2,J) + CONJG(C1NML(I,M1,J,K))*C2NML(I,M2,J,K)799:                     IMML(M1,M2,J) = IMML(M1,M2,J) + CONJG(C1NML(I,M1,J,K))*C2NML(I,M2,J,K)
801:                 ENDDO800:                 ENDDO
802:             ENDDO801:             ENDDO
803:         ENDDO802:         ENDDO
804:     ENDDO803:     ENDDO
805: ENDDO804: ENDDO
806: 805: 
807: END SUBROUTINE DOTHARMONICCOEFFSPERM806: END SUBROUTINE DOTHARMONICCOEFFSPERM
808: 807: 
809: SUBROUTINE CALCSIMILARITY(C1NML, C2NML, N, L, NPERMGROUP, NORM, MAXOVER)808: SUBROUTINE CALCSIMILARITY(C1NML, C2NML, N, L, NPERMGROUP, NORM, MAXOVER)
810: 809: 
811: IMPLICIT NONE810: IMPLICIT NONE
812: 811: 
813: INTEGER, INTENT(IN) :: N, L, NPERMGROUP812: INTEGER, INTENT(IN) :: N, L, NPERMGROUP
814: COMPLEX(KIND=REAL64), INTENT(IN) :: C1NML(0:N,-L:L,0:L,NPERMGROUP), C2NML(0:N,-L:L,0:L,NPERMGROUP)813: COMPLEX*16, INTENT(IN) :: C1NML(0:N,-L:L,0:L,NPERMGROUP), C2NML(0:N,-L:L,0:L,NPERMGROUP)
815: DOUBLE PRECISION, INTENT(OUT) :: NORM, MAXOVER814: DOUBLE PRECISION, INTENT(OUT) :: NORM, MAXOVER
816: 815: 
817: COMPLEX(KIND=REAL64) IMML(-L:L,-L:L,0:L), ILMM(0:L,0:2*L,0:2*L)816: COMPLEX*16 IMML(-L:L,-L:L,0:L), ILMM(0:L,0:2*L,0:2*L)
818: DOUBLE PRECISION OVERLAP(2*L+2,2*L+2,2*L+2)817: DOUBLE PRECISION OVERLAP(2*L+2,2*L+2,2*L+2)
819: 818: 
820: INTEGER J,M1,M2819: INTEGER J,M1,M2
821: 820: 
822: CALL DOTHARMONICCOEFFSPERM(C1NML, C2NML, N, L, IMML, NPERMGROUP)821: CALL DOTHARMONICCOEFFSPERM(C1NML, C2NML, N, L, IMML, NPERMGROUP)
823: 822: 
824: ! Calculated average overlap823: ! Calculated average overlap
825: DO J=0,L824: DO J=0,L
826:     DO M2=-J,J825:     DO M2=-J,J
827:         DO M1=-J,J826:         DO M1=-J,J
833: ! Calculate max overlap832: ! Calculate max overlap
834: CALL CALCOVERLAP(IMML, OVERLAP, L, ILMM)833: CALL CALCOVERLAP(IMML, OVERLAP, L, ILMM)
835: MAXOVER = MAXVAL(OVERLAP)834: MAXOVER = MAXVAL(OVERLAP)
836: 835: 
837: END SUBROUTINE CALCSIMILARITY836: END SUBROUTINE CALCSIMILARITY
838: 837: 
839: SUBROUTINE CALCSIMILARITIES(C1NMLLIST,N1LIST,C2NMLLIST,N2LIST,N,L,NPERMGROUP,NORMS,MAXOVERS,SYM)838: SUBROUTINE CALCSIMILARITIES(C1NMLLIST,N1LIST,C2NMLLIST,N2LIST,N,L,NPERMGROUP,NORMS,MAXOVERS,SYM)
840: 839: 
841: IMPLICIT NONE840: IMPLICIT NONE
842: INTEGER, INTENT(IN) :: N1LIST, N2LIST, N, L, NPERMGROUP841: INTEGER, INTENT(IN) :: N1LIST, N2LIST, N, L, NPERMGROUP
843: COMPLEX(KIND=REAL64), INTENT(IN) :: C1NMLLIST(0:N,-L:L,0:L,NPERMGROUP,N1LIST), &842: COMPLEX*16, INTENT(IN) :: C1NMLLIST(0:N,-L:L,0:L,NPERMGROUP,N1LIST), &
844:     & C2NMLLIST(0:N,-L:L,0:L,NPERMGROUP,N2LIST)843:     & C2NMLLIST(0:N,-L:L,0:L,NPERMGROUP,N2LIST)
845: LOGICAL, INTENT(IN) :: SYM844: LOGICAL, INTENT(IN) :: SYM
846: DOUBLE PRECISION, INTENT(OUT) :: NORMS(N1LIST,N2LIST), MAXOVERS(N1LIST,N2LIST)845: DOUBLE PRECISION, INTENT(OUT) :: NORMS(N1LIST,N2LIST), MAXOVERS(N1LIST,N2LIST)
847: 846: 
848: INTEGER I1, I2847: INTEGER I1, I2
849: 848: 
850: IF (SYM) THEN849: IF (SYM) THEN
851:     ! if C1NMLLIST == C2NMLLIST then only need to calculate half the values850:     ! if C1NMLLIST == C2NMLLIST then only need to calculate half the values
852:     DO I1=1,N1LIST851:     DO I1=1,N1LIST
853:         DO I2=I1,N1LIST852:         DO I2=I1,N1LIST
870: END SUBROUTINE CALCSIMILARITIES869: END SUBROUTINE CALCSIMILARITIES
871: 870: 
872: SUBROUTINE CALCOVERLAPMATRICES(COORDSLIST,NATOMS,NLIST,N,L,HWIDTH,KWIDTH,NORMS,MAXOVERS)871: SUBROUTINE CALCOVERLAPMATRICES(COORDSLIST,NATOMS,NLIST,N,L,HWIDTH,KWIDTH,NORMS,MAXOVERS)
873: 872: 
874: IMPLICIT NONE873: IMPLICIT NONE
875: 874: 
876: INTEGER, INTENT(IN) :: NATOMS, NLIST, N, L875: INTEGER, INTENT(IN) :: NATOMS, NLIST, N, L
877: DOUBLE PRECISION, INTENT(IN) :: COORDSLIST(3*NATOMS, NLIST), HWIDTH, KWIDTH876: DOUBLE PRECISION, INTENT(IN) :: COORDSLIST(3*NATOMS, NLIST), HWIDTH, KWIDTH
878: DOUBLE PRECISION, INTENT(OUT) :: NORMS(NLIST,NLIST), MAXOVERS(NLIST,NLIST)877: DOUBLE PRECISION, INTENT(OUT) :: NORMS(NLIST,NLIST), MAXOVERS(NLIST,NLIST)
879: 878: 
880: COMPLEX(KIND=REAL64) CNMLLIST(0:N,-L:L,0:L,1:NPERMGROUP, NLIST)879: COMPLEX*16 CNMLLIST(0:N,-L:L,0:L,1:NPERMGROUP, NLIST)
881: 880: 
882: CALL HARMONICCOEFFSMULTI(COORDSLIST,NATOMS,NLIST,CNMLLIST,N,L,HWIDTH,KWIDTH,NPERMGROUP)881: CALL HARMONICCOEFFSMULTI(COORDSLIST,NATOMS,NLIST,CNMLLIST,N,L,HWIDTH,KWIDTH,NPERMGROUP)
883: CALL CALCSIMILARITIES(CNMLLIST,NLIST,CNMLLIST,NLIST,N,L,NPERMGROUP,NORMS,MAXOVERS,.TRUE.)882: CALL CALCSIMILARITIES(CNMLLIST,NLIST,CNMLLIST,NLIST,N,L,NPERMGROUP,NORMS,MAXOVERS,.TRUE.)
884: 883: 
885: END SUBROUTINE CALCOVERLAPMATRICES884: END SUBROUTINE CALCOVERLAPMATRICES
886: 885: 
887: SUBROUTINE FOURIERCOEFFS(COORDSB, COORDSA, NATOMS, L, KWIDTH, IMML, YMLB, YMLA)886: SUBROUTINE FOURIERCOEFFS(COORDSB, COORDSA, NATOMS, L, KWIDTH, IMML, YMLB, YMLA)
888: !887: !
889: ! Calculates S03 Coefficients of the overlap integral of two structures888: ! Calculates S03 Coefficients of the overlap integral of two structures
890: ! does this calculation by direct calculation of the overlap between every pair889: ! does this calculation by direct calculation of the overlap between every pair
891: ! of atoms, slower than the Harmonic basis, but slightly more accurate.890: ! of atoms, slower than the Harmonic basis, but slightly more accurate.
892: !891: !
893: 892: 
894: IMPLICIT NONE893: IMPLICIT NONE
895: INTEGER, INTENT(IN) :: NATOMS, L894: INTEGER, INTENT(IN) :: NATOMS, L
896: DOUBLE PRECISION, INTENT(IN) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS), KWIDTH895: DOUBLE PRECISION, INTENT(IN) :: COORDSA(3*NATOMS), COORDSB(3*NATOMS), KWIDTH
897: COMPLEX(KIND=REAL64), INTENT(OUT) :: IMML(-L:L,-L:L,0:L)896: COMPLEX*16, INTENT(OUT) :: IMML(-L:L,-L:L,0:L)
898: 897: 
899: COMPLEX(KIND=REAL64), INTENT(OUT) ::  YMLA(-L:L,0:L,NATOMS), YMLB(-L:L,0:L,NATOMS)898: COMPLEX*16, INTENT(OUT) ::  YMLA(-L:L,0:L,NATOMS), YMLB(-L:L,0:L,NATOMS)
900: DOUBLE PRECISION RA(NATOMS), RB(NATOMS), IL(0:L), R1R2, EXPRA(NATOMS), EXPRB(NATOMS), FACT, TMP899: DOUBLE PRECISION RA(NATOMS), RB(NATOMS), IL(0:L), R1R2, EXPRA(NATOMS), EXPRB(NATOMS), FACT, TMP
901: 900: 
902: INTEGER IA,IB,I,J,K,M1,M2,INDM1,INDM2901: INTEGER IA,IB,I,J,K,M1,M2,INDM1,INDM2
903: 902: 
904: YMLA = CMPLX(0.D0, 0.D0, REAL64)903: YMLA = CMPLX(0.D0,0.D0,8)
905: YMLB = CMPLX(0.D0, 0.D0, REAL64)904: YMLB = CMPLX(0.D0,0.D0,8)
906: ! Precalculate some values905: ! Precalculate some values
907: DO I=1,NATOMS906: DO I=1,NATOMS
908:     CALL RYML(COORDSA(3*I-2:3*I), RA(I), YMLA(:,:,I), L)907:     CALL RYML(COORDSA(3*I-2:3*I), RA(I), YMLA(:,:,I), L)
909:     CALL RYML(COORDSB(3*I-2:3*I), RB(I), YMLB(:,:,I), L)908:     CALL RYML(COORDSB(3*I-2:3*I), RB(I), YMLB(:,:,I), L)
910:     EXPRA(I) = EXP(-0.25D0 * RA(I)**2 / KWIDTH**2)909:     EXPRA(I) = EXP(-0.25D0 * RA(I)**2 / KWIDTH**2)
911:     EXPRB(I) = EXP(-0.25D0 * RB(I)**2 / KWIDTH**2)910:     EXPRB(I) = EXP(-0.25D0 * RB(I)**2 / KWIDTH**2)
912: ENDDO911: ENDDO
913: 912: 
914: FACT = 4.D0 * PI**2.5 * KWIDTH**3913: FACT = 4.D0 * PI**2.5 * KWIDTH**3
915: 914: 
916: IMML = CMPLX(0.D0, 0.D0, REAL64)915: IMML = CMPLX(0.D0,0.D0,8)
917: DO IA=1,NATOMS916: DO IA=1,NATOMS
918:     DO IB=1,NATOMS917:     DO IB=1,NATOMS
919:         ! Don't calculate cross terms for points separated by 4 kwidths to speed up calculation918:         ! Don't calculate cross terms for points separated by 4 kwidths to speed up calculation
920:         IF (ABS(RA(IA)-RB(IB)).LT.(4*KWIDTH)) THEN919:         IF (ABS(RA(IA)-RB(IB)).LT.(4*KWIDTH)) THEN
921:             R1R2 = 0.5D0 * RA(IA)*RB(IB)/KWIDTH**2920:             R1R2 = 0.5D0 * RA(IA)*RB(IB)/KWIDTH**2
922:             CALL SPHI(L, R1R2, K, IL)921:             CALL SPHI(L, R1R2, K, IL)
923:             TMP = FACT*EXPRA(IA)*EXPRB(IB)!*SQRT(PI/2/R1R2)922:             TMP = FACT*EXPRA(IA)*EXPRB(IB)!*SQRT(PI/2/R1R2)
924:             DO J=0,L923:             DO J=0,L
925:                 DO M2=-L,L924:                 DO M2=-L,L
926:                     DO M1=-L,L925:                     DO M1=-L,L
935: END SUBROUTINE FOURIERCOEFFS934: END SUBROUTINE FOURIERCOEFFS
936: 935: 
937: SUBROUTINE CALCOVERLAP(IMML, OVERLAP, L, ILMM)936: SUBROUTINE CALCOVERLAP(IMML, OVERLAP, L, ILMM)
938: ! Converts an array of SO(3) Fourier Coefficients to a discrete937: ! Converts an array of SO(3) Fourier Coefficients to a discrete
939: ! overlap array using a fast discrete SO(3) Fourier Transform (DSOFT)938: ! overlap array using a fast discrete SO(3) Fourier Transform (DSOFT)
940: 939: 
941: USE DSOFT, ONLY : ISOFT940: USE DSOFT, ONLY : ISOFT
942: 941: 
943: IMPLICIT NONE942: IMPLICIT NONE
944: INTEGER, INTENT(IN) :: L943: INTEGER, INTENT(IN) :: L
945: COMPLEX(KIND=REAL64), INTENT(IN) :: IMML(-L:L,-L:L,0:L)944: COMPLEX*16, INTENT(IN) :: IMML(-L:L,-L:L,0:L)
946: DOUBLE PRECISION, INTENT(OUT) :: OVERLAP(2*L+2,2*L+2,2*L+2)945: DOUBLE PRECISION, INTENT(OUT) :: OVERLAP(2*L+2,2*L+2,2*L+2)
947: 946: 
948: COMPLEX(KIND=REAL64), INTENT(OUT) :: ILMM(0:L,0:2*L,0:2*L)947: COMPLEX*16, INTENT(OUT) :: ILMM(0:L,0:2*L,0:2*L)
949: COMPLEX(KIND=REAL64) FROT(2*L+2,2*L+2,2*L+2)948: COMPLEX*16 FROT(2*L+2,2*L+2,2*L+2)
950: INTEGER I,J,M1,M2, NJ949: INTEGER I,J,M1,M2, NJ
951: INTEGER(KIND=INT64) BW950: INTEGER*8 BW
952: 951: 
953: ! Convert array into format usable by DSOFT:952: ! Convert array into format usable by DSOFT:
954: BW = INT(L+1,8)953: BW = INT(L+1,8)
955: NJ = 2*L + 1954: NJ = 2*L + 1
956: 955: 
957: ILMM = CMPLX(0.D0, 0.D0, REAL64)956: ILMM = CMPLX(0.D0, 0.D0, 8)
958: DO J=0,L957: DO J=0,L
959:     ILMM(J,0,0) = IMML(0,0,J)958:     ILMM(J,0,0) = IMML(0,0,J)
960:     DO M2=1,J959:     DO M2=1,J
961:         ILMM(J,0,M2) = IMML(0,M2,J)960:         ILMM(J,0,M2) = IMML(0,M2,J)
962:         ILMM(J,0,NJ-M2) = IMML(0,-M2,J)961:         ILMM(J,0,NJ-M2) = IMML(0,-M2,J)
963:         ILMM(J,M2,0) = IMML(M2,0,J)962:         ILMM(J,M2,0) = IMML(M2,0,J)
964:         ILMM(J,NJ-M2,0) = IMML(-M2,0,J)963:         ILMM(J,NJ-M2,0) = IMML(-M2,0,J)
965:         DO M1=1,J964:         DO M1=1,J
966:             ILMM(J,M1,M2) = IMML(M1,M2,J)965:             ILMM(J,M1,M2) = IMML(M1,M2,J)
967:             ILMM(J,NJ-M1,M2) = IMML(-M1,M2,J)966:             ILMM(J,NJ-M1,M2) = IMML(-M1,M2,J)


r33314/fastutils.f90 2017-09-15 20:30:18.513834088 +0100 r33313/fastutils.f90 2017-09-15 20:30:19.837851793 +0100
 56: !         FINDPEAK 56: !         FINDPEAK
 57: !         FINDPEAKS 57: !         FINDPEAKS
 58: !     FFT subroutines 58: !     FFT subroutines
 59: !         FFT3D 59: !         FFT3D
 60: !         IFFT3D 60: !         IFFT3D
 61: !         FFT1D 61: !         FFT1D
 62: !         IFFT1D 62: !         IFFT1D
 63: !*********************************************************************** 63: !***********************************************************************
 64: USE COMMONS, ONLY : PERMGROUP, NPERMSIZE, NPERMGROUP, NATOMS, BESTPERM, NSETS, SETS, MYUNIT 64: USE COMMONS, ONLY : PERMGROUP, NPERMSIZE, NPERMGROUP, NATOMS, BESTPERM, NSETS, SETS, MYUNIT
 65: USE FFTW3 65: USE FFTW3
 66: USE PREC, ONLY: INT32, INT64, REAL64 
 67:  66: 
 68: IMPLICIT NONE 67: IMPLICIT NONE
 69:  68: 
 70: ! Variables and arrays needed for peakfinding 69: ! Variables and arrays needed for peakfinding
 71: INTEGER, PARAMETER :: DEFAULTWIDTH=2 70: INTEGER, PARAMETER :: DEFAULTWIDTH=2
 72: DOUBLE PRECISION, PARAMETER :: DEFAULTTOL=1.D-6 71: DOUBLE PRECISION, PARAMETER :: DEFAULTTOL=1.D-6
 73: INTEGER, SAVE :: FSIZE, FSHAPE(3) 72: INTEGER, SAVE :: FSIZE, FSHAPE(3)
 74: DOUBLE PRECISION, SAVE, ALLOCATABLE :: FSPACE(:,:,:),FSPACECOPY(:,:,:),GAUSARRAY(:,:,:),FVEC(:),FJAC(:,:) 73: DOUBLE PRECISION, SAVE, ALLOCATABLE :: FSPACE(:,:,:),FSPACECOPY(:,:,:),GAUSARRAY(:,:,:),FVEC(:),FJAC(:,:)
 75:  74: 
 76: ! Stuff for permutational alignment 75: ! Stuff for permutational alignment
595: 594: 
596: !***********************************************************************595: !***********************************************************************
597: ! FFT subroutines596: ! FFT subroutines
598: !***********************************************************************597: !***********************************************************************
599:     598:     
600: SUBROUTINE FFT3D(NX, NY, NZ, IN, OUT)599: SUBROUTINE FFT3D(NX, NY, NZ, IN, OUT)
601: ! calculates forward FFT in 3D600: ! calculates forward FFT in 3D
602: IMPLICIT NONE601: IMPLICIT NONE
603: 602: 
604: INTEGER, INTENT(IN) :: NX, NY, NZ603: INTEGER, INTENT(IN) :: NX, NY, NZ
605: COMPLEX(KIND=REAL64), INTENT(IN) :: IN(NX, NY, NZ)604: COMPLEX*16, INTENT(IN) :: IN(NX, NY, NZ)
606: COMPLEX(KIND=REAL64), INTENT(OUT) :: OUT(NX, NY, NZ)605: COMPLEX*16, INTENT(OUT) :: OUT(NX, NY, NZ)
607: 606: 
608: !INCLUDE "fftw3.f90"607: !INCLUDE "fftw3.f90"
609: INTEGER(KIND=INT64) PLAN_FORWARD608: INTEGER*8 PLAN_FORWARD
610: 609: 
611: CALL DFFTW_PLAN_DFT_3D_(PLAN_FORWARD, NX, NY, NZ, IN, OUT, FFTW_FORWARD, FFTW_ESTIMATE )610: CALL DFFTW_PLAN_DFT_3D_(PLAN_FORWARD, NX, NY, NZ, IN, OUT, FFTW_FORWARD, FFTW_ESTIMATE )
612: CALL DFFTW_EXECUTE_(PLAN_FORWARD)611: CALL DFFTW_EXECUTE_(PLAN_FORWARD)
613: !CALL DFFTW_DESTROY_PLAN(PLAN_FORWARD)612: !CALL DFFTW_DESTROY_PLAN(PLAN_FORWARD)
614: 613: 
615: END SUBROUTINE FFT3D614: END SUBROUTINE FFT3D
616: 615: 
617: !***********************************************************************616: !***********************************************************************
618: 617: 
619: SUBROUTINE IFFT3D(NX, NY, NZ, IN, OUT)618: SUBROUTINE IFFT3D(NX, NY, NZ, IN, OUT)
620: 619: 
621: ! calculates UNNORMALISED inverse fourier transform so,620: ! calculates UNNORMALISED inverse fourier transform so,
622: ! IN == IFFT3D(NX,NY,NZ, FFT3D(NX,NY,NZ, IN))/(NX*NY*NZ)621: ! IN == IFFT3D(NX,NY,NZ, FFT3D(NX,NY,NZ, IN))/(NX*NY*NZ)
623: 622: 
624: IMPLICIT NONE623: IMPLICIT NONE
625: 624: 
626: INTEGER, INTENT(IN) :: NX, NY, NZ625: INTEGER, INTENT(IN) :: NX, NY, NZ
627: COMPLEX(KIND=REAL64), INTENT(IN) :: IN(NX, NY, NZ)626: COMPLEX*16, INTENT(IN) :: IN(NX, NY, NZ)
628: COMPLEX(KIND=REAL64), INTENT(OUT) :: OUT(NX, NY, NZ)627: COMPLEX*16, INTENT(OUT) :: OUT(NX, NY, NZ)
629: 628: 
630: !INCLUDE "fftw3.f90"629: !INCLUDE "fftw3.f90"
631: INTEGER(KIND=INT64) PLAN_BACKWARD630: INTEGER*8 PLAN_BACKWARD
632: 631: 
633: CALL DFFTW_PLAN_DFT_3D_(PLAN_BACKWARD,NX,NY,NZ,IN,OUT,FFTW_BACKWARD,FFTW_ESTIMATE)632: CALL DFFTW_PLAN_DFT_3D_(PLAN_BACKWARD,NX,NY,NZ,IN,OUT,FFTW_BACKWARD,FFTW_ESTIMATE)
634: CALL DFFTW_EXECUTE_(PLAN_BACKWARD)633: CALL DFFTW_EXECUTE_(PLAN_BACKWARD)
635: CALL DFFTW_DESTROY_PLAN_(PLAN_BACKWARD)634: CALL DFFTW_DESTROY_PLAN_(PLAN_BACKWARD)
636: 635: 
637: END SUBROUTINE IFFT3D636: END SUBROUTINE IFFT3D
638: 637: 
639: SUBROUTINE FFT1D(N, IN, OUT)638: SUBROUTINE FFT1D(N, IN, OUT)
640: ! calculates forward FFT in 1D639: ! calculates forward FFT in 1D
641: 640: 
642: IMPLICIT NONE641: IMPLICIT NONE
643: 642: 
644: INTEGER(KIND=INT32), INTENT(IN) :: N643: INTEGER*4, INTENT(IN) :: N
645: COMPLEX(KIND=REAL64), INTENT(IN) :: IN(N)644: COMPLEX*16, INTENT(IN) :: IN(N)
646: COMPLEX(KIND=REAL64), INTENT(OUT) :: OUT(N)645: COMPLEX*16, INTENT(OUT) :: OUT(N)
647: 646: 
648: !INCLUDE "fftw3.f90"647: !INCLUDE "fftw3.f90"
649: INTEGER(KIND=INT64) PLAN_FORWARD648: INTEGER*8 PLAN_FORWARD
650: 649: 
651: CALL DFFTW_PLAN_DFT_1D_(PLAN_FORWARD, N, IN, OUT, FFTW_FORWARD, FFTW_ESTIMATE )650: CALL DFFTW_PLAN_DFT_1D_(PLAN_FORWARD, N, IN, OUT, FFTW_FORWARD, FFTW_ESTIMATE )
652: CALL DFFTW_EXECUTE_(PLAN_FORWARD)651: CALL DFFTW_EXECUTE_(PLAN_FORWARD)
653: CALL DFFTW_DESTROY_PLAN_(PLAN_FORWARD)652: CALL DFFTW_DESTROY_PLAN_(PLAN_FORWARD)
654: 653: 
655: END SUBROUTINE FFT1D654: END SUBROUTINE FFT1D
656: 655: 
657: !***********************************************************************656: !***********************************************************************
658: 657: 
659: SUBROUTINE IFFT1D(N, IN, OUT)658: SUBROUTINE IFFT1D(N, IN, OUT)
660: 659: 
661: ! calculates UNNORMALISED inverse fourier transform so,660: ! calculates UNNORMALISED inverse fourier transform so,
662: ! IN == IFFT1D(N, FFT1D(N, IN))/N661: ! IN == IFFT1D(N, FFT1D(N, IN))/N
663: 662: 
664: IMPLICIT NONE663: IMPLICIT NONE
665: 664: 
666: INTEGER(KIND=INT32), INTENT(IN) :: N665: INTEGER*4, INTENT(IN) :: N
667: COMPLEX(KIND=REAL64), INTENT(IN) :: IN(N)666: COMPLEX*16, INTENT(IN) :: IN(N)
668: COMPLEX(KIND=REAL64), INTENT(OUT) :: OUT(N)667: COMPLEX*16, INTENT(OUT) :: OUT(N)
669: 668: 
670: !INCLUDE "fftw3.f90"669: !INCLUDE "fftw3.f90"
671: INTEGER(KIND=INT64) PLAN_BACKWARD670: INTEGER*8 PLAN_BACKWARD
672: 671: 
673: CALL DFFTW_PLAN_DFT_1D_(PLAN_BACKWARD, N, IN, OUT, FFTW_BACKWARD, FFTW_ESTIMATE )672: CALL DFFTW_PLAN_DFT_1D_(PLAN_BACKWARD, N, IN, OUT, FFTW_BACKWARD, FFTW_ESTIMATE )
674: CALL DFFTW_EXECUTE_(PLAN_BACKWARD)673: CALL DFFTW_EXECUTE_(PLAN_BACKWARD)
675: CALL DFFTW_DESTROY_PLAN_(PLAN_BACKWARD)674: CALL DFFTW_DESTROY_PLAN_(PLAN_BACKWARD)
676: 675: 
677: END SUBROUTINE IFFT1D676: END SUBROUTINE IFFT1D
678: 677: 
679: SUBROUTINE ARGSORT(A,A2,ARGS,N)678: SUBROUTINE ARGSORT(A,A2,ARGS,N)
680: 679: 
681: IMPLICIT NONE680: IMPLICIT NONE
815: !  Reference:814: !  Reference:
816: !815: !
817: !    Shanjie Zhang, Jianming Jin,816: !    Shanjie Zhang, Jianming Jin,
818: !    Computation of Special Functions,817: !    Computation of Special Functions,
819: !    Wiley, 1996,818: !    Wiley, 1996,
820: !    ISBN: 0-471-11963-6,819: !    ISBN: 0-471-11963-6,
821: !    LC: QA351.C45.820: !    LC: QA351.C45.
822: !821: !
823: !  Parameters:822: !  Parameters:
824: !823: !
825: !    Input, integer ( kind = INT32 ) N, the order of the Bessel function.824: !    Input, integer ( kind = 4 ) N, the order of the Bessel function.
826: !825: !
827: !    Input, real ( kind = REAL64 ) X, the absolute value of the argument.826: !    Input, real ( kind = 8 ) X, the absolute value of the argument.
828: !827: !
829: !    Output, real ( kind = REAL64 ) ENVJ, the value.828: !    Output, real ( kind = 8 ) ENVJ, the value.
830: !829: !
831:   USE PREC, ONLY: INT32, REAL64 
832:  
833:   implicit none830:   implicit none
834: 831: 
835:   real ( kind = REAL64 ) envj832:   real ( kind = 8 ) envj
836:   real ( kind = REAL64 ) logten833:   real ( kind = 8 ) logten
837:   integer ( kind = INT32 ) n834:   integer ( kind = 4 ) n
838:   real ( kind = REAL64 ) n_r8835:   real ( kind = 8 ) n_r8
839:   real ( kind = REAL64 ) r8_gamma_log836:   real ( kind = 8 ) r8_gamma_log
840:   real ( kind = REAL64 ) x837:   real ( kind = 8 ) x
841: !838: !
842: !  Original code839: !  Original code
843: !840: !
844:   if ( .true. ) then841:   if ( .true. ) then
845: 842: 
846:     envj = 0.5D+00 * log10 ( 6.28D+00 * n ) &843:     envj = 0.5D+00 * log10 ( 6.28D+00 * n ) &
847:       - n * log10 ( 1.36D+00 * x / n )844:       - n * log10 ( 1.36D+00 * x / n )
848: !845: !
849: !  Modification suggested by Vincent Lafage.846: !  Modification suggested by Vincent Lafage.
850: !847: !
890: !  Reference:887: !  Reference:
891: !888: !
892: !    Shanjie Zhang, Jianming Jin,889: !    Shanjie Zhang, Jianming Jin,
893: !    Computation of Special Functions,890: !    Computation of Special Functions,
894: !    Wiley, 1996,891: !    Wiley, 1996,
895: !    ISBN: 0-471-11963-6,892: !    ISBN: 0-471-11963-6,
896: !    LC: QA351.C45.893: !    LC: QA351.C45.
897: !894: !
898: !  Parameters:895: !  Parameters:
899: !896: !
900: !    Input, real ( kind = REAL64 ) X, the argument.897: !    Input, real ( kind = 8 ) X, the argument.
901: !898: !
902: !    Input, integer ( kind = INT32 ) MP, the negative logarithm of the 899: !    Input, integer ( kind = 4 ) MP, the negative logarithm of the 
903: !    desired magnitude.900: !    desired magnitude.
904: !901: !
905: !    Output, integer ( kind = INT32 ) MSTA1, the starting point.902: !    Output, integer ( kind = 4 ) MSTA1, the starting point.
906: !903: !
907:   USE PREC, ONLY: INT32, REAL64 
908:  
909:   implicit none904:   implicit none
910: 905: 
911:   real ( kind = REAL64 ) a0906:   real ( kind = 8 ) a0
912:   real ( kind = REAL64 ) envj907:   real ( kind = 8 ) envj
913:   real ( kind = REAL64 ) f908:   real ( kind = 8 ) f
914:   real ( kind = REAL64 ) f0909:   real ( kind = 8 ) f0
915:   real ( kind = REAL64 ) f1910:   real ( kind = 8 ) f1
916:   integer ( kind = INT32 ) it911:   integer ( kind = 4 ) it
917:   integer ( kind = INT32 ) mp912:   integer ( kind = 4 ) mp
918:   integer ( kind = INT32 ) msta1913:   integer ( kind = 4 ) msta1
919:   integer ( kind = INT32 ) n0914:   integer ( kind = 4 ) n0
920:   integer ( kind = INT32 ) n1915:   integer ( kind = 4 ) n1
921:   integer ( kind = INT32 ) nn916:   integer ( kind = 4 ) nn
922:   real ( kind = REAL64 ) x917:   real ( kind = 8 ) x
923: 918: 
924:   a0 = abs ( x )919:   a0 = abs ( x )
925:   n0 = int ( 1.1D+00 * a0 ) + 1920:   n0 = int ( 1.1D+00 * a0 ) + 1
926:   f0 = envj ( n0, a0 ) - mp921:   f0 = envj ( n0, a0 ) - mp
927:   n1 = n0 + 5922:   n1 = n0 + 5
928:   f1 = envj ( n1, a0 ) - mp923:   f1 = envj ( n1, a0 ) - mp
929:   do it = 1, 20       924:   do it = 1, 20       
930:     nn = n1 - ( n1 - n0 ) / ( 1.0D+00 - f0 / f1 )                  925:     nn = n1 - ( n1 - n0 ) / ( 1.0D+00 - f0 / f1 )                  
931:     f = envj ( nn, a0 ) - mp926:     f = envj ( nn, a0 ) - mp
932:     if ( abs ( nn - n1 ) < 1 ) then927:     if ( abs ( nn - n1 ) < 1 ) then
973: !  Reference:968: !  Reference:
974: !969: !
975: !    Shanjie Zhang, Jianming Jin,970: !    Shanjie Zhang, Jianming Jin,
976: !    Computation of Special Functions,971: !    Computation of Special Functions,
977: !    Wiley, 1996,972: !    Wiley, 1996,
978: !    ISBN: 0-471-11963-6,973: !    ISBN: 0-471-11963-6,
979: !    LC: QA351.C45.974: !    LC: QA351.C45.
980: !975: !
981: !  Parameters:976: !  Parameters:
982: !977: !
983: !    Input, real ( kind = REAL64 ) X, the argument of Jn(x).978: !    Input, real ( kind = 8 ) X, the argument of Jn(x).
984: !979: !
985: !    Input, integer ( kind = INT32 ) N, the order of Jn(x).980: !    Input, integer ( kind = 4 ) N, the order of Jn(x).
986: !981: !
987: !    Input, integer ( kind = INT32 ) MP, the number of significant digits.982: !    Input, integer ( kind = 4 ) MP, the number of significant digits.
988: !983: !
989: !    Output, integer ( kind = INT32 ) MSTA2, the starting point.984: !    Output, integer ( kind = 4 ) MSTA2, the starting point.
990: !985: !
991:   USE PREC, ONLY: INT32, REAL64 
992:  
993:   implicit none986:   implicit none
994: 987: 
995:   real ( kind = REAL64 ) a0988:   real ( kind = 8 ) a0
996:   real ( kind = REAL64 ) ejn989:   real ( kind = 8 ) ejn
997:   real ( kind = REAL64 ) envj990:   real ( kind = 8 ) envj
998:   real ( kind = REAL64 ) f991:   real ( kind = 8 ) f
999:   real ( kind = REAL64 ) f0992:   real ( kind = 8 ) f0
1000:   real ( kind = REAL64 ) f1993:   real ( kind = 8 ) f1
1001:   real ( kind = REAL64 ) hmp994:   real ( kind = 8 ) hmp
1002:   integer ( kind = INT32 ) it995:   integer ( kind = 4 ) it
1003:   integer ( kind = INT32 ) mp996:   integer ( kind = 4 ) mp
1004:   integer ( kind = INT32 ) msta2997:   integer ( kind = 4 ) msta2
1005:   integer ( kind = INT32 ) n998:   integer ( kind = 4 ) n
1006:   integer ( kind = INT32 ) n0999:   integer ( kind = 4 ) n0
1007:   integer ( kind = INT32 ) n11000:   integer ( kind = 4 ) n1
1008:   integer ( kind = INT32 ) nn1001:   integer ( kind = 4 ) nn
1009:   real ( kind = REAL64 ) obj1002:   real ( kind = 8 ) obj
1010:   real ( kind = REAL64 ) x1003:   real ( kind = 8 ) x
1011: 1004: 
1012:   a0 = abs ( x )1005:   a0 = abs ( x )
1013:   hmp = 0.5D+00 * mp1006:   hmp = 0.5D+00 * mp
1014:   ejn = envj ( n, a0 )1007:   ejn = envj ( n, a0 )
1015: 1008: 
1016:   if ( ejn <= hmp ) then1009:   if ( ejn <= hmp ) then
1017:     obj = mp1010:     obj = mp
1018: !1011: !
1019: !  Original code:1012: !  Original code:
1020: !1013: !
1072: !  Reference:1065: !  Reference:
1073: !1066: !
1074: !    Shanjie Zhang, Jianming Jin,1067: !    Shanjie Zhang, Jianming Jin,
1075: !    Computation of Special Functions,1068: !    Computation of Special Functions,
1076: !    Wiley, 1996,1069: !    Wiley, 1996,
1077: !    ISBN: 0-471-11963-6,1070: !    ISBN: 0-471-11963-6,
1078: !    LC: QA351.C45.1071: !    LC: QA351.C45.
1079: !1072: !
1080: !  Parameters:1073: !  Parameters:
1081: !1074: !
1082: !    Input, integer ( kind = INT32 ) N, the order of In(X).1075: !    Input, integer ( kind = 4 ) N, the order of In(X).
1083: !1076: !
1084: !    Input, real ( kind = REAL64 ) X, the argument.1077: !    Input, real ( kind = 8 ) X, the argument.
1085: !1078: !
1086: !    Output, integer ( kind = INT32 ) NM, the highest order computed.1079: !    Output, integer ( kind = 4 ) NM, the highest order computed.
1087: !1080: !
1088: !    Output, real ( kind = REAL64 ) SI(0:N), DI(0:N), the values and derivatives1081: !    Output, real ( kind = 8 ) SI(0:N), DI(0:N), the values and derivatives
1089: !    of the function of orders 0 through N.1082: !    of the function of orders 0 through N.
1090: !1083: !
1091:   USE PREC, ONLY: INT32, REAL64 
1092:  
1093:   implicit none1084:   implicit none
1094: 1085: 
1095:   integer ( kind = INT32 ), intent(in) :: n1086:   integer ( kind = 4 ), intent(in) :: n
1096: 1087: 
1097:   real ( kind = REAL64 ) cs1088:   real ( kind = 8 ) cs
1098:   real ( kind = REAL64 ) f1089:   real ( kind = 8 ) f
1099:   real ( kind = REAL64 ) f01090:   real ( kind = 8 ) f0
1100:   real ( kind = REAL64 ) f11091:   real ( kind = 8 ) f1
1101:   integer ( kind = INT32 ) k1092:   integer ( kind = 4 ) k
1102:   integer ( kind = INT32 ) m1093:   integer ( kind = 4 ) m
1103:   integer ( kind = INT32 ) msta11094:   integer ( kind = 4 ) msta1
1104:   integer ( kind = INT32 ) msta21095:   integer ( kind = 4 ) msta2
1105:   integer ( kind = INT32 ), intent(out) :: nm1096:   integer ( kind = 4 ), intent(out) :: nm
1106:   real ( kind = REAL64 ), intent(out) :: si(0:n)1097:   real ( kind = 8 ), intent(out) :: si(0:n)
1107:   real ( kind = REAL64 ) si01098:   real ( kind = 8 ) si0
1108:   real ( kind = REAL64 ), intent(in) :: x1099:   real ( kind = 8 ), intent(in) :: x
1109: 1100: 
1110:   nm = n1101:   nm = n
1111: 1102: 
1112:   if ( abs ( x ) < 1.0D-100 ) then1103:   if ( abs ( x ) < 1.0D-100 ) then
1113:     do k = 0, n1104:     do k = 0, n
1114:       si(k) = 0.0D+001105:       si(k) = 0.0D+00
1115:     end do1106:     end do
1116:     si(0) = 1.0D+001107:     si(0) = 1.0D+00
1117:     return1108:     return
1118:   end if1109:   end if
1173: !  Reference:1164: !  Reference:
1174: !1165: !
1175: !    Shanjie Zhang, Jianming Jin,1166: !    Shanjie Zhang, Jianming Jin,
1176: !    Computation of Special Functions,1167: !    Computation of Special Functions,
1177: !    Wiley, 1996,1168: !    Wiley, 1996,
1178: !    ISBN: 0-471-11963-6,1169: !    ISBN: 0-471-11963-6,
1179: !    LC: QA351.C45.1170: !    LC: QA351.C45.
1180: !1171: !
1181: !  Parameters:1172: !  Parameters:
1182: !1173: !
1183: !    Input, real ( kind = REAL64 ) A, B, parameters.1174: !    Input, real ( kind = 8 ) A, B, parameters.
1184: !1175: !
1185: !    Input, real ( kind = REAL64 ) X, the argument.1176: !    Input, real ( kind = 8 ) X, the argument.
1186: !1177: !
1187: !    Output, real ( kind = REAL64 ) HG, the value of M(a,b,x).1178: !    Output, real ( kind = 8 ) HG, the value of M(a,b,x).
1188: !1179: !
1189:   USE PREC, ONLY: INT32, REAL64 
1190:  
1191:   implicit none1180:   implicit none
1192: 1181: 
1193:   real ( kind = REAL64 ), intent(in) :: ain1182:   real ( kind = 8 ), intent(in) :: ain
1194:   real ( kind = REAL64 ), intent(in) :: bin1183:   real ( kind = 8 ), intent(in) :: bin
1195:   real ( kind = REAL64 ), intent(in) :: xin1184:   real ( kind = 8 ), intent(in) :: xin
1196:   real ( kind = REAL64 ), intent(out) :: hg1185:   real ( kind = 8 ), intent(out) :: hg
1197: 1186: 
1198:   real ( kind = REAL64 ) a1187:   real ( kind = 8 ) a
1199:   real ( kind = REAL64 ) b1188:   real ( kind = 8 ) b
1200:   real ( kind = REAL64 ) x1189:   real ( kind = 8 ) x
1201: 1190: 
1202:   real ( kind = REAL64 ) a01191:   real ( kind = 8 ) a0
1203:   real ( kind = REAL64 ) a11192:   real ( kind = 8 ) a1
1204:   real ( kind = REAL64 ) aa1193:   real ( kind = 8 ) aa
1205: 1194: 
1206:   real ( kind = REAL64 ) hg11195:   real ( kind = 8 ) hg1
1207:   real ( kind = REAL64 ) hg21196:   real ( kind = 8 ) hg2
1208:   integer ( kind = INT32 ) i1197:   integer ( kind = 4 ) i
1209:   integer ( kind = INT32 ) j1198:   integer ( kind = 4 ) j
1210:   integer ( kind = INT32 ) k1199:   integer ( kind = 4 ) k
1211:   integer ( kind = INT32 ) la1200:   integer ( kind = 4 ) la
1212:   integer ( kind = INT32 ) m1201:   integer ( kind = 4 ) m
1213:   integer ( kind = INT32 ) n1202:   integer ( kind = 4 ) n
1214:   integer ( kind = INT32 ) nl1203:   integer ( kind = 4 ) nl
1215:   real ( kind = REAL64 ) pi1204:   real ( kind = 8 ) pi
1216:   real ( kind = REAL64 ) r1205:   real ( kind = 8 ) r
1217:   real ( kind = REAL64 ) r11206:   real ( kind = 8 ) r1
1218:   real ( kind = REAL64 ) r21207:   real ( kind = 8 ) r2
1219:   real ( kind = REAL64 ) rg1208:   real ( kind = 8 ) rg
1220:   real ( kind = REAL64 ) sum11209:   real ( kind = 8 ) sum1
1221:   real ( kind = REAL64 ) sum21210:   real ( kind = 8 ) sum2
1222:   real ( kind = REAL64 ) ta1211:   real ( kind = 8 ) ta
1223:   real ( kind = REAL64 ) tb1212:   real ( kind = 8 ) tb
1224:   real ( kind = REAL64 ) tba1213:   real ( kind = 8 ) tba
1225:   real ( kind = REAL64 ) x01214:   real ( kind = 8 ) x0
1226:   real ( kind = REAL64 ) xg1215:   real ( kind = 8 ) xg
1227:   real ( kind = REAL64 ) y01216:   real ( kind = 8 ) y0
1228:   real ( kind = REAL64 ) y11217:   real ( kind = 8 ) y1
1229: 1218: 
1230:   a=ain1219:   a=ain
1231:   b=bin1220:   b=bin
1232:   x=xin1221:   x=xin
1233:   pi = 3.141592653589793D+001222:   pi = 3.141592653589793D+00
1234:   a0 = a1223:   a0 = a
1235:   a1 = a1224:   a1 = a
1236:   x0 = x1225:   x0 = x
1237:   hg = 0.0D+001226:   hg = 0.0D+00
1238: 1227: 
1373: !  Reference:1362: !  Reference:
1374: !1363: !
1375: !    Shanjie Zhang, Jianming Jin,1364: !    Shanjie Zhang, Jianming Jin,
1376: !    Computation of Special Functions,1365: !    Computation of Special Functions,
1377: !    Wiley, 1996,1366: !    Wiley, 1996,
1378: !    ISBN: 0-471-11963-6,1367: !    ISBN: 0-471-11963-6,
1379: !    LC: QA351.C451368: !    LC: QA351.C45
1380: !1369: !
1381: !  Parameters:1370: !  Parameters:
1382: !1371: !
1383: !    Input, real ( kind = REAL64 ) X, the argument.1372: !    Input, real ( kind = 8 ) X, the argument.
1384: !    X must not be 0, or any negative integer.1373: !    X must not be 0, or any negative integer.
1385: !1374: !
1386: !    Output, real ( kind = REAL64 ) GA, the value of the Gamma function.1375: !    Output, real ( kind = 8 ) GA, the value of the Gamma function.
1387: !1376: !
1388:   USE PREC, ONLY: INT32, REAL64 
1389:  
1390:   implicit none1377:   implicit none
1391: 1378: 
1392:   real ( kind = REAL64 ), intent(in) :: x1379:   real ( kind = 8 ), intent(in) :: x
1393:   real ( kind = REAL64 ), intent(out) :: ga1380:   real ( kind = 8 ), intent(out) :: ga
1394: 1381: 
1395:   real ( kind = REAL64 ), dimension ( 26 ) :: g = (/ &1382:   real ( kind = 8 ), dimension ( 26 ) :: g = (/ &
1396:     1.0D+00, &1383:     1.0D+00, &
1397:     0.5772156649015329D+00, &1384:     0.5772156649015329D+00, &
1398:    -0.6558780715202538D+00, &1385:    -0.6558780715202538D+00, &
1399:    -0.420026350340952D-01, &1386:    -0.420026350340952D-01, &
1400:     0.1665386113822915D+00, &1387:     0.1665386113822915D+00, &
1401:    -0.421977345555443D-01, &1388:    -0.421977345555443D-01, &
1402:    -0.96219715278770D-02, &1389:    -0.96219715278770D-02, &
1403:     0.72189432466630D-02, &1390:     0.72189432466630D-02, &
1404:    -0.11651675918591D-02, &1391:    -0.11651675918591D-02, &
1405:    -0.2152416741149D-03, &1392:    -0.2152416741149D-03, &
1413:    -0.11812746D-08, &1400:    -0.11812746D-08, &
1414:     0.1043427D-09, & 1401:     0.1043427D-09, & 
1415:     0.77823D-11, &1402:     0.77823D-11, &
1416:    -0.36968D-11, &1403:    -0.36968D-11, &
1417:     0.51D-12, &1404:     0.51D-12, &
1418:    -0.206D-13, &1405:    -0.206D-13, &
1419:    -0.54D-14, &1406:    -0.54D-14, &
1420:     0.14D-14, &1407:     0.14D-14, &
1421:     0.1D-15 /)1408:     0.1D-15 /)
1422: 1409: 
1423:   real ( kind = REAL64 ) gr1410:   real ( kind = 8 ) gr
1424:   integer ( kind = INT32 ) k1411:   integer ( kind = 4 ) k
1425:   integer ( kind = INT32 ) m1412:   integer ( kind = 4 ) m
1426:   integer ( kind = INT32 ) m11413:   integer ( kind = 4 ) m1
1427:   real ( kind = REAL64 ), parameter :: pi = 3.141592653589793D+001414:   real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00
1428:   real ( kind = REAL64 ) r1415:   real ( kind = 8 ) r
1429:   real ( kind = REAL64 ) z1416:   real ( kind = 8 ) z
1430: 1417: 
1431:   if ( x == aint ( x ) ) then1418:   if ( x == aint ( x ) ) then
1432: 1419: 
1433:     if ( 0.0D+00 < x ) then1420:     if ( 0.0D+00 < x ) then
1434:       ga = 1.0D+001421:       ga = 1.0D+00
1435:       m1 = int ( x ) - 11422:       m1 = int ( x ) - 1
1436:       do k = 2, m11423:       do k = 2, m1
1437:         ga = ga * k1424:         ga = ga * k
1438:       end do1425:       end do
1439:     else1426:     else
1526: !    Jorge More, Burton Garbow, Kenneth Hillstrom,1513: !    Jorge More, Burton Garbow, Kenneth Hillstrom,
1527: !    User Guide for MINPACK-1,1514: !    User Guide for MINPACK-1,
1528: !    Technical Report ANL-80-74,1515: !    Technical Report ANL-80-74,
1529: !    Argonne National Laboratory, 1980.1516: !    Argonne National Laboratory, 1980.
1530: !1517: !
1531: !  Parameters:1518: !  Parameters:
1532: !1519: !
1533: !    Input, external FCN, the name of the user-supplied subroutine which1520: !    Input, external FCN, the name of the user-supplied subroutine which
1534: !    calculates the functions and the jacobian.  FCN should have the form:1521: !    calculates the functions and the jacobian.  FCN should have the form:
1535: !      subroutine fcn ( m, n, x, fvec, fjac, ldfjac, iflag )1522: !      subroutine fcn ( m, n, x, fvec, fjac, ldfjac, iflag )
1536: !      integer ( kind = INT32 ) ldfjac1523: !      integer ( kind = 4 ) ldfjac
1537: !      integer ( kind = INT32 ) n1524: !      integer ( kind = 4 ) n
1538: !      real ( kind = REAL64 ) fjac(ldfjac,n)1525: !      real ( kind = 8 ) fjac(ldfjac,n)
1539: !      real ( kind = REAL64 ) fvec(m)1526: !      real ( kind = 8 ) fvec(m)
1540: !      integer ( kind = INT32 ) iflag1527: !      integer ( kind = 4 ) iflag
1541: !      real ( kind = REAL64 ) x(n)1528: !      real ( kind = 8 ) x(n)
1542: !1529: !
1543: !    If IFLAG = 0 on input, then FCN is only being called to allow the user1530: !    If IFLAG = 0 on input, then FCN is only being called to allow the user
1544: !    to print out the current iterate.1531: !    to print out the current iterate.
1545: !    If IFLAG = 1 on input, FCN should calculate the functions at X and1532: !    If IFLAG = 1 on input, FCN should calculate the functions at X and
1546: !    return this vector in FVEC.1533: !    return this vector in FVEC.
1547: !    If IFLAG = 2 on input, FCN should calculate the jacobian at X and1534: !    If IFLAG = 2 on input, FCN should calculate the jacobian at X and
1548: !    return this matrix in FJAC.1535: !    return this matrix in FJAC.
1549: !    To terminate the algorithm, FCN may set IFLAG negative on return.1536: !    To terminate the algorithm, FCN may set IFLAG negative on return.
1550: !1537: !
1551: !    Input, integer ( kind = INT32 ) M, is the number of functions.1538: !    Input, integer ( kind = 4 ) M, is the number of functions.
1552: !1539: !
1553: !    Input, integer ( kind = INT32 ) N, is the number of variables.  1540: !    Input, integer ( kind = 4 ) N, is the number of variables.  
1554: !    N must not exceed M.1541: !    N must not exceed M.
1555: !1542: !
1556: !    Input/output, real ( kind = REAL64 ) X(N).  On input, X must contain an initial1543: !    Input/output, real ( kind = 8 ) X(N).  On input, X must contain an initial
1557: !    estimate of the solution vector.  On output X contains the final1544: !    estimate of the solution vector.  On output X contains the final
1558: !    estimate of the solution vector.1545: !    estimate of the solution vector.
1559: !1546: !
1560: !    Output, real ( kind = REAL64 ) FVEC(M), the functions evaluated at the output X.1547: !    Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X.
1561: !1548: !
1562: !    Output, real ( kind = REAL64 ) FJAC(LDFJAC,N), an M by N array.  The upper1549: !    Output, real ( kind = 8 ) FJAC(LDFJAC,N), an M by N array.  The upper
1563: !    N by N submatrix of FJAC contains an upper triangular matrix R with1550: !    N by N submatrix of FJAC contains an upper triangular matrix R with
1564: !    diagonal elements of nonincreasing magnitude such that1551: !    diagonal elements of nonincreasing magnitude such that
1565: !      P' * ( JAC' * JAC ) * P = R' * R,1552: !      P' * ( JAC' * JAC ) * P = R' * R,
1566: !    where P is a permutation matrix and JAC is the final calculated jacobian.1553: !    where P is a permutation matrix and JAC is the final calculated jacobian.
1567: !    Column J of P is column IPVT(J) of the identity matrix.  The lower1554: !    Column J of P is column IPVT(J) of the identity matrix.  The lower
1568: !    trapezoidal part of FJAC contains information generated during1555: !    trapezoidal part of FJAC contains information generated during
1569: !    the computation of R.1556: !    the computation of R.
1570: !1557: !
1571: !    Input, integer ( kind = INT32 ) LDFJAC, the leading dimension of FJAC.1558: !    Input, integer ( kind = 4 ) LDFJAC, the leading dimension of FJAC.
1572: !    LDFJAC must be at least M.1559: !    LDFJAC must be at least M.
1573: !1560: !
1574: !    Input, real ( kind = REAL64 ) FTOL.  Termination occurs when both the actual1561: !    Input, real ( kind = 8 ) FTOL.  Termination occurs when both the actual
1575: !    and predicted relative reductions in the sum of squares are at most FTOL.1562: !    and predicted relative reductions in the sum of squares are at most FTOL.
1576: !    Therefore, FTOL measures the relative error desired in the sum of1563: !    Therefore, FTOL measures the relative error desired in the sum of
1577: !    squares.  FTOL should be nonnegative.1564: !    squares.  FTOL should be nonnegative.
1578: !1565: !
1579: !    Input, real ( kind = REAL64 ) XTOL.  Termination occurs when the relative error1566: !    Input, real ( kind = 8 ) XTOL.  Termination occurs when the relative error
1580: !    between two consecutive iterates is at most XTOL.  XTOL should be1567: !    between two consecutive iterates is at most XTOL.  XTOL should be
1581: !    nonnegative.1568: !    nonnegative.
1582: !1569: !
1583: !    Input, real ( kind = REAL64 ) GTOL.  Termination occurs when the cosine of the1570: !    Input, real ( kind = 8 ) GTOL.  Termination occurs when the cosine of the
1584: !    angle between FVEC and any column of the jacobian is at most GTOL in1571: !    angle between FVEC and any column of the jacobian is at most GTOL in
1585: !    absolute value.  Therefore, GTOL measures the orthogonality desired1572: !    absolute value.  Therefore, GTOL measures the orthogonality desired
1586: !    between the function vector and the columns of the jacobian.  GTOL should1573: !    between the function vector and the columns of the jacobian.  GTOL should
1587: !    be nonnegative.1574: !    be nonnegative.
1588: !1575: !
1589: !    Input, integer ( kind = INT32 ) MAXFEV.  Termination occurs when the number of1576: !    Input, integer ( kind = 4 ) MAXFEV.  Termination occurs when the number of
1590: !    calls to FCN with IFLAG = 1 is at least MAXFEV by the end of an iteration.1577: !    calls to FCN with IFLAG = 1 is at least MAXFEV by the end of an iteration.
1591: !1578: !
1592: !    Input/output, real ( kind = REAL64 ) DIAG(N).  If MODE = 1, then DIAG is set1579: !    Input/output, real ( kind = 8 ) DIAG(N).  If MODE = 1, then DIAG is set
1593: !    internally.  If MODE = 2, then DIAG must contain positive entries that1580: !    internally.  If MODE = 2, then DIAG must contain positive entries that
1594: !    serve as multiplicative scale factors for the variables.1581: !    serve as multiplicative scale factors for the variables.
1595: !1582: !
1596: !    Input, integer ( kind = INT32 ) MODE, scaling option.1583: !    Input, integer ( kind = 4 ) MODE, scaling option.
1597: !    1, variables will be scaled internally.1584: !    1, variables will be scaled internally.
1598: !    2, scaling is specified by the input DIAG vector.1585: !    2, scaling is specified by the input DIAG vector.
1599: !1586: !
1600: !    Input, real ( kind = REAL64 ) FACTOR, determines the initial step bound.  This1587: !    Input, real ( kind = 8 ) FACTOR, determines the initial step bound.  This
1601: !    bound is set to the product of FACTOR and the euclidean norm of DIAG*X if1588: !    bound is set to the product of FACTOR and the euclidean norm of DIAG*X if
1602: !    nonzero, or else to FACTOR itself.  In most cases, FACTOR should lie1589: !    nonzero, or else to FACTOR itself.  In most cases, FACTOR should lie
1603: !    in the interval (0.1, 100) with 100 the recommended value.1590: !    in the interval (0.1, 100) with 100 the recommended value.
1604: !1591: !
1605: !    Input, integer ( kind = INT32 ) NPRINT, enables controlled printing of iterates1592: !    Input, integer ( kind = 4 ) NPRINT, enables controlled printing of iterates
1606: !    if it is positive.  In this case, FCN is called with IFLAG = 0 at the1593: !    if it is positive.  In this case, FCN is called with IFLAG = 0 at the
1607: !    beginning of the first iteration and every NPRINT iterations thereafter1594: !    beginning of the first iteration and every NPRINT iterations thereafter
1608: !    and immediately prior to return, with X and FVEC available1595: !    and immediately prior to return, with X and FVEC available
1609: !    for printing.  If NPRINT is not positive, no special calls1596: !    for printing.  If NPRINT is not positive, no special calls
1610: !    of FCN with IFLAG = 0 are made.1597: !    of FCN with IFLAG = 0 are made.
1611: !1598: !
1612: !    Output, integer ( kind = INT32 ) INFO, error flag.  If the user has terminated1599: !    Output, integer ( kind = 4 ) INFO, error flag.  If the user has terminated
1613: !    execution, INFO is set to the (negative) value of IFLAG. See description1600: !    execution, INFO is set to the (negative) value of IFLAG. See description
1614: !    of FCN.  Otherwise, INFO is set as follows:1601: !    of FCN.  Otherwise, INFO is set as follows:
1615: !    0, improper input parameters.1602: !    0, improper input parameters.
1616: !    1, both actual and predicted relative reductions in the sum of1603: !    1, both actual and predicted relative reductions in the sum of
1617: !       squares are at most FTOL.1604: !       squares are at most FTOL.
1618: !    2, relative error between two consecutive iterates is at most XTOL.1605: !    2, relative error between two consecutive iterates is at most XTOL.
1619: !    3, conditions for INFO = 1 and INFO = 2 both hold.1606: !    3, conditions for INFO = 1 and INFO = 2 both hold.
1620: !    4, the cosine of the angle between FVEC and any column of the jacobian1607: !    4, the cosine of the angle between FVEC and any column of the jacobian
1621: !       is at most GTOL in absolute value.1608: !       is at most GTOL in absolute value.
1622: !    5, number of calls to FCN with IFLAG = 1 has reached MAXFEV.1609: !    5, number of calls to FCN with IFLAG = 1 has reached MAXFEV.
1623: !    6, FTOL is too small.  No further reduction in the sum of squares1610: !    6, FTOL is too small.  No further reduction in the sum of squares
1624: !       is possible.1611: !       is possible.
1625: !    7, XTOL is too small.  No further improvement in the approximate1612: !    7, XTOL is too small.  No further improvement in the approximate
1626: !       solution X is possible.1613: !       solution X is possible.
1627: !    8, GTOL is too small.  FVEC is orthogonal to the columns of the1614: !    8, GTOL is too small.  FVEC is orthogonal to the columns of the
1628: !       jacobian to machine precision.1615: !       jacobian to machine precision.
1629: !1616: !
1630: !    Output, integer ( kind = INT32 ) NFEV, the number of calls to FCN with1617: !    Output, integer ( kind = 4 ) NFEV, the number of calls to FCN with
1631: !    IFLAG = 1.1618: !    IFLAG = 1.
1632: !1619: !
1633: !    Output, integer ( kind = INT32 ) NJEV, the number of calls to FCN with1620: !    Output, integer ( kind = 4 ) NJEV, the number of calls to FCN with
1634: !    IFLAG = 2.1621: !    IFLAG = 2.
1635: !1622: !
1636: !    Output, integer ( kind = INT32 ) IPVT(N), defines a permutation matrix P1623: !    Output, integer ( kind = 4 ) IPVT(N), defines a permutation matrix P
1637: !    such that JAC*P = Q*R, where JAC is the final calculated jacobian, Q is1624: !    such that JAC*P = Q*R, where JAC is the final calculated jacobian, Q is
1638: !    orthogonal (not stored), and R is upper triangular with diagonal1625: !    orthogonal (not stored), and R is upper triangular with diagonal
1639: !    elements of nonincreasing magnitude.  Column J of P is column1626: !    elements of nonincreasing magnitude.  Column J of P is column
1640: !    IPVT(J) of the identity matrix.1627: !    IPVT(J) of the identity matrix.
1641: !1628: !
1642: !    Output, real ( kind = REAL64 ) QTF(N), contains the first N elements of Q'*FVEC.1629: !    Output, real ( kind = 8 ) QTF(N), contains the first N elements of Q'*FVEC.
1643: !1630: !
1644:   USE PREC, ONLY: INT32, REAL64 
1645:  
1646:   implicit none1631:   implicit none
1647: 1632: 
1648:   integer ( kind = INT32 ), INTENT(IN) :: ldfjac1633:   integer ( kind = 4 ), INTENT(IN) :: ldfjac
1649:   integer ( kind = INT32 ), INTENT(IN) ::  m1634:   integer ( kind = 4 ), INTENT(IN) ::  m
1650:   integer ( kind = INT32 ), INTENT(IN) ::  n1635:   integer ( kind = 4 ), INTENT(IN) ::  n
1651: 1636: 
1652:   real ( kind = REAL64 ) actred1637:   real ( kind = 8 ) actred
1653:   real ( kind = REAL64 ) delta1638:   real ( kind = 8 ) delta
1654:   real ( kind = REAL64 ), INTENT(INOUT) :: diag(n)1639:   real ( kind = 8 ), INTENT(INOUT) :: diag(n)
1655:   real ( kind = REAL64 ) dirder1640:   real ( kind = 8 ) dirder
1656:   real ( kind = REAL64 ) enorm1641:   real ( kind = 8 ) enorm
1657:   real ( kind = REAL64 ) epsmch1642:   real ( kind = 8 ) epsmch
1658:   real ( kind = REAL64 ), INTENT(IN) :: factor1643:   real ( kind = 8 ), INTENT(IN) :: factor
1659:   external  fcn1644:   external  fcn
1660:   real ( kind = REAL64 ), INTENT(OUT) :: fjac(ldfjac,n)1645:   real ( kind = 8 ), INTENT(OUT) :: fjac(ldfjac,n)
1661:   real ( kind = REAL64 ) fnorm1646:   real ( kind = 8 ) fnorm
1662:   real ( kind = REAL64 ) fnorm11647:   real ( kind = 8 ) fnorm1
1663:   real ( kind = REAL64 ), INTENT(IN) :: ftol1648:   real ( kind = 8 ), INTENT(IN) :: ftol
1664:   real ( kind = REAL64 ), INTENT(OUT) :: fvec(m)1649:   real ( kind = 8 ), INTENT(OUT) :: fvec(m)
1665:   real ( kind = REAL64 ) gnorm1650:   real ( kind = 8 ) gnorm
1666:   real ( kind = REAL64 ), INTENT(IN) :: gtol1651:   real ( kind = 8 ), INTENT(IN) :: gtol
1667:   integer ( kind = INT32 ) i1652:   integer ( kind = 4 ) i
1668:   integer ( kind = INT32 ) iflag1653:   integer ( kind = 4 ) iflag
1669:   integer ( kind = INT32 ), INTENT(OUT) :: info1654:   integer ( kind = 4 ), INTENT(OUT) :: info
1670:   integer ( kind = INT32 ) ipvt(n)1655:   integer ( kind = 4 ) ipvt(n)
1671:   integer ( kind = INT32 ) iter1656:   integer ( kind = 4 ) iter
1672:   integer ( kind = INT32 ) j1657:   integer ( kind = 4 ) j
1673:   integer ( kind = INT32 ) l1658:   integer ( kind = 4 ) l
1674:   integer ( kind = INT32 ), INTENT(IN) :: maxfev1659:   integer ( kind = 4 ), INTENT(IN) :: maxfev
1675:   integer ( kind = INT32 ), INTENT(IN) :: mode1660:   integer ( kind = 4 ), INTENT(IN) :: mode
1676:   integer ( kind = INT32 ), INTENT(OUT) :: nfev1661:   integer ( kind = 4 ), INTENT(OUT) :: nfev
1677:   integer ( kind = INT32 ), INTENT(OUT) :: njev1662:   integer ( kind = 4 ), INTENT(OUT) :: njev
1678:   integer ( kind = INT32 ), INTENT(IN) :: nprint1663:   integer ( kind = 4 ), INTENT(IN) :: nprint
1679:   real ( kind = REAL64 ) par1664:   real ( kind = 8 ) par
1680:   logical pivot1665:   logical pivot
1681:   real ( kind = REAL64 ) pnorm1666:   real ( kind = 8 ) pnorm
1682:   real ( kind = REAL64 ) prered1667:   real ( kind = 8 ) prered
1683:   real ( kind = REAL64 ), INTENT(OUT) :: qtf(n)1668:   real ( kind = 8 ), INTENT(OUT) :: qtf(n)
1684:   real ( kind = REAL64 ) ratio1669:   real ( kind = 8 ) ratio
1685:   real ( kind = REAL64 ) sum21670:   real ( kind = 8 ) sum2
1686:   real ( kind = REAL64 ) temp1671:   real ( kind = 8 ) temp
1687:   real ( kind = REAL64 ) temp11672:   real ( kind = 8 ) temp1
1688:   real ( kind = REAL64 ) temp21673:   real ( kind = 8 ) temp2
1689:   real ( kind = REAL64 ) wa1(n)1674:   real ( kind = 8 ) wa1(n)
1690:   real ( kind = REAL64 ) wa2(n)1675:   real ( kind = 8 ) wa2(n)
1691:   real ( kind = REAL64 ) wa3(n)1676:   real ( kind = 8 ) wa3(n)
1692:   real ( kind = REAL64 ) wa4(m)1677:   real ( kind = 8 ) wa4(m)
1693:   real ( kind = REAL64 ) xnorm1678:   real ( kind = 8 ) xnorm
1694:   real ( kind = REAL64 ), INTENT(INOUT) ::  x(n)1679:   real ( kind = 8 ), INTENT(INOUT) ::  x(n)
1695:   real ( kind = REAL64 ), INTENT(IN) :: xtol1680:   real ( kind = 8 ), INTENT(IN) :: xtol
1696: 1681: 
1697:   epsmch = epsilon ( epsmch )1682:   epsmch = epsilon ( epsmch )
1698: 1683: 
1699:   info = 01684:   info = 0
1700:   iflag = 01685:   iflag = 0
1701:   nfev = 01686:   nfev = 0
1702:   njev = 01687:   njev = 0
1703: !1688: !
1704: !  Check the input parameters for errors.1689: !  Check the input parameters for errors.
1705: !1690: !
2057: !    Jorge More, Burton Garbow, Kenneth Hillstrom,2042: !    Jorge More, Burton Garbow, Kenneth Hillstrom,
2058: !    User Guide for MINPACK-1,2043: !    User Guide for MINPACK-1,
2059: !    Technical Report ANL-80-74,2044: !    Technical Report ANL-80-74,
2060: !    Argonne National Laboratory, 1980.2045: !    Argonne National Laboratory, 1980.
2061: !2046: !
2062: !  Parameters:2047: !  Parameters:
2063: !2048: !
2064: !    Input, external FCN, the name of the user-supplied subroutine which2049: !    Input, external FCN, the name of the user-supplied subroutine which
2065: !    calculates the functions and the jacobian.  FCN should have the form:2050: !    calculates the functions and the jacobian.  FCN should have the form:
2066: !      subroutine fcn ( m, n, x, fvec, fjac, ldfjac, iflag )2051: !      subroutine fcn ( m, n, x, fvec, fjac, ldfjac, iflag )
2067: !      integer ( kind = INT32 ) ldfjac2052: !      integer ( kind = 4 ) ldfjac
2068: !      integer ( kind = INT32 ) n2053: !      integer ( kind = 4 ) n
2069: !      real ( kind = REAL64 ) fjac(ldfjac,n)2054: !      real ( kind = 8 ) fjac(ldfjac,n)
2070: !      real ( kind = REAL64 ) fvec(m)2055: !      real ( kind = 8 ) fvec(m)
2071: !      integer ( kind = INT32 ) iflag2056: !      integer ( kind = 4 ) iflag
2072: !      real ( kind = REAL64 ) x(n)2057: !      real ( kind = 8 ) x(n)
2073: !2058: !
2074: !    If IFLAG = 0 on input, then FCN is only being called to allow the user2059: !    If IFLAG = 0 on input, then FCN is only being called to allow the user
2075: !    to print out the current iterate.2060: !    to print out the current iterate.
2076: !    If IFLAG = 1 on input, FCN should calculate the functions at X and2061: !    If IFLAG = 1 on input, FCN should calculate the functions at X and
2077: !    return this vector in FVEC.2062: !    return this vector in FVEC.
2078: !    If IFLAG = 2 on input, FCN should calculate the jacobian at X and2063: !    If IFLAG = 2 on input, FCN should calculate the jacobian at X and
2079: !    return this matrix in FJAC.2064: !    return this matrix in FJAC.
2080: !    To terminate the algorithm, FCN may set IFLAG negative on return.2065: !    To terminate the algorithm, FCN may set IFLAG negative on return.
2081: !2066: !
2082: !    Input, integer ( kind = INT32 ) M, the number of functions.2067: !    Input, integer ( kind = 4 ) M, the number of functions.
2083: !2068: !
2084: !    Input, integer ( kind = INT32 ) N, is the number of variables.  2069: !    Input, integer ( kind = 4 ) N, is the number of variables.  
2085: !    N must not exceed M.2070: !    N must not exceed M.
2086: !2071: !
2087: !    Input/output, real ( kind = REAL64 ) X(N).  On input, X must contain an initial2072: !    Input/output, real ( kind = 8 ) X(N).  On input, X must contain an initial
2088: !    estimate of the solution vector.  On output X contains the final2073: !    estimate of the solution vector.  On output X contains the final
2089: !    estimate of the solution vector.2074: !    estimate of the solution vector.
2090: !2075: !
2091: !    Output, real ( kind = REAL64 ) FVEC(M), the functions evaluated at the output X.2076: !    Output, real ( kind = 8 ) FVEC(M), the functions evaluated at the output X.
2092: !2077: !
2093: !    Output, real ( kind = REAL64 ) FJAC(LDFJAC,N), an M by N array.  The upper2078: !    Output, real ( kind = 8 ) FJAC(LDFJAC,N), an M by N array.  The upper
2094: !    N by N submatrix contains an upper triangular matrix R with2079: !    N by N submatrix contains an upper triangular matrix R with
2095: !    diagonal elements of nonincreasing magnitude such that2080: !    diagonal elements of nonincreasing magnitude such that
2096: !      P' * ( JAC' * JAC ) * P = R' * R,2081: !      P' * ( JAC' * JAC ) * P = R' * R,
2097: !    where P is a permutation matrix and JAC is the final calculated2082: !    where P is a permutation matrix and JAC is the final calculated
2098: !    jacobian.  Column J of P is column IPVT(J) of the identity matrix.2083: !    jacobian.  Column J of P is column IPVT(J) of the identity matrix.
2099: !    The lower trapezoidal part of FJAC contains information generated during2084: !    The lower trapezoidal part of FJAC contains information generated during
2100: !    the computation of R.2085: !    the computation of R.
2101: !2086: !
2102: !    Input, integer ( kind = INT32 ) LDFJAC, is the leading dimension of FJAC,2087: !    Input, integer ( kind = 4 ) LDFJAC, is the leading dimension of FJAC,
2103: !    which must be no less than M.2088: !    which must be no less than M.
2104: !2089: !
2105: !    Input, real ( kind = REAL64 ) TOL.  Termination occurs when the algorithm2090: !    Input, real ( kind = 8 ) TOL.  Termination occurs when the algorithm
2106: !    estimates either that the relative error in the sum of squares is at2091: !    estimates either that the relative error in the sum of squares is at
2107: !    most TOL or that the relative error between X and the solution is at2092: !    most TOL or that the relative error between X and the solution is at
2108: !    most TOL.2093: !    most TOL.
2109: !2094: !
2110: !    Output, integer ( kind = INT32 ) INFO, error flag.  If the user has terminated2095: !    Output, integer ( kind = 4 ) INFO, error flag.  If the user has terminated
2111: !    execution, INFO is set to the (negative) value of IFLAG. See description2096: !    execution, INFO is set to the (negative) value of IFLAG. See description
2112: !    of FCN.  Otherwise, INFO is set as follows:2097: !    of FCN.  Otherwise, INFO is set as follows:
2113: !    0, improper input parameters.2098: !    0, improper input parameters.
2114: !    1, algorithm estimates that the relative error in the sum of squares2099: !    1, algorithm estimates that the relative error in the sum of squares
2115: !       is at most TOL.2100: !       is at most TOL.
2116: !    2, algorithm estimates that the relative error between X and the2101: !    2, algorithm estimates that the relative error between X and the
2117: !       solution is at most TOL.2102: !       solution is at most TOL.
2118: !    3, conditions for INFO = 1 and INFO = 2 both hold.2103: !    3, conditions for INFO = 1 and INFO = 2 both hold.
2119: !    4, FVEC is orthogonal to the columns of the jacobian to machine precision.2104: !    4, FVEC is orthogonal to the columns of the jacobian to machine precision.
2120: !    5, number of calls to FCN with IFLAG = 1 has reached 100*(N+1).2105: !    5, number of calls to FCN with IFLAG = 1 has reached 100*(N+1).
2121: !    6, TOL is too small.  No further reduction in the sum of squares is2106: !    6, TOL is too small.  No further reduction in the sum of squares is
2122: !       possible.2107: !       possible.
2123: !    7, TOL is too small.  No further improvement in the approximate2108: !    7, TOL is too small.  No further improvement in the approximate
2124: !       solution X is possible.2109: !       solution X is possible.
2125: !2110: !
2126:   USE PREC, ONLY: INT32, REAL64 
2127:  
2128:   implicit none2111:   implicit none
2129: 2112: 
2130:   integer ( kind = INT32 ), INTENT(IN) ::  ldfjac2113:   integer ( kind = 4 ), INTENT(IN) ::  ldfjac
2131:   integer ( kind = INT32 ), INTENT(IN) ::  m2114:   integer ( kind = 4 ), INTENT(IN) ::  m
2132:   integer ( kind = INT32 ), INTENT(IN) ::  n2115:   integer ( kind = 4 ), INTENT(IN) ::  n
2133: 2116: 
2134:   real ( kind = REAL64 ) diag(n)2117:   real ( kind = 8 ) diag(n)
2135:   real ( kind = REAL64 ) factor2118:   real ( kind = 8 ) factor
2136:   external fcn2119:   external fcn
2137:   real ( kind = REAL64 ), INTENT(OUT) ::  fjac(ldfjac,n)2120:   real ( kind = 8 ), INTENT(OUT) ::  fjac(ldfjac,n)
2138:   real ( kind = REAL64 ) ftol2121:   real ( kind = 8 ) ftol
2139:   real ( kind = REAL64 ), INTENT(OUT) ::  fvec(m)2122:   real ( kind = 8 ), INTENT(OUT) ::  fvec(m)
2140:   real ( kind = REAL64 ) gtol2123:   real ( kind = 8 ) gtol
2141:   integer ( kind = INT32 ), INTENT(OUT) ::  info2124:   integer ( kind = 4 ), INTENT(OUT) ::  info
2142:   integer ( kind = INT32 ) ipvt(n)2125:   integer ( kind = 4 ) ipvt(n)
2143:   integer ( kind = INT32 ) maxfev2126:   integer ( kind = 4 ) maxfev
2144:   integer ( kind = INT32 ) mode2127:   integer ( kind = 4 ) mode
2145:   integer ( kind = INT32 ) nfev2128:   integer ( kind = 4 ) nfev
2146:   integer ( kind = INT32 ) njev2129:   integer ( kind = 4 ) njev
2147:   integer ( kind = INT32 ) nprint2130:   integer ( kind = 4 ) nprint
2148:   integer ( kind = INT32 ) iflag2131:   integer ( kind = 4 ) iflag
2149:   real ( kind = REAL64 ) qtf(n)2132:   real ( kind = 8 ) qtf(n)
2150:   real ( kind = REAL64 ), INTENT(IN) ::  tol2133:   real ( kind = 8 ), INTENT(IN) ::  tol
2151:   real ( kind = REAL64 ) x(n)2134:   real ( kind = 8 ) x(n)
2152:   real ( kind = REAL64 ) xtol2135:   real ( kind = 8 ) xtol
2153: 2136: 
2154:   info = 02137:   info = 0
2155: 2138: 
2156:   if ( n <= 0 ) then2139:   if ( n <= 0 ) then
2157:     return2140:     return
2158:   else if ( m < n ) then2141:   else if ( m < n ) then
2159:     return2142:     return
2160:   else if ( ldfjac < m ) then2143:   else if ( ldfjac < m ) then
2161:     return2144:     return
2162:   else if ( tol < 0.0D+00 ) then2145:   else if ( tol < 0.0D+00 ) then
2211: !2194: !
2212: !  Reference:2195: !  Reference:
2213: !2196: !
2214: !    Jorge More, Burton Garbow, Kenneth Hillstrom,2197: !    Jorge More, Burton Garbow, Kenneth Hillstrom,
2215: !    User Guide for MINPACK-1,2198: !    User Guide for MINPACK-1,
2216: !    Technical Report ANL-80-74,2199: !    Technical Report ANL-80-74,
2217: !    Argonne National Laboratory, 1980.2200: !    Argonne National Laboratory, 1980.
2218: !2201: !
2219: !  Parameters:2202: !  Parameters:
2220: !2203: !
2221: !    Input, integer ( kind = INT32 ) N, is the length of the vector.2204: !    Input, integer ( kind = 4 ) N, is the length of the vector.
2222: !2205: !
2223: !    Input, real ( kind = REAL64 ) X(N), the vector whose norm is desired.2206: !    Input, real ( kind = 8 ) X(N), the vector whose norm is desired.
2224: !2207: !
2225: !    Output, real ( kind = REAL64 ) ENORM, the Euclidean norm of the vector.2208: !    Output, real ( kind = 8 ) ENORM, the Euclidean norm of the vector.
2226: !2209: !
2227:   USE PREC, ONLY: INT32, REAL64 
2228:  
2229:   implicit none2210:   implicit none
2230: 2211: 
2231:   integer ( kind = INT32 ) n2212:   integer ( kind = 4 ) n
2232:   real ( kind = REAL64 ) x(n)2213:   real ( kind = 8 ) x(n)
2233:   real ( kind = REAL64 ) enorm2214:   real ( kind = 8 ) enorm
2234: 2215: 
2235:   enorm = sqrt ( sum ( x(1:n) ** 2 ))2216:   enorm = sqrt ( sum ( x(1:n) ** 2 ))
2236: 2217: 
2237:   return2218:   return
2238: end2219: end
2239: 2220: 
2240: function enorm2 ( n, x )2221: function enorm2 ( n, x )
2241: 2222: 
2242: !*****************************************************************************802223: !*****************************************************************************80
2243: !2224: !
2275: !2256: !
2276: !  Reference:2257: !  Reference:
2277: !2258: !
2278: !    Jorge More, Burton Garbow, Kenneth Hillstrom,2259: !    Jorge More, Burton Garbow, Kenneth Hillstrom,
2279: !    User Guide for MINPACK-12260: !    User Guide for MINPACK-1
2280: !    Argonne National Laboratory,2261: !    Argonne National Laboratory,
2281: !    Argonne, Illinois.2262: !    Argonne, Illinois.
2282: !2263: !
2283: !  Parameters:2264: !  Parameters:
2284: !2265: !
2285: !    Input, integer ( kind = INT32 ) N, is the length of the vector.2266: !    Input, integer ( kind = 4 ) N, is the length of the vector.
2286: !2267: !
2287: !    Input, real ( kind = REAL64 ) X(N), the vector whose norm is desired.2268: !    Input, real ( kind = 8 ) X(N), the vector whose norm is desired.
2288: !2269: !
2289: !    Output, real ( kind = REAL64 ) ENORM2, the Euclidean norm of the vector.2270: !    Output, real ( kind = 8 ) ENORM2, the Euclidean norm of the vector.
2290: !2271: !
2291:   USE PREC, ONLY: INT32, REAL64 
2292:  
2293:   implicit none2272:   implicit none
2294: 2273: 
2295:   integer ( kind = INT32 ) n2274:   integer ( kind = 4 ) n
2296: 2275: 
2297:   real ( kind = REAL64 ) agiant2276:   real ( kind = 8 ) agiant
2298:   real ( kind = REAL64 ) enorm22277:   real ( kind = 8 ) enorm2
2299:   integer ( kind = INT32 ) i2278:   integer ( kind = 4 ) i
2300:   real ( kind = REAL64 ) rdwarf2279:   real ( kind = 8 ) rdwarf
2301:   real ( kind = REAL64 ) rgiant2280:   real ( kind = 8 ) rgiant
2302:   real ( kind = REAL64 ) s12281:   real ( kind = 8 ) s1
2303:   real ( kind = REAL64 ) s22282:   real ( kind = 8 ) s2
2304:   real ( kind = REAL64 ) s32283:   real ( kind = 8 ) s3
2305:   real ( kind = REAL64 ) x(n)2284:   real ( kind = 8 ) x(n)
2306:   real ( kind = REAL64 ) xabs2285:   real ( kind = 8 ) xabs
2307:   real ( kind = REAL64 ) x1max2286:   real ( kind = 8 ) x1max
2308:   real ( kind = REAL64 ) x3max2287:   real ( kind = 8 ) x3max
2309: 2288: 
2310:   rdwarf = sqrt ( tiny ( rdwarf ) )2289:   rdwarf = sqrt ( tiny ( rdwarf ) )
2311:   rgiant = sqrt ( huge ( rgiant ) )2290:   rgiant = sqrt ( huge ( rgiant ) )
2312: 2291: 
2313:   s1 = 0.0D+002292:   s1 = 0.0D+00
2314:   s2 = 0.0D+002293:   s2 = 0.0D+00
2315:   s3 = 0.0D+002294:   s3 = 0.0D+00
2316:   x1max = 0.0D+002295:   x1max = 0.0D+00
2317:   x3max = 0.0D+002296:   x3max = 0.0D+00
2318:   agiant = rgiant / real ( n, kind = 8 )2297:   agiant = rgiant / real ( n, kind = 8 )
2429: !2408: !
2430: !  Reference:2409: !  Reference:
2431: !2410: !
2432: !    Jorge More, Burton Garbow, Kenneth Hillstrom,2411: !    Jorge More, Burton Garbow, Kenneth Hillstrom,
2433: !    User Guide for MINPACK-1,2412: !    User Guide for MINPACK-1,
2434: !    Technical Report ANL-80-74,2413: !    Technical Report ANL-80-74,
2435: !    Argonne National Laboratory, 1980.2414: !    Argonne National Laboratory, 1980.
2436: !2415: !
2437: !  Parameters:2416: !  Parameters:
2438: !2417: !
2439: !    Input, integer ( kind = INT32 ) N, the order of R.2418: !    Input, integer ( kind = 4 ) N, the order of R.
2440: !2419: !
2441: !    Input/output, real ( kind = REAL64 ) R(LDR,N),the N by N matrix.  The full2420: !    Input/output, real ( kind = 8 ) R(LDR,N),the N by N matrix.  The full
2442: !    upper triangle must contain the full upper triangle of the matrix R.2421: !    upper triangle must contain the full upper triangle of the matrix R.
2443: !    On output the full upper triangle is unaltered, and the strict lower2422: !    On output the full upper triangle is unaltered, and the strict lower
2444: !    triangle contains the strict upper triangle (transposed) of the upper2423: !    triangle contains the strict upper triangle (transposed) of the upper
2445: !    triangular matrix S.2424: !    triangular matrix S.
2446: !2425: !
2447: !    Input, integer ( kind = INT32 ) LDR, the leading dimension of R.  LDR must be2426: !    Input, integer ( kind = 4 ) LDR, the leading dimension of R.  LDR must be
2448: !    no less than N.2427: !    no less than N.
2449: !2428: !
2450: !    Input, integer ( kind = INT32 ) IPVT(N), defines the permutation matrix P 2429: !    Input, integer ( kind = 4 ) IPVT(N), defines the permutation matrix P 
2451: !    such that A*P = Q*R.  Column J of P is column IPVT(J) of the 2430: !    such that A*P = Q*R.  Column J of P is column IPVT(J) of the 
2452: !    identity matrix.2431: !    identity matrix.
2453: !2432: !
2454: !    Input, real ( kind = REAL64 ) DIAG(N), the diagonal elements of the matrix D.2433: !    Input, real ( kind = 8 ) DIAG(N), the diagonal elements of the matrix D.
2455: !2434: !
2456: !    Input, real ( kind = REAL64 ) QTB(N), the first N elements of the vector Q'*B.2435: !    Input, real ( kind = 8 ) QTB(N), the first N elements of the vector Q'*B.
2457: !2436: !
2458: !    Input, real ( kind = REAL64 ) DELTA, an upper bound on the euclidean norm2437: !    Input, real ( kind = 8 ) DELTA, an upper bound on the euclidean norm
2459: !    of D*X.  DELTA should be positive.2438: !    of D*X.  DELTA should be positive.
2460: !2439: !
2461: !    Input/output, real ( kind = REAL64 ) PAR.  On input an initial estimate of the2440: !    Input/output, real ( kind = 8 ) PAR.  On input an initial estimate of the
2462: !    Levenberg-Marquardt parameter.  On output the final estimate.2441: !    Levenberg-Marquardt parameter.  On output the final estimate.
2463: !    PAR should be nonnegative.2442: !    PAR should be nonnegative.
2464: !2443: !
2465: !    Output, real ( kind = REAL64 ) X(N), the least squares solution of the system2444: !    Output, real ( kind = 8 ) X(N), the least squares solution of the system
2466: !    A*X = B, sqrt(PAR)*D*X = 0, for the output value of PAR.2445: !    A*X = B, sqrt(PAR)*D*X = 0, for the output value of PAR.
2467: !2446: !
2468: !    Output, real ( kind = REAL64 ) SDIAG(N), the diagonal elements of the upper2447: !    Output, real ( kind = 8 ) SDIAG(N), the diagonal elements of the upper
2469: !    triangular matrix S.2448: !    triangular matrix S.
2470: !2449: !
2471:   USE PREC, ONLY: INT32, REAL64 
2472:  
2473:   implicit none2450:   implicit none
2474: 2451: 
2475:   integer ( kind = INT32 ) ldr2452:   integer ( kind = 4 ) ldr
2476:   integer ( kind = INT32 ) n2453:   integer ( kind = 4 ) n
2477: 2454: 
2478:   real ( kind = REAL64 ) delta2455:   real ( kind = 8 ) delta
2479:   real ( kind = REAL64 ) diag(n)2456:   real ( kind = 8 ) diag(n)
2480:   real ( kind = REAL64 ) dwarf2457:   real ( kind = 8 ) dwarf
2481:   real ( kind = REAL64 ) dxnorm2458:   real ( kind = 8 ) dxnorm
2482:   real ( kind = REAL64 ) enorm2459:   real ( kind = 8 ) enorm
2483:   real ( kind = REAL64 ) gnorm2460:   real ( kind = 8 ) gnorm
2484:   real ( kind = REAL64 ) fp2461:   real ( kind = 8 ) fp
2485:   integer ( kind = INT32 ) i2462:   integer ( kind = 4 ) i
2486:   integer ( kind = INT32 ) ipvt(n)2463:   integer ( kind = 4 ) ipvt(n)
2487:   integer ( kind = INT32 ) iter2464:   integer ( kind = 4 ) iter
2488:   integer ( kind = INT32 ) j2465:   integer ( kind = 4 ) j
2489:   integer ( kind = INT32 ) k2466:   integer ( kind = 4 ) k
2490:   integer ( kind = INT32 ) l2467:   integer ( kind = 4 ) l
2491:   integer ( kind = INT32 ) nsing2468:   integer ( kind = 4 ) nsing
2492:   real ( kind = REAL64 ) par2469:   real ( kind = 8 ) par
2493:   real ( kind = REAL64 ) parc2470:   real ( kind = 8 ) parc
2494:   real ( kind = REAL64 ) parl2471:   real ( kind = 8 ) parl
2495:   real ( kind = REAL64 ) paru2472:   real ( kind = 8 ) paru
2496:   real ( kind = REAL64 ) qnorm2473:   real ( kind = 8 ) qnorm
2497:   real ( kind = REAL64 ) qtb(n)2474:   real ( kind = 8 ) qtb(n)
2498:   real ( kind = REAL64 ) r(ldr,n)2475:   real ( kind = 8 ) r(ldr,n)
2499:   real ( kind = REAL64 ) sdiag(n)2476:   real ( kind = 8 ) sdiag(n)
2500:   real ( kind = REAL64 ) sum22477:   real ( kind = 8 ) sum2
2501:   real ( kind = REAL64 ) temp2478:   real ( kind = 8 ) temp
2502:   real ( kind = REAL64 ) wa1(n)2479:   real ( kind = 8 ) wa1(n)
2503:   real ( kind = REAL64 ) wa2(n)2480:   real ( kind = 8 ) wa2(n)
2504:   real ( kind = REAL64 ) x(n)2481:   real ( kind = 8 ) x(n)
2505: !2482: !
2506: !  DWARF is the smallest positive magnitude.2483: !  DWARF is the smallest positive magnitude.
2507: !2484: !
2508:   dwarf = tiny ( dwarf )2485:   dwarf = tiny ( dwarf )
2509: !2486: !
2510: !  Compute and store in X the Gauss-Newton direction.2487: !  Compute and store in X the Gauss-Newton direction.
2511: !2488: !
2512: !  If the jacobian is rank-deficient, obtain a least squares solution.2489: !  If the jacobian is rank-deficient, obtain a least squares solution.
2513: !2490: !
2514:   nsing = n2491:   nsing = n
2731: !2708: !
2732: !  Reference:2709: !  Reference:
2733: !2710: !
2734: !    Jorge More, Burton Garbow, Kenneth Hillstrom,2711: !    Jorge More, Burton Garbow, Kenneth Hillstrom,
2735: !    User Guide for MINPACK-1,2712: !    User Guide for MINPACK-1,
2736: !    Technical Report ANL-80-74,2713: !    Technical Report ANL-80-74,
2737: !    Argonne National Laboratory, 1980.2714: !    Argonne National Laboratory, 1980.
2738: !2715: !
2739: !  Parameters:2716: !  Parameters:
2740: !2717: !
2741: !    Input, integer ( kind = INT32 ) N, the order of R.2718: !    Input, integer ( kind = 4 ) N, the order of R.
2742: !2719: !
2743: !    Input/output, real ( kind = REAL64 ) R(LDR,N), the N by N matrix.2720: !    Input/output, real ( kind = 8 ) R(LDR,N), the N by N matrix.
2744: !    On input the full upper triangle must contain the full upper triangle2721: !    On input the full upper triangle must contain the full upper triangle
2745: !    of the matrix R.  On output the full upper triangle is unaltered, and2722: !    of the matrix R.  On output the full upper triangle is unaltered, and
2746: !    the strict lower triangle contains the strict upper triangle2723: !    the strict lower triangle contains the strict upper triangle
2747: !    (transposed) of the upper triangular matrix S.2724: !    (transposed) of the upper triangular matrix S.
2748: !2725: !
2749: !    Input, integer ( kind = INT32 ) LDR, the leading dimension of R, which must be2726: !    Input, integer ( kind = 4 ) LDR, the leading dimension of R, which must be
2750: !    at least N.2727: !    at least N.
2751: !2728: !
2752: !    Input, integer ( kind = INT32 ) IPVT(N), defines the permutation matrix P such 2729: !    Input, integer ( kind = 4 ) IPVT(N), defines the permutation matrix P such 
2753: !    that A*P = Q*R.  Column J of P is column IPVT(J) of the identity matrix.2730: !    that A*P = Q*R.  Column J of P is column IPVT(J) of the identity matrix.
2754: !2731: !
2755: !    Input, real ( kind = REAL64 ) DIAG(N), the diagonal elements of the matrix D.2732: !    Input, real ( kind = 8 ) DIAG(N), the diagonal elements of the matrix D.
2756: !2733: !
2757: !    Input, real ( kind = REAL64 ) QTB(N), the first N elements of the vector Q'*B.2734: !    Input, real ( kind = 8 ) QTB(N), the first N elements of the vector Q'*B.
2758: !2735: !
2759: !    Output, real ( kind = REAL64 ) X(N), the least squares solution.2736: !    Output, real ( kind = 8 ) X(N), the least squares solution.
2760: !2737: !
2761: !    Output, real ( kind = REAL64 ) SDIAG(N), the diagonal elements of the upper2738: !    Output, real ( kind = 8 ) SDIAG(N), the diagonal elements of the upper
2762: !    triangular matrix S.2739: !    triangular matrix S.
2763: !2740: !
2764:   USE PREC, ONLY: INT32, REAL64 
2765:  
2766:   implicit none2741:   implicit none
2767: 2742: 
2768:   integer ( kind = INT32 ) ldr2743:   integer ( kind = 4 ) ldr
2769:   integer ( kind = INT32 ) n2744:   integer ( kind = 4 ) n
2770: 2745: 
2771:   real ( kind = REAL64 ) c2746:   real ( kind = 8 ) c
2772:   real ( kind = REAL64 ) cotan2747:   real ( kind = 8 ) cotan
2773:   real ( kind = REAL64 ) diag(n)2748:   real ( kind = 8 ) diag(n)
2774:   integer ( kind = INT32 ) i2749:   integer ( kind = 4 ) i
2775:   integer ( kind = INT32 ) ipvt(n)2750:   integer ( kind = 4 ) ipvt(n)
2776:   integer ( kind = INT32 ) j2751:   integer ( kind = 4 ) j
2777:   integer ( kind = INT32 ) k2752:   integer ( kind = 4 ) k
2778:   integer ( kind = INT32 ) l2753:   integer ( kind = 4 ) l
2779:   integer ( kind = INT32 ) nsing2754:   integer ( kind = 4 ) nsing
2780:   real ( kind = REAL64 ) qtb(n)2755:   real ( kind = 8 ) qtb(n)
2781:   real ( kind = REAL64 ) qtbpj2756:   real ( kind = 8 ) qtbpj
2782:   real ( kind = REAL64 ) r(ldr,n)2757:   real ( kind = 8 ) r(ldr,n)
2783:   real ( kind = REAL64 ) s2758:   real ( kind = 8 ) s
2784:   real ( kind = REAL64 ) sdiag(n)2759:   real ( kind = 8 ) sdiag(n)
2785:   real ( kind = REAL64 ) sum22760:   real ( kind = 8 ) sum2
2786:   real ( kind = REAL64 ) t2761:   real ( kind = 8 ) t
2787:   real ( kind = REAL64 ) temp2762:   real ( kind = 8 ) temp
2788:   real ( kind = REAL64 ) wa(n)2763:   real ( kind = 8 ) wa(n)
2789:   real ( kind = REAL64 ) x(n)2764:   real ( kind = 8 ) x(n)
2790: !2765: !
2791: !  Copy R and Q'*B to preserve input and initialize S.2766: !  Copy R and Q'*B to preserve input and initialize S.
2792: !2767: !
2793: !  In particular, save the diagonal elements of R in X.2768: !  In particular, save the diagonal elements of R in X.
2794: !2769: !
2795:   do j = 1, n2770:   do j = 1, n
2796:     r(j:n,j) = r(j,j:n)2771:     r(j:n,j) = r(j,j:n)
2797:     x(j) = r(j,j)2772:     x(j) = r(j,j)
2798:   end do2773:   end do
2799: 2774: 
2934: !2909: !
2935: !  Reference:2910: !  Reference:
2936: !2911: !
2937: !    Jorge More, Burton Garbow, Kenneth Hillstrom,2912: !    Jorge More, Burton Garbow, Kenneth Hillstrom,
2938: !    User Guide for MINPACK-1,2913: !    User Guide for MINPACK-1,
2939: !    Technical Report ANL-80-74,2914: !    Technical Report ANL-80-74,
2940: !    Argonne National Laboratory, 1980.2915: !    Argonne National Laboratory, 1980.
2941: !2916: !
2942: !  Parameters:2917: !  Parameters:
2943: !2918: !
2944: !    Input, integer ( kind = INT32 ) M, the number of rows of A.2919: !    Input, integer ( kind = 4 ) M, the number of rows of A.
2945: !2920: !
2946: !    Input, integer ( kind = INT32 ) N, the number of columns of A.2921: !    Input, integer ( kind = 4 ) N, the number of columns of A.
2947: !2922: !
2948: !    Input/output, real ( kind = REAL64 ) A(LDA,N), the M by N array.2923: !    Input/output, real ( kind = 8 ) A(LDA,N), the M by N array.
2949: !    On input, A contains the matrix for which the QR factorization is to2924: !    On input, A contains the matrix for which the QR factorization is to
2950: !    be computed.  On output, the strict upper trapezoidal part of A contains2925: !    be computed.  On output, the strict upper trapezoidal part of A contains
2951: !    the strict upper trapezoidal part of R, and the lower trapezoidal2926: !    the strict upper trapezoidal part of R, and the lower trapezoidal
2952: !    part of A contains a factored form of Q (the non-trivial elements of2927: !    part of A contains a factored form of Q (the non-trivial elements of
2953: !    the U vectors described above).2928: !    the U vectors described above).
2954: !2929: !
2955: !    Input, integer ( kind = INT32 ) LDA, the leading dimension of A, which must2930: !    Input, integer ( kind = 4 ) LDA, the leading dimension of A, which must
2956: !    be no less than M.2931: !    be no less than M.
2957: !2932: !
2958: !    Input, logical PIVOT, is TRUE if column pivoting is to be carried out.2933: !    Input, logical PIVOT, is TRUE if column pivoting is to be carried out.
2959: !2934: !
2960: !    Output, integer ( kind = INT32 ) IPVT(LIPVT), defines the permutation matrix P 2935: !    Output, integer ( kind = 4 ) IPVT(LIPVT), defines the permutation matrix P 
2961: !    such that A*P = Q*R.  Column J of P is column IPVT(J) of the identity 2936: !    such that A*P = Q*R.  Column J of P is column IPVT(J) of the identity 
2962: !    matrix.  If PIVOT is false, IPVT is not referenced.2937: !    matrix.  If PIVOT is false, IPVT is not referenced.
2963: !2938: !
2964: !    Input, integer ( kind = INT32 ) LIPVT, the dimension of IPVT, which should 2939: !    Input, integer ( kind = 4 ) LIPVT, the dimension of IPVT, which should 
2965: !    be N if pivoting is used.2940: !    be N if pivoting is used.
2966: !2941: !
2967: !    Output, real ( kind = REAL64 ) RDIAG(N), contains the diagonal elements of R.2942: !    Output, real ( kind = 8 ) RDIAG(N), contains the diagonal elements of R.
2968: !2943: !
2969: !    Output, real ( kind = REAL64 ) ACNORM(N), the norms of the corresponding2944: !    Output, real ( kind = 8 ) ACNORM(N), the norms of the corresponding
2970: !    columns of the input matrix A.  If this information is not needed,2945: !    columns of the input matrix A.  If this information is not needed,
2971: !    then ACNORM can coincide with RDIAG.2946: !    then ACNORM can coincide with RDIAG.
2972: !2947: !
2973:   USE PREC, ONLY: INT32, REAL64 
2974:  
2975:   implicit none2948:   implicit none
2976: 2949: 
2977:   integer ( kind = INT32 ) lda2950:   integer ( kind = 4 ) lda
2978:   integer ( kind = INT32 ) lipvt2951:   integer ( kind = 4 ) lipvt
2979:   integer ( kind = INT32 ) m2952:   integer ( kind = 4 ) m
2980:   integer ( kind = INT32 ) n2953:   integer ( kind = 4 ) n
2981: 2954: 
2982:   real ( kind = REAL64 ) a(lda,n)2955:   real ( kind = 8 ) a(lda,n)
2983:   real ( kind = REAL64 ) acnorm(n)2956:   real ( kind = 8 ) acnorm(n)
2984:   real ( kind = REAL64 ) ajnorm2957:   real ( kind = 8 ) ajnorm
2985:   real ( kind = REAL64 ) enorm2958:   real ( kind = 8 ) enorm
2986:   real ( kind = REAL64 ) epsmch2959:   real ( kind = 8 ) epsmch
2987:   integer ( kind = INT32 ) i2960:   integer ( kind = 4 ) i
2988:   integer ( kind = INT32 ) i4_temp2961:   integer ( kind = 4 ) i4_temp
2989:   integer ( kind = INT32 ) ipvt(lipvt)2962:   integer ( kind = 4 ) ipvt(lipvt)
2990:   integer ( kind = INT32 ) j2963:   integer ( kind = 4 ) j
2991:   integer ( kind = INT32 ) k2964:   integer ( kind = 4 ) k
2992:   integer ( kind = INT32 ) kmax2965:   integer ( kind = 4 ) kmax
2993:   integer ( kind = INT32 ) minmn2966:   integer ( kind = 4 ) minmn
2994:   logical pivot2967:   logical pivot
2995:   real ( kind = REAL64 ) r8_temp(m)2968:   real ( kind = 8 ) r8_temp(m)
2996:   real ( kind = REAL64 ) rdiag(n)2969:   real ( kind = 8 ) rdiag(n)
2997:   real ( kind = REAL64 ) temp2970:   real ( kind = 8 ) temp
2998:   real ( kind = REAL64 ) wa(n)2971:   real ( kind = 8 ) wa(n)
2999: 2972: 
3000:   epsmch = epsilon ( epsmch )2973:   epsmch = epsilon ( epsmch )
3001: !2974: !
3002: !  Compute the initial column norms and initialize several arrays.2975: !  Compute the initial column norms and initialize several arrays.
3003: !2976: !
3004:   do j = 1, n2977:   do j = 1, n
3005:     acnorm(j) = enorm ( m, a(1:m,j) )2978:     acnorm(j) = enorm ( m, a(1:m,j) )
3006:   end do2979:   end do
3007: 2980: 
3008:   rdiag(1:n) = acnorm(1:n)2981:   rdiag(1:n) = acnorm(1:n)


r33314/fftw3.f90 2017-09-15 20:30:18.729836978 +0100 r33313/fftw3.f90 2017-09-15 20:30:20.061854790 +0100
  1: MODULE FFTW3  1: MODULE FFTW3
  2:   2: 
  3:   USE PREC, ONLY: INT32  3:   integer ( kind = 4 ), parameter :: fftw_r2hc = 0
  4:   4:   integer ( kind = 4 ), parameter :: fftw_hc2r = 1
  5:   integer ( kind = INT32 ), parameter :: fftw_r2hc = 0  5:   integer ( kind = 4 ), parameter :: fftw_dht = 2
  6:   integer ( kind = INT32 ), parameter :: fftw_hc2r = 1  6:   integer ( kind = 4 ), parameter :: fftw_redft00 = 3
  7:   integer ( kind = INT32 ), parameter :: fftw_dht = 2  7:   integer ( kind = 4 ), parameter :: fftw_redft01 = 4
  8:   integer ( kind = INT32 ), parameter :: fftw_redft00 = 3  8:   integer ( kind = 4 ), parameter :: fftw_redft10 = 5
  9:   integer ( kind = INT32 ), parameter :: fftw_redft01 = 4  9:   integer ( kind = 4 ), parameter :: fftw_redft11 = 6
 10:   integer ( kind = INT32 ), parameter :: fftw_redft10 = 5 10:   integer ( kind = 4 ), parameter :: fftw_rodft00 = 7
 11:   integer ( kind = INT32 ), parameter :: fftw_redft11 = 6 11:   integer ( kind = 4 ), parameter :: fftw_rodft01 = 8
 12:   integer ( kind = INT32 ), parameter :: fftw_rodft00 = 7 12:   integer ( kind = 4 ), parameter :: fftw_rodft10 = 9
 13:   integer ( kind = INT32 ), parameter :: fftw_rodft01 = 8 13:   integer ( kind = 4 ), parameter :: fftw_rodft11 = 10
 14:   integer ( kind = INT32 ), parameter :: fftw_rodft10 = 9 14:   integer ( kind = 4 ), parameter :: fftw_forward = -1
 15:   integer ( kind = INT32 ), parameter :: fftw_rodft11 = 10 15:   integer ( kind = 4 ), parameter :: fftw_backward = +1
 16:   integer ( kind = INT32 ), parameter :: fftw_forward = -1 16:   integer ( kind = 4 ), parameter :: fftw_measure = 0
 17:   integer ( kind = INT32 ), parameter :: fftw_backward = +1 17:   integer ( kind = 4 ), parameter :: fftw_destroy_input = 1
 18:   integer ( kind = INT32 ), parameter :: fftw_measure = 0 18:   integer ( kind = 4 ), parameter :: fftw_unaligned = 2
 19:   integer ( kind = INT32 ), parameter :: fftw_destroy_input = 1 19:   integer ( kind = 4 ), parameter :: fftw_conserve_memory = 4
 20:   integer ( kind = INT32 ), parameter :: fftw_unaligned = 2 20:   integer ( kind = 4 ), parameter :: fftw_exhaustive = 8
 21:   integer ( kind = INT32 ), parameter :: fftw_conserve_memory = 4 21:   integer ( kind = 4 ), parameter :: fftw_preserve_input = 16
 22:   integer ( kind = INT32 ), parameter :: fftw_exhaustive = 8 22:   integer ( kind = 4 ), parameter :: fftw_patient = 32
 23:   integer ( kind = INT32 ), parameter :: fftw_preserve_input = 16 23:   integer ( kind = 4 ), parameter :: fftw_estimate = 64
 24:   integer ( kind = INT32 ), parameter :: fftw_patient = 32 24:   integer ( kind = 4 ), parameter :: fftw_estimate_patient = 128
 25:   integer ( kind = INT32 ), parameter :: fftw_estimate = 64 25:   integer ( kind = 4 ), parameter :: fftw_believe_pcost = 256
 26:   integer ( kind = INT32 ), parameter :: fftw_estimate_patient = 128 26:   integer ( kind = 4 ), parameter :: fftw_dft_r2hc_icky = 512
 27:   integer ( kind = INT32 ), parameter :: fftw_believe_pcost = 256 27:   integer ( kind = 4 ), parameter :: fftw_nonthreaded_icky = 1024
 28:   integer ( kind = INT32 ), parameter :: fftw_dft_r2hc_icky = 512 28:   integer ( kind = 4 ), parameter :: fftw_no_buffering = 2048
 29:   integer ( kind = INT32 ), parameter :: fftw_nonthreaded_icky = 1024 29:   integer ( kind = 4 ), parameter :: fftw_no_indirect_op = 4096
 30:   integer ( kind = INT32 ), parameter :: fftw_no_buffering = 2048 30:   integer ( kind = 4 ), parameter :: fftw_allow_large_generic = 8192
 31:   integer ( kind = INT32 ), parameter :: fftw_no_indirect_op = 4096 31:   integer ( kind = 4 ), parameter :: fftw_no_rank_splits = 16384
 32:   integer ( kind = INT32 ), parameter :: fftw_allow_large_generic = 8192 32:   integer ( kind = 4 ), parameter :: fftw_no_vrank_splits = 32768
 33:   integer ( kind = INT32 ), parameter :: fftw_no_rank_splits = 16384 33:   integer ( kind = 4 ), parameter :: fftw_no_vrecurse = 65536
 34:   integer ( kind = INT32 ), parameter :: fftw_no_vrank_splits = 32768 34:   integer ( kind = 4 ), parameter :: fftw_no_simd = 131072
 35:   integer ( kind = INT32 ), parameter :: fftw_no_vrecurse = 65536 
 36:   integer ( kind = INT32 ), parameter :: fftw_no_simd = 131072 
 37:  35: 
 38: END MODULE FFTW3 36: END MODULE FFTW3


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0