Code covered by the BSD License  

Highlights from
block diagonal representation

from block diagonal representation by David Holdaway
Converts a matrix into a cell array of different sized matrices that make up the blocks.

blkdiagconv( matform,blocksizes )
function [ cellform,blocksizes ] = blkdiagconv( matform,blocksizes )
%blkdiagconv: Converts a full / sparse block diagonal matrix into a cell
%array of blocks of the appropriate sizes.
%if 'blocksizes' given just takes the blocks, this should be a k by 2 matrix
%where k is the number of blocks, else work out block sizes

%to go back to a full block diagonal matrix use 
%matform = blkdiag(cellform{}); 

if iscell(matform)==1
    error('input is cell, must be full or sparse matrix')
end

if nargin==2
    
    numblocks  = size(blocksizes,1);
    cellform(numblocks) = [];
    countpar1 = 0; countpar2 = 0; cellform{1} = []; 
    
    for k = 1:numblocks
    
        cellform{k} = matform(countpar1+(1:blocksizes(k,1))...
            ,countpar2+(1:blocksizes(k,2)));
        countpar1= countpar1 + blocksizes(k,1);
        countpar2= countpar2 + blocksizes(k,2);
        
    end
   return 
end

%% meth1

sz = size(matform); 

%calculate max and min position of non zeros entries for all rows / columns
rngud = zeros(sz(1),2);

for k =1:sz(1)
   if isempty(find(matform(k,:),1))       
%empty rnglr is a full row of zeros, empty rngud is a full column of zeros
%just treat these as being of the same length as the previous row / column
       if k ==1
       rngud(k,:) =[1,1];    
       else
       rngud(k,:) = rngud(k-1,:);
       end
   else
    rngud(k,1) =  find(matform(k,:),1);
   rngud(k,2) =  find(matform(k,:),1,'last');      
   end    


end
rnglr = zeros(sz(2),2);
for k =1:sz(2)
   if isempty( find(matform(:,k),1))
       if k ==1
       rnglr(k,:) =[1,1];    
       else
       rnglr(k,:) = rnglr(k-1,:);
       end
   else  
   rnglr(k,1) =  find(matform(:,k),1);  
   rnglr(k,2) =  find(matform(:,k),1,'last'); 
   end
end

%when no rows past row 'k' have a first element before the last in all the
%previous rows this denotes a new block top, likewise when all columns past
% column 'j' have first elements after the last element for columns <=j a
% new block side is created
numblocksud=0;
for k = 1:sz(1)

    if all(rngud(k+1:sz(1),1)>max(rngud(1:k,2)))
        numblocksud = numblocksud+1;
        bldiveud(numblocksud) = k; %#ok<AGROW>    
    end
end
numblockslr=0;
for j = 1:sz(2)
    if all(rnglr(j+1:sz(2),1)>max(rnglr(1:j,2)))
        numblockslr = numblockslr+1;
        bldivelr(numblockslr) = j; %#ok<AGROW>    
    end
end
%intersections between these vertical and horizonal lines make blocks, the
% first row line intersecting the first column line makes the first box, 
%the second row line intersecting the first makes the top of the next and
% intersecting with the second makes the bottom etc etc
verttp = [1,1]; numblocks = min(numblocksud,numblockslr);
for k = 1:numblocks
    
    vertbot = [bldiveud(k),bldivelr(k)];
    cellform{k} = matform(verttp(1):vertbot(1),verttp(2):vertbot(2)); %#ok<AGROW>
    blocksizes(k,:)= [length(verttp(1):vertbot(1)),length(verttp(2):vertbot(2))];
    verttp = vertbot+1;
    
end

return


Contact us