jan/meanas
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

store ZtAZ instead of AZU

Browse Source 
U is small (~number of modes)^2
This commit is contained in:
Jan Petykiewicz 2022-11-22 14:44:43 -08:00
parent bec0137c99
commit 03c15c8486
1 changed files with 8 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
meanas/fdfd/bloch.py
View File

@ -550,15 +550,16 @@ def eigsolve(
break break
for i in range(max_iters): for i in range(max_iters):
ZtZ = Z.conj().T @ Z Zt = Z.conj().T
ZtZ = Zt @ Z
U = numpy.linalg.inv(ZtZ) U = numpy.linalg.inv(ZtZ)
AZ = scipy_op @ Z AZ = scipy_op @ Z
AZU = AZ @ U ZtAZ = Zt @ AZ
ZtAZU = Z.conj().T @ AZU ZtAZU = ZtAZ @ U
E_signed = real(trace(ZtAZU)) E_signed = real(trace(ZtAZU))
sgn = numpy.sign(E_signed) sgn = numpy.sign(E_signed)
E = numpy.abs(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): 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)) 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_scale = num_modes / norm(D)
D *= d_scale D *= d_scale
ZtAZ = Z.conj().T @ AZ
AD = scipy_op @ D AD = scipy_op @ D
DtD = D.conj().T @ D DtD = D.conj().T @ D
DtAD = D.conj().T @ AD DtAD = D.conj().T @ AD
symZtD = _symmetrize(Z.conj().T @ D) symZtD = _symmetrize(Zt @ D)
symZtAD = _symmetrize(Z.conj().T @ AD) symZtAD = _symmetrize(Zt @ AD)
Qi_memo: List[Optional[float]] = [None, None] Qi_memo: List[Optional[float]] = [None, None]
@ -673,7 +672,7 @@ def eigsolve(
''' '''
Recover eigenvectors from Z Recover eigenvectors from Z
''' '''
U = numpy.linalg.inv(Z.conj().T @ Z) U = numpy.linalg.inv(ZtZ)
Y = Z @ scipy.linalg.sqrtm(U) Y = Z @ scipy.linalg.sqrtm(U)
W = Y.conj().T @ (scipy_op @ Y) W = Y.conj().T @ (scipy_op @ Y)