mldivide.m

几个关于多小波的程序

function [Q,method] = mldivide(P1,P2)

% MLDIVIDE -- left division of matrix polynomials
%
%        Q = P1 \ P2
%        Q = mldivide(P1,P2)
%
% This function is not meant to be called by the user. It is called by
% Matlab if an expression of the form P1\P2 is encountered, where at
% least one of the two operands is a matrix polynomial. P1 must be
% square, either of a size compatible with P2 or a scalar.
%
% This routine knows how to handle 3 special cases at the moment,
% otherwise it gives up:
%
% 1. P1 has a nonsingular first or last coefficient, in which case
%    we can do long division; if it fails, we know that division is
%    not possible. The cases P1 = monomial and P1 = scalar fall into
%    this category
%
% 2. P1 is a projection factor, which has a known inverse
%
% 3. P1 is upper or lower triangular, with an invertible diagonal; this
%    includes the case P1 = diagonal matrix.
%
% If none of these cases apply, division may or may not be possible.
% As a last resort, you could try Q = longinv(P1) * P2.
%
% For symbolic matrices, this routine may take a long time and
% return a very complicated result.
%
% The METHOD output parameter is strictly for testing purposes. It 
% indicates which of the special cases was applied.

% Copyright (c) 2004 by Fritz Keinert ([e-mail]),
% Dept. of Mathematics, Iowa State University, Ames, IA 50011.
% This software may be freely used and distributed for non-commercial
% purposes, provided this copyright statement is preserved, and
% appropriate credit for its use is given.
%
% Last update: Feb 20, 2004

% make sure P1, P2 are both matrix polynomials

if ~isa(P1,'mpoly')
    P1 = mpoly(P1);
end
if ~isa(P2,'mpoly')
    P2 = mpoly(P2);
end
P1 = trim(P1);
P2 = trim(P2);
P2max = get(P2,'max');

% check that P1 is square and compatible with P2

[n1,n2] = size(P1);
[m1,m2] = size(P2);
if (n1 ~= n2 | (n1 ~= 1 & n1 ~= m1))
    error('incompatible sizes in division');
end

% Case 1: check if either the first or last coefficient of P1 is
% non-singular. If yes, attempt long division. 
%
% If the first coefficient is nonsingular, we divide left to right. 
% If the last coefficient is nonsingular, but not the first, we would
% have to do the division right to left. To save on programming
% effort, I instead reverse the order of the coefficients before and
% after the division.
%
% If both the first and last coefficients of P1 are nonsingular 
% (which includes the cases P1 = scalar and P1 = monomial), then 
% degree of Q = degree of P2 - degree of P1. 
%
% Otherwise, the degree of Q could be larger. The worst case I can
% think of that could be useful is when P1 is unit lower or upper
% triangular, in which case the degree of P1inv could be as large as
% (n-1)*(degree of P1), n = size of P1. If the division does not
% terminate at that point, we give up.

% Algorithm:
%
% 1. Pad P2 with P2extra zero matrices, where P2extra=0 if the first
% and last coefficient of P1 are both non-singular, and 
% P2extra=n*(length(P1)-1) otherwise. 
%
% The earliest possibility for the division to terminate
% successfully is when we reach the end of the original P2, which is
% after 
%
%        kmin = length(P2) - length(P1) + 1
%
% division steps. The last possibility we consider is when we reach
% the end of (P2 + P2extra zeros), which happens after
%
%        kmax = kmin + P2extra
%
% steps.
%
% We start with (kmin-1) division steps without checking the
% remainder, then we start checking.

method = 1;

first = ~issingular(P1.coef(:,:,1));
last  = ~issingular(P1.coef(:,:,end));

if (last & ~first)
    P1 = reverse(P1);
    P2 = reverse(P2);
end
    
if (~first & size(P1.coef,3) == 1)
% we are trying to divide by a single singular matrix
    error('left division failed');
end

if (first | last)
    P1length = size(P1.coef,3);
    P2length = size(P2.coef,3);
    if (first & last)
	if (P2length < P1length)
	    error('left division failed');
	end
	P2extra = 0;
    else
	P2extra = n1 * (P1length - 1);
    end
    Qlengthmin = P2length - P1length + 1;
    Qlength = Qlengthmin + P2extra;
    % pad P2 with P2extra zeros
    P2 = trim(P2,[P2.min,P2max+P2extra]);
    % do division
    for i = 1:Qlength
	Qcoef(:,:,i) = P1.coef(:,:,1) \ P2.coef(:,:,i);
	for j = 2:P1length
	    P2.coef(:,:,i+j-1) = P2.coef(:,:,i+j-1) - P1.coef(:,:,j)*Qcoef(:,:,i);
	end
	if (i >= Qlengthmin) % check remainder
	    remainder = mpoly(P2.coef(:,:,i+1:i+P1length-1));
	    if iszero(remainder) % success!
		Q = mpoly(Qcoef,P2.min-P1.min);
		Q = trim(Q);
		if (last & ~first)
		    Q = reverse(Q);
		end
		[type,m,r] = match_type(P1,P2);
		Q = set(Q,'type',type,'m',m,'r',r);
		return
	    end
	end
    end
% If we get to this point, long division failed.
% If the first and last coefficients do not contain symbolic variables,
% we know that division is not possible
    if (~issymbolic(P1.coef(:,:,1)) & ~issymbolic(P1.coef(:,:,end)))
	error('left division failed');
    end
end


% case 2: check if P1 is a generalized projection factor
% That means we try to factor P1 as
%
%                              s   k
%        P1 = P1  [ A + (I-A) z ] z
%               0
%       2
% with A  = A. If this succeeds, then
%
%          -1               -s    -1  -k
%        P1   = [A + (I-A) z  ] P1   z
%                                 0

method = 2;

P10 = moment(P1);
if ~issingular(P10)
    F = P1;
    for i = 1:size(F.coef,3)
	F.coef(:,:,i) = P10 \ F.coef(:,:,i);
    end
    A = F.coef(:,:,1);
    if iszero(A^2 - A)
	if iszero(F.coef(:,:,end) - eye(n1) + A)
	    k = F.min;
	    s = F.min + size(F.coef,3) - 1 - k;
	    middle = mpoly(F.coef(:,:,2:end-1));
	    if iszero(middle) % success!
		P1inv = F;
		P1inv.coef(:,:,1) = eye(n1) - A;
		P1inv.coef(:,:,end) = A;
		P1inv.min = - k - s;
		for i = 1:size(P1inv.coef,3)
		    P1inv.coef(:,:,i) = P1inv.coef(:,:,i) / P10;
		end
		Q = trim(P1inv * P2);
		[type,m,r] = match_type(P1,P2);
		Q = set(Q,'type',type,'m',m,'r',r);
		return
	    end
	end
    end
end

% case 3: check if the matrix is upper or lower triangular
% If it is, we attempt forward or back substitution.

method = 3;

if iszero(triu(P1,1)) % matrix is lower triangular
    Q = mpoly(zeros(m1,m2),0);
    if (isa(P1.coef,'sym') | isa(P2.coef,'sym'))
	Q = sym(Q);
    end
% Q(1,:) = P1(1,1) \ P2(1,:);
    Q = subsasgn(Q,substruct('()',{1,':'}),...
		 subsref(P1,substruct('()',{1,1})) \ ...
		 subsref(P2,substruct('()',{1,':'}))); 
    for j = 2:m1
% Q(j,:) = P1(j,j) \ (P2(j,:) - P1(j,1:j-1)*Q(1:j-1,:));
	Q = subsasgn(Q,substruct('()',{j,':'}),...
		     subsref(P1,substruct('()',{j,j})) \ ...
		     (subsref(P2,substruct('()',{j,':'})) - ...
		      subsref(P1,substruct('()',{j,1:j-1})) * ...
		      subsref(Q,substruct('()',{1:j-1,':'}))));
    end
    Q = trim(Q);
    [type,m,r] = match_type(P1,P2);
    Q = set(Q,'type',type,'m',m,'r',r);
    return
elseif iszero(tril(P1,-1)) % matrix is upper triangular
    Q = mpoly(zeros(m1,m2),0);
    % Q(m1,:) = P1(m1,m1) \ P2(m1,:);
    Q = subsasgn(Q,substruct('()',{m1,':'}),...
	subsref(P1,substruct('()',{m1,m1})) \ ...
	subsref(P2,substruct('()',{m1,':'}))); 
    for j = m1-1:-1:1
	% Q(j,:) = P1(j,j) \ (P2(j,:) - P1(j,j+1:m1)*Q(j+1:m1,:));
	Q = subsasgn(Q,substruct('()',{j,':'}),...
	    subsref(P1,substruct('()',{j,j})) \ ...
	    (subsref(P2,substruct('()',{j,':'})) - ...
	    subsref(P1,substruct('()',{j,j+1:m1})) * ...
	    subsref(Q,substruct('()',{j+1:m1,':'}))));
    end
    Q = trim(Q);
    [type,m,r] = match_type(P1,P2);
    Q = set(Q,'type',type,'m',m,'r',r);
    return
end

% if we get to this point, division has failed

error('left division failed');