|
5164
|
1 ## Copyright (C) 2000-2001 Paul Kienzle |
|
|
2 ## |
|
|
3 ## This program is free software; you can redistribute it and/or modify |
|
|
4 ## it under the terms of the GNU General Public License as published by |
|
|
5 ## the Free Software Foundation; either version 2 of the License, or |
|
|
6 ## (at your option) any later version. |
|
|
7 ## |
|
|
8 ## This program is distributed in the hope that it will be useful, |
|
|
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
11 ## GNU General Public License for more details. |
|
|
12 ## |
|
|
13 ## You should have received a copy of the GNU General Public License |
|
|
14 ## along with this program; if not, write to the Free Software |
|
5307
|
15 ## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA |
|
|
16 ## 02110-1301 USA |
|
5164
|
17 |
|
|
18 ## -*- texinfo -*- |
|
|
19 ## @deftypefn {function File} {[@var{b}, @var{c}]} = spdiags (@var{a}) |
|
|
20 ## @deftypefnx {function File} {@var{b}} = spdiags (@var{a}, @var{c}) |
|
|
21 ## @deftypefnx {function File} {@var{b}} = spdiags (@var{v}, @var{c}, @var{a}) |
|
|
22 ## @deftypefnx {function File} {@var{b}} = spdiags (@var{v}, @var{c}, @var{m}, @var{n}) |
|
|
23 ## A generalization of the function @code{spdiag}. Called with a single |
|
|
24 ## input argument, the non-zero diagonals @var{c} of @var{A} are extracted. |
|
|
25 ## With two arguments the diagonals to extract are given by the vector |
|
|
26 ## @var{c}. |
|
|
27 ## |
|
|
28 ## The other two forms of @code{spdiags} modify the input matrix by |
|
|
29 ## replacing the diagonals. They use the columns of @var{v} to replace |
|
|
30 ## the columns represented by the vector @var{c}. If the sparse matrix |
|
|
31 ## @var{a} is defined then the diagonals of this matrix are replaced. |
|
|
32 ## Otherwise a matrix of @var{m} by @var{n} is created with the |
|
|
33 ## diagonals given by @var{v}. |
|
|
34 ## |
|
|
35 ## Negative values of @var{c} representive diagonals below the main |
|
|
36 ## diagonal, and positive values of @var{c} diagonals above the main |
|
|
37 ## diagonal. |
|
|
38 ## |
|
|
39 ## For example |
|
|
40 ## |
|
|
41 ## @example |
|
|
42 ## @group |
|
|
43 ## spdiags (reshape (1:12, 4, 3), [-1 0 1], 5, 4) |
|
|
44 ## @result{} 5 10 0 0 |
|
|
45 ## 1 6 11 0 |
|
|
46 ## 0 2 7 12 |
|
|
47 ## 0 0 3 8 |
|
|
48 ## 0 0 0 4 |
|
|
49 ## @end group |
|
|
50 ## @end example |
|
|
51 ## |
|
|
52 ## @end deftypefn |
|
|
53 |
|
|
54 function [A, c] = spdiags(v,c,m,n) |
|
|
55 |
|
|
56 wfi = warn_fortran_indexing; |
|
|
57 unwind_protect |
|
|
58 warn_fortran_indexing = 0; |
|
|
59 |
|
|
60 if nargin == 1 || nargin == 2 |
|
|
61 ## extract nonzero diagonals of v into A,c |
|
|
62 [i,j,v,nr,nc] = spfind(v); |
|
|
63 if nargin == 1 |
|
|
64 c = unique(j-i); # c contains the active diagonals |
|
|
65 endif |
|
|
66 ## FIXME: we can do this without a loop if we are clever |
|
|
67 offset = max(min(c,nc-nr),0); |
|
|
68 A = zeros(min(nr,nc),length(c)); |
|
|
69 for k=1:length(c) |
|
|
70 idx = find(j-i == c(k)); |
|
|
71 A(j(idx)-offset(k),k) = v(idx); |
|
|
72 end |
|
|
73 elseif nargin == 3 |
|
|
74 ## Replace specific diagonals c of m with v,c |
|
|
75 [nr,nc] = size(m); |
|
|
76 B = spdiags(m,c); |
|
|
77 A = m - spdiags(B,c,nr,nc) + spdiags(v,c,nr,nc); |
|
|
78 else |
|
|
79 ## Create new matrix of size mxn using v,c |
|
|
80 [j,i,v] = find(v); |
|
|
81 offset = max(min(c(:),n-m),0); |
|
|
82 j+=offset(i); |
|
|
83 i=j-c(:)(i); |
|
|
84 idx = i>0 & i<=m & j>0 & j<=n; |
|
|
85 A = sparse(i(idx),j(idx),v(idx),m,n); |
|
|
86 |
|
|
87 endif |
|
|
88 |
|
|
89 unwind_protect_cleanup |
|
|
90 warn_fortran_indexing = wfi; |
|
|
91 end_unwind_protect |
|
|
92 |
|
|
93 endfunction |
