function [A] = lap_bound_mat(A, b, varargin)

% If matrix dimensions are passed in, then use them to force
% boundary consitions on the outer rows and columns (ie. the
% surrounding rectangle)  Otherwise, boundary conditions are
% assumed only for non-zero entries
%
% If passed in args do not match dimensions for b, they are
% ignored.

sz = length(b);
m = floor(sqrt(sz));
n = m;
if nargin > 3
  % Rectangular matrix
  arg1 = varargin{1};
  arg2 = varargin{2};
  rows = arg1(1);
  cols = arg2(1);
  if ((rows*cols) == sz)
    m = rows;
    n = cols;
  end

  % now reshape b into a matrix, and force the boundary columns to
  % non-zero
  c = reshape(b, n ,m)';
  c(1, 1:n) = 1;
  c(m, 1:n) = 1;
  c(1:m, 1) = 1;
  c(1:m, n) = 1;
  b = reshape(c', m*n, 1);
end

% Create an identity matrix with ones in the boundary rows, and
% zeros elsewhere
i = find(b);
I = speye(sz);
if length(i) < 1
  A = I;
else
  % This takes FOREVER on large matrices, even sparse ones!
  %A(i,:) = I(i,:);

  % Use a different method; extract diagonals and operate on them only
  % Extend A on left and right so index errors don't occur
  A = [sparse(sz,n) A sparse(sz,n)];

  [B, d] = spdiags(A);
  B(i,:) = 0;
  B(i,3) = 1;
  A = spdiags(B, d, sz, (sz+2*n));

  % Reduce A back to original size
  A = A(:, ((n+1):(sz + n)));
end