move solve_general_linear to matrices and name gaussj_solve

This method is now used by default in solve instead of inversion
methods. It appears to be faster and leads to simpler results.
here are some examples:

>>> from sympy import *
>>> M=randMatrix(10)
>>> M.det()
-800548998337359151
>>> rhs=randMatrix(10,1)
>>> from time import time
>>> t=time();sol=M.solve(rhs);time()-t
0.24900007247924805
>>> t=time();sol=M.solve(rhs, method="GE");time()-t
0.33899998664855957

ok, your was a little faster

Now try introducing a symbol:

>>> M[0,0] = Symbol('x')
>>> t=time();sol=M.solve(rhs);time()-t
1.2119998931884766
>>> t=time();sol=M.solve(rhs, method="GE");time()-t

<stopped after more than a minute>

Let's try a smaller system:

>>> M=randMatrix(3)
>>> M.det()
-14532
>>> rhs=randMatrix(3,1)
>>> M[0,0]=Symbol('x')
>>> t=time();sol=M.solve(rhs);time()-t
0.031999826431274414
>>> t=time();gsol=M.solve(rhs, method="GE");time()-t
0.35899996757507324
>>> [count_ops(i) for i in (sol, gsol)]
[99, 169]
>>> sol.simplify() == gsol.simplify()
True
>>>

So the simplicity of the final result and the time needed to
get there is much better.