This is the mail archive of the
gcc-bugs@gcc.gnu.org
mailing list for the GCC project.
Mysterious core dump in g77
- To: egcs-bugs at egcs dot cygnus dot com
- Subject: Mysterious core dump in g77
- From: Frank Kuehnel <kuehnel at pa dot msu dot edu>
- Date: Tue, 21 Dec 1999 12:35:47 -0500
- Organization: Michigan State University
Hi,
I had to include a rather large an uninteresting program code for
this bug to appear, it seems to be very context sensitive. The only
two interesting lines are 1030 and 1031, where it appears that a simple
double precision multiplication produces a core dump! Compiled with DEC
Fortran, the code works fine, compiled with g77 I get the core dump.
This is my configuration:
g77 version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release) (from
FSF-g77 version 0.5.24-19981002)
Driving: f77 -v -c -xf77-version /dev/null -xnone
Reading specs from /usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/specs
gcc version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release)
/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/cpp -lang-c -v -undef
-D__GNUC__=2 -D__GNUC_MINOR__=91 -D__ELF__ -D__unix__ -D__i386__
-D__i386__ -D__linux__ -D__unix -D__i386 -D__linux -Asystem(posix)
-D_LANGUAGE_FORTRAN -traditional -Asystem(unix) -Acpu(i386)
-Amachine(i386) -Di386 -D__i386 -D__i386__ /dev/null /dev/null
GNU CPP version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release) (i386
Linux/ELF)
#include "..." search starts here:
#include <...> search starts here:
/usr/local/include
/usr/i386-redhat-linux/include
/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/include
/usr/include
End of search list.
/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/f771 -fnull-version
-quiet -dumpbase g77-version.f -version -fversion -o /tmp/ccG8c728.s
/dev/null
GNU F77 version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release)
(i386-redhat-linux) compiled by GNU C version egcs-2.91.66
19990314/Linux (egcs-1.1.2 release).
GNU Fortran Front End version 0.5.24-19981002
as -V -Qy -o /tmp/ccIrG6E9.o /tmp/ccG8c728.s
GNU assembler version 2.9.1 (i386-redhat-linux), using BFD version
2.9.1.0.23
ld -m elf_i386 -dynamic-linker /lib/ld-linux.so.2 -o /tmp/ccSsg7wg
/tmp/ccIrG6E9.o /usr/lib/crt1.o /usr/lib/crti.o
/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/crtbegin.o
-L/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66
-L/usr/i386-redhat-linux/lib -lg2c -lm -lgcc -lc -lgcc
/usr/lib/gcc-lib/i386-redhat-linux/egcs-2.91.66/crtend.o /usr/lib/crtn.o
/tmp/ccSsg7wg
__G77_LIBF77_VERSION__: 0.5.24
@(#)LIBF77 VERSION 19970919
__G77_LIBI77_VERSION__: 0.5.24-19981021
@(#) LIBI77 VERSION pjw,dmg-mods 19980617
__G77_LIBU77_VERSION__: 0.5.24-19990305
@(#) LIBU77 VERSION 19980709
To compile the program, you'll need the lapack-2.0-10 package:
f77 -Wall -Wimplicit flaw.f -llapack
Here is the program code:
PROGRAM MinAction
c ...
c ... program for minimizing the Euclidean short-range action
INTEGER NP, MAXSTP
DOUBLE PRECISION PI, LOG8, LOG4
PARAMETER (NP = 40, MAXSTP = 100)
PARAMETER (PI = 3.14159265358979d0)
PARAMETER (LOG8 = 2.07944154154168d0)
PARAMETER (LOG4 = 1.38629436111989d0)
c ... NP number of angular momentum basis-wavefunctions \psi_j(r)
c ...
c ... **********************************************************
c ... * *
c ... * for all matrices, only the upper right half is used *
c ... * *
c ... **********************************************************
c ...
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
c ... MVmn is the potential matrix v_{mn} which describes V
c ... MVjk is the interaction potential matrix of the nonlocal
c ... potential V (potential operator)
c ... MdEvpq is the derivative matrix of {\partial E\over \partial
v_{p,q}}
c ...
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
c ... MdVvpq is the derivative matrix of
c ... {\partial \over \partial v_{p,q}} \int V^2/2 d{\bf r}
c ... Mdovpq is the derivative matrix of
c ... {\partial \over \partial v_{p,q}} {1\over 2\sqrt{2} E}
c ... \int {\partial V\over\partial x} \psi_t^\ast\psi_b d{\bf r}
c ... ck is the matrix of eigenvectors, it physically shares the
same
c ... memory as MVjk, because of the diagonalization routine!
c ... Here ck(j,k) corresponds to the j component of the k
eigenvector
c ... c is the component vector of an initial guess wavefunction
EQUIVALENCE (MVjk, ck)
c ...
DOUBLE PRECISION eigenv(NP)
c ... eigenv contains the eigenvalues of the diagonalized matrix
MVjk
c ...
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION psqrt(NP), pbinom(NP,NP,NP)
COMMON /precalc/ psqrt, pbinom
c ... frequently used numbers are stored in arrays!
c ... psqrt(j) = DSQRT(j), j=1, NP used in ovlap(), dovpq()
c ... pbinom(j+1,k+1,m+1) = 2^(-k-m) \sqrt{ k+m \choose m}\sqrt{ k+m
\choose j},
c ... for j=0, NP-1; k=j+1, NP-1; m=0, NP-1-(k-j)
c ...
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ... FuncR is the value of the functional R = omega^2 d - \ln
(omega overl)^2
c ... ExpmR is \exp(-R)
c ... omega is the frequency (measured in terms of v)
c ... d is the weight integral {1/4E^2} \int V^2/2 d{\bf r}
c ... overl is the overlap integral 1/2E \int dV/dx \psi_t^\ast
\psi_b d{\bf r}
c ...
DOUBLE PRECISION Mxi(NP,NP), Mg(NP,NP), Mh(NP,NP)
DOUBLE PRECISION Mpcom(NP, NP), Mxicom(NP, NP)
COMMON /gradmethod/ Mxi, Mg, Mh, Mpcom, Mxicom
c ... notation is similar to Numerical Recipes (NR), chapter 10.6
c ... Mxi contains the gradient of R at the point Vmn
c ... Mg, Mh are downhill gradient and conjugate downhill gradient
(NR) ch. 10.6
c ... Mpcom, Mxicom are temporary vectors (matrices) for the line
minimization
c ...
DOUBLE PRECISION f, gg, dirW0, dilnG0, x0min
DOUBLE PRECISION dx0min, R0min, dR0min
DOUBLE PRECISION dzero, Gzero, omega0, lambda, qualit
INTEGER steps, iter
c ...
c ... precompute factors
CALL precom()
c ...
c ... prepare various matrices! clear the diagonal, and even
off-diagonal lines
CALL premat()
c ...
c ... test the gradient
c CALL tgrad()
c ...
c ... all right here we go!
GOTO 200
c ...
c ... first excercise: compare x-dependence of numerics to
variational method
PRINT '( "#",A7, 5(1X,A14) )', 'x0','W0',
& '-ln G0**2', 'direct W0', 'direct-lnG0**2'
omega = 1.0d0
x0min = 0.05d0
DO steps=1, MAXSTP
CALL initc(x0min)
c ... x0min corresponds to x0 in the notes
CALL gMVmn()
CALL EvalR()
dirW0 = PI*(1+DEXP(-x0min**2))
dilnG0= 2.0d0*x0min**2-DLOG(x0min**2)+LOG4
WRITE(*,10) x0min,d,-DLOG(overl**2),dirW0,dilnG0
x0min = x0min + 0.05d0
END DO
GO TO 300
10 FORMAT( F8.2, 4(1X,E14.6) )
c ... second excercise: compare omega-dependence in numerical to
direct variational method
100 PRINT '( "#",A5, 4(1X,A14) )', 'omega','x0-min','R0 min',
& 'direct x0-min', 'direct R0 min'
omega0 = 0.5d0
omega = omega0
DO steps=1, MAXSTP
c ... x0min, dxomin corresponds to x0 in the notes
CALL doptxp(dR0min, dx0min)
CALL optxps(R0min, x0min)
WRITE(*,20) omega,x0min,R0min,dx0min,dR0min
omega = omega0 + 0.05d0 * DBLE(steps)
END DO
GO TO 300
20 FORMAT( F6.3, 4(1X,E14.6) )
c ... now lets do the real minimization method
c ... first search for the optimal x0 distance
c ... then let it all go, and do real minimization
200 PRINT '( "#",A5,8(1X,A12),1X,A4,2(1X,A10))','omega','x0-min',
& 'W0/omega^2','(G0/omega)^2','R0min','W/omega^2',
& '(G/omega)^2','Rmin','g-metric','iter','lambda',
& 'quality'
omega0 = 0.4d0
omega = omega0
CALL initc(0.5d0)
CALL gMVmn()
c ... write startup potential
CALL wrMat(MVmn,"MVmn",0.0d0)
DO steps=1, MAXSTP
c ... x0min, corresponds to x0 in the notes
CALL optxps(R0min, x0min)
dzero = d
Gzero = overl
R0min = FuncR
c ... now do the real minimization with given setup
CALL initc(0.5d0)
CALL gMVmn()
c ... call conjugate gradient method
CALL cgradm(iter, DBLE(1.0e-8), f, gg)
c ... write optimal potential configuratio in file MVmn<omega>
CALL wrMat(MVmn,"MVmn",omega)
c ... found optimal potential, is it really optimal?
CALL veri(lambda,qualit)
WRITE(*,30) omega,x0min,dzero,Gzero**2,R0min,d,
& overl**2,f,gg,iter,lambda,qualit
omega = omega0 + 0.036d0*DBLE(steps)
END DO
30 FORMAT( F6.3, 8(1X,E12.5), 1X, I4, 2(1X,E10.4) )
300 STOP 'successfully terminated'
END
c ...
SUBROUTINE veri(lambda, qualit)
DOUBLE PRECISION lambda, qualit
c ...
INTEGER NP,MAXK
DOUBLE PRECISION SQRT2
PARAMETER (NP = 40, MAXK = 5)
PARAMETER (SQRT2 = 1.41421356237309d0)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
DOUBLE PRECISION Mxi(NP,NP), Mg(NP,NP), Mh(NP,NP)
DOUBLE PRECISION Mpcom(NP, NP), Mxicom(NP, NP)
COMMON /gradmethod/ Mxi, Mg, Mh, Mpcom, Mxicom
c ...
DOUBLE PRECISION psqrt(NP), pbinom(NP,NP,NP)
COMMON /precalc/ psqrt, pbinom
c ...
INTEGER j, k, n, p, q, switch
DOUBLE PRECISION pred,preov,sumk,maxvmn
DOUBLE PRECISION sumn, ovarr(NP-1)
c ...
c ... be sure to have the right data
CALL Evalal()
CALL gradR()
c ... obtain the reference potential
c ... calculate contribution from logarithmic correction
DO k=1, NP-1
switch = 1
sumn = DCMPLX(0.)
DO n=0, NP-2
IF (switch.EQ.1) THEN
sumn = sumn + psqrt(n+1)
& * ( ck(n+1,k)*ck(n+2,NP) + ck(n+1,NP)*ck(n+2,k) )
ELSE
sumn = sumn - psqrt(n+1)
& * ( ck(n+1,k)*ck(n+2,NP) + ck(n+1,NP)*ck(n+2,k) )
ENDIF
switch = -switch
END DO
ovarr(k) = sumn / (eigenv(NP) - eigenv(k))
END DO
c ... some useful prefactors
pred = omega/(2.0d0*eigenv(NP))
preov = 2.0d0*SQRT2*eigenv(NP)/(overl*omega)
c ...
c ... copy omega/(2 E)*MVmn into Mpcom
c ... 2\sqrt{2}E/G * ... into Mxicom
c ... the full reference potential into Mg
DO p=1, NP
switch = 0
DO q=p+1, NP
IF (switch.EQ.0) THEN
switch = 1
Mpcom(p,q) = pred * MVmn(p,q)
sumk = 0.0d0
DO k=1, NP-1
sumk = sumk + (ck(p,k)*ck(q,NP) +
& ck(q,k)*ck(p,NP))*ovarr(k)
END DO
Mxicom(p,q) = sumk*preov
Mg(p,q) = Mpcom(p,q) - Mxicom(p,q)
ELSE
switch = 0
ENDIF
END DO
END DO
c ... write reference potential and overlap correction
CALL wrMat(Mpcom,"Mref",omega)
CALL wrMat(Mxicom,"Movl",omega)
c ...
c ... copy eigenvector ck(j,NP) into c to generate potential
DO j=1, NP
c(j) = ck(j,NP)
END DO
CALL gMVmn()
c ... write eigenvector potential
CALL wrMat(MVmn,"Mpsi",omega)
c ... find lambda and compare both matrices
c ... find the largest entry in MVmn
c maxvmn = 0.0d0
c p = 0
c q = 0
c DO j=1, NP
c switch = 0
c DO k=j+1, NP
c IF (switch.EQ.0) THEN
c switch = 1
c IF (DBLE(MVmn(j,k)).GT.maxvmn) THEN
c maxvmn = DBLE(MVmn(j,k))
c p = j
c q = k
c ENDIF
c ELSE
c switch = 0
c ENDIF
c END DO
c END DO
p = 1
q = 2
c ... something went wrong?
IF (p.EQ.0) THEN
STOP 'verify: cant find lambda!'
ENDIF
c ... all right get lambda
lambda = DBLE(Mg(p,q)/MVmn(p,q))
c ... compare both matrices!
qualit = 0.0d0
DO p=1, NP
switch = 0
DO q=p+1, NP
IF (switch.EQ.0) THEN
switch = 1
qualit = qualit + (DBLE((Mg(p,q)
& - lambda * MVmn(p,q))))**2
ELSE
switch = 0
ENDIF
END DO
END DO
c ...
RETURN
20 FORMAT( <MAXK>('(',(E12.5,'+',E12.5,'i'),"),") )
END
c ...
SUBROUTINE wrMat(mat,prefix,omega)
DOUBLE PRECISION omega
CHARACTER prefix*4
c ...
INTEGER NP
PARAMETER (NP = 40)
c ...
DOUBLE PRECISION mat(NP,NP)
c ...
INTEGER j,k,MAXK,MAXJ, switch
PARAMETER (MAXK=15,MAXJ=15)
CHARACTER FNAME*10
c ...
c ... symmetrize mat
DO j=1,NP
switch = 0
DO k=j+1, NP
IF (switch.EQ.0) THEN
switch = 1
mat(k,j) = mat(j,k)
ELSE
switch = 0
ENDIF
END DO
END DO
c ...
IF (omega.LT.10.0d0) THEN
WRITE(FNAME, 20) prefix,omega
ELSE
WRITE(FNAME, 30) prefix,omega
ENDIF
c ...
c PRINT '( "write file: ",A10 )',FNAME
c ...
OPEN (UNIT=1, FILE=FNAME, FORM='FORMATTED')
WRITE (1, 10) ((DBLE(mat(j,k)), k=1, MAXK), j=1, MAXJ)
CLOSE (UNIT=1)
RETURN
10 FORMAT ( <MAXK>(E12.5,1X) )
20 FORMAT ( A4,F5.3 )
30 FORMAT ( A4,F5.2 )
END
c ...
SUBROUTINE tgrad()
c ... test the gradient at a given point
INTEGER NP,MAXK
DOUBLE PRECISION PI
PARAMETER (NP = 40,MAXK = 5)
PARAMETER (PI = 3.14159265358979d0)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION Mxi(NP,NP), Mg(NP,NP), Mh(NP,NP)
DOUBLE PRECISION Mpcom(NP, NP), Mxicom(NP, NP)
COMMON /gradmethod/ Mxi, Mg, Mh, Mpcom, Mxicom
c ...
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
DOUBLE PRECISION psqrt(NP), pbinom(NP,NP,NP)
COMMON /precalc/ psqrt, pbinom
c ...
DOUBLE PRECISION FuncR0
INTEGER j,k,switch,p,q
c ...
CALL initc(1.5d0)
CALL gMVmn()
c ...
PRINT '( "MVmn" )'
WRITE(*,20) ((MVmn(p,q), q=1,MAXK), p=1,MAXK)
c ...
DO j=1, NP
switch = 0
DO k=j+1, NP
IF (switch.EQ.0) THEN
switch = 1
Mxicom(j,k) = DABS(0.01d0*MVmn(j,k))
c ... only the upper right half is used
ELSE
switch = 0
ENDIF
END DO
END DO
c ...
CALL gMVjk()
omega = 1.0d0
CALL Evalal()
FuncR0 = FuncR
CALL gradR()
c ...
DO j=1, MAXK
switch = 0
DO k=j+1, MAXK
IF (switch.EQ.0) THEN
switch = 1
CALL gMVmn()
MVmn(j,k) = MVmn(j,k)+Mxicom(j,k)
CALL EvalR()
c IF (DABS((FuncR-FuncR0)/FuncR0).LT.
c & DBLE(1e-5)) THEN
c Mg(j,k) = 0.0d0
c ELSE
Mg(j,k) = (FuncR-FuncR0)/Mxicom(j,k)
c ENDIF
ELSE
switch = 0
ENDIF
END DO
END DO
PRINT '( "Mxi" )'
WRITE(*,20) ((Mxi(p,q), q=1,MAXK), p=1,MAXK)
PRINT '( "numerical" )'
WRITE(*,20) ((Mg(p,q), q=1,MAXK), p=1,MAXK)
STOP
20 FORMAT( <MAXK>(E12.5) )
RETURN
END
c ...
SUBROUTINE premat()
c ... prepare matrices
INTEGER NP
PARAMETER (NP = 40)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION Mxi(NP,NP), Mg(NP,NP), Mh(NP,NP)
DOUBLE PRECISION Mpcom(NP, NP), Mxicom(NP, NP)
COMMON /gradmethod/ Mxi, Mg, Mh, Mpcom, Mxicom
c ...
INTEGER j, k, switch
c ...
DO j=1, NP
switch = 0
DO k=j, NP
IF (switch.EQ.0) THEN
switch = 1
c ... k >= j
MVmn(j,k) = 0.0d0
c ... can be omitted, since each time MVjk is calculated
c ... the even off diagonal elements are deleted!
c ... diagonalization routine uses this as output array!
c MVjk(j,k) = 0.0d0
MdEvpq(j,k)= 0.0d0
MdVvpq(j,k)= 0.0d0
Mxi(j,k) = 0.0d0
Mg(j,k) = 0.0d0
Mh(j,k) = 0.0d0
Mpcom(j,k) = 0.0d0
Mxicom(j,k)= 0.0d0
ELSE
switch = 0
ENDIF
END DO
END DO
RETURN
END
c ...
SUBROUTINE precom()
c ... precompute factors
INTEGER NP
PARAMETER (NP = 40)
DOUBLE PRECISION psqrt(NP), pbinom(NP,NP,NP)
COMMON /precalc/ psqrt, pbinom
c ...
DOUBLE PRECISION p2j, p2k, p2km
DOUBLE PRECISION bjkm, bjk
INTEGER j, k, m
c ...
DO j=1, NP
psqrt(j) = DSQRT(DBLE(j))
END DO
c ...
p2j = 0.5d0
DO j=0, NP-1
p2k = p2j
bjk = DBLE(j+1)
DO k=j+1, NP-1
p2km = p2k
bjkm = DSQRT(bjk)
pbinom(j+1,k+1,1) = p2km * bjkm
DO m=1, NP-1-(k-j)
p2km = p2km * 0.5d0
bjkm = bjkm * DBLE(k+m) / DSQRT(DBLE(m*(k+m-j)))
pbinom(j+1,k+1,m+1) = p2km * bjkm
END DO
p2k = p2k * 0.5d0
bjk = bjk * DBLE(k+1) / DBLE(k-j+1)
END DO
p2j = p2j * 0.5d0
END DO
RETURN
END
c ...
SUBROUTINE cgradm(iter, ftol, fret, gg)
c ... here we start the minimization procedure
DOUBLE PRECISION ftol, fret
INTEGER iter
c ...
INTEGER NP, ITMAX
DOUBLE PRECISION EPS
PARAMETER (NP = 40, ITMAX = 100, EPS=1.e-10)
c ...
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
DOUBLE PRECISION Mxi(NP,NP), Mg(NP,NP), Mh(NP,NP)
DOUBLE PRECISION Mpcom(NP, NP), Mxicom(NP, NP)
COMMON /gradmethod/ Mxi, Mg, Mh, Mpcom, Mxicom
c ...
DOUBLE PRECISION gam, dgg, gg, fp
c ... gam, dgg, gg, fp are defined similarily to NR chapter 10.6
c ...
INTEGER j,k,switch,its
c ...
c ... setup initial point
CALL Evalal()
fp = FuncR
CALL gradR()
DO j=1, NP
switch = 0
DO k=j+1, NP
IF (switch.EQ.0) THEN
switch = 1
Mg(j,k) = -Mxi(j,k)
Mh(j,k) = Mg(j,k)
Mxi(j,k) = Mh(j,k)
ELSE
switch = 0
ENDIF
END DO
END DO
c ... lets minimize here
DO its=1, ITMAX
c PRINT '( "cgradm: iter", I4, 2X, "value", E16.8 )', its, fp
iter = its
c ... minimize R along the direction Mxi starting at the point MVmn
c ... returns the new minimum position MVmn, the distance vector Mxi
c ... and the values d, overl, FuncR at this very point
c ... PRINT '( "cgradm: call linemin" )'
CALL linmin(fret)
c ... PRINT '( "cgradm: back from linemin", F8.2 )', fret
IF (2.0d0*DABS(fret-fp).LT.
& ftol*(DABS(fret)+DABS(fp)+EPS)) RETURN
fp=fret
c ... calculate all eigenvalues and eigenvectors
c ... PRINT '( "cgradm: calculate eigenvalues and eigenvectors" )'
CALL Evalal()
c ... having evaluated R, return grad R at a given point MVmn
c ... PRINT '( "cgradm: calculate gradient" )'
CALL gradR()
gg = 0.0d0
dgg = 0.0d0
DO j=1, NP
switch = 0
DO k=j+1, NP
IF (switch.EQ.0) THEN
switch = 1
c ... this is the gradient metric, shows how well we`re doing
gg = gg + Mg(j,k)**2
dgg = dgg + (Mxi(j,k)+Mg(j,k))*Mxi(j,k)
ELSE
switch = 0
ENDIF
END DO
END DO
c ... very unlikely to happen!
c PRINT '( "cgradm: gradient metric", 1X, E16.8 )',gg
IF (gg.EQ.0.0d0) RETURN
gam = dgg / gg
DO j=1, NP
switch = 0
DO k=j+1, NP
IF (switch.EQ.0) THEN
switch = 1
Mg(j,k) = -Mxi(j,k)
Mh(j,k) = Mg(j,k) + gam * Mh(j,k)
Mxi(j,k) = Mh(j,k)
ELSE
switch = 0
ENDIF
END DO
END DO
END DO
PAUSE 'cgradm maximum iterations exceeded'
RETURN
END
c ...
SUBROUTINE optxps(fret, xmin)
DOUBLE PRECISION fret, xmin
c ... minimize R along by varying the initial x-pos setup
c ... returns the new minimum position MVmn,
c ... and the values d, overl, FuncR at this very point
c ...
DOUBLE PRECISION TOL
PARAMETER (TOL=1.e-4)
c ...
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
DOUBLE PRECISION fx1dim, brent
EXTERNAL fx1dim, brent
DOUBLE PRECISION ax,bx,fa,fb,fx,xx
c ...
ax = 0.001d0
xx = 0.05d0
CALL mnbrak(ax,xx,bx,fa,fx,fb,fx1dim)
fret=brent(ax,xx,bx,fx1dim,DBLE(TOL),xmin)
fret = FuncR
RETURN
END
c ...
SUBROUTINE doptxp(fret, xmin)
DOUBLE PRECISION fret, xmin
c ... minimize R along by varying the initial x-pos setup
c ... in direct variational method
c ...
DOUBLE PRECISION LOG8, LOG4, PI
PARAMETER (PI = 3.14159265358979d0)
PARAMETER (LOG8 = 2.07944154154168d0)
PARAMETER (LOG4 = 1.38629436111989d0)
DOUBLE PRECISION TOL
PARAMETER (TOL=1.e-4)
c ...
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
DOUBLE PRECISION dfx1di, brent,dirlnG
EXTERNAL dfx1di, brent
DOUBLE PRECISION ax,bx,fa,fb,fx,xx
c ...
ax = 0.001d0
xx = 0.05d0
CALL mnbrak(ax,xx,bx,fa,fx,fb,dfx1di)
fret=brent(ax,xx,bx,dfx1di,DBLE(TOL),xmin)
c ... return minimum
dirlnG= 2.0d0*xmin**2 - DLOG((omega*xmin)**2)+LOG4
fret= PI*(1+DEXP(-xmin**2))*omega**2 + dirlnG
RETURN
END
c ...
SUBROUTINE linmin(fret)
DOUBLE PRECISION fret
c ... minimize R along the direction Mxi starting at the point MVmn
c ... returns the new minimum position MVmn, the distance vector Mxi
c ... and the values d, overl, FuncR at this very point
c ...
INTEGER NP
DOUBLE PRECISION TOL
PARAMETER (NP = 40, TOL=1.e-5)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION Mxi(NP,NP), Mg(NP,NP), Mh(NP,NP)
DOUBLE PRECISION Mpcom(NP, NP), Mxicom(NP, NP)
COMMON /gradmethod/ Mxi, Mg, Mh, Mpcom, Mxicom
c ...
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
DOUBLE PRECISION f1dim, brent
EXTERNAL f1dim, brent
c ...
DOUBLE PRECISION ax,bx,fa,fb,fx,xmin,xx
INTEGER j, k, switch
c ... set up the common block
DO j=1, NP
switch = 0
DO k=j+1, NP
IF (switch.EQ.0) THEN
switch = 1
Mpcom(j,k) = MVmn(j,k)
Mxicom(j,k) = Mxi(j,k)
ELSE
switch = 0
ENDIF
END DO
END DO
ax = 0.0d0
xx = 0.002d0
c ... PRINT '( "linemin: call mnbrak" )'
CALL mnbrak(ax,xx,bx,fa,fx,fb,f1dim)
c ... PRINT '( "linemin: call brent" )'
fret=brent(ax,xx,bx,f1dim,DBLE(TOL),xmin)
c ... PRINT '( "linemin: return minimum" )'
c ... construct the vector results to return
DO j=1, NP
switch = 0
DO k=j+1, NP
IF (switch.EQ.0) THEN
switch = 1
Mxi(j,k) = xmin*Mxicom(j,k)
MVmn(j,k) = Mpcom(j,k) + Mxi(j,k)
ELSE
switch = 0
ENDIF
END DO
END DO
CALL EvalR()
c ... return minimum for R
fret = FuncR
RETURN
END
c ...
SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func)
DOUBLE PRECISION ax,bx,cx,fa,fb,fc,func,GOLD,GLIMIT,TINY
EXTERNAL func
PARAMETER (GLIMIT=100.0d0, TINY=1.0e-20)
PARAMETER (GOLD=1.61803398874989d0)
c ... Given a function func, and given distinct initial points ax
and bx,
c ... this routine searches in the uphill direction (defined by the
function
c ... as evaluated at the initial points) and returns new points ax,
bx, cx
c ... that bracket a maximum of the function. Also returned are the
function
c ... values at the three points, fa, fb, and fc. Parameters: GOLD
is the default
c ... ratio by which successive intervals are magnified; GLIMIT is
the maximum
c ... magnifcation allowed for a parabolic-fit step.
DOUBLE PRECISION dum,fu,q,r,u,ulim
fa=func(ax)
c ... PRINT '( "mnbrak: function value ",E16.8,1X,"at ",E16.8 )',
c & fa,ax
fb=func(bx)
c ... PRINT '( "mnbrak: function value ",E16.8,1X,"at ",E16.8 )',
c & fb,bx
IF (fa.GT.fb) THEN
c ... Switch roles of a and b so that we can go uphill in the
direction from a to b.
dum=ax
ax=bx
bx=dum
dum=fb
fb=fa
fa=dum
ENDIF
cx=bx+GOLD*(bx-ax)
c ... First guess for c.
fc=func(cx)
1 IF (fc.GE.fb) THEN
c ... "do while": keep returning here until we bracket.
r=(bx-ax)*(fb-fc)
c ... Compute u by parabolic extrapolation from a; b; c. TINY is
used to prevent
c ... any possible division by zero.
q=(bx-cx)*(fb-fa)
u=bx-((bx-cx)*q-(bx-ax)*r)/
& (2.0d0*DSIGN(DMAX1(DABS(q-r),DBLE(TINY)),q-r))
ulim=bx+GLIMIT*(cx-bx)
c ... We won't go farther than this. Test various possibilities:
IF ((bx-u)*(u-cx).GT.0.0d0) THEN
c ... Parabolic u is between b and c: try it.
fu=func(u)
IF (fu.GT.fc) THEN
c ... Got a maximum between b and c.
ax=bx
fa=fb
bx=u
fb=fu
RETURN
ELSE IF (fu.LT.fb) THEN
c ... Got a maximum between between a and u.
cx=u
fc=fu
RETURN
ENDIF
u=cx+GOLD*(cx-bx)
c ... Parabolic fit was no use. Use default magnifcation.
fu=func(u)
ELSE IF ((cx-u)*(u-ulim).GT.0.0d0) THEN
c ... Parabolic fit is between c and its allowed limit.
fu=func(u)
IF (fu.GT.fc) THEN
bx=cx
cx=u
u=cx+GOLD*(cx-bx)
fb=fc
fc=fu
fu=func(u)
ENDIF
ELSE IF ((u-ulim)*(ulim-cx).GE.0.0d0) THEN
c ... Limit parabolic u to maximum allowed value.
u=ulim
fu=func(u)
ELSE
c ... Reject parabolic u, use default magnifcation.
u=cx+GOLD*(cx-bx)
fu=func(u)
ENDIF
ax=bx
c ... Eliminate oldest point and continue.
bx=cx
cx=u
fa=fb
fb=fc
fc=fu
GOTO 1
ENDIF
c ... PRINT '( "mnbrak: ax bx cx",3(1X,E16.8) )',ax,bx,cx
c ... PRINT '( "mnbrak: fa fb fc",3(1X,E16.8) )',fa,fb,fc
RETURN
END
c ...
FUNCTION brent(ax,bx,cx,f,tol,xmax)
INTEGER ITMAX
DOUBLE PRECISION brent,ax,bx,cx,tol,xmax,f,CGOLD,ZEPS
EXTERNAL f
PARAMETER (CGOLD=0.38196601125011d0)
PARAMETER (ITMAX=50, ZEPS=1.0e-10)
c ... Given a function f, and given a bracketing triplet of
abscissas ax, bx, cx
c ... (such that bx is between ax and cx, and f(bx) is greater than
both f(ax) and f(cx)),
c ... this routine isolates the maximum to a fractional precision of
about tol using Brent's method.
c ... The abscissa of the maximum is returned as xmax, and the
maximum function value is
c ... returned as brent, the returned function value. Parameters:
Maximum allowed number of
c ... iterations;
c ... golden ratio; and a small number that protects against trying
to achieve fractional accuracy
c ... for a maximum that happens to be exactly zero.
INTEGER iter
DOUBLE PRECISION a,b,d,e,etemp,fu,fv,fw,fx,p,q,r
DOUBLE PRECISION tol1,tol2,u,v,w,x,xm
a=DMIN1(ax,cx)
c a and b must be in ascending order, though the input abscissas
need not be.
b=DMAX1(ax,cx)
v=bx
c ... Initializations...
w=v
x=v
e=0.0d0
c ,,, This will be the distance moved on the step before last.
fx=f(x)
fv=fx
fw=fx
DO iter=1, ITMAX
c ... Main program loop.
xm=0.5d0*(a+b)
tol1=tol*DABS(x)+DBLE(ZEPS)
tol2=2.0d0*tol1
IF (DABS(x-xm).LE.(tol2-0.5d0*(b-a))) GOTO 3
c ... Test for done here.
IF (DABS(e).GT.tol1) THEN
c ... Construct a trial parabolic-fit.
r=(x-w)*(fx-fv)
q=(x-v)*(fx-fw)
p=(x-v)*q-(x-w)*r
q=2.*(q-r)
IF (q.GT.0.0d0) p=-p
q=DABS(q)
etemp=e
e=d
IF (DABS(p).GE.DABS(0.5d0*q*etemp).OR.p.LE.
& q*(a-x).OR.p.GE.q*(b-x)) GOTO 1
c ... The above conditions determine the acceptability of the
parabolic fit. Here it is o.k.:
d=p/q
c ... Take the parabolic step.
u=x+d
IF (u-a.LT.tol2 .OR. b-u.LT.tol2) d=DSIGN(tol1,xm-x)
GOTO 2
c ... Skip over the golden section step.
ENDIF
1 IF (x.GE.xm) THEN
c ... We arrive here for a golden section step, which we take into
the larger of the two segments.
e=a-x
ELSE
e=b-x
ENDIF
d=CGOLD * e
c ... Take the golden section step.
2 IF (DABS(d).GE.tol1) THEN
c ... Arrive here with d computed either from parabolic fit, or else
from golden section.
u=x+d
ELSE
u=x+DSIGN(tol1,d)
ENDIF
fu=f(u)
c ... This is the one function evaluation per iteration,
IF (fu.GT.fx) THEN
c ... and now we have to decide what to do with our function
evaluation. Housekeeping follows:
IF (u.GE.x) THEN
a=x
ELSE
b=x
ENDIF
v=w
fv=fw
w=x
fw=fx
x=u
fx=fu
ELSE
IF (u.LT.x) THEN
a=u
ELSE
b=u
ENDIF
IF (fu.GT.fw .OR. w.EQ.x) THEN
v=w
fv=fw
w=u
fw=fu
ELSE IF (fu.GT.fv .OR. v.EQ.x .OR. v.EQ.w) THEN
v=u
fv=fu
ENDIF
ENDIF
c ... Done with housekeeping. Back for another iteration.
END DO
PAUSE 'brent exceed maximum iterations'
3 xmax=x
c ... Arrive here ready to exit with best values.
brent=fx
RETURN
END
c ...
FUNCTION dfx1di(x)
DOUBLE PRECISION dfx1di, x
c ... Used by doptxp as the function passed to mnbrak and brent.
DOUBLE PRECISION PI, LOG8, LOG4
PARAMETER (PI = 3.14159265358979d0)
PARAMETER (LOG8 = 2.07944154154168d0)
PARAMETER (LOG4 = 1.38629436111989d0)
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
c DOUBLE PRECISION dirlnG
c ...
c ... dirlnG= 2.0d0*x**2 - DLOG((omega*x)**2)+LOG4
c ... dfx1di= PI*(1+DEXP(-x**2))*omega**2 + dirlnG
dfx1di = 0.25d0*(omega*x)**2*
& DEXP(-PI*(1+DEXP(-x**2))*omega**2-2.0d0*x**2)
RETURN
END
c ...
FUNCTION fx1dim(x)
DOUBLE PRECISION fx1dim, x
c ... Used by optxps as the function passed to mnbrak and brent.
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
CALL initc(x)
CALL gMVmn()
c ... compute R with minimal computational work
CALL EvalR()
fx1dim = ExpmR
RETURN
END
c ...
FUNCTION f1dim(x)
DOUBLE PRECISION f1dim, x
c ... Used by linmin as the function passed to mnbrak and brent.
INTEGER NP
PARAMETER (NP = 40)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
DOUBLE PRECISION Mxi(NP,NP), Mg(NP,NP), Mh(NP,NP)
DOUBLE PRECISION Mpcom(NP, NP), Mxicom(NP, NP)
COMMON /gradmethod/ Mxi, Mg, Mh, Mpcom, Mxicom
c ...
INTEGER j, k, switch
c ...
c ... PRINT '( "f1dim: called with x=",E20.12 )',x
DO j=1, NP
switch = 0
DO k=j+1, NP
IF (switch.EQ.0) THEN
switch = 1
MVmn(j,k) = Mpcom(j,k) + x * Mxicom(j,k)
ELSE
switch = 0
ENDIF
END DO
END DO
c ... PRINT '( "f1dim: evaluate R" )'
c ... compute R with minimal computational work
CALL EvalR()
f1dim = ExpmR
RETURN
END
c ...
SUBROUTINE EvalR()
c ... evaluate R with minimal computation a given point MVmn
INTEGER NP
DOUBLE PRECISION PI, SQRT2
PARAMETER (NP = 40, PI = 3.14159265358979d0,
& SQRT2 = 1.41421356237309d0)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
DOUBLE PRECISION weight, ovlap
EXTERNAL weight, ovlap
c ...
c ... obtain the interaction matrix MVjk
CALL gMVjk()
c ... calculate largest eigenvalue with eigenvector
CALL diagV()
c ... PRINT '( "EvalR: eigenvalue ",E20.12 )', eigenv(NP)
c ... calculate \int V^2/2 d{\bf r}
d = weight() / (2.0d0*eigenv(NP))**2
c ... PRINT '( "EvalR: weight ",E20.12 )', d
c ... calculate overlap sqrt{2} \sum_j (-1)^j \sqrt{j+1} Re{c_j
c_{j+1}^\ast}
overl = ovlap() * SQRT2
c ... PRINT '( "EvalR: overl ",E20.12 )', overl
c ... the functional value R, to be minimized
FuncR = d*omega**2 - DLOG((overl*omega)**2)
c ... the exponential of the functional to be maximized
ExpmR = (overl*omega)**2 * DEXP(-d*omega**2)
RETURN
END
c ...
SUBROUTINE Evalal()
c ... evaluate R and all eigenvalues and vectors at a given point
MVmn
INTEGER NP
DOUBLE PRECISION PI, SQRT2
PARAMETER (NP = 40, PI = 3.14159265358979d0,
& SQRT2 = 1.41421356237309d0)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
DOUBLE PRECISION weight, ovlap
EXTERNAL weight, ovlap
c ...
c ... obtain the interaction matrix MVjk
CALL gMVjk()
c ... calculate all its eigenvectors and eigenvalues
c PRINT '( "Evalal: call diagVa()" )'
CALL diagVa()
c ... calculate \int V^2/2 d{\bf r}
c PRINT '( "Evalal: call weight()" )'
d = weight() / (4.0d0*eigenv(NP)*eigenv(NP))
c ... calculate overlap \sqrt{2} \sum_j (-1)^j \sqrt{j+1} Re{c_j
c_{j+1}^\ast}
overl = SQRT2 * ovlap()
PRINT '( "Here it goes well: ",E14.5 )', overl
PRINT '( "Now it produces a core dump: ",E14.5 )',
& overl*overl
PRINT '( "got through!" )'
c ... the functional value R, to be minimized
FuncR = d*omega**2 - DLOG((overl*omega)**2)
c ... the exponential of the functional to be maximized
ExpmR = (overl*omega)**2 * DEXP(-d*omega**2)
c PRINT '( "exit Evalal" )'
RETURN
END
c ...
SUBROUTINE gradR()
c ... having R evaluated, return grad R at a given point MVmn
INTEGER NP
DOUBLE PRECISION SQRT2
PARAMETER (NP = 40)
PARAMETER (SQRT2 = 1.41421356237309d0)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION d, overl, omega, FuncR, ExpmR
COMMON /minimize/ d, overl, omega, FuncR, ExpmR
c ...
DOUBLE PRECISION Mxi(NP,NP), Mg(NP,NP), Mh(NP,NP)
DOUBLE PRECISION Mpcom(NP, NP), Mxicom(NP, NP)
COMMON /gradmethod/ Mxi, Mg, Mh, Mpcom, Mxicom
c ...
DOUBLE PRECISION pred,preE,preov
INTEGER j,k,switch
c ... here we calculate its gradient at point MVmn
CALL dEvpq()
CALL dVvpq()
CALL dovpq()
c ... finally the gradient
pred = (omega/(2.0d0*eigenv(NP)))**2
preE = 2.0d0 * d * omega**2 / eigenv(NP)
preov = 2.0d0 * SQRT2 / overl
DO j=1, NP
switch = 0
DO k=j+1, NP
IF (switch.EQ.0) THEN
switch = 1
Mxi(j,k) = -preE * MdEvpq(j,k) + pred * MdVvpq(j,k)
& -preov * Mdovpq(j,k)
ELSE
switch = 0
ENDIF
END DO
END DO
RETURN
END
c ...
SUBROUTINE diagV()
c ... given the interaction potential MVjk, calculate the largest
eigenvalue
c ... with corresponding normalized eigenvector
INTEGER NP
PARAMETER (NP = 40)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION work(8*NP),Z(NP,1)
INTEGER iwork(5*NP), ifail(NP)
INTEGER j,info,M
c PRINT '( "diagonalize DSYEVX" )'
c ... call lapack routine to obtain largest eigenvalue and
eigenvector
CALL DSYEVX('V','I','U',NP,MVjk,NP,0.0d0,0.0d0,NP,NP,
& DBLE(1.e-6),M,eigenv,Z,NP,work,8*NP,iwork,ifail,info)
c ... info is zero ?
IF (info.NE.0) THEN
IF (info.GT.0) THEN
PAUSE 'eigenvector failed to converge in zheevx'
ELSE
STOP 'argument error in zheevx'
ENDIF
ENDIF
c ... trimm the solution cut off too small values < 10^-13
eigenv(NP) = eigenv(1)
DO j=1, NP
ck(j,NP) = Z(j,1)
END DO
RETURN
END
c ...
SUBROUTINE diagVa()
c ... given the interaction potential MVjk, calculate all
eigenvectors, eigenvalues
c ... returns eigenvalues and orthonormal eigenvectors in ascending
order
INTEGER NP
PARAMETER (NP = 40)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION work(1+5*NP+3*NP**2+2*NP*8)
INTEGER iwork(2+5*NP)
INTEGER info
c PRINT '( "diagonalize DSYEVD" )'
c ... call lapack routine to diagonolize
CALL DSYEVD('V','U',NP,MVjk,NP,eigenv,work,1+5*NP+3*NP**2+2*NP*8,
& iwork,2+5*NP,info)
c ... info is zero ?
IF (info.NE.0) THEN
IF (info.GT.0) THEN
PAUSE 'eigenvector failed to converge in zheevd'
ELSE
STOP 'argument error in zheevd'
ENDIF
ENDIF
c PRINT '( "DSYEVD: successful" )'
RETURN
END
c ...
SUBROUTINE initc(x)
c ... given the initial x position, calculate the coefficient vector
c
DOUBLE PRECISION x
c ...
INTEGER NP
PARAMETER (NP = 40)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
INTEGER j
DOUBLE PRECISION pref
DOUBLE PRECISION alp0, alpha
pref = DEXP(-x**2/4.0d0)
alp0 = x/DSQRT(2.0d0)
alpha = 1.0d0
DO j=1, NP
c(j) = pref * alpha
alpha = alpha * alp0 / DSQRT(DBLE(j))
END DO
RETURN
END
c ...
SUBROUTINE gMVmn()
c ... given the vector c, calculate MVmn(NP,NP)
c ...
INTEGER NP
PARAMETER (NP = 40)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
INTEGER j,k,switch
DO j=1, NP
switch = 0
DO k=j+1, NP
IF (switch.EQ.0) THEN
switch = 1
MVmn(j,k) = c(j)*c(k)*2.0d0
ELSE
switch = 0
ENDIF
END DO
END DO
RETURN
END
c ...
SUBROUTINE gMVjk()
c ... given the matrix MVmn, calculate the interaction matrix
MVjk(NP,NP)
c ... each time, the even off diagonal elements have to be zeroed,
c ... since diagonalization routine uses this as output array!
c ...
INTEGER NP
DOUBLE PRECISION PI, D4PI
PARAMETER (NP = 40, PI = 3.14159265358979d0)
PARAMETER (D4PI = 1.0d0/(4.0d0*PI))
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION psqrt(NP), pbinom(NP,NP,NP)
COMMON /precalc/ psqrt, pbinom
c ...
DOUBLE PRECISION summ
INTEGER j,k,m,switch
c ...
DO k=0, NP-1
switch = 0
DO j=k+1, NP-1
IF (switch.EQ.0) THEN
switch = 1
summ = MVmn(1,j-k+1) * pbinom(k+1,j+1,1)
DO m=1, NP-1-(j-k)
summ = summ + MVmn(m+1,m+j-k+1)*pbinom(k+1,j+1,m+1)
END DO
MVjk(k+1,j+1) = summ * D4PI
ELSE
switch = 0
MVjk(k+1, j+1) = 0.0d0
ENDIF
END DO
END DO
c ... clear diagonal
DO j=1, NP
MVjk(j,j) = 0.0d0
END DO
RETURN
END
c ...
FUNCTION ovlap()
DOUBLE PRECISION ovlap
c ... calculate the overlap integral
c ... \sum_j (-1)^j \sqrt{j+1} Re{c_j c_{j+1}^\ast}
c ...
INTEGER NP
PARAMETER (NP = 40)
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION psqrt(NP), pbinom(NP,NP, NP)
COMMON /precalc/ psqrt, pbinom
c ...
INTEGER j,s
DOUBLE PRECISION sumj
sumj = 0.0d0
s = 1
DO j=0, NP-2
IF (s.EQ.1) THEN
sumj = sumj + psqrt(j+1)
& * ck(j+1,NP)*ck(j+2,NP)
ELSE
sumj = sumj - psqrt(j+1)
& * ck(j+1,NP)*ck(j+2,NP)
ENDIF
s = -s
END DO
ovlap = sumj
RETURN
END
c ...
FUNCTION weight()
DOUBLE PRECISION weight
c ... calculate the weight of MVmn
c ...
INTEGER NP
DOUBLE PRECISION PI, D4PI
PARAMETER (NP = 40, PI = 3.14159265358979d0)
PARAMETER (D4PI = 1.0d0/(4.0d0*PI))
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION psqrt(NP), pbinom(NP,NP,NP)
COMMON /precalc/ psqrt, pbinom
c ...
INTEGER m,n,j,switch
DOUBLE PRECISION sum
c ...
sum = 0.0d0
DO m=0, NP-1
switch = 0
DO n=m+1, NP-1
IF (switch.EQ.0) THEN
switch = 1
c ... inner loop sum over j = 0 .. NP-1
DO j=0, NP-1-(n-m)
sum = sum + MVmn(m+1,n+1)*
& MVmn(j+1,j+n-m+1) * pbinom(m+1,n+1,j+1)
END DO
ELSE
switch = 0
ENDIF
END DO
END DO
weight=sum * D4PI
RETURN
END
c ...
SUBROUTINE dEvpq()
c ... given the eigenvector c_j^{(0)}, calculate derivative of E to
v_{p,q}
c ... and store the values in MdEvpq(NP,NP)
c ...
INTEGER NP
DOUBLE PRECISION PI, D2PI
PARAMETER (NP = 40, PI = 3.14159265358979d0)
PARAMETER (D2PI = 1.0d0/(2.0d0*PI))
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION psqrt(NP), pbinom(NP,NP,NP)
COMMON /precalc/ psqrt, pbinom
c ...
DOUBLE PRECISION sumj
INTEGER j,p,q,switch
c ...
DO p=0, NP-1
switch = 0
DO q=p+1, NP-1
IF (switch.EQ.0) THEN
switch = 1
c ... inner loop sum over j = 0 .. NP-1
sumj = ck(q-p+1,NP)*ck(1,NP) * pbinom(p+1,q+1,1)
DO j=1, NP-1-(q-p)
sumj = sumj + ck(q-p+j+1,NP)*ck(j+1,NP)*
& pbinom(p+1,q+1,j+1)
END DO
MdEvpq(p+1, q+1) = sumj * D2PI
ELSE
switch = 0
ENDIF
END DO
END DO
RETURN
END
c ...
SUBROUTINE dVvpq()
c ... given the matrix V_mn, calculate derivative of \int V^2/2 to
v_{p,q}
c ... and store the values in MdVvpq(NP,NP)
c ...
INTEGER NP
DOUBLE PRECISION PI, D2PI
PARAMETER (NP = 40, PI = 3.14159265358979d0)
PARAMETER (D2PI = 1.0d0/(2.0d0*PI))
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION psqrt(NP), pbinom(NP,NP,NP)
COMMON /precalc/ psqrt, pbinom
c ...
DOUBLE PRECISION sumj
INTEGER j,p,q,switch
c ...
DO p=0, NP-1
switch = 0
DO q=p+1, NP-1
IF (switch.EQ.0) THEN
switch = 1
c ... inner loop sum over j = 0 .. NP-1
sumj = MVmn(1,q-p+1) * pbinom(p+1,q+1,1)
DO j=1, NP-1-(q-p)
sumj = sumj + MVmn(j+1,q-p+j+1) *
& pbinom(p+1,q+1,j+1)
END DO
MdVvpq(p+1, q+1) = sumj * D2PI
ELSE
switch = 0
ENDIF
END DO
END DO
RETURN
END
c ...
SUBROUTINE dovpq()
c ... given all eigenvectors c_j^{(k)} and its eigenvalues E^{(k)}
c ... calculate the derivative of 1/(2\sqrt{2}E) \int \partial V/
\partial x
c ... \psi_t^\ast \psi_b d{\bf r} to v_{p,q} and store the values
in Mdovpq(NP,NP)
c ...
INTEGER NP
DOUBLE PRECISION PI, D4PI
PARAMETER (NP = 40, PI = 3.14159265358979d0)
PARAMETER (D4PI = 1.0d0/(4.0d0*PI))
DOUBLE PRECISION MVmn(NP, NP), MVjk(NP, NP), MdEvpq(NP,NP)
DOUBLE PRECISION MdVvpq(NP,NP), Mdovpq(NP,NP), ck(NP, NP), c(NP)
EQUIVALENCE (MVjk, ck)
DOUBLE PRECISION eigenv(NP)
COMMON /matrix/ MVmn, MVjk, MdEvpq, MdVvpq, Mdovpq, c, eigenv
c ...
DOUBLE PRECISION psqrt(NP), pbinom(NP,NP,NP)
COMMON /precalc/ psqrt, pbinom
c ...
DOUBLE PRECISION sumj, sumn, ovarr(NP-1)
DOUBLE PRECISION sumk
INTEGER j,k,p,q,n,switch
c ...
DO k=1, NP-1
switch = 1
sumn = 0.0d0
DO n=0, NP-2
IF (switch.EQ.1) THEN
sumn = sumn + psqrt(n+1)
& * ( ck(n+1,k)*ck(n+2,NP) + ck(n+1,NP)*ck(n+2,k) )
ELSE
sumn = sumn - psqrt(n+1)
& * ( ck(n+1,k)*ck(n+2,NP) + ck(n+1,NP)*ck(n+2,k) )
ENDIF
switch = -switch
END DO
ovarr(k) = sumn / (eigenv(NP) - eigenv(k))
END DO
c ...
DO p=0, NP-1
switch = 0
DO q=p+1, NP-1
IF (switch.EQ.0) THEN
switch = 1
sumk = 0.0d0
DO k=1, NP-1
c ... inner loop sum over j = 0 .. NP-1
sumj = (ck(q-p+1,k)*ck(1,NP)
& + ck(1,k)*ck(q-p+1,NP)) * pbinom(p+1,q+1,1)
DO j=1, NP-1-(q-p)
sumj = sumj + (ck(q-p+j+1,k) * ck(j+1,NP)
& + ck(j+1,k) * ck(q-p+j+1,NP)) * pbinom(p+1,q+1,j+1)
END DO
sumk = sumk + ovarr(k) * sumj
END DO
Mdovpq(p+1, q+1) = sumk * D4PI
ELSE
switch = 0
ENDIF
END DO
END DO
RETURN
END