chainmult.m

Matlab中优化多个矩阵相乘以提高计算效率的源代码。

         q.Lvar=node(1:i);
         q.Rvar=node(i+1:end);
         if ~isempty(field),
            y=orig;
            y=setfield(y,field{:},q);
         else
            y=q;
         end
         TreeStack('push',y);
      end
      stop=1; % forked new trees on stack
   end
else
   y=node;
   [y.Lvar,stop]=SplitOneLevel(node.Lvar,[idx 1],stop);
   if (stop~=1),
      stop=0; % Must reset to 0, since it could be -1 at this point
      [y.Rvar,stop]=SplitOneLevel(node.Rvar,[idx 2],stop);
      if (stop~=1),
         stop=-1;
         y=node; % we're done with this tree
      end
   end
end


% ==================================================================================
% genFwdStruct:
%   Generate structure describing forward-chain only.
%   Use for non-optimized case.
% ==================================================================================
function list = genFwdStruct(N)
if N==0,
    list={};
elseif N==1,
    list = {'A'};
else
    c = char('A'+(0:N-1));
    list.Lvar = c(1);
    list.Rvar = c(2);
    for i=3:N,
        y.Lvar = list;
        y.Rvar = c(i);
        list=y;
    end	
end
list={list};


% ==================================================================================
% TreeStack:
% Implement a stack that can hold structure trees.
% ==================================================================================
function [t,ok]=TreeStack(op,push)

persistent trees;
if ~iscell(trees), trees={}; end

ok = 1;
switch op
case 'pop'
   % pop   
   if isempty(trees),
      ok=0;
      t=[];
      trees={};
   else
      t=trees{end};
      trees=trees(1:end-1);
   end
case 'push'
   trees(end+1)={push};
   t=trees;
case 'query'
   t=trees;
case 'reset'
   trees={};
   t=trees;
otherwise,
	error('Unknown operation requested.');
end


% ==================================================================================
% convertCell2Str:
% Print the chain multiply permutation expression.
%
% chain is a vector cell array of 3-element vector indices
% describing the chain multiplication.
%
% varNames is an optional cell string array containing one
% string per variable name.  These variable names are used
% during construction of the MATLAB string expressions.
% If omitted, variable names are 'A', B', ...

% Ex:
%  {[4 5 -1] [3 -1 -2] [2 -2 -1] [1 -1 0]} % A * (B * (C*(D*E)))
%  {[3 4 -1] [-1 5 -2] [2 -2 -1] [1 -1 0]} % A * (B * ((C*D)*E))
%  {[2 3 -1] [-1 4 -2] [-2 5 -1] [1 -1 0]} % A * (((B*C)*D) * E)
%  {[3 4 -1] [2 -1 -2] [-2 5 -1] [1 -1 0]} % A * ((B*(C*D)) * E)
%  {[2 3 -1] [4 5 -2] [-1 -2 -3] [1 -3 0]} % A * ((B*C) * (D*E))
% ==================================================================================
function str = convertCell2Str(chain,varNames)

if nargin<2 || isempty(varNames),
    % Generate as many default variable names as required
    %
    % xxx Find max non-temporary (positive) index specified
    % across all input vectors:
	varNames = cellstr(char(double('A') + (0:25))');
end

str = '';
temp_strs = {};

for i=1:length(chain),
   x1_idx  = chain{i}(1);  % x1 matrix index
   x2_idx  = chain{i}(2);  % x2 matrix index
   
   if length(chain{i})==3,    % skip single-input case
      x12_idx = chain{i}(3);  % x12 result index
      
      % Get size vectors of x1 and x2
      % Each size vector specifies [rows cols]
      ch1 = getMatrixStr(x1_idx, temp_strs, varNames);
      ch2 = getMatrixStr(x2_idx, temp_strs, varNames);
      
      % Multiply these together:
      
      if x12_idx == 0,
         % Last operation:
         str = ['Y = ' ch1 '*' ch2 ';'];
      else
         temp_strs{-x12_idx} = ['(' ch1 '*' ch2 ')'];
      end
   else
      % 2-input case:
      ch1 = getMatrixStr(x1_idx, temp_strs, varNames);
      str = ['Y = ' ch1 ';'];
   end
end


% -----------------------------------------------------------
function x_str = getMatrixStr(x_idx, temp_strs, varNames)
% getMatrixStr Return string for i'th matrix chain operation
% temp_strs is a list of string expressions describing the
% contents of each generated temp.

if x_idx<0,
   % Get temporary-variable string expression:
   x_str = temp_strs{-x_idx};
   
elseif x_idx>0,
   % Return variable name of i'th input matrix
   % 1st = A, 2nd = B, etc.
   if x_idx>26,
      error('Too many input variables for prettyPrintVars.');
  end
   
   % Get name for x_idx'th variable:
   % x_str = char(double('A') + x_idx-1);
   x_str = varNames{x_idx};
   
else
   error(['Input index to matrix multiplier pair cannot be ' ...
         'specified as the output matrix (index 0).']);
end


% ==================================================================================
% convertCell2Struct:
% ==================================================================================
function c = convertCell2Struct(chain)

c={}; temp={};
for i=1:length(chain),
   if length(chain{i})==3,    % skip single-input case
      x1_idx  = chain{i}(1);  % x1 matrix index
      x2_idx  = chain{i}(2);  % x2 matrix index
      x12_idx = chain{i}(3);  % x12 result index
      
      % Get references to input matrices
      % Use letter variable names, or
      % structure entry from temp storage:
      p.Lvar=getMatrixRef(x1_idx, temp);
      p.Rvar=getMatrixRef(x2_idx, temp);
      
      if x12_idx == 0,
         % output area specified for destination
         % we are done:
         c=p;
      else
         % A temporary storage area has been specified
         tempIdx = -x12_idx; % create positive temp index
         temp{tempIdx} = p;
      end
   end
end


% ----------------------------------------------------
function x_ref = getMatrixRef(x_idx, temp_refs)
% getMatrixRef
%
% Positive indices are taken to be input matrices,
% negative indices are temp matrices.
% The 0'th index is the output matrix, and is not available
% for use in intermediate results (at least in the enumerated
% specifications).

%  If index is negative, return size-vector from temp_sizes
%  Otherwise, return size-vector from matrix_sizes
if x_idx<0,
   x_ref = temp_refs{-x_idx};
elseif x_idx>0,
   x_ref = char('A'-1+x_idx);
else
   error(['Input index to matrix multiplier pair cannot be ' ...
         'specified as the output matrix (index 0).']);
end


% ==================================================================================
% convertStruct2Cell:
% Converts cell-array of structures in P (representing multiple enumerations
%   of the chain multiplication problem) into a cell-array of cell-arrays
%   of "program" vectors.
% ==================================================================================
function y=convertStruct2Cell(pc)

% If this is one structure, wrap it in a cell-array:
single_struct = ~iscell(pc);
if single_struct,
    pc={pc};
end

% Is this a single-variable problem?
if length(pc)==1,
    if ~isstruct(pc{1}),
        y={{[1 0]}};
        return
    end	
end	

for i=1:length(pc),
    p=pc{i};
    
    localRecursion(p,0);
    c=localRecursion(p,1);
    
    % Overwrite last result index with the output index (0):
    c{end}(3)=0;
    y{i}=c;
end

% Convert back from nested cells if needed:
if single_struct,
    y=y{1};
end


% ------------------------------------------------
function c=localRecursion(p,init)

persistent Gtemp chainResult;
if nargin>1,
   if init==0,
      Gtemp=[]; chainResult={};
   else
      c=chainResult;
      return
   end
end

freeList=[];

% Get depth order:
% [1 2] -> Left is deeper than Right
% [2 1] -> Right is deeper than Left
DepthOrder = getDepthOrder(p);
children = {p.Lvar p.Rvar};
children = children(DepthOrder);

for i=1:2,
   % Get "left" contribution
   % If .Lvar is a string, use it directly
   % If .Lvar is a struct, it's the child
   if isstruct(children{i}),
      % Get child computation's temp name
      c=localRecursion(children{i});
      % c(3) is the destination
      result{i}=c(3);  % numeric, neg. temp index
      % Consumed the temp variable
      % record index for free'ing after use
      freeList(end+1)=c(3);
   else
      result{i}=children{i};
   end
   % Convert to # if a character:
   if ischar(result{i}),
      result{i} = result{i}-'A'+1;
   end
end

[T,Gtemp]=allocTemp(Gtemp);
c = [result{DepthOrder} T];

chainResult{end+1}=c;
Gtemp=freeTemp(Gtemp,freeList);


% ----------------------------------------------------
function depthOrder= getDepthOrder(p)

% Determine ordering of depth of subtree p
% 
% [1 2] -> Left is deeper than Right
% [2 1] -> Right is deeper than Left

%depthOrder = [2 1];
%depthOrder = [1 2];

cntL=descendDepth(p.Lvar,0);
cntR=descendDepth(p.Rvar,0);

if cntL>=cntR,
	depthOrder = [1 2];
else
    depthOrder = [2 1];
end


% ----------------------------------------------------
function cnt = descendDepth(p,cnt)

if isstruct(p),
	cntL = descendDepth(p.Lvar,cnt+1);
	cntR = descendDepth(p.Rvar,cnt+1);
    cnt = max(cntL,cntR);
end

% ----------------------------------------------------
function Gtemp=freeTemp(Gtemp,idx)
% Return temp variable - no longer being used:

% Temp index is negative, and so are the entries
% in the global temp list:
if isempty(Gtemp),
   error('No temps in list - cannot free.');
end
for j=1:length(idx),
   i=find(Gtemp==idx(j));
   if isempty(i),
      error('Temp index not found - cannot free');
   end
   Gtemp(i)=[];
end


% ----------------------------------------------------
function [tnew,Gtemp]=allocTemp(Gtemp)
% Determine index of next temp to use
% Maintain cache of all active temp #'s in Gtemp

if ~isempty(Gtemp),
   
   % Get +ve # of all temps currently in use
   posTemp=-Gtemp;
   maxTemp=max(posTemp);
   
   % Determine all temps possible in this range
   tempList=1:maxTemp;
   
   % Find which ones were not used ("holes" in the temp list)
   unused = setdiff(tempList,posTemp);
   if ~isempty(unused),
      % Get lowest-numbered "hole" to use as a temp
      nextTemp=min(unused);
   else
      nextTemp=maxTemp+1;
   end
else
   nextTemp=1;  % no temps used yet
end

% Use -ve # for temp:
tnew=-nextTemp;

% Put into global list:
Gtemp(end+1)=tnew;


% ==================================================================================
% convertStruct2Str:
% ==================================================================================
function c=convertStruct2Str(p)
L=p.Lvar;
if isstruct(L),
   L=convertStruct2Str(L);
end
R=p.Rvar;
if isstruct(R),
   R=convertStruct2Str(R);
end
c=['(' L '*' R ')'];


% [EOF] chainmult.m