newton_potential

0001 classdef newton_potential < TTeMPS_op_laplace
0002     % Class for the Newton potential derived from the TTeMPS_op_laplace class.
0003 
0004     %   TTeMPS Toolbox.
0005     %   Michael Steinlechner, 2013-2016
0006     %   Questions and contact: michael.steinlechner@epfl.ch
0007     %   BSD 2-clause license, see LICENSE.txt
0008 
0009     properties
0010         potential
0011         laplacerank
0012     end
0013 
0014     methods( Access = public )
0015 
0016         function A = newton_potential( n, d )
0017             
0018             one = ones(n,1);
0019             q = linspace( -1, 1, n)';
0020             h = abs(q(2) - q(1));
0021             L = -spdiags( [one, -2*one, one], [-1 0 1], n, n) / (h^2);            
0022 
0023             % superclass constructor
0024             A = A@TTeMPS_op_laplace( L, d );
0025             [A.V_L, A.E_L] = eig(full(A.L0));
0026             A.E_L = diag( A.E_L );
0027 
0028             % create potential
0029 
0030             a = d*min(q.^2);
0031             b = d*max(q.^2);
0032 
0033             k = 10;
0034             %  http://www.mis.mpg.de/scicomp/EXP_SUM/1_sqrtx/1_sqrtxk10_2E4
0035             %  0.0122353280385689973405283041685276401722   {omega[1]}
0036             %  0.0170048221567127200615394955196535420328   {omega[2]}
0037             %  0.0292580238127349109583679896348651361393   {omega[3]}
0038             %  0.0513674405142589591607559396796434114663   {omega[4]}
0039             %  0.0876522649954556273282240720645663856203   {omega[5]}
0040             %  0.1452687747918532464746523003018552344656   {omega[6]}
0041             %  0.2347776097665007749814821205736059539504   {omega[7]}
0042             %  0.3716699993620031342879666408363092955369   {omega[8]}
0043             %  0.5822367252037624550361008535226403637353   {omega[9]}
0044             %  0.9561416473224761119619093119315067497155   {omega[10]}
0045             %  0.0000278884968834587275257159520245177180   {alpha[1]}
0046             %  0.0003161367116046333453815662397571456532   {alpha[2]}
0047             %  0.0014174384495922332182950292714732065669   {alpha[3]}
0048             %  0.0052598987752533055729865546674972609509   {alpha[4]}
0049             %  0.0176506185630681001889793731510214236380   {alpha[5]}
0050             %  0.0548279527792261968593219480933020903990   {alpha[6]}
0051             %  0.1597705718538013465668734189306654513985   {alpha[7]}
0052             %  0.4411611832132119288894626929486975086547   {alpha[8]}
0053             %  1.1661865464705391523319785718193486445671   {alpha[9]}
0054             %  3.0303805868035630160465393467816852535179   {alpha[10]}
0055 
0056             omega = [0.0122353280385689973405283041685276401722; 0.0170048221567127200615394955196535420328;
0057                     0.0292580238127349109583679896348651361393; 0.0513674405142589591607559396796434114663;
0058                     0.0876522649954556273282240720645663856203; 0.1452687747918532464746523003018552344656;
0059                     0.2347776097665007749814821205736059539504; 0.3716699993620031342879666408363092955369;
0060                     0.5822367252037624550361008535226403637353; 0.9561416473224761119619093119315067497155];
0061             alpha = [0.0000278884968834587275257159520245177180; 0.0003161367116046333453815662397571456532;
0062                     0.0014174384495922332182950292714732065669; 0.0052598987752533055729865546674972609509;
0063                     0.0176506185630681001889793731510214236380; 0.0548279527792261968593219480933020903990;
0064                     0.1597705718538013465668734189306654513985; 0.4411611832132119288894626929486975086547;
0065                     1.1661865464705391523319785718193486445671; 3.0303805868035630160465393467816852535179];
0066 
0067             alpha = alpha/a;
0068             omega = omega/sqrt(a);
0069 
0070             C = cell(1, d);
0071 
0072             C = exp(-alpha(1)*q.^2);
0073             C = reshape( C, [1 n 1] );
0074             A.potential = omega(1)*TTeMPS( repmat( {C}, [1, d]));
0075             for i = 2:k
0076                 C = exp(-alpha(i)*q.^2);
0077                 C = reshape( C, [1 n 1] );
0078                 A.potential = A.potential + omega(i)*TTeMPS( repmat( {C}, [1, d]));
0079             end
0080 
0081             A.laplacerank = A.rank;
0082             A.rank(2:end-1) = A.rank(2:end-1) + A.potential.rank(2:end-1);
0083             
0084         end
0085 
0086 
0087         
0088         function y = apply(A, x, idx)
0089             
0090             if ~exist( 'idx', 'var' )
0091                 % apply laplace part
0092                 A.rank = A.laplacerank;
0093                 y = apply@TTeMPS_op_laplace( A, x );
0094                 A.rank(2:end-1) = A.rank(2:end-1) + A.potential.rank(2:end-1);
0095                 
0096                 % apply henon-heiles part
0097 
0098                 y = y + hadamard( A.potential, x );
0099                 %y = hadamard( A.potential, x );
0100 
0101             else
0102                 % apply laplace part
0103                 A.rank = A.laplacerank;
0104                 l = apply@TTeMPS_op_laplace( A, x, idx );
0105                 A.rank(2:end-1) = A.rank(2:end-1) + A.potential.rank(2:end-1);
0106                 
0107                 % apply henon-heiles part
0108 
0109                 h = hadamard( A.potential, x, idx);
0110 
0111                 if idx == 1
0112                     y = zeros( 1, size(l,2), size(l,3)+size(h,3), size(h,4));
0113                     y( 1, :, 1:size(l,3), : ) = l;
0114                     y( 1, :, size(l,3)+1:end, : ) = h;
0115                 elseif idx == A.order
0116                     y = zeros( size(l,1)+size(h,1), size(l,2), 1, size(h,4) );
0117                     y( 1:size(l,1), :, 1, :) = l;
0118                     y( size(l,1)+1:end, :, 1, :) = h;
0119                 else
0120                     y = zeros( size(l,1)+size(h,1), size(l,2), size(l,3)+size(h,3), size(h,4) );
0121                     y( 1:size(l,1), :, 1:size(l,3), : ) = l;
0122                     y( size(l,1)+1:end, :, size(l,3)+1:end, : ) = h;
0123                 end
0124 
0125             end
0126         end
0127         
0128        
0129         
0130         
0131 
0132         function P = constr_precond( A, k )
0133 
0134             d = A.order;
0135             ev = eig(full(A.L0));
0136 
0137             lmin = d*min(ev);
0138             lmax = d*max(ev);
0139 
0140             R = lmax/lmin;
0141 
0142             %  http://www.mis.mpg.de/scicomp/EXP_SUM/1_x/1_xk07_2E2
0143             %  0.0133615547183825570028305575534521842940   {omega[1]}
0144             %  0.0429728469424360175410925952177443321034   {omega[2]}
0145             %  0.1143029399081515586560726591147663100401   {omega[3]}
0146             %  0.2838881266934189482611071431161775535656   {omega[4]}
0147             %  0.6622322841999484042811198458711174907876   {omega[5]}
0148             %  1.4847175320092703810050463464342840325116   {omega[6]}
0149             %  3.4859753729916252771962870138366952232900   {omega[7]}
0150             %  0.0050213411684266507485648978019454613531   {alpha[1]}
0151             %  0.0312546410994290844202411500801774835168   {alpha[2]}
0152             %  0.1045970270084145620410366606112262388706   {alpha[3]}
0153             %  0.2920522758702768403556507270657505159761   {alpha[4]}
0154             %  0.7407504784499061527671195936939341208927   {alpha[5]}
0155             %  1.7609744335543204401530945069076494746696   {alpha[6]}
0156             %  4.0759036969145123916954953635638503328664   {alpha[7]}
0157             
0158             if k == 3
0159                 [omega, alpha] = load_coefficients( R );
0160 
0161             elseif k == 7
0162                 omega = [0.0133615547183825570028305575534521842940 0.0429728469424360175410925952177443321034 0.1143029399081515586560726591147663100401,...
0163                          0.2838881266934189482611071431161775535656 0.6622322841999484042811198458711174907876 1.4847175320092703810050463464342840325116,...
0164                          3.4859753729916252771962870138366952232900];
0165                 alpha = [0.0050213411684266507485648978019454613531 0.0312546410994290844202411500801774835168 0.1045970270084145620410366606112262388706,...
0166                          0.2920522758702768403556507270657505159761 0.7407504784499061527671195936939341208927 1.7609744335543204401530945069076494746696,...
0167                          4.0759036969145123916954953635638503328664];
0168             else
0169                 error('Unknown rank specified. Choose either k=3 or k=7');
0170             end
0171 
0172             omega = omega/lmin;
0173             alpha = alpha/lmin;
0174 
0175             E = reshape( expm( -alpha(1) * A.L0), [1, A.size_row(d), A.size_col(d), 1]);
0176             P = omega(1)*TTeMPS_op( repmat({E},1,d) );
0177             for i = 2:k
0178                 E = reshape( expm( -alpha(i) * A.L0), [1, A.size_row(d), A.size_col(d), 1]);
0179                 P = P + omega(i)*TTeMPS_op( repmat({E},1,d) );
0180             end
0181 
0182         end
0183 
0184         function expB = constr_precond_inner( A, X, mu )
0185             
0186             %V = reshape( V, [1, 1, size(L, 1), tmp.rank(i), tmp.rank(i+1)] );
0187             %V = permute( V, [1, 4, 3, 2, 5]);
0188             %tmp.U{i} = reshape( V, [tmp.rank(i), size(L, 1), tmp.rank{i+1}]);
0189 
0190             n = size(A.L0, 1);
0191             sz = [X.rank(mu), X.size(mu), X.rank(mu+1)];
0192 
0193             B1 = zeros( X.rank(mu) );
0194             % calculate B1 part:
0195             for i = 1:mu-1
0196                 % apply L to the i'th core
0197                 tmp = X;
0198                 Xi = matricize( tmp.U{i}, 2 );
0199                 Xi = A.L0*Xi;
0200                 tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
0201                 B1 = B1 + innerprod( X, tmp, 'LR', mu-1);
0202             end
0203 
0204             % calculate B2 part:
0205             B2 = A.L0;
0206 
0207             B3 = zeros( X.rank(mu+1) );
0208             % calculate B3 part:
0209             for i = mu+1:A.order
0210                 tmp = X;
0211                 Xi = matricize( tmp.U{i}, 2 );
0212                 Xi = A.L0*Xi;
0213                 tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
0214                 B3 = B3 + innerprod( X, tmp, 'RL', mu+1);
0215             end
0216             
0217             [V1,e1] = eig(B1);
0218             e1 = diag(e1);
0219             [V3,e3] = eig(B3);
0220             e3 = diag(e3);
0221 
0222             lmin = min(e1) + min(A.E_L) + min(e3);
0223             lmax = max(e1) + max(A.E_L) + max(e3);
0224 
0225             R = lmax/lmin;
0226             
0227             [omega, alpha] = load_coefficients( R );
0228 
0229             k = length(omega);
0230             omega = omega/lmin;
0231             alpha = alpha/lmin;
0232 
0233             expB = cell(3,k);
0234             
0235             for i = 1:k
0236                 expB{1,i} = omega(i) * V1*diag( exp( -alpha(i)*e1 ))*V1';    % include omega in first part
0237                 expB{2,i} = A.V_L*diag( exp( -alpha(i)*A.E_L ))*A.V_L';
0238                 expB{3,i} = V3*diag( exp( -alpha(i)*e3 ))*V3';
0239             end
0240         end
0241 
0242         function expB = constr_precond_outer( A, X, mu1, mu2 )
0243             
0244             %V = reshape( V, [1, 1, size(L, 1), tmp.rank(i), tmp.rank(i+1)] );
0245             %V = permute( V, [1, 4, 3, 2, 5]);
0246             %tmp.U{i} = reshape( V, [tmp.rank(i), size(L, 1), tmp.rank{i+1}]);
0247 
0248             n = size(A.L0, 1);
0249 
0250             B1 = zeros( X.rank(mu1) );
0251             % calculate B1 part:
0252             for i = 1:mu1-1
0253                 % apply L to the i'th core
0254                 tmp = X;
0255                 Xi = matricize( tmp.U{i}, 2 );
0256                 Xi = A.L0*Xi;
0257                 tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
0258                 B1 = B1 + innerprod( X, tmp, 'LR', mu1-1);
0259             end
0260 
0261             % calculate B2 part:
0262             %D = spdiags( ones(n,1), 0, n, n);
0263             %B2 = kron(A.L, D) + kron(D, A.L);
0264 
0265             B3 = zeros( X.rank(mu2+1) );
0266             % calculate B3 part:
0267             for i = mu2+1:A.order
0268                 tmp = X;
0269                 Xi = matricize( tmp.U{i}, 2 );
0270                 Xi = A.L0*Xi;
0271                 tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
0272                 B3 = B3 + innerprod( X, tmp, 'RL', mu2+1);
0273             end
0274             
0275             [V1,e1] = eig(B1);
0276             e1 = diag(e1);
0277             [V3,e3] = eig(B3);
0278             e3 = diag(e3);
0279 
0280             lmin = min(e1) + 2*min(A.E_L) + min(e3);
0281             lmax = max(e1) + 2*max(A.E_L) + max(e3);
0282 
0283             R = lmax/lmin;
0284             
0285             [omega, alpha] = load_coefficients( R );
0286 
0287             k = 3;
0288             omega = omega/lmin;
0289             alpha = alpha/lmin;
0290 
0291             expB = cell(4,k);
0292             
0293             for i = 1:k
0294                 expB{1,i} = omega(i) * V1*diag( exp( -alpha(i)*e1 ))*V1';    % include omega in first part
0295                 expB{2,i} = A.V_L*diag( exp( -alpha(i)*A.E_L ))*A.V_L';
0296                 expB{3,i} = A.V_L*diag( exp( -alpha(i)*A.E_L ))*A.V_L';
0297                 expB{4,i} = V3*diag( exp( -alpha(i)*e3 ))*V3';
0298             end
0299         end
0300 
0301     end
0302         
0303 end