Mercurial > hg > octave-nkf

comparison scripts/sparse/spdiags.m @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children 4c8a2e4e0717
comparison
equal deleted inserted replaced
5163:9f3299378193 5164:57077d0ddc8e
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