--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ssm/blockMatMult.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,96 @@
+% multiplies block defined matrix stored inside cell array
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: cell_mult multiplies block defined matrix stored inside cell
+% array
+%
+% CALL: [cell3] = ssm.blockMatMult(cell1,cell2,isnotempty_1,isnotempty_2)
+%
+% INPUTS:
+%       cell1 - cell array of matrices representing a matrix by blocs.
+%               blocs may be empty
+%       cell2 - cell array of matrices representing a matrix by blocs.
+%               blocs may be empty
+%
+% OPTIONNAL INPUTS:
+%       isnotempty_1 - logical array tells ~isequal(cell1{ii,jj},[])
+%       isnotempty_2 - logical array tells ~isequal(cell2{ii,jj},[])
+%
+% OUTPUTS:
+%       cell3 - cell array of matrices representing a matrix by blocs.
+%               blocs may be empty
+%
+% NOTE : function is private to the ssm class
+%
+% VERSION: '$Id: blockMatMult.m,v 1.7 2011/04/08 08:56:24 hewitson Exp $'
+%
+% TO DO :
+% check ME in case of mixed symbolic and double
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function c = blockMatMult(varargin)
+  a=varargin{1};
+  b =varargin{2};
+  n1 = size(a,1);
+  n2 = size(a,2);
+  n3 = size(b,2);
+  if nargin == 4
+    sizes1 = varargin{3};
+    sizes3 = varargin{4};
+    if numel(sizes1)~=n1
+      error('Incompatible user input (3rd argument does not correctly indicate final line heights)');
+    elseif numel(sizes3)~=n3
+      error('Incompatible user output (4th argument does not correctly indicate final column widths)');
+    end
+  end
+  c = cell(n1,n3);
+  if isempty(a)||isempty(b)
+    % if the content matrix is empty
+    if n2==0 && n1>0 && n3>0
+      if nargin == 4
+
+        c = ssm.blockMatFillDiag(c,sizes1,sizes3);
+      else
+        error('cannot build proper matrices')
+      end
+    end
+  else
+    % if matrices are not empty
+    isnotempty_a = not(cellfun(@isempty, a));
+    isnotempty_b = not(cellfun(@isempty, b));
+    % extra check for the extended diagonal
+    for ii=1:max(n1, n2)
+      isnotempty_a(min(ii,n1),min(ii,n2)) = true;
+    end
+    for ii=1:max(n2, n3)
+      isnotempty_b(min(ii,n2), min(ii,n3)) = true;
+    end
+    isempty_c = true(n1,n3);
+    for ii=1:n1
+      for kk=1:n2
+        if isnotempty_a(ii,kk)
+          for jj=1:n3
+            if isnotempty_b(kk,jj)
+              if isempty_c(ii,jj)
+                % note that the matrix multiplication preserves a non-empty block diagonal
+                if  (isempty(a{ii,kk}) || isempty(b{kk,jj}))
+                  % exception for the case where you have one or more empty matrix
+                  c{ii,jj} = zeros(size(a{ii,kk},1), size(b{kk,jj},2));
+                else
+                  c{ii,jj} = a{ii,kk}*b{kk,jj};
+                end
+                isempty_c(ii,jj) = false;
+              else
+                if (isempty(a{ii,kk}) || isempty(b{kk,jj}))
+                  % exception for the case where you have one or more empty matrix
+                else
+                  c{ii,jj} = c{ii,jj} + a{ii,kk}*b{kk,jj};
+                end
+              end
+            end
+          end
+        end
+      end
+    end
+  end
+end