docs: def krylovMethod(self,tol=1e-8): 
          """
          We obtain ``pi`` by using the :func:``gmres`` solver for the system of linear equations. 
          It searches in Krylov subspace for a vector with minimal residual. The result is stored in the class attribute ``pi``.   
  
          Example
          -------
          >>> P = np.array([[0.5,0.5],[0.6,0.4]])
          >>> mc = markovChain(P)
          >>> mc.krylovMethod()
          >>> print(mc.pi) 
          [ 0.54545455  0.45454545]
          
          Parameters
          ----------
          tol : float, optional(default=1e-8)
              Tolerance level for the precision of the end result. A lower tolerance leads to more accurate estimate of ``pi``.        
          
          Remarks
          -------
          For large state spaces, this method may not always give a solution. 
          Code due to http://stackoverflow.com/questions/21308848/
          """            
          P       = self.getIrreducibleTransitionMatrix()
          
          #if P consists of one element, then set self.pi = 1.0
          if P.shape == (1, 1):
              self.pi = np.array([1.0]) 
              return
              
          size    = P.shape[0]
          dP      = P - eye(size)
          #Replace the first equation by the normalizing condition.
          A       = vstack([np.ones(size), dP.T[1:,:]]).tocsr()
          rhs     = np.zeros((size,))
          rhs[0]  = 1
                  
          pi, info = gmres(A, rhs, tol=tol)
          if info != 0:
              raise RuntimeError("gmres did not converge")
          self.pi = pi