Repository: sympy/sympy

Issue:
Incorrect circular references in solutions from linsolve in underdetermined systems
**Edit**: See comment below pinpointing the problem and providing a possible solution

I am sorry that this is hardly a minimal working example, but this is a real world problem and though I've been playing around with it I can't reproduce it with simpler systems:

```python
V, f, J = symbols('V f J')
freevars_gen = numbered_symbols('x')
free_vars = [next(freevars_gen) for i in range(17)]
m = Matrix([
[     0,      0,      0,  2*I*J,     0, -2*I*J,  2*I*J,  2*I*J,  2*I*f,          0,           0,           0,   0,   0,      0,      0,          0,                   2*I*V**2*f**3/J**4],
[-2*I*J,      0,      0,      0, 2*I*J,      0, -2*I*J, -2*I*J,      0,          0,           0,           0,   0,   0,      0,      0,          0,               4*I*V**2*f**3/(3*J**4)],
[ 2*I*J,  2*I*J,      0,      0,     0, -2*I*J,  2*I*J,      0, -2*I*f,          0,           0,           0,   0,   0,  2*I*f,      0,          0,               8*I*V**2*f**3/(3*J**4)],
[-2*I*J, -2*I*J, -2*I*J,      0, 2*I*J,      0,      0,      0,      0, 2*I*V**2/f, -4*I*V**2/f, -2*I*V**2/f, 2*V, 2*V,      0, -2*I*f, 2*I*V**2/f, 2*I*V**2*f*(-2*V**2 + f**2)/(3*J**4)],
[     0,      0, -2*I*J, -2*I*J, 2*I*J,      0,      0, -2*I*J,      0,          0,           0,           0,   0,   0,      0,  2*I*f,          0,                  -4*I*V**2*f**3/J**4],
[     0, -2*I*J,      0, -2*I*J, 2*I*J,      0, -2*I*J,      0,      0,          0,           0,           0,   0,   0,      0,      0,          0,              -8*I*V**2*f**3/(3*J**4)],
[     0,  2*I*J,  2*I*J,  2*I*J,     0, -2*I*J,      0,      0,      0, 2*I*V**2/f, -4*I*V**2/f, -2*I*V**2/f, 2*V, 2*V, -2*I*f,      0, 2*I*V**2/f,   4*I*V**2*f*(-V**2 + f**2)/(3*J**4)],
[ 2*I*J,      0,  2*I*J,      0,     0, -2*I*J,      0,  2*I*J,      0,          0,           0,           0,   0,   0,      0,      0,          0,                                    0]])

solve_linear_system(m,*free_vars)
dict(zip(free_vars, [simplify(result) for result in list(linsolve(m, *free_vars))[0]]))

Out[43]: 
{x1: x7 + f*x14/(2*J) - f*x15/J - 5*V**2*f**3/(6*J**5),
 x8: x14 - x15 - 2*V**2*f**2/J**4,
 x3: x4 - x6 - x7 - f*x14/(2*J) + f*x15/J + 13*V**2*f**3/(6*J**5),
 x5: x4 + f*x14/(2*J) - 5*V**2*f**3/(6*J**5),
 x2: x6 + f*x14/(2*J) - V**2*f**3/(6*J**5),
 x9: 2*x10 + x11 - x16 + I*f*x12/V + I*f*x13/V + f**2*x14/V**2 - 2*V**2*f**2/(3*J**4) - 4*f**4/(3*J**4),
 x0: x4 - x6 - x7 - 2*V**2*f**3/(3*J**5)}

Out[44]: 
{x13: x13,
 x12: x12,
 x4: x6 - f*x14/(2*J) + 5*V**2*f**3/(6*J**5),
 x14: x14,
 x6: x7,
 x3: -x5 + x6 - x7 - f*x14/J + f*x15/J + 3*V**2*f**3/J**5,
 x10: x10,
 x11: x11,
 x1: x5 + f*x14/(2*J) - f*x15/J - 5*V**2*f**3/(6*J**5),
 x8: x14 - x15 - 2*V**2*f**2/J**4,
 x16: x16,
 x5: x6,
 x2: x7 + f*x14/(2*J) - V**2*f**3/(6*J**5),
 x9: 2*x10 + x11 - x16 + I*f*x12/V + I*f*x13/V + f**2*x14/V**2 - 2*V**2*f**2/(3*J**4) - 4*f**4/(3*J**4),
 x7: x5,
 x0: -x5 + x6 - x7 - f*x14/(2*J) + V**2*f**3/(6*J**5),
 x15: x15}
```

`linsolve` will only get this far with the changes in #10119 merged. Now observe in the `linsolve` answer we have a number of free variables denoted by terms in the dict like `x13: x13`. If we wanted to arrive at a result of the form in `solve_linear_system` we could just delete these terms from the dict, so that's fine.

But look at the x5, x6, x7 terms. We have circular references `x5: x6`, `x6: x7`, `x7: x5`, whereas one might expect something like `x5: x7`, `x6: x7`, `x7: x7` If you look at say the solution for x3 you find x5, x6, x7 all appear despite clearly not being independent. For simpler undetermined matrices this does not seem to happen, so I'm not sure what's the particular feature of this one causing the problem.

@aktech @hargup 


Discussion:

Patch:
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index 80a6a9c1a5..3a248551eb 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -4336,11 +4336,8 @@ def gauss_jordan_solve(self, b, freevar=False):
         if not v[rank:, 0].is_zero:
             raise ValueError("Linear system has no solution")
 
-        # Get index of free symbols (free parameter)
-        free_var_index = []
-        for r in range(col):
-            if r not in pivots:
-                free_var_index.append(r)  # non-pivots columns are free variables
+        # Get index of free symbols (free parameters)
+        free_var_index = permutation[len(pivots):]  # non-pivots columns are free variables
 
         # Free parameters
         dummygen = numbered_symbols("tau", Dummy)


Test patch:
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index db313f87ed..ec6f518db7 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -976,6 +976,11 @@ def test_linsolve():
     b = Matrix([0, 0, 1])
     assert linsolve((A, b), (x, y, z)) == EmptySet()
 
+    # Issue #10121 - Assignment of free variables
+    a, b, c, d, e = symbols('a, b, c, d, e')
+    Augmatrix = Matrix([[0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0]])
+    assert linsolve(Augmatrix, a, b, c, d, e) == FiniteSet((a, 0, c, 0, e))
+
 
 def test_issue_9556():
     x = Symbol('x')