forward_solver

[V] = forward_solver(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  = UNUSED
V   = The approximated nodal potential distribution (USED FOR PCG SOLN)

0001 function [V] = forward_solver(E,I,tol,pp,V, compat_param);
0002 %[V] = forward_solver(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 %pp  = UNUSED
0011 %V   = The approximated nodal potential distribution (USED FOR PCG SOLN)
0012 
0013 % (c) N. Polydorides 2003 % Copying permitted under terms of GNU GPL
0014 % $Id: forward_solver.m 3082 2012-06-07 17:56:34Z aadler $
0015 
0016 % The EIDORS3D V2 was [V] = forward_solver(vtx,E,I,tol,pp,V);
0017 % but vtx not used, so it's forward_solver(E,I,tol,pp,V, compat_param);
0018 warning('EIDORS:deprecated','FORWARD_SOLVER is deprecated as of 07-Jun-2012. Use LEFT_DIVIDE instead.');
0019 
0020 if size(E,2) == 3; % E is actually vtx
0021   if nargin>=2; E= I;             end
0022   if nargin>=3; I= tol;           end
0023   if nargin>=4; tol= pp;          end
0024   if nargin>=5; pp= V;            end
0025   if nargin>=6; V =compat_param;  end
0026 end 
0027 if ~exist('tol'); tol = 1e-8; end
0028 
0029 [n_nodes,n_stims] = size(I);
0030 
0031 try
0032 % V= E\I;
0033 % This takes MUCH longer when you have  more vectors in I,
0034 %  even if they are repeated. There must be some way to simplify
0035 %  this to speed it up. Matlab's sparse operators really should
0036 %  do this for you.
0037 
0038 % TODO, Q&R should be cached somewhere
0039    [Q,R] = qr(I,0);
0040    rnotzeros = any(R~=0,2);
0041    Q= Q(:,rnotzeros);
0042    R= R(rnotzeros,:);
0043    V= (E \ Q)*R;
0044 
0045 % TODO: Iteratively refine
0046 %  From GH Scott: "once we have
0047 %   computed the approximate solution x, we perform one step
0048 %   of iterative refinement by computing the residual: r = Ax - b
0049 %   and then recalling the solve routine to solve
0050 %   Adx = r for the correction dx.
0051 % However, we don't want to repeat the '\', so we implement
0052 %   the underlying algorithm:
0053 %   If A is sparse, then MATLAB software uses CHOLMOD (after 7.2) to compute X.
0054 %    The computations result in  P'*A*P = R'*R
0055 %   where P is a permutation matrix generated by amd, and R is
0056 %   an upper triangular matrix. In this case, X = P*(R\(R'\(P'*B)))
0057 %
0058 % See also:
0059 % http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf
0060 % especially page 15 where it discusses the value of iterative refinement
0061 %  without extra precision bits.  ALso, we need to enable
0062 
0063 
0064 catch 
0065    [lasterr_str,lasterr_id]= lasterr;
0066    if ~strcmp(lasterr_id , 'MATLAB:nomem')
0067       error(lasterr_str); % rethrow error
0068    end
0069 
0070    eidors_msg('Memory exhausted for inverse. Trying PCG',2);
0071 
0072    if nargin < 5
0073       sz= [size(E,1),n_stims];
0074       V = eidors_obj('get-cache', sz, 'forward_solver_V');
0075       if isempty(V); V= zeros(sz); end
0076    end
0077 
0078    if isreal(E)
0079       U = cholinc(E,tol*100); L = U'; 
0080       cgsolver = @pcg;
0081    else %Complex
0082       [L,U] = luinc(E,tol/10);
0083       cgsolver = @bicgstab;
0084    end
0085 
0086    for i=1:n_stims
0087       [V(:,i),flag] = feval( cgsolver, E,I(:,i), ...
0088                tol*norm(I(:,i)),n_nodes,L,U,V(:,i));
0089    end 
0090       eidors_obj('set-cache', sz, 'forward_solver_V', V);
0091 end
0092 
0093 
0094 %%% OLD CODE
0095    % Cholesky solver. Gives poor results matching others
0096    % so we no longer use it
0097    if 0 
0098        %Permute the rows and columns to make the factors sparser
0099        E = E(pp,pp);
0100        In = I(pp,:);
0101        rr(pp)=1:max(size(pp));  % this should be done only Once!
0102                                 % actually much better just to do the
0103                                 % renumbering when the mesh is generated!
0104        U = chol(E);
0105        q_c =  U' \ In;  
0106        Vn = U \ q_c;
0107        %De-permute the result for Cholesky
0108        V = Vn(rr,:);
0109    end