From octave-sources-request at bevo dot che dot wisc dot edu  Mon Jan 24 03:26:47 2000
Subject: powm.m ( mat-by-mat power of commutating matrices)
From: Rolf Fabian <fabian at tu-cottbus dot de>
To: "'octave-sources UWISC'" <octave-sources at bevo dot che dot wisc dot edu>
Date: Mon, 24 Jan 2000 10:25:40 +0100

##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 at tu-cottbus dot 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



##USAGE   C = powm ( B, A  {,CCHK=1} )
##
##	power function C = B^A
##	arbitrary (numerical) input types
##
##	square matrix A exponent possible if
##	* base B scalar  -or-
##	* base B square matrix commutating with exponent A  >>NEW<<
##CCHK	flag for commutator matrix check
##	ONLY relevant at equally sized  SQUARE MAT B,A  input.
##	     If you KNOW a priori, that 'commutate(A,B)==TRUE',
##	     you may set CCHK=0 to speed-up calculation.
##	     Especially useful if called within iterations.
##	WARN result is definitely INCORRECT in general, if
##	     CCHK is set to 0 , but 'commutate(A,B)==FALSE'
##	>>>  DON'T USE THIS OPTION  if you're in doubt <<<
##NOTE	syntax identical to C/C++ pow scalar-function
##SEE TOO commutate              (author -written)
##	expm, logm, operator ^ (octave built-in)
##AUTHOR  (c) Rolf Fabian 000120V <fabian at tu-cottbus dot de>

#used-scripts
#commutate.m	author-written

#FIXME	There's still the problem, that 
#	for e.g.  B=ones(2)  i.e. non-diagonalizable
#	powm (B,A) hangs     e.g. for A=B
#	even octave:>  B^(non-integer scalar A)  FAILS

function C=powm(B,A,CCHK)

	if (nargin>3||nargin<2)
		usage ("C = powm ( B, A  {,CCHK=1} )"); C=[]; return;
	else
		if (nargin <3) CCHK=1; end 
		tag  ="powm: ";
	   	merr0="%s square matrix  %s required .. quit.\n";
		merr1=\
		"%s  B^A  undefined  for A not commutating with B  .. quit.\n";
	   	merr2="%s A,B square matrices,  but sizes mismatch .. quit.\n";

		[ar,ac]=size(A);
		[br,bc]=size(B);

		C=[]; 				#init output

		if (ar!=ac)
			fprintf(stderr,merr0,tag,'A'); return;
		end
		if (br!=bc)
			fprintf(stderr,merr0,tag,'B'); return;
		end

		if	( (BS=all([br,bc]==1)) )#scalar base
			if (B==0)
			    if (any(any(A)))
				 C=zeros(ar);	#overwrite octave NaN-output
			    else C=  eye(ar);	#shortcut
			    end
			else
			    C=B^A;
			end
		elseif	(!BS&& all([ar,ac]==1)) #MAT base + scalar exponent 
			if (any(any(B)))
				C=B^A;
			else			#here	scalar A
				if (A) C=zeros(br);
				else   C=eye(br);
				end
			end
		elseif	(ar==br)		#equally sized A,B
			if     ( !any(any(B)) )	#zero MAT base
			    if (  any(any(A)) )
				 C=zeros(ar);	#zeros(N)^MAT(N)  := zeros(N)
			    else C=  eye(ar);	#zeros(N)^zeros(N) =\
			    end			#  MAT(N)^zeros(N):=   eye(N)
			else
			    #call 'commutate' ONLY if REQUESTED by nonzero CCHK 
			    #REM  if called, 'commutate' checks for equal sizes too
			    if ( !CCHK || commutate(A,B) )
				 C = expm(A*logm(B));
			    else
			     	 fprintf(stderr,merr1,tag);
			    end
			end
		else
			fprintf(stderr,merr2,tag);
		end
	end
end