Skip to content

Commit

Permalink
Changed back to use vector in stead of sparse matrix to save memory
Browse files Browse the repository at this point in the history
  • Loading branch information
tristan van Leeuwen committed Oct 25, 2015
1 parent 7442d97 commit 36e69f6
Showing 1 changed file with 113 additions and 103 deletions.
216 changes: 113 additions & 103 deletions opDiag.m
Original file line number Diff line number Diff line change
@@ -1,105 +1,115 @@
classdef opDiag < opSpot
%OPDIAG Diagonal operator.
%
% opDiag(D) creates an operator for multiplication by the
% diagonal matrix that has a vector D on its diagonal.
%
% See also diag.

% Copyright 2009, Ewout van den Berg and Michael P. Friedlander
% See the file COPYING.txt for full copyright information.
% Use the command 'spot.gpl' to locate this file.

% http://www.cs.ubc.ca/labs/scl/spot

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Properties
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
properties( SetAccess = private )
diag % Diagonal entries
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Methods - Public
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
methods

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% opDiag. Constructor
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function op = opDiag(d)
if nargin ~= 1
error('Invalid number of arguments.');
end

% Vectorize d and get size
d = d(:);
n = length(d);

% Construct operator
op = op@opSpot('Diag',n,n);
op.cflag = ~isreal(d);
op.diag = spdiags(d,0,n,n);
op.sweepflag = true;
end % function opDiag

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% double
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A = double(op)
A = full(op.diag);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% transpose
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function opOut = transpose(op)
opOut = op;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conj
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function opOut = conj(op)
opOut = opDiag(conj(spdiags(op.diag)));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ctranpose
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function opOut = ctranspose(op)
opOut = opDiag(conj(spdiags(op.diag)));
end

end % methods - public

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Methods - protected
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
methods( Access = protected )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% multiply
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = multiply(op,x,mode)
if mode == 1
y = op.diag*x;
else
y = op.diag'*x;
end
end % function multiply

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% divide
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = divide(op,b,mode)
if mode == 1
x = op.diag\b;
else
x = op.diag'\b;
end
end % function divide

end % methods - protected

%OPDIAG Diagonal operator.
%
% opDiag(D) creates an operator for multiplication by the
% diagonal matrix that has a vector D on its diagonal.
%
% See also diag.

% Copyright 2009, Ewout van den Berg and Michael P. Friedlander
% See the file COPYING.txt for full copyright information.
% Use the command 'spot.gpl' to locate this file.

% http://www.cs.ubc.ca/labs/scl/spot

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Properties
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
properties( SetAccess = private )
diag % Diagonal entries
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Methods - Public
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
methods

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% opDiag. Constructor
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function op = opDiag(d)
if nargin ~= 1
error('Invalid number of arguments.');
end

% Vectorize d and get size
d = d(:);
n = length(d);

% Construct operator
op = op@opSpot('Diag',n,n);
op.cflag = ~isreal(d);
op.diag = d(:);
op.sweepflag = true;
end % function opDiag

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% double
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function A = double(op)
A = diag(op.diag);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% transpose
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function opOut = transpose(op)
opOut = op;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conj
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function opOut = conj(op)
opOut = opDiag(conj(op.diag));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ctranpose
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function opOut = ctranspose(op)
opOut = opDiag(conj(op.diag));
end

end % methods - public

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Methods - protected
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
methods( Access = protected )

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% multiply
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = multiply(op,x,mode)
y = zeros(op.n,size(x,2));
if mode == 1
for k = 1:size(x,2) % for sweepflag, can also be achieved with bsxfun, but that seems slower
y(:,k) = op.diag.*x(:,k);
end
else
for k = 1:size(x,2) % for sweepflag, can also be achieved with bsxfun, but that seems slower
y(:,k) = conj(op.diag).*x(:,k);
end
end
end % function multiply

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% divide
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = divide(op,b,mode)
x = zeros(op.n,size(b,2));
if mode == 1
for k = 1:size(x,2) % for sweepflag, can also be achieved with bsxfun, but that seems slower
x(:,k) = b(:,k)./op.diag;
end
else
for k = 1:size(x,2) % for sweepflag, can also be achieved with bsxfun, but that seems slower
x = b(:,k)./conj(op.diag);
end
end
end % function divide

end % methods - protected

end % classdef

0 comments on commit 36e69f6

Please sign in to comment.