spot/opFactorization.m at master · mpf/spot · GitHub

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
85
86
87
88
89
90
91
92
classdef opFactorization < opSpot
%OPFACTORIZATION  Operator representing an inverse by way of a factorization.
%                 Useful to avoid multiple factorizations as in opInverse
%                 and save time by performing only forward and backsolves.

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Properties
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  properties( SetAccess = protected )
    A             % Input matrix as operator
    Ainv          % Inverse as a factorization operator
    rNorm         % Residual norm (if iterative refinement is performed)
  end

  properties( SetAccess = public )
    nitref = 3     % Default max number of iterative refinement steps
    itref_tol = 1.0e-8       % Default iterative refinement tolerance
    force_itref = false  % Force nitref steps of iterative refinement
  end

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

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % opLDL. Constructor
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function op = opFactorization(name, m, n)
      % Construct operator
      op = op@opSpot(name, m, n);
      op.sweepflag    = true;
    end % function opFactorization

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Setters
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function op = set.nitref(op, val)
      op.nitref = max(0, round(val));
    end

    function op = sef.itref_tol(op, val)
      op.itref_tol = max(0, val);
    end

    function op = set.force_itref(op, val)
      if val ~= false && val ~= true
        op.force_itref = false;
      else
        op.force_itref = val;
      end
    end

  end % methods - public

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

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % multiply
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function y = multiply(op, x, ~)
       y = op.Ainv * x;
       % Perform iterative refinement if necessary / requested
       if op.nitref > 0
          r = x - op.A * y;
          rNorm = norm(r, 'inf'); %#ok<*PROP>
          xNorm = norm(x, 'inf');
          nit = 0;
          while nit < op.nitref && (rNorm >= op.itref_tol * xNorm || op.force_itref)
             dy = op.Ainv * r;
             y = y + dy;
             r = x - op.A * y;
             rNorm = norm(r, 'inf');
             nit = nit + 1;
          end
          op.rNorm = rNorm;
       end
    end % function multiply

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % divide
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function x = divide(op, b, ~)
       x = op.A * b;
    end % function divide

  end % methods - protected

end % classdef