octave_packages/m/sparse/spdiags.m

   1 ## Copyright (C) 2000-2012 Paul Kienzle
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {[@var{b}, @var{c}] =} spdiags (@var{A})
  21 ## @deftypefnx {Function File} {@var{b} =} spdiags (@var{A}, @var{c})
  22 ## @deftypefnx {Function File} {@var{b} =} spdiags (@var{v}, @var{c}, @var{A})
  23 ## @deftypefnx {Function File} {@var{b} =} spdiags (@var{v}, @var{c}, @var{m}, @var{n})
  24 ## A generalization of the function @code{diag}.  Called with a single
  25 ## input argument, the non-zero diagonals @var{c} of @var{A} are extracted.
  26 ## With two arguments the diagonals to extract are given by the vector
  27 ## @var{c}.
  28 ##
  29 ## The other two forms of @code{spdiags} modify the input matrix by
  30 ## replacing the diagonals.  They use the columns of @var{v} to replace
  31 ## the columns represented by the vector @var{c}.  If the sparse matrix
  32 ## @var{A} is defined then the diagonals of this matrix are replaced.
  33 ## Otherwise a matrix of @var{m} by @var{n} is created with the
  34 ## diagonals given by @var{v}.
  35 ##
  36 ## Negative values of @var{c} represent diagonals below the main
  37 ## diagonal, and positive values of @var{c} diagonals above the main
  38 ## diagonal.
  39 ##
  40 ## For example:
  41 ##
  42 ## @example
  43 ## @group
  44 ## spdiags (reshape (1:12, 4, 3), [-1 0 1], 5, 4)
  45 ##    @result{} 5 10  0  0
  46 ##       1  6 11  0
  47 ##       0  2  7 12
  48 ##       0  0  3  8
  49 ##       0  0  0  4
  50 ## @end group
  51 ## @end example
  52 ##
  53 ## @end deftypefn
  54
  55 function [A, c] = spdiags (v, c, m, n)
  56
  57   if (nargin == 1 || nargin == 2)
  58     ## extract nonzero diagonals of v into A,c
  59     [nr, nc] = size (v);
  60     [i, j, v] = find (v);
  61
  62     if (nargin == 1)
  63       ## c contains the active diagonals
  64       c = unique (j-i);
  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     endfor
  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 = 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   endif
  87
  88 endfunction
  89
  90 %!test
  91 %assert(spdiags(zeros(1,0),1,1,1),0)
  92
  93 %!test
  94 %assert(spdiags(zeros(0,1),1,1,1),0)