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.m 4832 2015-03-29 21:13:53Z bgrychtol-ipa $ 0016 0017 fwd_model= inv_model.fwd_model; 0018 0019 img_bkgnd= calc_jacobian_bkgnd( inv_model ); 0020 J = calc_jacobian(img_bkgnd); 0021 0022 % The one_step reconstruction matrix is cached 0023 JtJ = eidors_cache(@calc_hessian, J, copt); 0024 0025 l_data1= length(data1); l1_0 = l_data1 ~=0; 0026 l_data2= length(data2); l2_0 = l_data2 ~=0; 0027 l_data= max( l_data1, l_data2 ); 0028 0029 0030 dva = calc_difference_data( data1, data2, fwd_model); 0031 0032 tol= 1e-4; 0033 maxiter= 50; 0034 if isfield(inv_model,'parameters'); 0035 tol= inv_model.parameters.term_tolerance; 0036 maxiter= inv_model.parameters.max_iterations; 0037 end 0038 sol = pcg(JtJ, J'*dva, tol, maxiter); 0039 0040 0041 % create a data structure to return 0042 img.name= 'solved by inv_solve_trunc_iterative'; 0043 img.elem_data = sol; 0044 img.fwd_model= fwd_model; 0045 0046 function JtJ = calc_hessian(J) 0047 JtJ = J'*J;