Mercurial > hg > octave-nkf

annotate scripts/sparse/spdiags.m @ 14953:711e64a11d36

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix some segfaults
author Max Brister <max@2bass.com>
date Fri, 15 Jun 2012 13:15:48 -0500
parents 1e77f6078692
children 1c89599167a6
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
14138
72c96de7a403 maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents: 12558
diff changeset
1 ## Copyright (C) 2000-2012 Paul Kienzle
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
2 ##
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
3 ## This file is part of Octave.
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
4 ##
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
5 ## Octave is free software; you can redistribute it and/or modify it
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
6 ## under the terms of the GNU General Public License as published by
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
7 ## the Free Software Foundation; either version 3 of the License, or (at
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
8 ## your option) any later version.
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
9 ##
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
10 ## Octave is distributed in the hope that it will be useful, but
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
13 ## General Public License for more details.
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
14 ##
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
15 ## You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
16 ## along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6547
diff changeset
17 ## <http://www.gnu.org/licenses/>.
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
18
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
19 ## -*- texinfo -*-
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 10821
diff changeset
20 ## @deftypefn {Function File} {[@var{b}, @var{c}] =} spdiags (@var{A})
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 10821
diff changeset
21 ## @deftypefnx {Function File} {@var{b} =} spdiags (@var{A}, @var{c})
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 10821
diff changeset
22 ## @deftypefnx {Function File} {@var{b} =} spdiags (@var{v}, @var{c}, @var{A})
7144
30c0533e39ae [project @ 2007-11-09 17:37:24 by jwe]
jwe
parents: 7017
diff changeset
23 ## @deftypefnx {Function File} {@var{b} =} spdiags (@var{v}, @var{c}, @var{m}, @var{n})
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
24 ## A generalization of the function @code{diag}. Called with a single
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
25 ## input argument, the non-zero diagonals @var{c} of @var{A} are extracted.
11587
c792872f8942 all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents: 11523
diff changeset
26 ## With two arguments the diagonals to extract are given by the vector
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
27 ## @var{c}.
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
28 ##
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
29 ## The other two forms of @code{spdiags} modify the input matrix by
9051
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
30 ## replacing the diagonals. They use the columns of @var{v} to replace
1bf0ce0930be Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents: 8920
diff changeset
31 ## the columns represented by the vector @var{c}. If the sparse matrix
11471
994e2a93a8e2 Use uppercase 'A' to refer to matrix inputs in m-files.
Rik <octave@nomad.inbox5.com>
parents: 10821
diff changeset
32 ## @var{A} is defined then the diagonals of this matrix are replaced.
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
33 ## Otherwise a matrix of @var{m} by @var{n} is created with the
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
34 ## diagonals given by @var{v}.
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
35 ##
8325
b93ac0586e4b spelling corrections
Brian Gough<bjg@network-theory.co.uk>
parents: 7515
diff changeset
36 ## Negative values of @var{c} represent diagonals below the main
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
37 ## diagonal, and positive values of @var{c} diagonals above the main
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
38 ## diagonal.
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
39 ##
10821
693e22af08ae Grammarcheck documentation of m-files
Rik <octave@nomad.inbox5.com>
parents: 10793
diff changeset
40 ## For example:
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
41 ##
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
42 ## @example
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
43 ## @group
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
44 ## spdiags (reshape (1:12, 4, 3), [-1 0 1], 5, 4)
14327
4d917a6a858b doc: Use Octave coding conventions in @example blocks of docstrings.
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
45 ## @result{} 5 10 0 0
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
46 ## 1 6 11 0
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
47 ## 0 2 7 12
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
48 ## 0 0 3 8
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
49 ## 0 0 0 4
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
50 ## @end group
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
51 ## @end example
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
52 ##
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
53 ## @end deftypefn
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
54
5568
e9cde940b271 [project @ 2005-12-08 02:28:22 by jwe]
jwe
parents: 5307
diff changeset
55 function [A, c] = spdiags (v, c, m, n)
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
56
6498
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
57 if (nargin == 1 || nargin == 2)
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
58 ## extract nonzero diagonals of v into A,c
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7505
diff changeset
59 [nr, nc] = size (v);
7505
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions
David Bateman <dbateman@free.fr>
parents: 7144
diff changeset
60 [i, j, v] = find (v);
f5005d9510f4 Remove dispatched sparse functions and treat in the generic versions of the functions
David Bateman <dbateman@free.fr>
parents: 7144
diff changeset
61
6498
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
62 if (nargin == 1)
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
63 ## c contains the active diagonals
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
64 c = unique (j-i);
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
65 endif
6498
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
66 ## FIXME: we can do this without a loop if we are clever
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
67 offset = max (min (c, nc-nr), 0);
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
68 A = zeros (min (nr, nc), length (c));
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
69 for k = 1:length (c)
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
70 idx = find (j-i == c(k));
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
71 A(j(idx)-offset(k),k) = v(idx);
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
72 endfor
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
73 elseif (nargin == 3)
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
74 ## Replace specific diagonals c of m with v,c
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
75 [nr, nc] = size (m);
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
76 B = spdiags (m, c);
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
77 A = m - spdiags (B, c, nr, nc) + spdiags (v, c, nr, nc);
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
78 else
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
79 ## Create new matrix of size mxn using v,c
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
80 [j, i, v] = find (v);
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
81 offset = max (min (c(:), n-m), 0);
14666
1e77f6078692 spdiags.m: correctly treat row vector arguments
David Bateman <dbateman@free.fr>
parents: 14363
diff changeset
82 j = j(:) + offset(i(:));
1e77f6078692 spdiags.m: correctly treat row vector arguments
David Bateman <dbateman@free.fr>
parents: 14363
diff changeset
83 i = j - c(:)(i(:));
6498
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
84 idx = i > 0 & i <= m & j > 0 & j <= n;
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
85 A = sparse (i(idx), j(idx), v(idx), m, n);
2c85044aa63f [project @ 2007-04-05 17:59:47 by jwe]
jwe
parents: 5568
diff changeset
86 endif
5164
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
87
57077d0ddc8e [project @ 2005-02-25 19:55:24 by jwe]
jwe
parents:
diff changeset
88 endfunction
12558
1f0243949b77 spdiags.m: Treat empty vector (1x0 or 0x1) the same as diag().
Marco Caliari <marco.caliari@univr.it>
parents: 11587
diff changeset
89
1f0243949b77 spdiags.m: Treat empty vector (1x0 or 0x1) the same as diag().
Marco Caliari <marco.caliari@univr.it>
parents: 11587
diff changeset
90
14363
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14327
diff changeset
91 %!assert (spdiags (zeros (1,0),1,1,1), sparse (0))
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14327
diff changeset
92 %!assert (spdiags (zeros (0,1),1,1,1), sparse (0))
14666
1e77f6078692 spdiags.m: correctly treat row vector arguments
David Bateman <dbateman@free.fr>
parents: 14363
diff changeset
93 %!assert (spdiags ([0.5 -1 0.5], 0:2, 1, 1), sparse(0.5))