|
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 |
|
|
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
|
|
16 |
|
|
17 ## -*- texinfo -*- |
|
|
18 ## @deftypefn {function File} {[@var{b}, @var{c}]} = spdiags (@var{a}) |
|
|
19 ## @deftypefnx {function File} {@var{b}} = spdiags (@var{a}, @var{c}) |
|
|
20 ## @deftypefnx {function File} {@var{b}} = spdiags (@var{v}, @var{c}, @var{a}) |
|
|
21 ## @deftypefnx {function File} {@var{b}} = spdiags (@var{v}, @var{c}, @var{m}, @var{n}) |
|
|
22 ## A generalization of the function @code{spdiag}. Called with a single |
|
|
23 ## input argument, the non-zero diagonals @var{c} of @var{A} are extracted. |
|
|
24 ## With two arguments the diagonals to extract are given by the vector |
|
|
25 ## @var{c}. |
|
|
26 ## |
|
|
27 ## The other two forms of @code{spdiags} modify the input matrix by |
|
|
28 ## replacing the diagonals. They use the columns of @var{v} to replace |
|
|
29 ## the columns represented by the vector @var{c}. If the sparse matrix |
|
|
30 ## @var{a} is defined then the diagonals of this matrix are replaced. |
|
|
31 ## Otherwise a matrix of @var{m} by @var{n} is created with the |
|
|
32 ## diagonals given by @var{v}. |
|
|
33 ## |
|
|
34 ## Negative values of @var{c} representive diagonals below the main |
|
|
35 ## diagonal, and positive values of @var{c} diagonals above the main |
|
|
36 ## diagonal. |
|
|
37 ## |
|
|
38 ## For example |
|
|
39 ## |
|
|
40 ## @example |
|
|
41 ## @group |
|
|
42 ## spdiags (reshape (1:12, 4, 3), [-1 0 1], 5, 4) |
|
|
43 ## @result{} 5 10 0 0 |
|
|
44 ## 1 6 11 0 |
|
|
45 ## 0 2 7 12 |
|
|
46 ## 0 0 3 8 |
|
|
47 ## 0 0 0 4 |
|
|
48 ## @end group |
|
|
49 ## @end example |
|
|
50 ## |
|
|
51 ## @end deftypefn |
|
|
52 |
|
|
53 function [A, c] = spdiags(v,c,m,n) |
|
|
54 |
|
|
55 wfi = warn_fortran_indexing; |
|
|
56 unwind_protect |
|
|
57 warn_fortran_indexing = 0; |
|
|
58 |
|
|
59 if nargin == 1 || nargin == 2 |
|
|
60 ## extract nonzero diagonals of v into A,c |
|
|
61 [i,j,v,nr,nc] = spfind(v); |
|
|
62 if nargin == 1 |
|
|
63 c = unique(j-i); # c contains the active diagonals |
|
|
64 endif |
|
|
65 ## FIXME: we can do this without a loop if we are clever |
|
|
66 offset = max(min(c,nc-nr),0); |
|
|
67 A = zeros(min(nr,nc),length(c)); |
|
|
68 for k=1:length(c) |
|
|
69 idx = find(j-i == c(k)); |
|
|
70 A(j(idx)-offset(k),k) = v(idx); |
|
|
71 end |
|
|
72 elseif nargin == 3 |
|
|
73 ## Replace specific diagonals c of m with v,c |
|
|
74 [nr,nc] = size(m); |
|
|
75 B = spdiags(m,c); |
|
|
76 A = m - spdiags(B,c,nr,nc) + spdiags(v,c,nr,nc); |
|
|
77 else |
|
|
78 ## Create new matrix of size mxn using v,c |
|
|
79 [j,i,v] = find(v); |
|
|
80 offset = max(min(c(:),n-m),0); |
|
|
81 j+=offset(i); |
|
|
82 i=j-c(:)(i); |
|
|
83 idx = i>0 & i<=m & j>0 & j<=n; |
|
|
84 A = sparse(i(idx),j(idx),v(idx),m,n); |
|
|
85 |
|
|
86 endif |
|
|
87 |
|
|
88 unwind_protect_cleanup |
|
|
89 warn_fortran_indexing = wfi; |
|
|
90 end_unwind_protect |
|
|
91 |
|
|
92 endfunction |
