function    [A] = transinplace(A);
%
%   computes the transpose of a CSC matrix A
%   in place 
%
%   Uses Tim Davis' bucket sort linear time
%   construction of CSC format from Triplet Format
%
%   Uses Joe Young's idea of compressing (j,i) -> j + m*(i-1) 
%
%   Warning:   Only good when A is square i.e. A.n = A.m
%              is required (should include a test and error flag) 
%
%   D.C. Sorensen
%   8 Sep 06
%
    n = A.n;
    m = A.m;
    nz = A.nz;
%
%   reconstruct column indices j so that
%   [A.i, j, A.x]  is in   [i,j,Aij] format
%   and then swap i <-> j to get transpose
%   and also compress (j,i) into one location 
%

    c = zeros (1,m) ;
    for j = 1:n,
	for (p = A.p(j) : A.p(j+1)-1),
		i = A.i(p);
		c(i) = c(i) + 1;

%               A.i(p) = j;
%      
		A.i(p) = j-1 + m*(i-1);
		% this records the (i,j) in one location
		% care taken so j = m doesn't fail
	end
    end

    A.p = cumsum ([1 c]) ;
    w = A.p ;
    for col = 1:n
	for k = (A.p (col)) : (A.p (col+1)-1)

            Ik = mod(A.i(k),m) ; 
	    Jk = (A.i(k) - Ik)/m + 1;  
            while (col  ~= Jk)
		% the kth triplet is out of place
		ik = mod(A.i (k),m) ;
	        jk = (A.i(k) - ik)/m + 1;  
	        ik = ik + 1;
		xk = A.x (k) ;
		% find out where it should go in the final result
		p = w (jk) ;
		w (jk) = w (jk) + 1 ;

		% swap

                Ik = mod(A.i(p),m); 
	        Jk = (A.i(p) - Ik)/m + 1;  

	        itemp = A.i(k)
		A.i (k) = A.i (p) ;
		A.x (k) = A.x (p) ;
%	A.i (p) = ik ;
		A.i (p) = itemp ;
		A.x (p) = xk ;

		% now the kth triplet might still be out of place, but the
		% pth triplet is in its final position.  The pth triplet will
		% not be considered again as a "p"th triplet, since it is
		% inside the sorted range for its column.  Thus, we must
		% eventually place the kth triplet in its final place.
		% The algorithm will progress to the final solution: QED.

	    end

	end
    end

    for k = 1:nz,
        ik = mod(A.i (k),m) ;
	jk = (A.i(k) - ik)/m + 1;  
	A.i(k) = ik + 1;
    end
    %
    %   If  nzcmax = (max no. nonzeros per column) << n  then
    %   sorting within columns essentially takes
    %   O(1) time  w.r.t    n both in space and time.  
    %
    %   So, it seems this is linear in n
    %
    for col = 1:n
	range = (A.p (col)) : (A.p (col+1)-1)
	[indices perm] = sort(A.i(range));
	A.i(range) = indices;
	xtemp = A.x(range);
	A.x(range) = xtemp(perm);
    end