code: def davidson(
      A: LinearOperator,
      neig: int,
      mode: str,
      M: Optional[LinearOperator] = None,
      max_niter: int = 1000,
      nguess: Optional[int] = None,
      v_init: str = "randn",
      max_addition: Optional[int] = None,
      min_eps: float = 1e-6,
      verbose: bool = False,
      **unused,
  ) -> Tuple[torch.Tensor, torch.Tensor]:
      """
      Using Davidson method for large sparse matrix eigendecomposition [2]_.
  
      Examples
      --------
      >>> import torch
      >>> import numpy as np
      >>> from deepchem.utils.differentiation_utils import LinearOperator
      >>> A = LinearOperator.m(torch.rand(2, 2))
      >>> neig = 2
      >>> mode = "lowest"
      >>> eigen_val, eigen_vec = davidson(A, neig, mode)
  
      Parameters
      ----------
      A: LinearOperator
          Linear operator to be diagonalized. Shape: ``(*BA, q, q)``.
      neig: int
          Number of eigenvalues and eigenvectors to be calculated.
      mode: str
          Mode of the eigenvalues to be calculated (``"lowest"``, ``"uppest"``)
      M: Optional[LinearOperator] (default None)
          The overlap matrix. If None, identity matrix is used. Shape: ``(*BM, q, q)``.
      max_niter: int
          Maximum number of iterations
      v_init: str
          Mode of the initial guess (``"randn"``, ``"rand"``, ``"eye"``)
      max_addition: int or None
          Maximum number of new guesses to be added to the collected vectors.
          If None, set to ``neig``.
      min_eps: float
          Minimum residual error to be stopped
      verbose: bool
          Option to be verbose
  
      Returns
      -------
      evals: torch.Tensor
          Eigenvalues of the linear operator.
      evecs: torch.Tensor
          Eigenvectors of the linear operator.
  
      References
      ----------
      .. [2] P. Arbenz, "Lecture Notes on Solving Large Scale Eigenvalue Problems"
             http://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter12.pdf
  
      """
      # TODO: optimize for large linear operator and strict min_eps
      # Ideas:
      # (1) use better strategy to get the estimate on eigenvalues
      # (2) use restart strategy
  
      if nguess is None:
          nguess = neig
      if max_addition is None:
          max_addition = neig
  
      # get the shape of the transformation
      na = A.shape[-1]
      if M is None:
          bcast_dims = A.shape[:-2]
      else:
          bcast_dims = get_bcasted_dims(A.shape[:-2], M.shape[:-2])
      dtype = A.dtype
      device = A.device
  
      prev_eigvalT = None
  
      # set up the initial guess
      V = _set_initial_v(
          v_init.lower(), dtype, device, bcast_dims, na, nguess, M=M
      )  # (*BAM, na, nguess)
  
      best_resid: Union[float, torch.Tensor] = float("inf")
      AV = A.mm(V)
      for i in range(max_niter):
          VT = V.transpose(-2, -1)  # (*BAM,nguess,na)
          # Can be optimized by saving AV from the previous iteration and only
          # operate AV for the new V. This works because the old V has already
          # been orthogonalized, so it will stay the same
          # AV = A.mm(V) # (*BAM,na,nguess)
          T = torch.matmul(VT, AV)  # (*BAM,nguess,nguess)
  
          # eigvals are sorted from the lowest
          # eval: (*BAM, nguess), evec: (*BAM, nguess, nguess)
          eigvalT, eigvecT = torch.linalg.eigh(T)
          eigvalT, eigvecT = _take_eigpairs(
              eigvalT, eigvecT, neig, mode
          )  # (*BAM, neig) and (*BAM, nguess, neig)
  
          # calculate the eigenvectors of A
          eigvecA = torch.matmul(V, eigvecT)  # (*BAM, na, neig)
  
          # calculate the residual
          AVs = torch.matmul(AV, eigvecT)  # (*BAM, na, neig)
          LVs = eigvalT.unsqueeze(-2) * eigvecA  # (*BAM, na, neig)
          if M is not None:
              LVs = M.mm(LVs)
          resid = AVs - LVs  # (*BAM, na, neig)
  
          # print information and check convergence
          max_resid = resid.abs().max()
          if prev_eigvalT is not None:
              deigval = eigvalT - prev_eigvalT
              max_deigval = deigval.abs().max()
              if verbose:
                  print(
                      "Iter %3d (guess size: %d): resid: %.3e, devals: %.3e"
                      % (i + 1, nguess, max_resid, max_deigval)
                  )  # type:ignore
  
          if max_resid < best_resid:
              best_resid = max_resid
              best_eigvals = eigvalT
              best_eigvecs = eigvecA
          if max_resid < min_eps:
              break
          if AV.shape[-1] == AV.shape[-2]:
              break
          prev_eigvalT = eigvalT
  
          # apply the preconditioner
          t = -resid  # (*BAM, na, neig)
  
          # orthogonalize t with the rest of the V
          t = to_fortran_order(t)
          Vnew = torch.cat((V, t), dim=-1)
          if Vnew.shape[-1] > Vnew.shape[-2]:
              Vnew = Vnew[..., : Vnew.shape[-2]]
          nadd = Vnew.shape[-1] - V.shape[-1]
          nguess = nguess + nadd
          if M is not None:
              MV_ = M.mm(Vnew)
              V, R = tallqr(Vnew, MV=MV_)
          else:
              V, R = tallqr(Vnew)
          AVnew = A.mm(V[..., -nadd:])  # (*BAM,na,nadd)
          AVnew = to_fortran_order(AVnew)
          AV = torch.cat((AV, AVnew), dim=-1)
  
      eigvals = best_eigvals  # (*BAM, neig)
      eigvecs = best_eigvecs  # (*BAM, na, neig)
      return eigvals, eigvecs
  
file_path: deepchem/deepchem/utils/differentiation_utils/symeig.py
imports: import torch
  from typing import Optional, Sequence, Tuple, Union
  from deepchem.utils.differentiation_utils import LinearOperator
  import functools
  from deepchem.utils.pytorch_utils import tallqr, to_fortran_order
  from deepchem.utils.differentiation_utils import get_bcasted_dims
label_desc: Use image node_id 2 for calling a global function with example usage: davidson(A, neig, mode, M, max_niter, nguess, v_init, max_addition, min_eps, verbose) and returns: eigvals,  eigvecs