forward_solver

PURPOSE ^

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

SYNOPSIS ^

function [V] = forward_solver(E,I,tol,pp,V, compat_param);

DESCRIPTION ^

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

This function solves the forward problem using matlab's \ sovler, or
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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [V] = forward_solver(E,I,tol,pp,V, compat_param);
0002 %[V] = forward_solver(E,I,tol,pp,V);
0003 %
0004 %This function solves the forward problem using matlab's \ sovler, or
0005 %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.html 2819 2011-09-07 16:43:11Z 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 if size(E,2) == 3; % E is actually vtx
0019   if nargin>=2; E= I;             end
0020   if nargin>=3; I= tol;           end
0021   if nargin>=4; tol= pp;          end
0022   if nargin>=5; pp= V;            end
0023   if nargin>=6; V =compat_param;  end
0024 end 
0025 
0026 [n_nodes,n_stims] = size(I);
0027 
0028 try
0029 % V= E\I;
0030 % This takes MUCH longer when you have  more vectors in I,
0031 %  even if they are repeated. There must be some way to simplify
0032 %  this to speed it up. Matlab's sparse operators really should
0033 %  do this for you.
0034 
0035 % TODO, Q&R should be cached somewhere
0036    [Q,R] = qr(I,0);
0037    rnotzeros = any(R~=0,2);
0038    Q= Q(:,rnotzeros);
0039    R= R(rnotzeros,:);
0040    V= (E \ Q)*R;
0041 
0042 % TODO: Iteratively refine
0043 %  From GH Scott: "once we have
0044 %   computed the approximate solution x, we perform one step
0045 %   of iterative refinement by computing the residual: r = Ax - b
0046 %   and then recalling the solve routine to solve
0047 %   Adx = r for the correction dx.
0048 % However, we don't want to repeat the '\', so we implement
0049 %   the underlying algorithm:
0050 %   If A is sparse, then MATLAB software uses CHOLMOD (after 7.2) to compute X.
0051 %    The computations result in  P'*A*P = R'*R
0052 %   where P is a permutation matrix generated by amd, and R is
0053 %   an upper triangular matrix. In this case, X = P*(R\(R'\(P'*B)))
0054 
0055 catch 
0056    [lasterr_str,lasterr_id]= lasterr;
0057    if ~strcmp(lasterr_id , 'MATLAB:nomem')
0058       error(lasterr_str); % rethrow error
0059    end
0060 
0061    eidors_msg('Memory exhausted for inverse. Trying PCG',2);
0062 
0063    if nargin < 5
0064       sz= [size(E,1),n_stims];
0065       V = eidors_obj('get-cache', sz, 'forward_solver_V');
0066       if isempty(V); V= zeros(sz); end
0067    end
0068 
0069    if isreal(E)
0070       U = cholinc(E,tol*100); L = U'; 
0071       cgsolver = @pcg;
0072    else %Complex
0073       [L,U] = luinc(E,tol/10);
0074       cgsolver = @bicgstab;
0075    end
0076 
0077    for i=1:n_stims
0078       [V(:,i),flag] = feval( cgsolver, E,I(:,i), ...
0079                tol*norm(I(:,i)),n_nodes,L,U,V(:,i));
0080    end 
0081       eidors_obj('set-cache', sz, 'forward_solver_V', V);
0082 end
0083 
0084 
0085 %%% OLD CODE
0086    % Cholesky solver. Gives poor results matching others
0087    % so we no longer use it
0088    if 0 
0089        %Permute the rows and columns to make the factors sparser
0090        E = E(pp,pp);
0091        In = I(pp,:);
0092        rr(pp)=1:max(size(pp));  % this should be done only Once!
0093                                 % actually much better just to do the
0094                                 % renumbering when the mesh is generated!
0095        U = chol(E);
0096        q_c =  U' \ In;  
0097        Vn = U \ q_c;
0098        %De-permute the result for Cholesky
0099        V = Vn(rr,:);
0100    end

Generated on Tue 09-Aug-2011 11:38:31 by m2html © 2005