function bd = mne_block_diag(A,n) Make or extract a sparse block diagonal matrix 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.
0001 function bd = mne_block_diag(A,n); 0002 % 0003 % function bd = mne_block_diag(A,n) 0004 % 0005 % Make or extract a sparse block diagonal matrix 0006 % 0007 % If A is not sparse, then returns a sparse block diagonal "bd", diagonalized from the 0008 % elements in "A". 0009 % "A" is ma x na, comprising bdn=(na/"n") blocks of submatrices. 0010 % Each submatrix is ma x "n", and these submatrices are 0011 % placed down the diagonal of the matrix. 0012 % 0013 % If A is already sparse, then the operation is reversed, yielding a block 0014 % row matrix, where each set of n columns corresponds to a block element 0015 % from the block diagonal. 0016 % 0017 % Routine uses NO for-loops for speed considerations. 0018 0019 0020 0021 % 0022 % Principal Investigators and Developers: 0023 % ** Richard M. Leahy, PhD, Signal & Image Processing Institute, 0024 % University of Southern California, Los Angeles, CA 0025 % ** John C. Mosher, PhD, Biophysics Group, 0026 % Los Alamos National Laboratory, Los Alamos, NM 0027 % ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory, 0028 % CNRS, Hopital de la Salpetriere, Paris, France 0029 % 0030 % Copyright (c) 2005 BrainStorm by the University of Southern California 0031 % This software distributed under the terms of the GNU General Public License 0032 % as published by the Free Software Foundation. Further details on the GPL 0033 % license can be found at http://www.gnu.org/copyleft/gpl.html . 0034 % 0035 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 0036 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 0037 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 0038 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 0039 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 0040 % 0041 % Author: John C. Mosher 1993 - 2004 0042 % 0043 % 0044 % Modifications for mne Matlab toolbox 0045 % 0046 % Matti Hamalainen 0047 % 2006 0048 % Revision 1.2 2006/04/23 15:29:40 msh 0049 % Added MGH to the copyright 0050 % 0051 % Revision 1.1 2006/04/18 20:44:46 msh 0052 % Added reading of forward solution. 0053 % Use length instead of size when appropriate 0054 % 0055 % 0056 0057 0058 if(~issparse(A)), % then make block sparse 0059 [ma,na] = size(A); 0060 bdn = na/n; % number of submatrices 0061 0062 if(bdn - fix(bdn)), 0063 error('Width of matrix must be even multiple of n'); 0064 end 0065 0066 tmp = reshape([1:(ma*bdn)]',ma,bdn); 0067 i = zeros(ma*n,bdn); 0068 for iblock = 1:n, 0069 i((iblock-1)*ma+[1:ma],:) = tmp; 0070 end 0071 0072 i = i(:); % row indices foreach sparse bd 0073 0074 0075 j = [1:na]; 0076 j = j(ones(ma,1),:); 0077 j = j(:); % column indices foreach sparse bd 0078 0079 bd = sparse(i,j,A(:)); 0080 0081 else % already is sparse, unblock it 0082 0083 [mA,na] = size(A); % matrix always has na columns 0084 % how many entries in the first column? 0085 bdn = na/n; % number of blocks 0086 ma = mA/bdn; % rows in first block 0087 0088 % blocks may themselves contain zero entries. Build indexing as above 0089 tmp = reshape([1:(ma*bdn)]',ma,bdn); 0090 i = zeros(ma*n,bdn); 0091 for iblock = 1:n, 0092 i((iblock-1)*ma+[1:ma],:) = tmp; 0093 end 0094 0095 i = i(:); % row indices foreach sparse bd 0096 0097 0098 j = [0:mA:(mA*(na-1))]; 0099 j = j(ones(ma,1),:); 0100 j = j(:); 0101 0102 i = i + j; 0103 0104 bd = full(A(i)); % column vector 0105 bd = reshape(bd,ma,na); % full matrix 0106 end 0107 0108 return