%
% Copyright (C) May/June 2005 by Henrik Pagenkopf
%
function P = product(A, B, zero, rec_depth) % approximates A* B via interval sectioning with "regula-falsi" for square matrices.
%
% - cond: A and B are square matrices and have same sizes > 0.
%
sa = size(A);
sb = size(B);
P = zeros(sa(1)); % assign output arguments in case of recursion.
if sa(1) >= 1 && sa(1)== sa(2) && sb(1)== sb(2) && sa(1)== sb(1) % => sa(2)== sb(2)
s = sa(1);
S = P; % deep-copy s.b.
[rA,cA,dA]= max_coord_dist(A);
[rB,cB,dB]= max_coord_dist(B);
if dA > zero || dB > zero
rsumA = zeros(s,1); % column vector/tuple for row sums of A.
% -> to be refined (e.g. via partial-vector sums)...
rsumB = zeros(s,1); % ...the same for B.
ravgA = zeros(s,1); % column vector/tuple for row averages of A.
ravgB = zeros(s,1); % ...for B.
csumA = zeros(1,s); % row vector/tuple for column sums of A.
csumB = zeros(1,s);
cavgA = zeros(1,s); % row vector/tuple for column averages of A.
cavgB = zeros(1,s);
for j = 1: s
for k = 1: s
rsumA(j)= rsumA(j)+ A(j,k);
rsumB(j)= rsumB(j)+ B(j,k);
end;
ravgA(j)= rsumA(j)/s; % the average scalar value of j-th row of A.
ravgB(j)= rsumB(j)/s; % ...for B.
end; % for
for k = 1: s
for j = 1: s
csumA(k)= csumA(k)+ A(j,k);
csumB(k)= csumB(k)+ B(j,k);
end;
cavgA(k)= csumA(k)/s; % the average scalar value of k-th column of A.
cavgB(k)= csumB(k)/s; % ...for B.
end; % for
C = zeros(s); % approximate A* B =approx.= C = A *([cavgB, cavgB,..., cavgB]')
% where ' denotes the symbolic matrix transposition.
for j = 1: s
for k = 1: s
C(j,k)= rsumA(j)* cavgB(k);
end;
end; % for
S = C;
if rec_depth > 0 % make use of the distributive laws in the ring of square matrices of constant size s.
S = correct_PR_column(A, B, C, cB);
S = correct_PR_row (A, B, S, rA);
DeltaA = zeros(s);
DeltaB = zeros(s);
for k = 1: s
DeltaA(:,k)= A(:,k)- ravgA;
DeltaB(k,:)= B(k,:)- cavgB;
end;
DeltaA(rA,:)= zeros(1,s); % re-correct after row correcture of C into S .
DeltaB(:,cB)= zeros(s,1); % re-correct after column correcture of S (former C).
S = S + product(tril(DeltaA, 0)/2, B, zero, rec_depth -1); % split factors Delta<.> into lower + upper triangular parts,
S = S + product(triu(DeltaA, 1)/2, B, zero, rec_depth -1); % the main diagonal only once, to avoid zero avg matrices;
S = S + product(A, tril(DeltaB, 0)/2, zero, rec_depth -1); % 2 ways of approximation in parallel, using linearity (* 0.5)
S = S + product(A, triu(DeltaB, 1)/2, zero, rec_depth -1); % of the participating matrices.
end % if(rec_depth)
end; % if(dA)
P = correct_PR(A, B, S, rA, cA, rB, cB);
else % error
disp('Dimension size error.');
end
% return.
