1 ## Copyright (C) 2010 VZLU Prague
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.
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.
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/>.
18 ## @deftypefn {Function File} mrdivide (@var{x}, @var{y})
19 ## Performs a left division with a block sparse matrix.
20 ## If @var{y} is a block sparse matrix, it must be either diagonal
21 ## or triangular, and @var{x} must be full.
22 ## If @var{y} is built-in sparse or full, @var{x} is converted
23 ## accordingly, then the built-in division is used.
26 function c = mrdivide (a, b)
27 if (isa (b, "blksparse"))
29 error ("blksparse: sparse / block sparse not implemented");
31 c = mrdivide_ms (a, b);
40 function y = mrdivide_ms (x, s)
43 if (bsiz(1) != bsiz(2) || siz(1) != siz(2))
44 error ("blksparse: can only divide by square matrices with square blocks");
49 if (xc != siz(2)*bsiz(2))
50 gripe_nonconformant (siz.*bsiz, [xr, xc]);
53 if (isempty (s) || isempty (x))
59 x = reshape (x, xr, bsiz(2), siz(2));
69 full_diag = length (d) == n;
71 isdiag = full_diag && ns == n; # block diagonal
72 islower = full_diag && all (si >= sj); # block upper triangular
73 isupper = full_diag && all (si <= sj); # block lower triangular
76 xx = num2cell (x, [1, 2]);
77 ss = num2cell (sv, [1, 2]);
78 yy = cellfun (@mldivide, ss, xx, "uniformoutput", false);
83 ## this is the dot version
84 y(:,:,1) = x(:,:,1) / sv(:,:,1);
87 xy = blkmm (y(:,:,si(k)), sv(:,:,k));
88 y(:,:,j) = (x(:,:,j) - sum (xy, 3)) / sv(:,:,d(j));
92 ## this is the dot version
93 y(:,:,n) = x(:,:,n) / sv(:,:,ns);
96 xy = blkmm (y(:,:,si(k)), sv(:,:,k));
97 y(:,:,j) = (x(:,:,j) - sum (xy, 3)) / sv(:,:,d(j));
100 error ("blksparse: mldivide: matrix must be block triangular or diagonal");
104 y = reshape (y, xr, bsiz(2)*siz(2));
107 function gripe_nonconformant (s1, s2, what = "arguments")
108 error ("Octave:nonconformant-args", ...
109 "nonconformant %s (op1 is %dx%d, op2 is %dx%d)", what, s1, s2);