@Carl Love 

1) Thanks for comment on zero entries.

2)  Transpose of M does not help.  

In Maple 12:

 convert(M^%T, list)[];

@Carl Love  Done in Maple 12 Classic.

# Return the angle of a point in interval [0, 2*Pi).

Angle:= proc(P::[realcons,realcons])

local arg:= evalf(arctan(P[2],P[1]));

     `if`(arg < 0, evalf(arg+2*Pi), arg)

end proc:

SortByAngle:= (A::Matrix)->

     Matrix(sort(convert(A^%T, listlist), (P,Q)-> Angle(P) < Angle(Q)))^%T:

M:= < -1,2,3,1,1,-1 ; 0,2,3,1,0,0 >:

SortByAngle(M);

@wingjammer   My example is of a general nature, you can write in the brackets the order of columns that you need.

Specify your problem more clearly. I do not understand what you want in the general case.

@Carl Love   Your code is compact and elegant, but conceptually more complex and less effective:

t := time():  SpM(1000):  time()-t;

t := time():  Am(1000):  time()-t;

Of course, this can be done in different ways. Here are 2 other ways:

restart;

f:=x+y+z:

V:=<seq(diff(f, var), var=[x,y,z])>;

restart;

f:=x+y+z:  V:=Vector():  var:=[x,y,z]:

for i to nops(var) do

V(i):=diff(f, var[i]):

od:

V;

@Carl Love   It is sufficient to find the exact roots of the polynomial (the first command in my answer)

Another example:

factor(3*x^2-7);

factor(3*x^2-7, {sqrt(3), sqrt(7)});

 -op(1,%)*3*``(expand(op(2,%)/sqrt(3)))*``(-expand(op(3,%)/sqrt(3)));

@Carl Love   In your example we can do exact symbolic factorization:

A := factor(x^4+1, sqrt(2));

B := ``(op(1, A)*op(2, A))*op(3, A);

@Markiyan Hirnyk  My program has the following differences from these:

1) The program works for systems of all sizes, and these tutors only if the every dimension of the matrix  no greater than 5.

2) My program in the during of Gauss - Jordan elimination , simulating manual solution produces some simplifying transformations. For example, solving the system  {19*x+20*y=39, x+y=2}, to get rid from fractions, useful to swap positions of equations.

@Axel Vogt   Yes, if we increase Digits, the error disappears:

Digits := 30:

KsT := 0.10e11: KaT := 4.6*0.10e6: KsE2 := 3.14*0.10e10: KaE2 := 4.21*0.10e6: KsDHT := 3*0.10e6: KaDHT := 3.5*0.10e6:

eq1 := E2 = fE2*(1+(KaE2+Ca)/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsE2*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq2 := T = fT*(1+KaT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq3 := DHT = fDHT*(1+KaDHT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsDHT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

E2, T, DHT, Ca, Cshbg := 1, 2, 3, 4, 5:

solve({eq1, eq2, eq3}, useassumptions) assuming fE2 > 0, fT > 0, fDHT > 0;

     {fDHT = 0.300792608884171123405264656623e-5, fE2 = 3.02872370161763976538808187091*10^(-10), fT    = 1.90199449982478958012781172005*10^(-10)}

@velimir74  I do not know the cause of this extra solution. In fact, this is a bug, because easy to verify that this is not a solution:

KsT := 0.10e11: KaT := 4.6*0.10e6: KsE2 := 3.14*0.10e10: KaE2 := 4.21*0.10e6: KsDHT := 3*0.10e6: KaDHT := 3.5*0.10e6:

eq1 := E2 = fE2*(1+(KaE2+Ca)/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsE2*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq2 := T = fT*(1+KaT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

eq3 := DHT = fDHT*(1+KaDHT*Ca/(1+KaE2*fE2+KaT*fT+KaDHT*fDHT)+KsDHT*Cshbg/(1+KsE2*fE2+KsT*fT+KsDHT*fDHT)):

E2, T, DHT, Ca, Cshbg:=1, 2, 3, 4, 5:

eval({eq1, eq2, eq3}, {fE2 = 7.15539742720000*10^10, fT = 3.30345789526199, fDHT =

4.72656319086592*10^14});

                               {1 = .7155397427e11, 2 = 3.303457895, 3 = .4726563191e15}

@wingjammer

UniqueColumns1:=proc(A::Matrix)

local L;

uses LinearAlgebra, ListTools;

L:=[Categorize((x,y)->Equal(x,y), [seq(A[..,i],i=1..ColumnDimension(A))])];

seq(L[i,1], i=1..nops(L));

end;

Your example:

Q := Matrix([[1,1,1,0,0,0,-2,-1], [0,0,0,1,1,1,0,-3]]);

UniqueColumns1(Q);

Addition: another way.

UniqueColumns2:=proc(A::Matrix)

uses LinearAlgebra, ListTools;

MakeUnique([seq(A[..,i], i=1..ColumnDimension(A))], 1, (x,y)->Equal(x,y));

end; 

@Carl Love   

sqrt(139.^2+105^2);

                      174.2010333

@as5987   If you are using the sequence then you should write

with(plots):  with(plottools):

d := seq(circle([1, 1/i], 1/i, color = blue), i=0.5..20, 0.1):

Example:

display(d[1], d[10], d[50]);

The number in brackets should be positive integer.  It shows the sequence number a member of the sequence.

@as5987  plottools[disk]  command creates plot structure of filled circle centered at  (1, 0)  and the radius  0.01 .  The same can be written shorter

Point1 := disk(eval([x,y], Sol[1]), 0.01, color = red):

@Carl Love   Thanks for the helpful comments.