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