function TestGenBlockSparse()

clc
close all
% Example 1. Single Block
m = 1;
n = 1000;
p = 0;
q = .003;
A = GenBlockSparse(m,n,p,q);
spy(A)
title(sprintf('Expected Number of Nonzeros = %6.0f',(m*m*(1-p))*(n*n*q)))

% Example 2
figure
m = 20;
n = 100;
p = .8;
q = .01;
A = GenBlockSparse(m,n,p,q);
spy(A)
title(sprintf('Expected Number of Nonzeros = %6.0f',(m*m*(1-p))*(n*n*q)))
shg


function A = GenBlockSparseCVL(m,n,p,q)

% B is an m-by-m 0-1 matrix with the property that if B(i,j)==1, then
%    the (i,j) block of A is nonzero...
B = rand(m,m)<=1-p;
% i and j have the property that blocks (i(k),j(k)) are nonzero for
% k=1:length(i)...
[i j] = find(B);
% Build up the index vectors that define A...
idx = []; jdx = [];
for k=1:length(i)
    % Generate a 0-1 matrix that tells where the nonzeros are in 
    % block (i(k),j(k))...
    C = rand(n,n)<=q;
    [i0,j0] = find(C);
    idx = [idx ; (i(k)-1)*n+i0];
    jdx = [jdx ; (j(k)-1)*n+j0];
end
A = sparse(idx,jdx,ones(length(idx),1),m*n,m*n);

        