1% Sparse operator for differentiation.
2%
3% diffx = sparse_diff(sz1,m,cyclic)
4%
5% Sparse operator for differentiation along dimension m for "collapsed" matrix
6% of the size sz1.
7%
8% Input:
9%   sz1: size of rhs
10%   m: dimension to differentiate
11%   cyclic: true if domain is cyclic along dimension m. False is the
12%   default value
13
14function diffx = sparse_diff(sz1,m,cyclic)
15
16if nargin == 2
17    cyclic = false;
18end
19
20n1 = prod(sz1);
21sz2 = sz1;
22
23if ~cyclic
24    sz2(m) = sz2(m)-1;
25end
26
27n2 = prod(sz2);
28
29n = length(sz1);
30
31vi = cell(1,n);
32for i=1:n
33    vi{i} = 1:sz2(i);
34end
35
36IJ = cell(n,1);
37
38[IJ{:}] = ndgrid(vi{:});
39
40for i=1:n
41    IJ{i} = IJ{i}(:);
42end
43
44L1 = [1:n2]';
45
46L2 = sub2ind(sz1,IJ{:});
47one = ones(size(L1));
48
49IJ{m} = IJ{m} + 1;
50
51if cyclic
52    IJ{m} = mod(IJ{m}-1,sz1(m))+1;
53end
54
55
56L2o = sub2ind(sz1,IJ{:});
57
58diffx = sparse( ...
59    [L1;     L1;  ],  ...
60    [L2;     L2o; ],  ...
61    [-one;   one  ], n2 , n1 );
62
63
64
65% Copyright (C) 2009,2012 Alexander Barth <a.barth@ulg.ac.be>
66%
67% This program is free software; you can redistribute it and/or modify
68% it under the terms of the GNU General Public License as published by
69% the Free Software Foundation; either version 2 of the License, or
70% (at your option) any later version.
71%
72% This program is distributed in the hope that it will be useful,
73% but WITHOUT ANY WARRANTY; without even the implied warranty of
74% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
75% GNU General Public License for more details.
76%
77% You should have received a copy of the GNU General Public License
78% along with this program; If not, see <http://www.gnu.org/licenses/>.
79
80
81% LocalWords:  diffx sz
82