mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-28 22:25:48 -04:00
Merge branch 'master' of https://github.com/nwchemgit/nwchem
This commit is contained in:
Mark Driver
commit
1c72592ee9
37 changed files with 1122 additions and 17819 deletions
|
|
@ -42,7 +42,7 @@ python noprint
|
|||
def printexp(z):
|
||||
print("\n Exponents:")
|
||||
for i in range(len(z)):
|
||||
print " %14.8f" % (z[i]*z[i]),
|
||||
print (" %14.8f" % (z[i]*z[i])),
|
||||
if ((i+1)%5) == 0:
|
||||
print("")
|
||||
print(" ")
|
||||
|
|
|
|||
|
|
@ -65,7 +65,7 @@ python noprint
|
|||
def printexp(z):
|
||||
print("\n Exponents:")
|
||||
for i in range(len(z)):
|
||||
print " %14.8f" % (z[i]*z[i]),
|
||||
print ( " %14.8f" % (z[i]*z[i])),
|
||||
if ((i+1)%5) == 0:
|
||||
print("")
|
||||
print(" ")
|
||||
|
|
|
|||
|
|
@ -8,7 +8,7 @@ def zerovector(n):
|
|||
return a
|
||||
|
||||
def zeromatrix(n,m):
|
||||
a = range(n)
|
||||
a = list(range(n))
|
||||
for i in range(n):
|
||||
a[i] = zerovector(m)
|
||||
return a
|
||||
|
|
@ -152,10 +152,17 @@ def quadratic_step(trust, g0, h0):
|
|||
if h0 > 0:
|
||||
delta2 = -g0/h0
|
||||
if abs(delta2) > trust:
|
||||
<<<<<<< HEAD
|
||||
print(" Step restriction: %f %f " % (delta2, trust))
|
||||
delta2 = abs(trust*delta2)/delta2
|
||||
else:
|
||||
print(" Negative curvature ")
|
||||
=======
|
||||
print (" Step restriction: %12.5f %12.5f " % (delta2, trust))
|
||||
delta2 = abs(trust*delta2)/delta2
|
||||
else:
|
||||
print (" Negative curvature ")
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
delta2 = -abs(trust*g0)/g0
|
||||
return delta2
|
||||
|
||||
|
|
@ -169,24 +176,42 @@ def linesearch(func, x0, s, lsgrad, eps):
|
|||
# bracketed the minimum or gone downhil with enough
|
||||
# energy difference to start fitting
|
||||
|
||||
<<<<<<< HEAD
|
||||
print(" Line search: step alpha grad hess value")
|
||||
print(" ---- --------- -------- -------- ----------------")
|
||||
=======
|
||||
print (" Line search: step alpha grad hess value")
|
||||
print (" ---- --------- -------- -------- ----------------")
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
|
||||
trust = 0.2
|
||||
alpha0 = 0.0
|
||||
f0 = func(x0)
|
||||
<<<<<<< HEAD
|
||||
print(" %9.2e %8.1e %16.8f" % (alpha0, lsgrad, f0))
|
||||
=======
|
||||
print (" %9.2e %8.1e %16.8f" % (alpha0 , lsgrad , f0))
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
|
||||
if lsgrad < 0:
|
||||
alpha1 = alpha0 + trust
|
||||
else:
|
||||
alpha1 = alpha0 - trust
|
||||
f1 = func(takestep(x0,s,alpha1))
|
||||
<<<<<<< HEAD
|
||||
print(" %9.2e %16.8f" % (alpha1, f1))
|
||||
|
||||
while f1 > f0:
|
||||
if trust < 0.00125:
|
||||
print(" system is too badly conditioned for initial step")
|
||||
=======
|
||||
print (" %9.2e %16.8f" % \
|
||||
(alpha1, f1))
|
||||
|
||||
while f1 > f0:
|
||||
if trust < 0.00125:
|
||||
print (" system is too badly conditioned for initial step")
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
return (alpha0,f0) # Cannot seem to find my way
|
||||
trust = trust * 0.5
|
||||
if lsgrad < 0:
|
||||
|
|
@ -194,7 +219,12 @@ def linesearch(func, x0, s, lsgrad, eps):
|
|||
else:
|
||||
alpha1 = alpha0 - trust
|
||||
f1 = func(takestep(x0,s,alpha1))
|
||||
<<<<<<< HEAD
|
||||
print(" %9.2e %16.8f" % (alpha1, f1))
|
||||
=======
|
||||
print (" %9.2e %16.8f" % \
|
||||
(alpha1, f1))
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
|
||||
g0 = lsgrad
|
||||
h0 = (f1-f0-alpha1*g0)/alpha1**2
|
||||
|
|
@ -214,20 +244,35 @@ def linesearch(func, x0, s, lsgrad, eps):
|
|||
if iter == 1:
|
||||
f2prev = f2
|
||||
|
||||
<<<<<<< HEAD
|
||||
print(" %9.2e %16.8f" % (alpha2, f2))
|
||||
|
||||
# Check for convergence or insufficient precision to proceed further
|
||||
if (abs(f0-f1)<(10*eps)) and (abs(f1-f2)<(10*eps)):
|
||||
print(" ")
|
||||
print(" Insufficient precision ... terminating LS")
|
||||
=======
|
||||
print (" %9.2e %16.8f" % \
|
||||
(alpha2, f2))
|
||||
|
||||
# Check for convergence or insufficient precision to proceed further
|
||||
if (abs(f0-f1)<(10*eps)) and (abs(f1-f2)<(10*eps)):
|
||||
print (" "),
|
||||
print (" Insufficient precision ... terminating LS")
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
break
|
||||
|
||||
if (f2-f2prev) > 0:
|
||||
# New point is higher than previous worst
|
||||
if nbackstep < 3:
|
||||
nbackstep = nbackstep + 1
|
||||
<<<<<<< HEAD
|
||||
print(" ")
|
||||
print(" Back stepping due to uphill step")
|
||||
=======
|
||||
print (" "),
|
||||
print (" Back stepping due to uphill step")
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
trust = max(0.01,0.2*abs(alpha2 - alpha0)) # Reduce trust radius
|
||||
alpha2 = alpha0 + 0.2*(alpha2 - alpha0)
|
||||
continue
|
||||
|
|
@ -247,11 +292,20 @@ def linesearch(func, x0, s, lsgrad, eps):
|
|||
alpha1, alpha2, f1, f2 = alpha2, alpha1, f2, f1
|
||||
|
||||
(f0, g0, h0) = quadfit(alpha0, f0, alpha1, f1, alpha2, f2)
|
||||
<<<<<<< HEAD
|
||||
print(" %4i %9.2e %8.1e %8.1e %16.8f" %(iter, alpha0, g0, h0, f0))
|
||||
|
||||
if (h0>0.0) and (abs(g0) < 0.03*abs(lsgrad)):
|
||||
print(" ")
|
||||
print(" gradient reduced 30-fold ... terminating LS")
|
||||
=======
|
||||
print (" %4i %9.2e %8.1e %8.1e %16.8f" % \
|
||||
(iter, alpha0, g0, h0, f0))
|
||||
|
||||
if (h0>0.0) and (abs(g0) < 0.03*abs(lsgrad)):
|
||||
print (" "),
|
||||
print (" gradient reduced 30-fold ... terminating LS")
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
break
|
||||
|
||||
# Determine the next step
|
||||
|
|
@ -259,8 +313,13 @@ def linesearch(func, x0, s, lsgrad, eps):
|
|||
alpha2 = alpha0 + delta
|
||||
df = g0*delta + 0.5*h0*delta*delta
|
||||
if abs(df) < 10.0*eps:
|
||||
<<<<<<< HEAD
|
||||
print(" ")
|
||||
print(" projected energy reduction < 10*eps ... terminating LS")
|
||||
=======
|
||||
print (" "),
|
||||
print (" projected energy reduction < 10*eps ... terminating LS")
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
break
|
||||
|
||||
|
||||
|
|
@ -399,7 +458,11 @@ def hessian_update_bfgs(hp, dx, g, gp):
|
|||
for j in range(n):
|
||||
h[i][j] = h[i][j] + dg[i]*dg[j]/dxdg - hdx[i]*hdx[j]/dxhdx
|
||||
else:
|
||||
<<<<<<< HEAD
|
||||
print(' BFGS not updating dxdg (%e), dgdg (%e), dxhdx (%f), dxdx(%e)' % (dxdg, dgdg, dxhdx, dxdx))
|
||||
=======
|
||||
print (' BFGS not updating dxdg (%e), dgdg (%e), dxhdx (%f), dxdx(%e)' % (dxdg, dgdg, dxhdx, dxdx))
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
|
||||
return h
|
||||
|
||||
|
|
@ -437,16 +500,27 @@ def quasinr(func, guess, tol, eps, printvar=None):
|
|||
(value,g,h) = numderiv(func, x, step, eps)
|
||||
gmax = max(map(abs,g))
|
||||
|
||||
<<<<<<< HEAD
|
||||
print(' ')
|
||||
print(' iter gmax value ')
|
||||
print(' ---- --------- ----------------')
|
||||
print("%4i %9.2e %16.8f" % (iter,gmax,value))
|
||||
=======
|
||||
print (' ')
|
||||
print (' iter gmax value ')
|
||||
print (' ---- --------- ----------------')
|
||||
print ("%4i %9.2e %16.8f" % (iter,gmax,value))
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
|
||||
if (printvar):
|
||||
printvar(x)
|
||||
|
||||
if gmax < tol:
|
||||
<<<<<<< HEAD
|
||||
print(' Converged!')
|
||||
=======
|
||||
print (' Converged!')
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
break
|
||||
|
||||
if iter == 0:
|
||||
|
|
@ -458,14 +532,22 @@ def quasinr(func, guess, tol, eps, printvar=None):
|
|||
(v,e) = jacobi(hessian)
|
||||
emax = max(map(abs,e))
|
||||
emin = emax*1e-4 # Control noise in small eigenvalues
|
||||
<<<<<<< HEAD
|
||||
print('\n Eigenvalues of the Hessian:')
|
||||
=======
|
||||
print ('\n Eigenvalues of the Hessian:')
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
printvector(e)
|
||||
|
||||
# Transform to spectral form, take step, transform back
|
||||
gs = mxv(v,g)
|
||||
for i in range(n):
|
||||
if e[i] < emin:
|
||||
<<<<<<< HEAD
|
||||
print(' Mode %d: small/negative eigenvalue (%f).' % (i, e[i]))
|
||||
=======
|
||||
print (' Mode %d: small/negative eigenvalue (%f).' % (i, e[i]))
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
s[i] = -gs[i]/emin
|
||||
else:
|
||||
s[i] = -gs[i]/e[i]
|
||||
|
|
@ -476,7 +558,12 @@ def quasinr(func, guess, tol, eps, printvar=None):
|
|||
for i in range(n):
|
||||
trust = max(abs(x[i]),abs(x[i]/sqrt(max(1e-4,abs(hessian[i][i])))))
|
||||
if abs(s[i]) > trust:
|
||||
<<<<<<< HEAD
|
||||
print(' restricting ', i, trust, abs(x[i]), abs(x[i]/sqrt(abs(hessian[i][i]))), s[i])
|
||||
=======
|
||||
print (' restricting ', i, trust, abs(x[i]), \
|
||||
abs(x[i]/sqrt(abs(hessian[i][i]))), s[i])
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
scale = min(scale,trust/abs(s[i]))
|
||||
if scale != 1.0:
|
||||
for i in range(n):
|
||||
|
|
@ -485,7 +572,11 @@ def quasinr(func, guess, tol, eps, printvar=None):
|
|||
(alpha,value) = linesearch(func, x, s, dot(s,g), eps)
|
||||
|
||||
if alpha == 0.0:
|
||||
<<<<<<< HEAD
|
||||
print(' Insufficient precision to proceed further')
|
||||
=======
|
||||
print (' Insufficient precision to proceed further')
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
break
|
||||
|
||||
for i in range(n):
|
||||
|
|
@ -510,6 +601,7 @@ def cgminold(func, dfunc, guess, tol):
|
|||
g = dfunc(x)
|
||||
gmax = max(map(abs,g))
|
||||
|
||||
<<<<<<< HEAD
|
||||
print(' ')
|
||||
print(' iter gmax value ')
|
||||
print(' ---- --------- ----------------')
|
||||
|
|
@ -517,6 +609,15 @@ def cgminold(func, dfunc, guess, tol):
|
|||
|
||||
if gmax < tol:
|
||||
print(' Converged!')
|
||||
=======
|
||||
print (' ')
|
||||
print (' iter gmax value ')
|
||||
print (' ---- --------- ----------------')
|
||||
print ("%4i %9.2e %16.8f" % (iter,gmax,value))
|
||||
|
||||
if gmax < tol:
|
||||
print (' Converged!')
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
break
|
||||
|
||||
if (iter == 0) or ((iter%20) == 0):
|
||||
|
|
@ -557,6 +658,7 @@ def cgmin(func, dfunc, guess, tol, precond=None, reset=None):
|
|||
g = dfunc(x)
|
||||
gmax = max(map(abs,g))
|
||||
|
||||
<<<<<<< HEAD
|
||||
print(' ')
|
||||
print(' iter gmax value ')
|
||||
print(' ---- --------- ----------------')
|
||||
|
|
@ -564,6 +666,15 @@ def cgmin(func, dfunc, guess, tol, precond=None, reset=None):
|
|||
|
||||
if gmax < tol:
|
||||
print(' Converged!')
|
||||
=======
|
||||
print (' ')
|
||||
print (' iter gmax value ')
|
||||
print (' ---- --------- ----------------')
|
||||
print ("%4i %9.2e %16.8f" % (iter,gmax,value))
|
||||
|
||||
if gmax < tol:
|
||||
print (' Converged!')
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
break
|
||||
|
||||
if precond:
|
||||
|
|
@ -631,16 +742,27 @@ def cgmin2(func, guess, tol, eps, printvar=None,reset=None):
|
|||
(value,g,hh) = numderiv(func, x, step, eps)
|
||||
gmax = max(map(abs,g))
|
||||
|
||||
<<<<<<< HEAD
|
||||
print(' ')
|
||||
print(' iter gmax value ')
|
||||
print(' ---- --------- ----------------')
|
||||
print("%4i %9.2e %16.8f" % (iter,gmax,value))
|
||||
=======
|
||||
print (' ')
|
||||
print (' iter gmax value ')
|
||||
print (' ---- --------- ----------------')
|
||||
print ("%4i %9.2e %16.8f" % (iter,gmax,value))
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
|
||||
if (printvar):
|
||||
printvar(x)
|
||||
|
||||
if gmax < tol:
|
||||
<<<<<<< HEAD
|
||||
print(' Converged!')
|
||||
=======
|
||||
print (' Converged!')
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
break
|
||||
|
||||
if (iter % reset) == 0:
|
||||
|
|
@ -671,12 +793,20 @@ def cgmin2(func, guess, tol, eps, printvar=None,reset=None):
|
|||
if alpha == 0.0:
|
||||
# LS failed, probably due to lack of precision.
|
||||
if beta != 0.0:
|
||||
<<<<<<< HEAD
|
||||
print("LS failed - trying preconditioned steepest descent direction")
|
||||
=======
|
||||
print ("LS failed - trying preconditioned steepest descent direction")
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
for i in range(n):
|
||||
s[i] = -g[i]
|
||||
(alpha,value) = linesearch(func, x, s, dot(s,g), eps)
|
||||
if alpha == 0.0:
|
||||
<<<<<<< HEAD
|
||||
print(" Insufficient precision to proceed further")
|
||||
=======
|
||||
print (" Insufficient precision to proceed further")
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
break
|
||||
|
||||
for i in range(n):
|
||||
|
|
@ -713,6 +843,7 @@ if __name__ == '__main__':
|
|||
return sum
|
||||
|
||||
|
||||
<<<<<<< HEAD
|
||||
print('\n\n TESTING QUASI-NR SOLVER \n\n')
|
||||
quasinr(f, [1.,0.5,0.3,-0.4], 1e-4, 1e-10)
|
||||
|
||||
|
|
@ -729,6 +860,24 @@ if __name__ == '__main__':
|
|||
cgminold(f, df, [1.,0.5,0.3,-0.4], 1e-4)
|
||||
|
||||
print('\n\n TESTING JACOBI EIGENSOLVER\n\n')
|
||||
=======
|
||||
print ('\n\n TESTING QUASI-NR SOLVER \n\n')
|
||||
quasinr(f, [1.,0.5,0.3,-0.4], 1e-4, 1e-10)
|
||||
|
||||
print ('\n\n TESTING GC WITH NUM. GRAD. AND DIAG. PRECOND.\n\n')
|
||||
cgmin2(f, [1.,0.5,0.3,-0.4], 1e-4, 1e-10, reset=20)
|
||||
|
||||
print ('\n\n TESTING GC WITH ANAL. GRAD. AND WITHOUT OPTIONAL PRECOND.\n\n')
|
||||
cgmin(f, df, [1.,0.5,0.3,-0.4], 1e-4)
|
||||
|
||||
print ('\n\n TESTING GC WITH ANAL. GRAD. AND WITH OPTIONAL PRECOND.\n\n')
|
||||
cgmin(f, df, [1.,0.5,0.3,-0.4], 1e-4, precond=precond)
|
||||
|
||||
print ('\n\n TESTING GC WITH ANAL. GRAD. AND NO PRECOND.\n\n')
|
||||
cgminold(f, df, [1.,0.5,0.3,-0.4], 1e-4)
|
||||
|
||||
print ('\n\n TESTING JACOBI EIGENSOLVER\n\n')
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
n = 50
|
||||
a = zeromatrix(n,n)
|
||||
for i in range(n):
|
||||
|
|
@ -737,7 +886,11 @@ if __name__ == '__main__':
|
|||
|
||||
(v,e)= jacobi(a)
|
||||
|
||||
<<<<<<< HEAD
|
||||
print(' eigenvalues')
|
||||
=======
|
||||
print (' eigenvalues')
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
printvector(e)
|
||||
#print ' v '
|
||||
#printmatrix(v)
|
||||
|
|
@ -748,4 +901,8 @@ if __name__ == '__main__':
|
|||
err = max(err,abs(ev[i][j] - e[i]*v[i][j]))
|
||||
err = err/(n*max(1.0,abs(e[i])))
|
||||
if err > 1e-12:
|
||||
<<<<<<< HEAD
|
||||
print(' Error in eigenvector ', i, err)
|
||||
=======
|
||||
print (' Error in eigenvector ', i, err)
|
||||
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue