code: def relaxation_vectors(A, R, k, alpha):
      
      # random n x R block in column ordering
      n = A.shape[0]
      x = np.random.rand(n * R) - 0.5
      x = np.reshape(x, (n, R), order='F')
      # for i in range(R):
      #     x[:,i] = x[:,i] - np.mean(x[:,i])
      b = np.zeros((n, 1))
  
      for r in range(0, R):
          jacobi(A, x[:, r], b, iterations=k, omega=alpha)
          # x[:,r] = x[:,r]/norm(x[:,r])
  
      return x
docs: Generate test vectors by relaxing on Ax=0 for some random vectors x.
  
      Parameters
      ----------
      A : csr_matrix
          Sparse NxN matrix
      alpha : scalar
          Weight for Jacobi
      R : integer
          Number of random vectors
      k : integer
          Number of relaxation passes
  
      Returns
      -------
      x : array
          Dense array N x k array of relaxation vectors