##USAGE	[ ret,C,DIM ] = commutate ( A , B {,TOL=1e6*eps} )
##
##	check for ZERO commutator matrix
##		C = A*B - B*A
##	of equally sized square-matrices A,B
##
##ret	FALSE=0    non-commutating matrices A,B
##	TRUE =1    elseif all( all( abs(C)<TOL))
##DIM	product C  square matrix dimension
##NOTE	* ret=TRUE scalar C=0 is returned instead of C=zeros(DIM)
##	* ret=TRUE holds ALWAYS if at least one of A,B is  scalar
##	* ret=TRUE holds ALWAYS if B  is an arbitrary  FUNCM( A )
##EXA	* commutate( A, logm( inv(expm(A)) + eye(DIM)  ) )   TRUE
##	* commutate( A, logm(  1./expm(A)  + eye(DIM)  ) )  FALSE
##	* commutate( A, logm( inv(expm(A)) + SCALAR!=0 ) )  FALSE
##	* commutate( A, logm( inv(expm(A)) + ones(DIM) ) )  FALSE
##	* commutate( A,sqrtm( inv(expm(A)) + eye (DIM) ) )   TRUE
##	  all tested /w arbitrary (i.e random+complex) A elements
#WARN	DON'T confuse with 'commutation_matrix'
##AUTHOR  (c) Rolf Fabian 000120V <fabian@tu-cottbus.de>

#MORE CHECKED EXA
#	commutate( diag(1:3),       diag(7:9)  )	TRUE
#	commutate( diag(1:3),flipud(diag(7:9)) ) 	FALSE
#	commutate( diag(1:3),fliplr(diag(7:9)) ) 	FALSE
#	commutate( diag(1:3),rot90 (diag(7:9)) ) 	FALSE
#	commutate( diag(1:3),rot90 (diag(7:9),2) ) 	TRUE

function [ret,C,DIM]=commutate(A,B,TOL)

	if (nargin>3||nargin<2)
		fprintf(stderr,"commutate argument error .. quit.\n");
		[ret,C]=[]; return;
	end

	#-- parameter --#
	TOLSTD=1e6*eps;	# tolerance by default
	TRUE =1;
	FALSE=0;
	WARN =1;	# warn at size incompatibility
	#-- editable ---#

	ret=FALSE;	#init output
	DIM=[];

	if (nargin<3) TOL=TOLSTD; end

	if ( any((sa=size(A)).-(sb=size(B)))||any(diff(sa))||any(diff(sb)) )	
		if (all(sa==1)||all(sb==1))	# A|B scalar, ret is always TRUE
			ret=1; C=0; DIM=max(max(sa,sb)); return;
		else
			if(WARN) fprintf(stderr,\
	  "commutate WARN: size incompatibility .. undefined commutator.\n");
			end
			C=[];
		end
	else
		if (all(all(abs( C=A*B.-B*A )<TOL))) ret=TRUE; C=0; end
		DIM=sa(1);
	end
end



