Merge pull request #102 from adc-connect/cvs-relax

Small fixes and modularisation of Lanczos

adcc/AdcMethod.py

@@ -40,17 +40,17 @@ def __init__(self, method):
                              + ",".join(self.available_methods) + " are known.")
 
         split = method.split("-")
-        self.base_method = split[-1]
+        self.__base_method = split[-1]
         split = split[:-1]
         self.is_core_valence_separated = "cvs" in split
 
         try:
-            if self.base_method == "adc2x":
+            if self.__base_method == "adc2x":
                 self.level = 2
             else:
-                self.level = int(self.base_method[-1])
+                self.level = int(self.__base_method[-1])
         except ValueError:
-            raise ValueError("Not a valid base method: " + self.base_method)
+            raise ValueError("Not a valid base method: " + self.__base_method)
 
     def at_level(self, newlevel):
         """
@@ -65,9 +65,9 @@ def at_level(self, newlevel):
     @property
     def name(self):
         if self.is_core_valence_separated:
-            return "cvs-" + self.base_method
+            return "cvs-" + self.__base_method
         else:
-            return self.base_method
+            return self.__base_method
 
     @property
     def property_method(self):
@@ -76,11 +76,19 @@ def property_method(self):
         for this ADC method. This only differs from the name property
         for the ADC(2)-x family of methods.
         """
-        if self.base_method == "adc2x":
+        if self.__base_method == "adc2x":
             return AdcMethod(self.name.replace("adc2x", "adc2")).name
         else:
             return self.name
 
+    @property
+    def base_method(self):
+        """
+        The base (full) method, i.e. with all approximations such as
+        CVS stripped off.
+        """
+        return AdcMethod(self.__base_method)
+
     def __eq__(self, other):
         return self.name == other.name
 

adcc/ElectronicTransition.py

@@ -92,6 +92,7 @@ def __init__(self, data, method=None, property_method=None):
         if isinstance(data, EigenSolverStateBase):
             self._excitation_vector = data.eigenvectors
             self._excitation_energy_uncorrected = data.eigenvalues
+            self.residual_norm = data.residual_norms
         else:
             if hasattr(data, "eigenvalues"):
                 self._excitation_energy_uncorrected = data.eigenvalues

adcc/solver/LanczosIterator.py

@@ -20,7 +20,9 @@
 ## along with adcc. If not, see <http://www.gnu.org/licenses/>.
 ##
 ## ---------------------------------------------------------------------
+import warnings
 import numpy as np
+import scipy.linalg as la
 
 from adcc import evaluate, lincomb
 from adcc.timings import Timer
@@ -248,3 +250,17 @@ def matrix_product(self):
                 vectors.append(r[ires])
             AV.append(lincomb(np.array(coefficients), vectors, evaluate=True))
         return AV
+
+    def check_orthogonality(self, tolerance=None):
+        if tolerance is None:
+            tolerance = self.n_problem * np.finfo(float).eps
+        orth = np.array([[SSi @ SSj for SSi in self.subspace]
+                         for SSj in self.subspace])
+        orth -= np.eye(len(self.subspace))
+        orth = np.max(np.abs(orth))
+        if orth > tolerance:
+            warnings.warn(la.LinAlgWarning(
+                "LanczosSubspace has lost orthogonality. "
+                "Expect inaccurate results."
+            ))
+        return orth

adcc/solver/lanczos.py

@@ -71,33 +71,70 @@ def default_print(state, identifier, file=sys.stdout):
         print("=== Restart ===", file=file)
 
 
-def amend_true_residuals(state, subspace, rvals, rvecs, epair_mask):
+def compute_true_residuals(subspace, rvals, rvecs, epair_mask):
     """
     Compute the true residuals and residual norms (and not the ones estimated from
-    the Lanczos subspace) and amend the `state` accordingly.
+    the Lanczos subspace).
     """
     V = subspace.subspace
     AV = subspace.matrix_product
 
     def form_residual(rval, rvec):
         coefficients = np.hstack((rvec, -rval * rvec))
         return lincomb(coefficients, AV + V, evaluate=True)
-    state.residuals = [form_residual(rvals[i], rvec)
-                       for i, rvec in enumerate(np.transpose(rvecs))
-                       if i in epair_mask]
-    state.eigenvectors = [lincomb(rvec, V, evaluate=True)
-                          for i, rvec in enumerate(np.transpose(rvecs))
-                          if i in epair_mask]
+    residuals = [form_residual(rvals[i], rvec)
+                 for i, rvec in enumerate(np.transpose(rvecs))
+                 if i in epair_mask]
+    eigenvectors = [lincomb(rvec, V, evaluate=True)
+                    for i, rvec in enumerate(np.transpose(rvecs))
+                    if i in epair_mask]
+    rnorms = np.array([np.sqrt(r @ r) for r in residuals])
+
+    # Note here the actual residual norm (and not the residual norm squared)
+    # is returned.
+    return eigenvectors, residuals, rnorms
 
-    rnorms = np.array([np.sqrt(r @ r) for r in state.residuals])
-    state.residual_norms = rnorms
+
+def amend_true_residuals(state, subspace, rvals, rvecs, epair_mask):
+    """
+    Compute the true residuals and residual norms (and not the ones estimated from
+    the Lanczos subspace) and amend the `state` accordingly.
+    """
+    res = compute_true_residuals(subspace, rvals, rvecs, epair_mask)
+    state.eigenvectors, state.residuals, state.residual_norms = res
 
     # TODO For consistency with the Davidson the residual norms are
-    #      squared again to give output in the same order of magnitude.
+    #      squared to give output in the same order of magnitude.
     state.residual_norms = state.residual_norms**2
+
     return state
 
 
+def check_convergence(subspace, rvals, rvecs, tol):
+    b = subspace.rayleigh_extension
+
+    # Norm of the residual vector block
+    norm_residual = np.sqrt(np.sum(subspace.residual[p] @ subspace.residual[p]
+                                   for p in range(subspace.n_block)))
+
+    # Minimal tolerance for convergence criterion
+    # same settings as in ARPACK are used:
+    #    norm(r) * norm(b^T * rvec) <= max(mintol, tol * abs(rval)
+    mintol = np.finfo(float).eps * np.max(np.abs(rvals))
+    eigenpair_error = []
+    eigenpair_converged = np.ones_like(rvals, dtype=bool)
+    for i, rval in enumerate(rvals):
+        lhs = norm_residual * np.linalg.norm(b.T @ rvecs[:, i])
+        rhs = max(mintol, tol * abs(rval))
+        if lhs > rhs:
+            eigenpair_converged[i] = False
+        if mintol < tol * abs(rval):
+            eigenpair_error.append(lhs / abs(rval))
+        else:
+            eigenpair_error.append(lhs * tol / mintol)
+    return eigenpair_converged, np.array(eigenpair_error)
+
+
 def lanczos_iterations(iterator, n_ep, min_subspace, max_subspace, conv_tol=1e-9,
                        which="LA", max_iter=100, callback=None,
                        debug_checks=False, state=None):
@@ -148,44 +185,18 @@ def callback(state, identifier):
         n_applies_offset = state.n_applies
 
     for subspace in iterator:
-        T = subspace.subspace_matrix
         b = subspace.rayleigh_extension
-        eps = np.finfo(float).eps
         with state.timer.record("rayleigh_ritz"):
-            rvals, rvecs = np.linalg.eigh(T)
+            rvals, rvecs = np.linalg.eigh(subspace.subspace_matrix)
 
         if debug_checks:
-            orth = np.array([[SSi @ SSj for SSi in subspace.subspace]
-                             for SSj in subspace.subspace])
-            orth -= np.eye(len(subspace.subspace))
-            state.subspace_orthogonality = np.max(np.abs(orth))
+            eps = np.finfo(float).eps
             orthotol = max(tol / 1000, subspace.n_problem * eps)
-            if state.subspace_orthogonality > orthotol:
-                warnings.warn(la.LinAlgWarning(
-                    "Subspace in lanczos has lost orthogonality. "
-                    "Expect inaccurate results."
-                ))
-
-        # Norm of the residual vector block
-        norm_residual = np.sqrt(np.sum(subspace.residual[p] @ subspace.residual[p]
-                                       for p in range(subspace.n_block)))
-
-        # Minimal tolerance for convergence criterion
-        # same settings as in ARPACK are used:
-        #    norm(r) * norm(b^T * rvec) <= max(mintol, tol * abs(rval)
-        mintol = eps * np.max(np.abs(rvals))
-        eigenpair_error = []
-        is_rval_converged = np.ones_like(rvals, dtype=bool)
-        for i, rval in enumerate(rvals):
-            lhs = norm_residual * np.linalg.norm(b.T @ rvecs[:, i])
-            rhs = max(mintol, tol * abs(rval))
-            if lhs > rhs:
-                is_rval_converged[i] = False
-            if mintol < tol * abs(rval):
-                eigenpair_error.append(lhs / abs(rval))
-            else:
-                eigenpair_error.append(lhs * tol / mintol)
-        eigenpair_error = np.array(eigenpair_error)
+            orth = subspace.check_orthogonality(orthotol)
+            state.subspace_orthogonality = orth
+
+        is_rval_converged, eigenpair_error = check_convergence(subspace, rvals,
+                                                               rvecs, tol)
 
         # Update state
         state.n_iter += 1