diff --git a/krypy/linsys.py b/krypy/linsys.py
index 361e1ec..926aa82 100644
--- a/krypy/linsys.py
+++ b/krypy/linsys.py
@@ -285,6 +285,7 @@ def __init__(
         maxiter=None,
         explicit_residual=False,
         store_arnoldi=False,
+        store_all_xk = False,
         dtype=None,
     ):
         r"""Init standard attributes and perform checks.
@@ -351,6 +352,7 @@ def __init__(
         self.flat_vecs, (self.x0,) = utils.shape_vecs(x0)
         self.explicit_residual = explicit_residual
         self.store_arnoldi = store_arnoldi
+        self.store_all_xk = store_all_xk
 
         # get initial guess
         self.x0 = self._get_initial_guess(self.x0)
@@ -364,6 +366,7 @@ def __init__(
         self.tol = tol
 
         self.xk = None
+        self.xk_all = [self.x0]
         """Approximate solution."""
 
         # find common dtype
@@ -653,6 +656,8 @@ def _solve(self):
 
             # update solution
             yk += alpha * p
+            if self.store_all_xk:
+                self.xk_all.append(self._get_xk(yk))
 
             # update residual
             self.Mlrk -= alpha * Ap
@@ -846,6 +851,9 @@ def _solve(self):
             yk = yk + y[0] * z
             y = [y[1], 0]
 
+            if self.store_all_xk:
+                self.xk_all.append(self._get_xk(yk))
+
             self._finalize_iteration(yk, numpy.abs(y[0]))
 
         # compute solution if not yet done
@@ -990,6 +998,9 @@ def _solve(self):
             self.R[k : k + 2, k] = G[k].apply(self.R[k : k + 2, k])
             y[k : k + 2] = G[k].apply(y[k : k + 2])
 
+            if self.store_all_xk:
+                self.xk_all.append(self._get_xk(y[: k + 1]))
+
             self._finalize_iteration(y[: k + 1], abs(y[k + 1, 0]))
 
         # compute solution if not yet done
diff --git a/krypy/utils.py b/krypy/utils.py
index 7ff5077..f51154f 100644
--- a/krypy/utils.py
+++ b/krypy/utils.py
@@ -16,7 +16,7 @@
 from scipy.sparse import isspmatrix
 
 # from scipy.sparse.linalg import LinearOperator, aslinearoperator
-from scipy.sparse.sputils import isintlike
+# from scipy.sparse.sputils import isintlike
 
 __all__ = [
     "ArgumentError",
@@ -1369,7 +1369,7 @@ class LinearOperator(object):
     """
 
     def __init__(self, shape, dtype, dot=None, dot_adj=None):
-        if len(shape) != 2 or not isintlike(shape[0]) or not isintlike(shape[1]):
+        if len(shape) != 2 or int(shape[0]) != shape[0] or int(shape[1]) != shape[1]:
             raise LinearOperatorError("shape must be (m,n) with m and n " "integer")
         self.shape = shape
         self.dtype = numpy.dtype(dtype)  # defaults to float64
@@ -1525,7 +1525,7 @@ def __init__(self, A, p):
             raise LinearOperatorError("LinearOperator expected as A")
         if A.shape[0] != A.shape[1]:
             raise LinearOperatorError("square LinearOperator expected as A")
-        if not isintlike(p):
+        if int(p) != p:
             raise LinearOperatorError("integer expected as p")
         self.args = (A, p)
         super(_PowerLinearOperator, self).__init__(