multisolve3x3(AA,b)

function xx = multisolve3x3(AA,b)
%     xx = multisolve3x3(AA,b)
%     
%     Vectorized solve of multiple 3x3 linear systems. 
%     
%     The m 3-by-3-coefficient matrices can be given in 2 forms:
%     (1) The argument matrix AA is 3m x 3 (i.e., the k-th linear system occupies rows 3k-2:3k).
%     (2) The argument AA is an 3x3xm array. (the k-th linear system is AA(:,:,k) ).
%     
%     Similarly, the result vectors b can be given in 2 forms (independent of AA):
%     (1) a 3m x 1 vector, where the k-th result is in rows 3k-2:3k,
%     (2) a 3xm matrix, whose k-th column contains the k-th result.
%     
%     Either way, the solution xx is given in size identical to that of b.

    if isvector(b)
	  bb=b;
	  m3 = numel(bb);
	  m = m3/3;
    elseif length(size(b)) && size(b,1)==3
	  bb=b(:);
	  m=size(b,2);
	  m3=3*m;
    else
	  error('improper input arg size.');
    end
    
    
    if ndims(AA)==3 && all(size(AA)== [3 3 m])
	  AdjAA = permute(cross( AA([2 3 1],:,:), AA([3 1 2],:,:) ,2) ,[2 1 3]) ; % 3x3 adjoints in a 3x3xm array
	  Dets = dot(squeeze(AdjAA(:,1,:)), squeeze(AA(1,:,:)));

    elseif ndims(AA)==2 && all(size(AA)==[m3,3] )
	  ii=reshape(1:m3 , 3, m) ;   % TODO: transpose, use columns
	  AdjAA = cross( AA(ii([2 3 1],:),:), AA(ii([3 1 2],:),:) ,2) .';
	  Dets = dot(AdjAA(:,1:3:end), AA(1:3:end,:).');
	  
    else
	  error('improper input arg size.');
    end
    
    	  invDets = 1 ./ Dets ;
	  invDets = invDets([1 1 1],:);
	  invDets = invDets(:);

	  Iidx = [ 1:3:m3 ; 2:3:m3; 3:3:m3; 1:3:m3 ; 2:3:m3; 3:3:m3; 1:3:m3 ; 2:3:m3; 3:3:m3] ;
	  % = [ 1 2 3 1 2 3 1 2 3 4 5 6 4 5 6 4 5 6 ...]
	  Jidx =  [ 1:m3 ; 1:m3;  1:m3] ;
	  % = [ 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 ...

	  sparseAdjAA = sparse(Iidx,Jidx,AdjAA(:),m3,m3);
	  xx = sparseAdjAA * bb .* invDets;

	  if length(size(b)) && size(b,1)==3
		xx= reshape(xx,3,[]);
	  end