docs: def advance(self):
          """Carry out one iteration of Arnoldi."""
          if self.iter >= self.maxiter:
              raise ArgumentError('Maximum number of iterations reached.')
          if self.invariant:
              raise ArgumentError('Krylov subspace was found to be invariant '
                                  'in the previous iteration.')
  
          N = self.V.shape[0]
          k = self.iter
  
          # the matrix-vector multiplication
          Av = self.A * self.V[:, [k]]
  
          if self.ortho == 'house':
              # Householder
              for j in range(k+1):
                  Av[j:] = self.houses[j].apply(Av[j:])
                  Av[j] *= numpy.conj(self.houses[j].alpha)
              if k+1 < N:
                  house = House(Av[k+1:])
                  self.houses.append(house)
                  Av[k+1:] = house.apply(Av[k+1:]) * numpy.conj(house.alpha)
                  self.H[:k+2, [k]] = Av[:k+2]
              else:
                  self.H[:k+1, [k]] = Av[:k+1]
              # next line is safe due to the multiplications with alpha
              self.H[k+1, k] = numpy.abs(self.H[k+1, k])
              if self.H[k+1, k] / numpy.linalg.norm(self.H[:k+2, :k+1], 2)\
                      <= 1e-14:
                  self.invariant = True
              else:
                  vnew = numpy.zeros((N, 1), dtype=self.dtype)
                  vnew[k+1] = 1
                  for j in range(k+1, -1, -1):
                      vnew[j:] = self.houses[j].apply(vnew[j:])
                  self.V[:, [k+1]] = vnew * self.houses[-1].alpha
          else:
              # determine vectors for orthogonalization
              start = 0
  
              # Lanczos?
              if self.ortho == 'lanczos':
                  start = k
                  if k > 0:
                      self.H[k-1, k] = self.H[k, k-1]
                      if self.M is not None \
                              and not isinstance(self.M, IdentityLinearOperator):
                          Av -= self.H[k, k-1] * self.P[:, [k-1]]
                      else:
                          Av -= self.H[k, k-1] * self.V[:, [k-1]]
  
              # (double) modified Gram-Schmidt
              for reortho in range(self.reorthos+1):
                  # orthogonalize
                  for j in range(start, k+1):
                      alpha = inner(self.V[:, [j]], Av, ip_B=self.ip_B)[0, 0]
                      if self.ortho == 'lanczos':
                          # check if alpha is real
                          if abs(alpha.imag) > 1e-10:
                              warnings.warn(
                                  'Iter {0}: abs(alpha.imag) = {1} > 1e-10. '
                                  'Is your operator self-adjoint in the '
                                  'provided inner product?'
                                  .format(self.iter, abs(alpha.imag)))
                          alpha = alpha.real
                      self.H[j, k] += alpha
                      if self.M is not None:
                          Av -= alpha * self.P[:, [j]]
                      else:
                          Av -= alpha * self.V[:, [j]]
              if self.M is not None:
                  MAv = self.M * Av
                  self.H[k+1, k] = norm(Av, MAv, ip_B=self.ip_B)
              else:
                  self.H[k+1, k] = norm(Av, ip_B=self.ip_B)
              if self.H[k+1, k] / numpy.linalg.norm(self.H[:k+2, :k+1], 2)\
                      <= 1e-14:
                  self.invariant = True
              else:
                  if self.M is not None:
                      self.P[:, [k+1]] = Av / self.H[k+1, k]
                      self.V[:, [k+1]] = MAv / self.H[k+1, k]
                  else:
                      self.V[:, [k+1]] = Av / self.H[k+1, k]
  
          # increase iteration counter
          self.iter += 1