blk_diag

(Toolbox/blk_diag.m in BrainStorm 2.0 (Alpha))

Function Synopsis

bd = blk_diag(A,n);

Help Text

BLK_DIAG - Make or extract a sparse block diagonal matrix
 function bd = blk_diag(A,n);
 If A is not sparse, then
 returns a sparse block diagonal "bd", diagonalized from the
 elements in "A".
 "A" is ma x na, comprising bdn=(na/"n") blocks of submatrices.
 Each submatrix is ma x "n", and these submatrices are
 placed down the diagonal of the matrix.

 If A is already sparse, then the operation is reversed, yielding a block
 row matrix, where each set of n columns corresponds to a block element
 from the block diagonal.

 Routine uses NO for-loops for speed considerations.

Cross-Reference Information

cross_mat                          C:\BrainStorm_2001\Toolbox\cross_mat.m
montreal_data_create               C:\BrainStorm_2001\Toolbox\montreal_data_create.m
os_meg                             C:\BrainStorm_2001\Toolbox\os_meg.m
simulator                          C:\BrainStorm_2001\Toolbox\simulator.m

Listing of function C:\BrainStorm_2001\Toolbox\blk_diag.m

function bd = blk_diag(A,n);
%BLK_DIAG - Make or extract a sparse block diagonal matrix
% function bd = blk_diag(A,n);
% If A is not sparse, then
% returns a sparse block diagonal "bd", diagonalized from the
% elements in "A".
% "A" is ma x na, comprising bdn=(na/"n") blocks of submatrices.
% Each submatrix is ma x "n", and these submatrices are
% placed down the diagonal of the matrix.
%
% If A is already sparse, then the operation is reversed, yielding a block
% row matrix, where each set of n columns corresponds to a block element
% from the block diagonal.
%
% Routine uses NO for-loops for speed considerations.

%<autobegin> ---------------------- 14-Jun-2004 17:09:48 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Utility - Numeric
%
% At Check-in: $Author: Mosher $  $Revision: 15 $  $Date: 6/14/04 3:37p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 14-Jun-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 14-Jun-2004 17:09:48 -----------------------

% ----------------------------- Script History ---------------------------------
% July 29, 1993 Author
% September 28, 1993 JCM Conversion to sparse
% July 27, 1995 JCM inverse block diagonal added
% JCM 19-May-2004  Comments cleaning
% ----------------------------- Script History ---------------------------------


if(~issparse(A)),        % then make block sparse
  [ma,na] = size(A);
  bdn = na/n;             % number of submatrices

  if(bdn - fix(bdn)),
    error('Width of matrix must be even multiple of n');
  end

if(0)
  i = [1:ma]';
  i = i(:,ones(1,n));
  i = i(:);             % row indices first submatrix
  
  ml = length(i);         % ma*n
  
  % ndx = [0:(bdn-1)]*ma;     % row offsets per submatrix
  ndx = [0:ma:(ma*(bdn-1))];     % row offsets per submatrix
  
  i = i(:,ones(1,bdn)) + ndx(ones(ml,1),:);
else
  tmp = reshape([1:(ma*bdn)]',ma,bdn);
  i = zeros(ma*n,bdn);
  for iblock = 1:n,
    i((iblock-1)*ma+[1:ma],:) = tmp;
  end
end
  
  i = i(:);             % row indices foreach sparse bd
  
  
  j = [1:na];
  j = j(ones(ma,1),:);
  j = j(:);             % column indices foreach sparse bd
  
  bd = sparse(i,j,A(:));

else                 % already is sparse, unblock it
  
  [mA,na] = size(A);        % matrix always has na columns
  % how many entries in the first column?
  bdn = na/n;            % number of blocks
  ma = mA/bdn;            % rows in first block
  
  % blocks may themselves contain zero entries.  Build indexing as above
if(0)
  i = [1:ma]';
  i = i(:,ones(1,n));
  i = i(:);             % row indices first submatrix
  
  ml = length(i);         % ma*n
  
  % ndx = [0:(bdn-1)]*ma;     % row offsets per submatrix
  ndx = [0:ma:(ma*(bdn-1))];     % row offsets per submatrix
  
  i = i(:,ones(1,bdn)) + ndx(ones(ml,1),:);
else
  tmp = reshape([1:(ma*bdn)]',ma,bdn);
  i = zeros(ma*n,bdn);
  for iblock = 1:n,
    i((iblock-1)*ma+[1:ma],:) = tmp;
  end  
end

  i = i(:);             % row indices foreach sparse bd


if(0)  
  j = [1:na];
  j = j(ones(ma,1),:);
  j = j(:);             % column indices foreach sparse bd
  
  % so now we have the complete two dimensional indexing. Convert to
  % one dimensional
  
  i = i + (j-1)*mA;
else
  j = [0:mA:(mA*(na-1))];
  j = j(ones(ma,1),:);
  j = j(:);
  
  i = i + j;
end
  
  bd = full(A(i));     % column vector
  bd = reshape(bd,ma,na);    % full matrix
end  

return