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