[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)
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