inv_solve_trunc_iterative

 INV_SOLVE_TRUNC_ITERATIVE using Morozov truncated iteration
 img= inv_solve_trunc_iterative( inv_model, data1, data2)
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data1      => differential data at earlier time
 data2      => differential data at later time

 both data1 and data2 may be matrices (MxT) each of
  M measurements at T times
 if either data1 or data2 is a vector, then it is expanded
  to be the same size matrix

0001 function img= inv_solve_trunc_iterative( inv_model, data1, data2)
0002 % INV_SOLVE_TRUNC_ITERATIVE using Morozov truncated iteration
0003 % img= inv_solve_trunc_iterative( inv_model, data1, data2)
0004 % img        => output image (or vector of images)
0005 % inv_model  => inverse model struct
0006 % data1      => differential data at earlier time
0007 % data2      => differential data at later time
0008 %
0009 % both data1 and data2 may be matrices (MxT) each of
0010 %  M measurements at T times
0011 % if either data1 or data2 is a vector, then it is expanded
0012 %  to be the same size matrix
0013 
0014 % (C) 2005 David Stephenson. License: GPL version 2 or version 3
0015 % $Id: inv_solve_trunc_iterative.html 2819 2011-09-07 16:43:11Z aadler $
0016 
0017 fwd_model= inv_model.fwd_model;
0018 
0019 img_bkgnd= calc_jacobian_bkgnd( inv_model );
0020 J = calc_jacobian( fwd_model, img_bkgnd);
0021 
0022 % The one_step reconstruction matrix is cached
0023 JtJ = eidors_obj('get-cache', inv_model, 'Hessian');
0024 if ~isempty(JtJ)
0025     eidors_msg('inv_solve_trunc_iterative: using cached value', 2);
0026 else
0027     JtJ= J'*J;
0028 
0029     eidors_obj('set-cache', inv_model, 'Hessian', JtJ);
0030     eidors_msg('inv_solve_trunc_iterative: setting cached value', 2);
0031 end
0032 
0033 l_data1= length(data1); l1_0 = l_data1 ~=0;
0034 l_data2= length(data2); l2_0 = l_data2 ~=0;
0035 l_data= max( l_data1, l_data2 );
0036 
0037 
0038 dva = calc_difference_data( data1, data2, fwd_model);
0039 
0040 tol= 1e-4;
0041 maxiter= 50;
0042 if isfield(inv_model,'parameters');
0043     tol=     inv_model.parameters.term_tolerance;
0044     maxiter= inv_model.parameters.max_iterations;
0045 end
0046 sol = pcg(JtJ, J'*dva, tol, maxiter);
0047 
0048 
0049 % create a data structure to return
0050 img.name= 'solved by inv_solve_trunc_iterative';
0051 img.elem_data = sol;
0052 img.fwd_model= fwd_model;