constr_precond

   TTeMPS Toolbox. 
   Michael Steinlechner, 2013-2016
   Questions and contact: michael.steinlechner@epfl.ch
   BSD 2-clause license, see LICENSE.txt

0001 function P = constr_precond( A, k )
0002 
0003     %   TTeMPS Toolbox.
0004     %   Michael Steinlechner, 2013-2016
0005     %   Questions and contact: michael.steinlechner@epfl.ch
0006     %   BSD 2-clause license, see LICENSE.txt
0007 
0008     d = A.order;
0009     ev = eig(A.L0);
0010 
0011     lmin = d*min(ev);
0012     lmax = d*max(ev);
0013 
0014     R = lmax/lmin
0015     
0016     if k == 3
0017         [omega, alpha] = load_coefficients( R );
0018 
0019     elseif k == 7
0020         omega = [0.0133615547183825570028305575534521842940 0.0429728469424360175410925952177443321034 0.1143029399081515586560726591147663100401,...
0021                  0.2838881266934189482611071431161775535656 0.6622322841999484042811198458711174907876 1.4847175320092703810050463464342840325116,...
0022                  3.4859753729916252771962870138366952232900];
0023         alpha = [0.0050213411684266507485648978019454613531 0.0312546410994290844202411500801774835168 0.1045970270084145620410366606112262388706,...
0024                  0.2920522758702768403556507270657505159761 0.7407504784499061527671195936939341208927 1.7609744335543204401530945069076494746696,...
0025                  4.0759036969145123916954953635638503328664];
0026     else
0027         error('Unknown rank specified. Choose either k=3 or k=7');
0028     end
0029 
0030     omega = omega/lmin;
0031     alpha = alpha/lmin;
0032 
0033     E = reshape( expm( -alpha(1) * A.L0), [1, A.size_row(1), A.size_col(1), 1]);
0034     P = omega(1)*TTeMPS_op( repmat({E},1,d) );
0035     for i = 2:k
0036         E = reshape( expm( -alpha(i) * A.L0), [1, A.size_row(1), A.size_col(1), 1]);
0037         P = P + omega(i)*TTeMPS_op( repmat({E},1,d) );
0038     end
0039 
0040 end