Remove czher

parent a749cfaf
Pipeline #105076713 passed with stage
in 3 minutes and 31 seconds
......@@ -12,7 +12,7 @@ plat = get_platform()
platform_id = os.getenv('CPU_ARCH')
if platform_id:
plat += '-' + platform_id
build_path = join(__path__[0], '..', 'build')
build_path = join(__path__[0], '..', 'build') # noqa
arch = '{}-{}.{}'.format(plat, *sys.version_info[0:2])
# If we are running the code from the source directory, then we will
......@@ -22,7 +22,7 @@ sys.path.insert(0, join(build_path, 'lib.' + arch))
if 'OMP_NUM_THREADS' not in os.environ:
os.environ['OMP_NUM_THREADS'] = '1'
from gpaw.broadcast_imports import broadcast_imports
from gpaw.broadcast_imports import broadcast_imports # noqa
with broadcast_imports:
import os
......
from __future__ import print_function, division
from functools import partial
import numpy as np
from scipy.spatial import Delaunay
from ase.utils import devnull
from ase.utils import convert_string_to_fd
from ase.utils import devnull, convert_string_to_fd
from ase.utils.timing import timer, Timer
from _gpaw import tetrahedron_weight
from scipy.spatial import Delaunay
from scipy.linalg.blas import zher
import gpaw.mpi as mpi
from gpaw.utilities.blas import gemm, rk, czher, mmm
from gpaw.utilities.blas import gemm, rk, mmm
from gpaw.utilities.progressbar import ProgressBar
from functools import partial
from _gpaw import tetrahedron_weight
def czher(alpha: float, x, A) -> None:
"""Hermetian rank-1 update of upper half of A.
class Integrator():
A += alpha * np.outer(x.conj(), x)
"""
AT = A.T
out = zher(alpha, x, 1, 1, 0, len(x), AT, 1)
assert out is AT
class Integrator:
def __init__(self, cell_cv, response='density', comm=mpi.world,
txt='-', timer=None, nblocks=1, eshift=0.0):
"""Baseclass for Brillouin zone integration and band summation.
......@@ -27,7 +34,7 @@ class Integrator():
comm: mpi.communicator
nblocks: block parallelization
"""
self.response = response
self.comm = comm
self.eshift = eshift
......@@ -148,7 +155,7 @@ class PointIntegrator(Integrator):
mydomain_t = self.distribute_domain(domain)
nbz = len(domain[0])
get_matrix_element, get_eigenvalues = integrand
prefactor = (2 * np.pi)**3 / self.vol / nbz
out_wxx /= prefactor
......@@ -201,13 +208,13 @@ class PointIntegrator(Integrator):
else:
for out_xx in out_wxx:
out_xx[iu] = out_xx[il].conj()
out_wxx *= prefactor
@timer('CHI_0 update')
def update(self, n_mG, deps_m, wd, chi0_wGG, timeordered=False, eta=None):
"""Update chi."""
omega_w = wd.get_data()
deps_m += self.eshift * np.sign(deps_m)
if timeordered:
......@@ -216,7 +223,7 @@ class PointIntegrator(Integrator):
else:
deps1_m = deps_m + 1j * eta
deps2_m = deps_m - 1j * eta
for omega, chi0_GG in zip(omega_w, chi0_wGG):
if self.response == 'density':
x_m = (1 / (omega + deps1_m) - 1 / (omega - deps2_m))
......@@ -226,7 +233,7 @@ class PointIntegrator(Integrator):
nx_mG = n_mG[:, self.Ga:self.Gb] * x_m[:, np.newaxis]
else:
nx_mG = n_mG * x_m[:, np.newaxis]
gemm(1.0, n_mG.conj(), np.ascontiguousarray(nx_mG.T),
1.0, chi0_GG)
......
from __future__ import print_function
from gpaw.utilities.blas import czher, axpy
from gpaw.response.integrators import czher
import numpy as np
from time import time
alpha = 0.5
x = np.random.rand(3) + 1j * np.random.rand(3)
a = np.random.rand(9).reshape(3,3) + np.random.rand(9).reshape(3,3) * 1j
a = np.random.rand(9).reshape(3, 3) + np.random.rand(9).reshape(3, 3) * 1j
# make a hermitian
for i in range(3):
for j in range(3):
a[i,j] = a[j,i].conj()
a[i,i] = np.real(a[i,i])
a[i, j] = a[j, i].conj()
a[i, i] = np.real(a[i, i])
b = alpha * np.outer(x.conj(), x) + a
czher(alpha, x, a)
for i in range(3):
for j in range(i,3):
a[j,i] = a[i,j].conj()
assert np.abs(b-a).sum() < 1e-14
for j in range(i, 3):
a[j, i] = a[i, j].conj()
assert np.abs(b - a).sum() < 1e-14
# testing speed
t_czher = 0
......@@ -33,7 +32,7 @@ for i in np.arange(1000):
t0 = time()
xx = np.outer(x.conj(), x)
axpy(alpha, xx, a)
a += alpha * xx
t_axpy += time() - t0
print('t_czher:', t_czher)
......
......@@ -121,27 +121,6 @@ def axpy(alpha, x, y):
assert z is y, (x, y, x.shape, y.shape)
def czher(alpha, x, a):
"""alpha x * x.conj() + a.
Performs the operation::
y <- alpha * x * x.conj() + a
where x is a N element vector and a is a N by N hermitian matrix, alpha
is a real scalar.
"""
assert isinstance(alpha, float)
assert is_contiguous(x, complex) and is_contiguous(a, complex)
assert x.flags.contiguous and a.flags.contiguous
assert x.ndim == 1 and a.ndim == 2
assert x.shape[0] == a.shape[0]
# Use zherk instead?
_gpaw.czher(alpha, x, a)
def rk(alpha, a, beta, c, trans='c'):
"""Rank-k update of a matrix.
......@@ -271,7 +250,7 @@ def _gemmdot(a, b, alpha=1.0, beta=1.0, out=None, trans='n'):
return out.reshape(outshape)
if 1:#not hasattr(_gpaw, 'mmm'):
if 1: # not hasattr(_gpaw, 'mmm'):
def gemm(alpha, a, b, beta, c, transa='n'): # noqa
if c.size == 0:
return
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment