left_divide

PURPOSE ^

[V] = left_divide(E,I,tol,pp,V);

SYNOPSIS ^

function [V] = left_divide(E,I,tol,pp,V)

DESCRIPTION ^

[V] = left_divide(E,I,tol,pp,V);

 Implements left division and overcomes many small inefficiencies
 of matlab's left division. Also uses conjugate gradients (for large problems). 

E   = The full rank system matrix 
I   = The currents matrix (RHS) 
tol = The tolerance in the forward solution, e.g. 1e-5

 pp,V are old options from previous solver

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [V] = left_divide(E,I,tol,pp,V)
0002 %[V] = left_divide(E,I,tol,pp,V);
0003 %
0004 % Implements left division and overcomes many small inefficiencies
0005 % of matlab's left division. Also uses conjugate gradients (for large problems).
0006 %
0007 %E   = The full rank system matrix
0008 %I   = The currents matrix (RHS)
0009 %tol = The tolerance in the forward solution, e.g. 1e-5
0010 %
0011 % pp,V are old options from previous solver
0012 
0013 % (c) N. Polydorides 2003 % Copying permitted under terms of GNU GPL
0014 % $Id: left_divide.m 5169 2016-01-28 20:24:42Z alistair_boyle $
0015 
0016 if ~exist('tol'); tol = 1e-8; end
0017 
0018 [n_nodes,n_stims] = size(I);
0019 
0020 try
0021 % V= E\I;
0022 % This takes MUCH longer when you have  more vectors in I,
0023 %  even if they are repeated. There must be some way to simplify
0024 %  this to speed it up. Matlab's sparse operators really should
0025 %  do this for you.
0026 
0027 % TODO, Q&R should be cached somewhere
0028    [Q,R] = qr(I,0);
0029    rnotzeros = any(R~=0,2);
0030    Q= Q(:,rnotzeros);
0031    R= R(rnotzeros,:);
0032    V= (E \ Q)*R;
0033 
0034 % TODO: Iteratively refine
0035 %  From GH Scott: "once we have
0036 %   computed the approximate solution x, we perform one step
0037 %   of iterative refinement by computing the residual: r = Ax - b
0038 %   and then recalling the solve routine to solve
0039 %   Adx = r for the correction dx.
0040 % However, we don't want to repeat the '\', so we implement
0041 %   the underlying algorithm:
0042 %   If A is sparse, then MATLAB software uses CHOLMOD (after 7.2) to compute X.
0043 %    The computations result in  P'*A*P = R'*R
0044 %   where P is a permutation matrix generated by amd, and R is
0045 %   an upper triangular matrix. In this case, X = P*(R\(R'\(P'*B)))
0046 %
0047 % See also:
0048 % http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf
0049 % especially page 15 where it discusses the value of iterative refinement
0050 %  without extra precision bits.  ALso, we need to enable
0051 
0052 
0053 catch 
0054    [lasterr_str,lasterr_id]= lasterr;
0055    if ~strcmp(lasterr_id , 'MATLAB:nomem')
0056       error(lasterr_str); % rethrow error
0057    end
0058 
0059    eidors_msg('Memory exhausted for inverse. Trying PCG',2);
0060 
0061    if nargin < 5
0062       sz= [size(E,1),n_stims];
0063       V = eidors_obj('get-cache', sz, 'left_divide_V');
0064       if isempty(V); V= zeros(sz); end
0065    end
0066 
0067    ver = eidors_obj('interpreter_version'); % Matlab2013 renamed cholinc -> ichol
0068    if isreal(E)
0069       opts.droptol = tol*100;
0070       opt.type = 'ict';
0071       if ver.isoctave || ver.ver < 7.012
0072          U = cholinc(E, opt.droptol);
0073       else
0074          U = ichol(E, opt);
0075       end
0076       L = U';
0077       cgsolver = @pcg;
0078    else %Complex
0079       opts.droptol = tol/10;
0080       if ver.isoctave || ver.ver < 7.012 % Matlab2007 introduced ilu, luinc has now been dropped
0081          [L,U] = luinc(E, opt.droptol);
0082       else
0083          [L,U] = ilu(E, opt);
0084       end
0085       cgsolver = @bicgstab;
0086    end
0087 
0088    for i=1:n_stims
0089       [V(:,i),flag] = feval( cgsolver, E,I(:,i), ...
0090                tol*norm(I(:,i)),n_nodes,L,U,V(:,i));
0091    end 
0092       eidors_obj('set-cache', sz, 'left_divide_V', V);
0093 end
0094 
0095 
0096 %%% OLD CODE
0097    % Cholesky solver. Gives poor results matching others
0098    % so we no longer use it
0099    if 0 
0100        %Permute the rows and columns to make the factors sparser
0101        E = E(pp,pp);
0102        In = I(pp,:);
0103        rr(pp)=1:max(size(pp));  % this should be done only Once!
0104                                 % actually much better just to do the
0105                                 % renumbering when the mesh is generated!
0106        U = chol(E);
0107        q_c =  U' \ In;  
0108        Vn = U \ q_c;
0109        %De-permute the result for Cholesky
0110        V = Vn(rr,:);
0111    end

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005