Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 97f6a6f0 authored by hhakim's avatar hhakim
Browse files

Fix matfaust.lazylinop.LazyLinearOp.mtimes in case the operands are reversed... 

Fix matfaust.lazylinop.LazyLinearOp.mtimes in case the operands are reversed like in M * L (where M a matrix and L a LazyLinearOp).
parent dbdb334e
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
wrapper/matlab/+matfaust/+lazylinop/@LazyLinearOp/LazyLinearOp.m
+ 50
− 45
View file @ 97f6a6f0
...@@ -117,54 +117,59 @@ classdef LazyLinearOp < handle % needed to use references on objects ...@@ -117,54 +117,59 @@ classdef LazyLinearOp < handle % needed to use references on objects
end end
function LM = mtimes(L, A) function LM = mtimes(L, A)
import matfaust.lazylinop.LazyLinearOp import matfaust.lazylinop.LazyLinearOp
if isscalar(L) && LazyLinearOp.isLazyLinearOp(A) if LazyLinearOp.isLazyLinearOp(A)
LM = mtimes(A, L); if isscalar(L)
return; LM = mtimes(A, L);
end return;
check_meth(L, 'mtimes'); elseif ~ LazyLinearOp.isLazyLinearOp(L)
op_is_scalar = all(size(A) == [1, 1]); LM = ctranspose(mtimes(A.', L.'));
if ~ op_is_scalar && ~ all(size(L, 2) == size(A, 1)) return;
error('Dimensions must agree') end
end end
if op_is_scalar check_meth(L, 'mtimes');
new_size = size(L); op_is_scalar = all(size(A) == [1, 1]);
else if ~ op_is_scalar && ~ all(size(L, 2) == size(A, 1))
new_size = [size(L, 1), size(A, 2)]; error('Dimensions must agree')
end end
if op_is_scalar
new_size = size(L);
else
new_size = [size(L, 1), size(A, 2)];
end
function l = mul_index_lambda(L, A, S) function l = mul_index_lambda(L, A, S)
% L and A must be LazyLinearOp % L and A must be LazyLinearOp
import matfaust.lazylinop.LazyLinearOp import matfaust.lazylinop.LazyLinearOp
Sr.type = '()'; Sr.type = '()';
Sr.subs = {S.subs{1}, ':'}; Sr.subs = {S.subs{1}, ':'};
Sc.type = '()'; Sc.type = '()';
Sc.subs = {':', S.subs{2}}; Sc.subs = {':', S.subs{2}};
L.lambdas{L.I}(Sr) * A.lambdas{L.I}(Sc); L.lambdas{L.I}(Sr) * A.lambdas{L.I}(Sc);
end end
if ~ LazyLinearOp.isLazyLinearOp(A) && ismatrix(A) && isnumeric(A) && any(size(A) ~= [1, 1]) if ~ LazyLinearOp.isLazyLinearOp(A) && ismatrix(A) && isnumeric(A) && any(size(A) ~= [1, 1])
% A is a dense matrix that is not limited to one element % A is a dense matrix that is not limited to one element
LM = L.lambdas{L.MUL}(A); LM = L.lambdas{L.MUL}(A);
else else
if isscalar(A) if isscalar(A)
LM_size = size(L); LM_size = size(L);
else else
if ~ LazyLinearOp.isLazyLinearOp(A) if ~ LazyLinearOp.isLazyLinearOp(A)
A = LazyLinearOp.create_from_op(A); A = LazyLinearOp.create_from_op(A);
end end
LM_size = [size(L, 1), size(A, 2)] LM_size = [size(L, 1), size(A, 2)]
end end
lambdas = {@(o) L * (A * o), ... %MUL lambdas = {@(o) L * (A * o), ... %MUL
@() A.' * L.', ... % T @() A.' * L.', ... % T
@() A' * L', ... % H @() A' * L', ... % H
@(S) mul_index_lambda(L, A, S)% I @(S) mul_index_lambda(L, A, S)% I
}; };
LM = LazyLinearOp(lambdas, LM_size); LM = LazyLinearOp(lambdas, LM_size);
end end
end end
function LP = plus(L, op) function LP = plus(L, op)
......
wrapper/matlab/+matfaust/+lazylinop/@LazyLinearOpKron/LazyLinearOpKron.m deleted 100644 → 0
+ 0
− 103
View file @ dbdb334e
% ================================
%> @brief This class implements a specialization of LazyLinearOp dedicated to the
%> Kronecker product of two linear operators.
%>
%> @b See @b also: matfaust.lazylinop.kron.
% ================================
classdef LazyLinearOpKron < matfaust.lazylinop.LazyLinearOp
properties (SetAccess = private, Hidden = true)
A;
B;
end
methods
%================================
%> Constructor for the A x B Kronecker product.
%================================
function LK = LazyLinearOpKron(A, B)
shape = [size(A, 1) * size(B, 1), size(A, 2) * size(B, 2)];
LK = LK@matfaust.lazylinop.LazyLinearOp(@() A, shape, A);
LK.A = A;
LK.B = B;
end
%================================
%> Returns the LazyLinearOpKron conjugate.
%================================
function CLK = conj(LK)
import matfaust.lazylinop.*
CLK = LazyLinearOpKron(conj(asLazyLinearOp(LK.A)), conj(asLazyLinearOp(LK.B)));
end
%================================
%> Returns the LazyLinearOpKron transpose.
%================================
function CLK = transpose(LK)
import matfaust.lazylinop.*
CLK = LazyLinearOpKron(transpose(asLazyLinearOp(LK.A)), transpose(asLazyLinearOp((LK.B))));
end
%================================
%> Returns the LazyLinearOpKron adjoint/transconjugate.
%================================
function CTLK = ctranspose(LK)
import matfaust.lazylinop.*
CTLK = LazyLinearOpKron(ctranspose(asLazyLinearOp(LK.A)), ctranspose(asLazyLinearOp((LK.B))));
end
%=============================================================
%> @brief Returns the LazyLinearOp for the multiplication self * op
%> or if op is a full matrix it returns the full matrix (self * op).
%>
%> @note this specialization is particularly optimized for multiplying the
%> operator by a vector.
%>
%> @param op: an object compatible with self for this binary operation.
%=============================================================
function LmK = mtimes(LK, op)
import matfaust.lazylinop.LazyLinearOp
if isscalar(LK) && LazyLinearOp.isLazyLinearOp(op)
LmK = mtimes(op, LK);
return;
end
check_meth(LK, 'mtimes');
op_is_scalar = all(size(op) == [1, 1]);
if ~ op_is_scalar && ~ all(size(LK, 2) == size(op, 1))
error('Dimensions must agree')
end
if op_is_scalar
new_size = LK.shape;
else
new_size = [size(LK, 1), size(op, 2)];
end
if ~ LazyLinearOp.isLazyLinearOp(op) && ismatrix(op) && isnumeric(op) && any(size(op) ~= [1, 1]) && any(ismember(methods(op), 'reshape')) && any(ismember(methods(op), 'mtimes')) %&& any(ismember(methods(op), 'subsref'))
% op is a dense matrix that is not limited to one element
LmK = zeros(new_size);
A = LK.A;
B = LK.B;
for j=1:new_size(2)
op_mat = reshape(op(:, j), [size(B, 2), size(A, 2)]);
LmK(:, j) = reshape(LazyLinearOp.eval_if_lazy(B) * op_mat * transpose(LazyLinearOp.eval_if_lazy(A)), [new_size(1), 1]);
end
else
LmK = LazyLinearOp(@() LK.eval() * LazyLinearOp.eval_if_lazy(op), new_size, LK.root_obj);
end
end
%================================
%> See LazyLinearOp.eval.
%================================
function E = eval(LK)
A = LK.A;
B = LK.B;
if any(ismember(methods(A), 'full'))
A = full(A);
end
if any(ismember(methods(B), 'full'))
B = full(B);
end
E = kron(A, B);
end
end
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment