octave_packages/linear-algebra-2.2.0/@blksparse/mtimes.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 ## -*- texinfo -*-
  18 ## @deftypefn {Function File} mtimes (@var{x}, @var{y})
  19 ## Multiplies a block sparse matrix with a full matrix, or two block sparse
  20 ## matrices. Multiplication of block sparse * sparse is not implemented.
  21 ## If one of arguments is a scalar, it's a scalar multiply.
  22 ## @end deftypefn
  23
  24 function c = mtimes (a, b)
  25   if (isa (a, "blksparse"))
  26     if (isa (b, "blksparse"))
  27       c = mtimes_ss (a, b);
  28     else
  29       c = mtimes_sm (a, b);
  30     endif
  31   elseif (isa (b, "blksparse"))
  32     c = mtimes_ms (a, b);
  33   else
  34     error ("blksparse: invalid arguments to mtimes");
  35   endif
  36 endfunction
  37
  38 function y = mtimes_sm (s, x)
  39   if (isscalar (x))
  40     y = s;
  41     y.sv *= x;
  42     return;
  43   elseif (issparse (x))
  44     error ("blksparse * sparse not implemented.");
  45   endif
  46   siz = s.siz;
  47   bsiz = s.bsiz;
  48   ## Check sizes.
  49   [xr, xc] = size (x);
  50   if (xr != siz(2)*bsiz(2))
  51     gripe_nonconformant (siz.*bsiz, [xr, xc]);
  52   endif
  53   ## Form blocks.
  54   x = reshape (x, bsiz(2), siz(2), xc);
  55   x = permute (x, [1, 3, 2]);
  56   ## Scatter.
  57   xs = x(:,:,s.j);
  58   ## Multiply.
  59   ys = blkmm (s.sv, xs);
  60   ## Gather.
  61   y = accumdim (s.i, ys, 3, siz(1));
  62   y = permute (y, [1, 3, 2]);
  63   ## Narrow blocks.
  64   y = reshape (y, bsiz(1)*siz(1), xc);
  65 endfunction
  66
  67 function y = mtimes_ms (x, s)
  68   if (isscalar (x))
  69     y = s;
  70     y.sv *= x;
  71     return;
  72   elseif (issparse (x))
  73     error ("sparse * blksparse not implemented.");
  74   endif
  75   siz = s.siz;
  76   bsiz = s.bsiz;
  77   ## Check sizes.
  78   [xr, xc] = size (x);
  79   if (xc != siz(1)*bsiz(1))
  80     gripe_nonconformant ([xr, xc], siz.*bsiz);
  81   endif
  82   ## Form blocks.
  83   x = reshape (x, xr, bsiz(2), siz(2));
  84   ## Scatter.
  85   xs = x(:,:,s.i);
  86   ## Multiply.
  87   ys = blkmm (xs, s.sv);
  88   ## Gather.
  89   y = accumdim (s.j, ys, 3, siz(2));
  90   ## Narrow blocks.
  91   y = reshape (y, xr, bsiz(2)*siz(2));
  92 endfunction
  93
  94 function s = mtimes_ss (s1, s2)
  95
  96   ## Conformance check.
  97   siz1 = s1.siz;
  98   bsiz1 = s1.bsiz;
  99   siz2 = s2.siz;
 100   bsiz2 = s2.bsiz;
 101   if (bsiz1(2) != bsiz2(1))
 102     gripe_nonconformant (bsiz1, bsiz2, "block sizes");
 103   elseif (siz1(2) != siz2(1))
 104     gripe_nonconformant (bsiz1.*siz1, bsiz2.*siz2);
 105   endif
 106
 107   ## Hardcore hacks, man!
 108   ss = sparse (s1.i, s1.j, 1:length (s1.i), "unique");
 109   ss = ss(:,s2.i);
 110   [i, j, k] = find (ss);
 111   sv = blkmm (s1.sv(:,:,k), s2.sv(:,:,j));
 112   j = s2.j(j);
 113
 114   s = blksparse (i, j, sv, siz1(1), siz2(2));
 115
 116 endfunction
 117
 118 function gripe_nonconformant (s1, s2, what = "arguments")
 119   error ("Octave:nonconformant-args", ...
 120   "nonconformant %s (op1 is %dx%d, op2 is %dx%d)", what, s1, s2);
 121 endfunction