octave_packages/linear-algebra-2.2.0/@blksparse/subsref.m

   1 ## Copyright (C) 2010 VZLU Prague
   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 3 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 Octave; see the file COPYING.  If not, see
  15 ## <http://www.gnu.org/licenses/>.
  16
  17 function ss = subsref (s, subs)
  18   if (length (subs) != 1)
  19     error ("blksparse: invalid index chain");
  20   endif
  21   if (strcmp (subs(1).type, "()"))
  22     ind = subs(1).subs;
  23     if (length (ind) == 2)
  24       idx = make_block_index (ind{1}, s.bsiz(1));
  25       jdx = make_block_index (ind{2}, s.bsiz(2));
  26       ## Use sparse indexing to solve it all.
  27       sn = sparse (s.i, s.j, 1:size (s.sv, 3), s.siz(1), s.siz (2));
  28       sn = sn(idx, jdx);
  29       [i, j, k] = find (sn);
  30       ss = s;
  31       ss.i = i;
  32       ss.j = j;
  33       ss.sv = s.sv(:,:,k);
  34       ss.siz = size (sn);
  35     else
  36       error ("blksparse: linear indexing is not supported");
  37     endif
  38   else
  39     error ("blksparse: only supports () indexing");
  40   endif
  41
  42 endfunction
  43
  44 function bi = make_block_index (i, bs)
  45   if (strcmp (i, ':'))
  46     bi = i;
  47   else
  48     if (rem (numel (i), bs) == 0)
  49       ba = reshape (i, bs, []);
  50       bi = ba(1,:);
  51       if (any (rem (bi, bs) != 1) || any ((ba != bsxfun (@plus, bi, [0:bs-1].'))(:)))
  52         error ("blksparse: index must preserve block structure");
  53       else
  54         bi = ceil (bi / bs);
  55       endif
  56     else
  57       error ("blksparse: index must preserve block structure");
  58     endif
  59   endif
  60 endfunction