forward_solver.sci

来自「用来实现三维阻抗及光学断层成像重建的matlab程序」· SCI 代码 · 共 98 行

function [V]=forward_solver(vtx,E,I,tol,pp,V)
[nargout,nargin] = argn(0)
//function [V] = forward_solver(vtx,E,I,tol,pp,V);
// 
//This function solves the forward problem using the Cholesky or LU method or 
//conjugate gradients. 
// 
// 
// 
//vtx = The vertices
//E   = The full rank system matrix 
//I   = The currents matrix (RHS) 
//pp  = the column permutation vector 
//V   = The approximated nodal potential distribution
//tol = The tolerance in the forward solution, e.g. 1e-5
 
 
 
 
// d: number of current patterns
[bla,d] = size(I);
 
if nargin<6 then
  V = zeros(size(E,1),d);
end
 
 
if isreal(E,0)==%t then
   
  if pp~=1:size(E,1) then
    //There is a colum permutation, Cholesky opted
    //Permute the columns to make the factors sparser
    E = E(pp,pp);
    In = I(pp,:);
    rr(1,pp) = 1:max(size(pp));
     
    //!! Unknown function cholinc ,the original calling sequence is used
    U = cholinc(E,tol);
    q_c = U'\In;
    Vn = U\q_c;
    //De-permute the result for Cholesky
    V = Vn(rr,:);
  else
     
    //Alternatively use pcg ********
     
    //!! Unknown function cholinc ,the original calling sequence is used
    K = cholinc(E,tol*100);
    for i = 1:d
       
      //! mtlb(relres) can be replaced by relres() or relres whether relres is an m-file or not
       
      //! mtlb(iter) can be replaced by iter() or iter whether iter is an m-file or not
       
      //! mtlb(resvec) can be replaced by resvec() or resvec whether resvec is an m-file or not
       
      //!! Unknown function pcg ,the original calling sequence is used
      [V(:,i),ones(8,1).*.[1 0 0;1 1 1;0 0 1;0 0 0],mtlb(relres),mtlb(iter),mtlb(resvec)]==pcg(E,I(:,i),0.000001*norm(I(:,i)),d^3,K',K,V(:,i));
    end
     
  end
   
end
//is real
 
if isreal(E,0)==%f then
   
  E = E(pp,pp);
  In = I(pp,:);
   
  //! rr = mtlb_i(rr,pp,1:max(size(pp))) may be replaced by rr(pp) = 1:max(size(pp))
  //!   if rr is not a  vector
  //!   or if rr and 1:max(size(pp)) are both row or column vectors
  rr = mtlb_i(rr,pp,1:max(size(pp)))
   
  //!! Unknown function luinc ,the original calling sequence is used
  [L,U] = luinc(E,tol);
  //Level of approximation, can be adjusted.
  q_c = L\In;
  Vn = U\q_c;
  //De-permute the result for Cholesky
  V = Vn(rr,:);
   
end
//is complex
 
 
 
 
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// This is part of the EIDORS suite.
// Copyright (c) N. Polydorides 2001
// Copying permitted under terms of GNU GPL
// See enclosed file gpl.html for details.
// EIDORS 3D version 1.0
// MATLAB version 5.3 R11
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%