-
Notifications
You must be signed in to change notification settings - Fork 29
/
opChol.m
84 lines (73 loc) · 2.78 KB
/
opChol.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
classdef opChol < opFactorization
%opCHOL Operator representing the Cholesky factorization of a
% symmetric and definite matrix with optional iterative
% refinement. Only the lower triangle of the input matrix
% is referenced.
%
% opChol(A) creates an operator for multiplication by the
% inverse of the matrix A implicitly represented by its Cholesky
% factorization. Optionally, iterative refinement is performed.
% Note that A is an explicit matrix.
%
% The following attributes may be changed by the user:
% * nitref : the maximum number of iterative refinement steps (3)
% * itref_tol : iterative refinement tolerance (1.0e-8)
% * force_itref : force iterative refinement (false)
%
% See also chol.
%
% Dominique Orban <dominique.orban@gerad.ca>, 2014.
%
% 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 )
L % Lower triangular Cholesky factor
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Methods - Public
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
methods
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% opChol. Constructor
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function op = opChol(A)
if nargin ~= 1
error('Invalid number of arguments.');
end
[m,n] = size(A);
% Construct operator
op = op@opFactorization('Chol', m, n);
B = A;
if ~issparse(A)
B = sparse(A);
end
op.A = opHermitian(B);
op.L = chol(B, 'lower');
op.Ainv = inv(op.L') * inv(op.L);
op.cflag = ~isreal(A);
end % function opChol
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% transpose
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function opOut = transpose(op)
opOut = inv(op.L.') * inv(op.L);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conj
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function opOut = conj(op)
opOut = inv(conj(op.L')) * inv(conj(op.L));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ctranpose
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function opOut = ctranspose(op)
opOut = inv(op.L') * inv(op.L);
end
end % methods - public
end % classdef