Easily-resolved efficiency issue in OPdiffusion

[rileyjmurray]

The OPdiffusion function performs a sequence of matrix-vector products with a linear operator $R^{-1} S$ (or its adjoint), where $R$ and $S$ are sparse matrices. The current implementation, here, use Matlab's matrix-division syntax with the for-loop for computing the matrix-vector products. It's much more efficient to factor $R$ before the for-loops and then use the factors there.

Here's modified code that works when $R$ is positive definite. If there are valid situations where $R$ is indefinite (I'm not sure if there are ...) then this can be modified to use LU decomposition instead of Cholesky.

function w = OPdiffusion(u,R,S,Tsteps,tflag)
n = size(R,1);
if strcmpi(tflag,'size')
    w(1) = n;
    w(2) = n;
    return;
end
[r,f,p] = chol(R);
assert(f == 0);
% r' r = p' R p
% p r'r p' = R
% inv(R) = p inv(r'r) p'
%        = p inv(r) inv(r') p'
if strcmpi(tflag,'notransp')
    w = u;
    for i = 1:Tsteps
        % w = R\(S*w);
        v = S*w;
        w = p * (r\((r')\(p' * v)));
    end
else
    w = u;
    for i = 1:Tsteps
        % w = S'*((R')\w);
        v = p * (r\((r')\(p' * w)));
        w = S' * v;
    end
end