From 03c15c84865939a4a11daa1f8ecbbca804fa4cc9 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Tue, 22 Nov 2022 14:44:43 -0800 Subject: [PATCH] store ZtAZ instead of AZU U is small (~number of modes)^2 --- meanas/fdfd/bloch.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/meanas/fdfd/bloch.py b/meanas/fdfd/bloch.py index c4c4b60..4f6d90c 100644 --- a/meanas/fdfd/bloch.py +++ b/meanas/fdfd/bloch.py @@ -550,15 +550,16 @@ def eigsolve( break for i in range(max_iters): - ZtZ = Z.conj().T @ Z + Zt = Z.conj().T + ZtZ = Zt @ Z U = numpy.linalg.inv(ZtZ) AZ = scipy_op @ Z - AZU = AZ @ U - ZtAZU = Z.conj().T @ AZU + ZtAZ = Zt @ AZ + ZtAZU = ZtAZ @ U E_signed = real(trace(ZtAZU)) sgn = numpy.sign(E_signed) E = numpy.abs(E_signed) - G = (AZU - Z @ U @ ZtAZU) * sgn + G = (AZ @ U - Z @ U @ ZtAZU) * sgn if i > 0 and abs(E - prev_E) < tolerance * 0.5 * (E + prev_E + 1e-7): logger.info('Optimization succeded: {} - 5e-8 < {} * {} / 2'.format(abs(E - prev_E), tolerance, E + prev_E)) @@ -577,14 +578,12 @@ def eigsolve( d_scale = num_modes / norm(D) D *= d_scale - ZtAZ = Z.conj().T @ AZ - AD = scipy_op @ D DtD = D.conj().T @ D DtAD = D.conj().T @ AD - symZtD = _symmetrize(Z.conj().T @ D) - symZtAD = _symmetrize(Z.conj().T @ AD) + symZtD = _symmetrize(Zt @ D) + symZtAD = _symmetrize(Zt @ AD) Qi_memo: List[Optional[float]] = [None, None] @@ -673,7 +672,7 @@ def eigsolve( ''' Recover eigenvectors from Z ''' - U = numpy.linalg.inv(Z.conj().T @ Z) + U = numpy.linalg.inv(ZtZ) Y = Z @ scipy.linalg.sqrtm(U) W = Y.conj().T @ (scipy_op @ Y)