%bandb.m extracts a block k diagonal banded matrix out of A (mu=semi-bandwidth)
%     B   = bandb(A,mu)         % Simple banded mu matrix
%[B,p,p2] = bandb(A, k, mL,mU)  % Diagonal with L / U bands  mL/mU blocks of size k
%[B,p,p2] = bandb(A,-k, mL,mU)  % Off-Diagonal with band mu blocks of size |k| 
%  where p is the magic permutation for turning a banded matrix to 
%              a block matrix of diagonal blocks (BMDB) and
%       p2 is the reverse permutation (Transpose-like)
%______________________________________
%ALSO  banda => banded arrow extraction
%      bandg => banded matrix with a gap
%e.g. A=rand(64); [B p1 p2]=bandb(A,-8,1,0);  [C p3 p4]=bandb(A,16,1,1);
%% compare the following
%  figure,spy(B),nx=8;mygrid;figure,spy(B(p1,p1)),mygrid;figure,spy(B(p2,p2))
%   figure,spy(C),nx=16;ny=16;mygrid;figure,spy(C(p3,p3)),figure,spy(C(p4,p4))

%Ex_1=rand(64); [B p0 p1]=bandb(Ex_1,16,3,9); close all; nx=16;ny=16; %---------
%    A=B(p0,p0); figure; subplot(121); spy(A) % Block size=16 lower=3/upper=9
%    title('[1] Ex1 : A block matrix of diagonal blocks (BMDB)'); mygrid
%    C=A(p1,p1); subplot(122);spy(C);title('[2] Permuted to block forms');mygrid
%    disp('Check that C is the same as B (the nice matrix)')
%Ex_2=rand(24); [B p p2]=bandb(Ex_2,-4,0,0); nx=4;ny=4; %-----------------------
%    figure; subplot(121); spy(B),
%    title('[1] Ex 2: A block matrix of diagonal blocks (BMDB)'); mygrid;nx=6;ny=6;
%    C=B(p2,p2); subplot(122);spy(C);title('[2] Permuted to block forms');mygrid
%    disp('Check that C is again the same as B (the nice matrix)')
function [B,p,p2]=bandb(A,k,m_L,m_U)
 if nargin<1,more on,help bandb; more off; return; end
  if nargin<2,m_U=1; m_L=1; k=1; end
  if nargin<3,m_L=k; m_U=k; k=1; end
  if k==0, help bandb; k, return; end
N = size(A,1); no = N/abs(k); p=[]; % no of blocks
  p2=[];
  for j=1:abs(k)
  for i=1:no
  p2=[p2 j+(i-1)*abs(k)]; 
  end
  end
  for i=1:no
  for j=1:abs(k)
  p=[p i+(j-1)*no]; 
  end
  end
if no==N
 B=tril(tril(A,m_U)',m_L)';
elseif k>0
 B=[];
 for block=1:no
 ii=((block-1)*k+1) : (block*k);
 B(ii,ii)=tril(tril(A(ii,ii),m_U)',m_L)';
 end
else
 B=[]; k=-k;
 for block=1:no
 ss=((block-2)*k+1) : ((block-1)*k); %Small
 ii=((block-1)*k+1) : (block*k);
 ll=(block*k+1) : ((block+1)*k); % Large
 B(ii,ii)=diag(diag(A(ii,ii)));
   if block<no
   B(ii,ll)=tril(tril(A(ii,ll),m_U)',m_L)';
   end
   if block>1
   B(ii,ss)=tril(tril(A(ii,ss),m_U)',m_L)';
   end
 end % block
end % if no

