anisotropicdiffusion

0001 classdef anisotropicdiffusion < TTeMPS_op_laplace
0002     % Class for anisotropic diffusion operator with tridiagonal diffusion matrix
0003     %
0004     %       [ 1 a 0     ...0 ]
0005     %       [ a 1 a 0    ..0 ]
0006     %   D = [ 0 a 1 a 0  ..0 ]
0007     %       [  ..   .. . .   ]
0008     %       [ 0 ... .. 0 a 1 ]
0009     %
0010     %
0011 
0012     %   TTeMPS Toolbox.
0013     %   Michael Steinlechner, 2013-2016
0014     %   Questions and contact: michael.steinlechner@epfl.ch
0015     %   BSD 2-clause license, see LICENSE.txt
0016 
0017     properties
0018         L
0019         D
0020         % precomputed spectral decomp of 1D Laplace:
0021 
0022     end
0023 
0024     methods
0025 
0026         function A = update_properties( A );
0027 
0028             A.rank = [1,  3*ones(1, length(A.U)-1), 1];  % the TT rank is always three for such Laplace-like tensors
0029             size_col_ = cellfun( @(y) size(y,1), A.U);
0030             A.size_col = size_col_ ./ (A.rank(1:end-1).*A.rank(2:end));
0031             A.size_row = cellfun( @(y) size(y,2), A.U);
0032             A.order = length( A.size_row );
0033         end
0034     end
0035 
0036 
0037     methods( Access = public )
0038 
0039         function A = anisotropicdiffusion( n, d, alpha )
0040 
0041             if ~exist('alpha', 'var')
0042                 alpha = 0.25;
0043             end
0044             
0045             one = ones(n,1);
0046             q = linspace( -10, 10, n)';
0047             h = abs(q(2) - q(1));
0048             L = -spdiags( [one, -2*one, one], [-1 0 1], n, n) / (h^2);
0049 
0050             % superclass constructor
0051             A = A@TTeMPS_op_laplace( L, d );
0052             % precompute eigenvalue information and exponential for use in local
0053             A = initialize_precond( A );
0054             % preconditioner:
0055             A.L = L;
0056 
0057             [A.V_L, A.E_L] = eig(full(A.L));
0058             A.E_L = diag(A.E_L);
0059 
0060             A.D = spdiags( [-one,one], [-1,1], n, n ) / (2*h); 
0061             I = speye( n, n );
0062 
0063             e1 = sparse( 1, 1, 1, 3, 1 );
0064             e2 = sparse( 2, 1, 1, 3, 1 );
0065             e3 = sparse( 3, 1, 1, 3, 1 );
0066 
0067             l_mid = sparse( 3, 1, 1, 9, 1 );                % e_3
0068             b_mid = sparse( 6, 1, 1, 9, 1 );                % e_6
0069             m_mid = sparse( [1;9], [1;1], [1;1], 9, 1 );    % e_1 + e_9
0070             c_mid = sparse( 2, 1, 1, 9, 1 );                % e_2
0071 
0072             A.U = cell( 1, d );
0073             A.U{1} = kron( A.L, e1 ) + kron( 2*alpha*A.D, e2 ) + kron( I, e3);
0074             A_mid = kron( A.L, l_mid ) + kron( 2*alpha*A.D, b_mid ) + kron( I, m_mid) + kron( A.D, c_mid );
0075             for i=2:d-1
0076                 A.U{i} = A_mid;
0077             end
0078             A.U{d} = kron( I, e1 ) + kron( A.D, e2 ) + kron( A.L, e3);
0079 
0080             A = update_properties( A );
0081            
0082         end
0083 
0084         function expB = constr_precond_inner( A, X, mu )
0085 
0086             n = size(A.L, 1);
0087             sz = [X.rank(mu), X.size(mu), X.rank(mu+1)]
0088 
0089             B1 = zeros( X.rank(mu) );
0090             % calculate B1 part:
0091             for i = 1:mu-1
0092                 % apply L to the i'th core
0093                 tmp = X;
0094                 Xi = matricize( tmp.U{i}, 2 );
0095                 Xi = A.L*Xi;
0096                 tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
0097                 B1 = B1 + innerprod( X, tmp, 'LR', mu-1);
0098             end
0099 
0100             B3 = zeros( X.rank(mu+1) );
0101             % calculate B3 part:
0102             for i = mu+1:A.order
0103                 tmp = X;
0104                 Xi = matricize( tmp.U{i}, 2 );
0105                 Xi = A.L*Xi;
0106                 tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
0107                 B3 = B3 + innerprod( X, tmp, 'RL', mu+1);
0108             end
0109             
0110             [V1,e1] = eig(B1);
0111             e1 = diag(e1);
0112             [V3,e3] = eig(B3);
0113             e3 = diag(e3);
0114 
0115             lmin = min(e1) + min(A.E_L) + min(e3);
0116             lmax = max(e1) + max(A.E_L) + max(e3);
0117 
0118             R = lmax/lmin
0119             
0120             [omega, alpha] = load_coefficients( R );
0121 
0122             k = 3;
0123             omega = omega/lmin;
0124             alpha = alpha/lmin;
0125 
0126             expB = cell(3,k);
0127             
0128             for i = 1:k
0129                 expB{1,i} = omega(i) * V1*diag( exp( -alpha(i)*e1 ))*V1';    % include omega in first part
0130                 expB{2,i} = A.V_L*diag( exp( -alpha(i)*A.E_L ))*A.V_L';
0131                 expB{3,i} = V3*diag( exp( -alpha(i)*e3 ))*V3';
0132             end
0133         end
0134 
0135         function expB = constr_precond_outer( A, X, mu1, mu2 )
0136             
0137             n = size(A.L, 1);
0138 
0139             B1 = zeros( X.rank(mu1) );
0140             % calculate B1 part:
0141             for i = 1:mu1-1
0142                 % apply L to the i'th core
0143                 tmp = X;
0144                 Xi = matricize( tmp.U{i}, 2 );
0145                 Xi = A.L*Xi;
0146                 tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
0147                 B1 = B1 + innerprod( X, tmp, 'LR', mu1-1);
0148             end
0149 
0150             B3 = zeros( X.rank(mu2+1) );
0151             % calculate B3 part:
0152             for i = mu2+1:A.order
0153                 tmp = X;
0154                 Xi = matricize( tmp.U{i}, 2 );
0155                 Xi = A.L*Xi;
0156                 tmp.U{i} = tensorize( Xi, 2, [X.rank(i), n, X.rank(i+1)] );
0157                 B3 = B3 + innerprod( X, tmp, 'RL', mu2+1);
0158             end
0159             
0160             [V1,e1] = eig(B1);
0161             e1 = diag(e1);
0162             [V3,e3] = eig(B3);
0163             e3 = diag(e3);
0164 
0165             lmin = min(e1) + 2*min(A.E_L) + min(e3);
0166             lmax = max(e1) + 2*max(A.E_L) + max(e3);
0167 
0168             R = lmax/lmin
0169             
0170             [omega, alpha] = load_coefficients( R );
0171 
0172             k = 3;
0173             omega = omega/lmin;
0174             alpha = alpha/lmin;
0175 
0176             expB = cell(4,k);
0177             
0178             for i = 1:k
0179                 expB{1,i} = omega(i) * V1*diag( exp( -alpha(i)*e1 ))*V1';    % include omega in first part
0180                 expB{2,i} = A.V_L*diag( exp( -alpha(i)*A.E_L ))*A.V_L';
0181                 expB{3,i} = A.V_L*diag( exp( -alpha(i)*A.E_L ))*A.V_L';
0182                 expB{4,i} = V3*diag( exp( -alpha(i)*e3 ))*V3';
0183             end
0184         end
0185 
0186         function P = constr_precond( A, k )
0187 
0188             d = A.order;
0189 
0190             lmin = d*min(A.E_L);
0191             lmax = d*max(A.E_L);
0192 
0193             R = lmax/lmin
0194 
0195             %  http://www.mis.mpg.de/scicomp/EXP_SUM/1_x/1_xk07_2E2
0196             %  0.0133615547183825570028305575534521842940   {omega[1]}
0197             %  0.0429728469424360175410925952177443321034   {omega[2]}
0198             %  0.1143029399081515586560726591147663100401   {omega[3]}
0199             %  0.2838881266934189482611071431161775535656   {omega[4]}
0200             %  0.6622322841999484042811198458711174907876   {omega[5]}
0201             %  1.4847175320092703810050463464342840325116   {omega[6]}
0202             %  3.4859753729916252771962870138366952232900   {omega[7]}
0203             %  0.0050213411684266507485648978019454613531   {alpha[1]}
0204             %  0.0312546410994290844202411500801774835168   {alpha[2]}
0205             %  0.1045970270084145620410366606112262388706   {alpha[3]}
0206             %  0.2920522758702768403556507270657505159761   {alpha[4]}
0207             %  0.7407504784499061527671195936939341208927   {alpha[5]}
0208             %  1.7609744335543204401530945069076494746696   {alpha[6]}
0209             %  4.0759036969145123916954953635638503328664   {alpha[7]}
0210             
0211             if k == 3
0212                 [omega, alpha] = load_coefficients( R );
0213 
0214             elseif k == 7
0215                 omega = [0.0133615547183825570028305575534521842940 0.0429728469424360175410925952177443321034 0.1143029399081515586560726591147663100401,...
0216                          0.2838881266934189482611071431161775535656 0.6622322841999484042811198458711174907876 1.4847175320092703810050463464342840325116,...
0217                          3.4859753729916252771962870138366952232900];
0218                 alpha = [0.0050213411684266507485648978019454613531 0.0312546410994290844202411500801774835168 0.1045970270084145620410366606112262388706,...
0219                          0.2920522758702768403556507270657505159761 0.7407504784499061527671195936939341208927 1.7609744335543204401530945069076494746696,...
0220                          4.0759036969145123916954953635638503328664];
0221             else
0222                 error('Unknown rank specified. Choose either k=3 or k=7');
0223             end
0224 
0225             omega = omega/lmin;
0226             alpha = alpha/lmin;
0227 
0228             % careful: all cores assumed to be of same size
0229             E = reshape( expm( -alpha(1) * A.L), [1, A.size_row(2), A.size_col(2), 1]);
0230             P = omega(1)*TTeMPS_op( repmat({E},1,d) );
0231             for i = 2:k
0232                 E = reshape( expm( -alpha(i) * A.L), [1, A.size_row(2), A.size_col(2), 1]);
0233                 P = P + omega(i)*TTeMPS_op( repmat({E},1,d) );
0234             end
0235 
0236         end
0237 
0238     end
0239         
0240 end