function [res1,res2] = spdiags(arg1,arg2,arg3,arg4)
%SPDIAGS Sparse matrix formed from diagonals.
%   SPDIAGS, which generalizes the function "diag", deals with three
%   matrices, in various combinations, as both input and output.
%
%   [B,d] = SPDIAGS(A) extracts all nonzero diagonals from the m-by-n
%   matrix A.  B is a min(m,n)-by-p matrix whose columns are the p
%   nonzero diagonals of A.  d is a vector of length p whose integer
%   components specify the diagonals in A.
%
%   B = SPDIAGS(A,d) extracts the diagonals specified by d.
%   A = SPDIAGS(B,d,A) replaces the diagonals of A specified by d with
%       the columns of B.  The output is sparse.
%   A = SPDIAGS(B,d,m,n) creates an m-by-n sparse matrix from the
%       columns of B and places them along the diagonals specified by d.
%
%   Roughly, A, B and d are related by
%       for k = 1:p
%           B(:,k) = diag(A,d(k))
%       end
%   Some elements of B, corresponding to positions "outside" of A, are
%   not actually used.  They are not referenced when B is input and
%   are set to zero when B is output.
%
%   Example: These commands generate a sparse tridiagonal representation
%   of the classic second difference operator on n points.
%       e = ones(n,1);
%       A = spdiags([e -2*e e], -1:1, n, n)
%
%   See also DIAG.

%   Rob Schreiber, RIACS, and Cleve Moler, 2/13/91, 6/1/92.
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.7 $  $Date: 1997/11/21 23:44:42 $


if nargin <= 2
   if size(arg1,1) ~= size(arg1,2) & min(size(arg1))==1,
      error('A cannot be a vector.');
   end
   % Extract diagonals
   A = arg1;
   if nargin == 1
      % Find all nonzero diagonals
      [i,j] = find(A);
      d = sort(j-i);
      d = d(find(diff([-inf; d])));
   else
      % Diagonals are specified
      d = arg2(:);
   end
   [m,n] = size(A);
   p = length(d);
   B = zeros(min(m,n),p);
   for k = 1:p
      if m >= n
         i = max(1,1+d(k)):min(n,m+d(k));
      else
         i = max(1,1-d(k)):min(m,n-d(k));
      end
      if ~isempty(i), B(i,k) = diag(A,d(k)); end
   end
   res1 = B;
   res2 = d;
end

if nargin >= 3
   % Create a sparse matrix from its diagonals
   B = arg1;
   d = arg2(:);
   p = length(d);
   moda = (nargin == 3);
   if moda
      % Modify a given matrix
      A = arg3;
   else
      % Start from scratch
      A = sparse(arg3,arg4);
   end
   % Process A in compact form
   [i,j,a] = find(A);
   a = [i j a];
   [m,n] = size(A);
   for k = 1:p
      if moda
         % Delete current d(k)-th diagonal
         i = find(a(:,2)-a(:,1) == d(k));
         a(i,:) = [];
      end 
      % Append new d(k)-th diagonal to compact form
      i = (max(1,1-d(k)):min(m,n-d(k)))';
      a = [a; i i+d(k) B(i+(m>=n)*d(k),k)];
   end
   res1 = sparse(a(:,1),a(:,2),full(a(:,3)),m,n);
end
