Home > manopt > manifolds > ttfixedrank > TTeMPS_1.1 > algorithms > eigenvalue > amen_eigenvalue.m

# amen_eigenvalue

## PURPOSE

AMEN_EIGENVALUE Calculate p smallest eigenvalues of a TTeMPS operator

## SYNOPSIS

function [X, C, evalue, residuums, micro_res, objective, elapsed_time] = amen_eigenvalue(A, prec, p, rr, opts)

## DESCRIPTION

```AMEN_EIGENVALUE Calculate p smallest eigenvalues of a TTeMPS operator

[X, C, evalue, residuums, micro_res, objective, elapsed_time] = amen_eigenvalue(A, PREC, P, RR, OPTS)
performs a rank-adaptive AMEn-type optimization scheme to compute the p smallest eigenvalues of A
using the algorithm described in [1].

A preconditioner can be given by passing the PREC argument:
[]      no preconditioner used
PR      with PR of type TTeMPS_op: PR is applied as a preconditioner
true    exponential sum preconditioner described in [1] (recommended setting)

RR defines the starting rank and should usually be set to ones(1,d) where d is the dimensionality.

Specify the wanted options using the OPTS argument. All options have
default values (denoted in brackets). Hence, you only have to specify those you want
to change.
The OPTS struct is organized as follows:
OPTS.maxiter        Maximum number of full sweeps to perform        [3]
OPTS.maxrank        Maximum rank during the iteration               [40]
OPTS.maxrankRes     Maximum rank of residual (0 = No limit)         [0]
OPTS.tol            Tolerance for the shift from one core
to the next                                     [1e-8]
OPTS.tolOP          Tolerance for truncation of residual            [1e-6]
OPTS.tolLOBPCG      Tolerance for inner LOBPCG solver               [1e-6]
OPTS.maxiterLOBPCG  Max. number of iterations for LOBPCG            [2000]
OPTS.verbose        Show iteration information (true/false)         [true]
OPTS.precInner      Precondition the LOBPCG (STRONGLY RECOMMENDED!) [true]

NOTE: To run amen_eigenvalue, Knyazev's implementation of LOBPCG is required. The corresponding
http://www.mathworks.com/matlabcentral/fileexchange/48-lobpcg-m

References:
[1] D. Kressner, M. Steinlechner, A. Uschmajew:
Low-rank tensor methods with subspace correction for symmetric eigenvalue problems
SIAM J. Sci. Comput., 36(5):A2346-A2368, 2014.```

## CROSS-REFERENCE INFORMATION

This function calls:
• cat CAT concatenation of two TT/MPS tensors.
• contract CONTRACT Contraction of two TT/MPS tensors.
• disp DISP Display TT/MPS tensor.
• innerprod INNERPROD Inner product between two TT/MPS tensors.
• norm NORM Norm of a TT/MPS tensor.
• orthogonalize ORTHOGONALIZE Orthogonalize tensor.
• truncate TRUNCATE Truncate TTeMPS tensor to prescribed rank.
• disp DISP Display TT/MPS block-mu tensor.
• innerprod INNERPROD Inner product between two TT/MPS tensors.
• norm NORM Norm of a TT/MPS block-mu tensor.
• orthogonalize ORTHOGONALIZE Orthogonalize TT/MPS Block-mu tensor.
• truncate ROUND Approximate TTeMPS tensor within a prescribed tolerance.
• apply APPLY Application of TT/MPS operator to a TT/MPS tensor
• contract CONTRACT Contraction of two TT/MPS tensors with inner TT/MPS operator.
• disp DISP Display TT/MPS operator.
• apply APPLY Application of TT/MPS Laplace-like operator to a TT/MPS tensor
• constr_precond_inner TTeMPS Toolbox.
• contract CONTRACT Contraction of two TT/MPS tensors with inner TT/MPS Laplace operator.
• disp DISP Display TT/MPS operator.
• TTeMPS_rand TTEMPS_RAND Create random TTeMPS tensor
• tensorprod_ttemps TENSORPROD_TTEMPS Tensor-times-Matrix product.
• trunc_singular REL_TRUNC_SINGULAR Helper routine to truncate singular values
• unfold UNFOLD Left/right-unfold a 3D array.
• orthogonalize Orthonormalizes a basis of tangent vectors in the Manopt framework.
This function is called by:

## SOURCE CODE

```0001 function [X, C, evalue, residuums, micro_res, objective, elapsed_time] = amen_eigenvalue(A, prec, p, rr, opts)
0002     %AMEN_EIGENVALUE Calculate p smallest eigenvalues of a TTeMPS operator
0003     %
0004     %   [X, C, evalue, residuums, micro_res, objective, elapsed_time] = amen_eigenvalue(A, PREC, P, RR, OPTS)
0005     %       performs a rank-adaptive AMEn-type optimization scheme to compute the p smallest eigenvalues of A
0006     %       using the algorithm described in [1].
0007     %
0008     %   A preconditioner can be given by passing the PREC argument:
0009     %       []      no preconditioner used
0010     %       PR      with PR of type TTeMPS_op: PR is applied as a preconditioner
0011     %       true    exponential sum preconditioner described in [1] (recommended setting)
0012     %
0013     %   RR defines the starting rank and should usually be set to ones(1,d) where d is the dimensionality.
0014     %
0015     %   Specify the wanted options using the OPTS argument. All options have
0016     %   default values (denoted in brackets). Hence, you only have to specify those you want
0017     %   to change.
0018     %   The OPTS struct is organized as follows:
0019     %       OPTS.maxiter        Maximum number of full sweeps to perform        [3]
0020     %       OPTS.maxrank        Maximum rank during the iteration               [40]
0021     %       OPTS.maxrankRes     Maximum rank of residual (0 = No limit)         [0]
0022     %       OPTS.tol            Tolerance for the shift from one core
0023     %                           to the next                                     [1e-8]
0024     %       OPTS.tolOP          Tolerance for truncation of residual            [1e-6]
0025     %       OPTS.tolLOBPCG      Tolerance for inner LOBPCG solver               [1e-6]
0026     %       OPTS.maxiterLOBPCG  Max. number of iterations for LOBPCG            [2000]
0027     %       OPTS.verbose        Show iteration information (true/false)         [true]
0028     %       OPTS.precInner      Precondition the LOBPCG (STRONGLY RECOMMENDED!) [true]
0029     %
0030     %
0031     %   NOTE: To run amen_eigenvalue, Knyazev's implementation of LOBPCG is required. The corresponding
0033     %               http://www.mathworks.com/matlabcentral/fileexchange/48-lobpcg-m
0034     %
0035     %    References:
0036     %    [1] D. Kressner, M. Steinlechner, A. Uschmajew:
0037     %        Low-rank tensor methods with subspace correction for symmetric eigenvalue problems
0038     %        SIAM J. Sci. Comput., 36(5):A2346-A2368, 2014.
0039
0040     %   TTeMPS Toolbox.
0041     %   Michael Steinlechner, 2013-2016
0042     %   Questions and contact: michael.steinlechner@epfl.ch
0044
0045
0046 nn = A.size_col;
0047 d = A.order;
0048
0049 if ~isempty(prec)
0050     if isa(prec,'TTeMPS_op')
0051         precondition = 1;
0052     else
0053         precondition = 2;
0054     end
0055 else
0056     precondition = 0;
0057 end
0058
0059 if ~isfield( opts, 'maxiter');       opts.maxiter = 3;           end
0060 if ~isfield( opts, 'maxrank');       opts.maxrank = 40;          end
0061 if ~isfield( opts, 'maxrankRes');    opts.maxrankRes = 0;        end
0062 if ~isfield( opts, 'tol');           opts.tol = 1e-8;            end
0063 if ~isfield( opts, 'tolOP');         opts.tolOP = 1e-6;          end
0064 if ~isfield( opts, 'tolLOBPCG');     opts.tolLOBPCG = 1e-6;      end
0065 if ~isfield( opts, 'maxiterLOBPCG'); opts.maxiterLOBPCG = 2000;  end
0066 if ~isfield( opts, 'verbose');       opts.verbose = true;        end
0067 if ~isfield( opts, 'precInner');     opts.precInner = true;     end
0068
0069 if p == 1
0070     tolLOBPCGmod = opts.tolLOBPCG;
0071 else
0072     tolLOBPCGmod = 0.01;
0073 end
0074
0075 X = TTeMPS_rand( rr, nn );
0076 % Left-orthogonalize the tensor:
0077 X = orthogonalize( X, X.order );
0078
0079 C = cell( 1, p );
0080 C{1} = X.U{d};
0081 for i = 2:p
0082     C{i} = rand(size(X.U{d}));
0083 end
0084
0085 X.U{d} = rand( X.rank(d), X.size(d), X.rank(d+1), p );
0086
0087 evalue = [];
0088 residuums = zeros(p,2*opts.maxiter);
0089 micro_res = [];
0090 resi_norm = zeros(p,1);
0091 objective = [];
0092 tic;
0093 elapsed_time = [];
0094
0095 for i = 1:opts.maxiter
0096
0097     disp(sprintf('Iteration %i:', i));
0098     % right-to-left sweep
0099     fprintf(1, 'RIGHT-TO-LEFT SWEEP. ---------------\n')
0100
0101     for mu = d:-1:2
0102
0103         sz = [X.rank(mu), X.size(mu), X.rank(mu+1)];
0104
0105         if opts.verbose
0106             fprintf(1,'Current core: %i. Current iterate (first eigenvalue): \n', mu);
0107             disp(X);
0108             fprintf(1,'Running LOBPCG: system size %i, tol %g ... ', prod(sz),  max(opts.tolLOBPCG, min( tolLOBPCGmod, sqrt(sum(residuums(:,2*(i-1)+1).^2))/sqrt(p)*tolLOBPCGmod ))) %opts.tolLOBPCG)
0109         else
0110             fprintf(1,'%i  ',mu);
0111         end
0112
0113         [left, right] = Afun_prepare( A, X, mu );
0114
0115         if opts.precInner
0116             if prod(sz) < 5*p
0117                 Amat = zeros( prod(sz) );
0118                 for i_mat = 1:prod(sz)
0119                     Amat(:,i_mat) = Afun_block_optim( A, [zeros(i_mat-1,1); 1; zeros(prod(sz)-i_mat,1)], sz, left, right, mu);
0120                 end
0121                 [vmat,emat] = eig(Amat);
0122                 [emat,sortind] = sort(diag(emat),'ascend');
0123                 vmat = vmat(:, sortind);
0124                 L = emat(1:p);
0125                 V = vmat(:, 1:p);
0126
0127                 disp('Reduced system too small for LOBPCG, using eig instead');
0128                 disp(['Current Eigenvalue approx: ', num2str( L(1:p)')]);
0129
0130             else
0131                 expB = constr_precond_inner( A, X, mu );
0132                 [V,L,failureFlag,Lhist ] = lobpcg( rand( prod(sz), p), ...
0133                                                    @(y) Afun_block_optim( A, y, sz, left, right, mu), [], ...
0134                                                    @(y) apply_local_precond( y, sz, expB), ...
0135                                                    max(opts.tolLOBPCG, min( tolLOBPCGmod, sqrt(sum(residuums(:,2*(i-1)+1).^2))/sqrt(p)*tolLOBPCGmod )), opts.maxiterLOBPCG, 0);
0136                 if opts.verbose
0137                     if failureFlag
0138                         fprintf(1,'NOT CONVERGED within %i steps!\n', opts.maxiterLOBPCG)
0139                     else
0140                         fprintf(1,'converged after %i steps!\n', size(Lhist,2));
0141                     end
0142                     disp(['Current Eigenvalue approx: ', num2str( L(1:p)')]);
0143                 end
0144             end
0145         else
0146             [V,L,failureFlag,Lhist ] = lobpcg( rand( prod(sz), p), ...
0147                                                @(y) Afun_block_optim( A, y, sz, left, right, mu), ...
0148                                                opts.tolLOBPCG, opts.maxiterLOBPCG, 0);
0149         end
0150
0151         evalue = [evalue, L(1:p)];
0152         objective = [objective, sum(evalue(:,end))];
0153         elapsed_time = [elapsed_time, toc];
0154
0155         X.U{mu} = reshape( V, [sz, p] );
0156         lamX = X;
0157         lamX.U{mu} = repmat( reshape(L, [1 1 1 p]), [X.rank(mu), X.size(mu), X.rank(mu+1), 1]).*lamX.U{mu};
0158
0159         res = apply(A, X) - lamX;
0160
0161         if opts.maxrankRes ~= 0
0162             res_sz = [size(res.U{mu},1), size(res.U{mu},2), size(res.U{mu},3), size(res.U{mu},4)];
0163             res.U{mu} = reshape( permute( res.U{mu}, [1 2 4 3]), [ res_sz(1), res_sz(2)*p, res_sz(3)]);
0164             res = truncate( res, [1 opts.maxrankRes*ones(1,d-1) 1] );
0165             res.U{mu} = ipermute( reshape( res.U{mu}, [res.rank(mu), res_sz(2), p, res.rank(mu+1)] ), [1 2 4 3]);
0166         end
0167
0168         % A not so nice way to calculate the residual norm
0169         for j = 1:p
0170             tmp = res;
0171             tmp.U{mu} = tmp.U{mu}(:,:,:,j);
0172             resi_norm(j) = norm(tmp);
0173             residuums(j,2*(i-1)+1) = resi_norm(j);
0174         end
0175         micro_res = [micro_res, resi_norm];
0176
0177         if precondition == 1
0178             res = apply(prec, res);
0179         end
0180
0181         res = contract( X, res, [mu-1, mu]);
0182         res_combined = unfold(res{1},'left') * reshape(res{2}, [size(res{2},1), size(res{2},2)*size(res{2},3)*p]);
0183
0184         if precondition == 2
0185             expBglobal = constr_precond_outer( A, X, mu-1, mu );
0186             res_size = [size(res{1},1), size(res{1},2), size(res{2},2), size(res{2},3), p];
0187             res_combined = reshape( res_combined, res_size );
0188             res_combined = apply_global_precond( res_combined, res_size, expBglobal );
0189         end
0190
0191         res_combined = reshape( res_combined, [size(res{1},1)*size(res{1},2), size(res{2},2)*size(res{2},3)*p]);
0192
0193         [uu,ss,vv] = svd( res_combined, 'econ');
0194         s = trunc_singular( diag(ss), opts.tolOP, true, opts.maxrankRes );
0195
0196         res{1} = reshape( uu(:,1:s), [size(res{1},1), size(res{1},2), s]);
0197         res{2} = reshape( ss(1:s,1:s)*vv(:,1:s)', [s, size(res{2},2), size(res{2},3), p]);
0198
0199         left = cat(3, X.U{mu-1}, res{1} );
0200         V = cat(1, X.U{mu}, res{2});
0201
0202         tmp = permute( V, [1, 4, 2, 3] );
0203         V = reshape( tmp, [size(tmp,1)*p, size(tmp,3)*size(tmp,4)] );
0204
0205         [U,S,V] = svd( V, 'econ' );
0206         s = trunc_singular( diag(S), opts.tol, true, opts.maxrank );
0207
0208         V = V(:,1:s)';
0209         X.U{mu} = reshape( V, [s, sz(2), sz(3)] );
0210
0211         W = U(:,1:s)*S(1:s,1:s);
0212         W = reshape( W, [size(tmp,1), p, s]);
0213         W = permute( W, [1, 3, 2]);
0214         for k = 1:p
0215             C{k} = tensorprod_ttemps( left, W(:,:,k)', 3);
0216         end
0217         X.U{mu-1} = C{1};
0218
0219         if opts.verbose
0220             disp( ['Augmented system of size (', num2str( [sz(1)*p, sz(2)*sz(3)]), '). Cut-off tol: ', num2str(opts.tol) ])
0221             disp( sprintf( 'Number of SVs: %i. Truncated to: %i => %g %%', length(diag(S)), s, s/length(diag(S))*100))
0222             disp(' ')
0223         end
0224
0225     end
0226
0227     % calculate current residuum
0228
0229
0230     fprintf(1, '---------------    finshed sweep.    ---------------\n')
0231     disp(['Current residuum: ', num2str(residuums(:,2*(i-1)+1).')])
0232     disp(' ')
0233     disp(' ')
0234     fprintf(1, '--------------- LEFT-TO-RIGHT SWEEP. ---------------\n')
0235     % left-to-right sweep
0236     for mu = 1:d-1
0237         sz = [X.rank(mu), X.size(mu), X.rank(mu+1)];
0238
0239         if opts.verbose
0240             fprintf(1,'Current core: %i. Current iterate (first eigenvalue): \n', mu);
0241             disp(X);
0242             fprintf(1,'Running LOBPCG: system size %i, tol %g ... ', prod(sz), max(opts.tolLOBPCG, min( tolLOBPCGmod, sqrt(sum(residuums(:,2*(i-1)+2).^2))/sqrt(p)*tolLOBPCGmod )))
0243         else
0244             fprintf(1,'%i  ',mu);
0245         end
0246
0247         [left, right] = Afun_prepare( A, X, mu );
0248         if opts.precInner
0249             expB = constr_precond_inner( A, X, mu );
0250             [U,L,failureFlag,Lhist ] = lobpcg( rand( prod(sz), p), ...
0251                                                @(y) Afun_block_optim( A, y, sz, left, right, mu), [], ...
0252                                                @(y) apply_local_precond( y, sz, expB), ...
0253                                                max(opts.tolLOBPCG, min( tolLOBPCGmod, sqrt(sum(residuums(:,2*(i-1)+2).^2))/sqrt(p)*tolLOBPCGmod )), ...
0254                                                opts.maxiterLOBPCG, 0);
0255         else
0256             [U,L,failureFlag,Lhist ] = lobpcg( rand( prod(sz), p), ...
0257                                                @(y) Afun_block_optim( A, y, sz, left, right, mu), ...
0258                                                opts.tolLOBPCG, opts.maxiterLOBPCG, 0);
0259         end
0260
0261         if opts.verbose
0262             if failureFlag
0263                 fprintf(1,'NOT CONVERGED within %i steps!\n', opts.maxiterLOBPCG)
0264             else
0265                 fprintf(1,'converged after %i steps!\n', size(Lhist,2));
0266             end
0267             disp(['Current Eigenvalue approx: ', num2str( L(1:p)')]);
0268         end
0269
0270         evalue = [evalue, L(1:p)];
0271         objective = [objective, sum(evalue(:,end))];
0272         elapsed_time = [elapsed_time, toc];
0273
0274         X.U{mu} = reshape( U, [sz, p] );
0275         lamX = X;
0276         lamX.U{mu} = repmat( reshape(L, [1 1 1 p]), [X.rank(mu), X.size(mu), X.rank(mu+1), 1]).*lamX.U{mu};
0277         res = apply(A, X) - lamX;
0278
0279         if opts.maxrankRes ~= 0
0280             res_sz = [size(res.U{mu},1), size(res.U{mu},2), size(res.U{mu},3), size(res.U{mu},4)];
0281             res.U{mu} = reshape( permute( res.U{mu}, [1 2 4 3]), [ res_sz(1), res_sz(2)*p, res_sz(3)]);
0282             res = truncate( res, [1 opts.maxrankRes*ones(1,d-1) 1] );
0283             res.U{mu} = ipermute( reshape( res.U{mu}, [ res.rank(mu), res_sz(2), p, res.rank(mu+1)] ), [1 2 4 3]);
0284         end
0285
0286         % ugly hack to calc the residual
0287         for j = 1:p
0288             tmp = res;
0289             tmp.U{mu} = tmp.U{mu}(:,:,:,j);
0290             resi_norm(j) = norm(tmp);
0291             residuums(j,2*(i-1)+2) = resi_norm(j);
0292         end
0293         micro_res = [micro_res, resi_norm];
0294
0295         if precondition == 1
0296             res = apply(prec, res);
0297         end
0298         res = contract( X, res, [mu, mu+1]);
0299         res_left = permute( res{1}, [1 2 4 3] );
0300         res_left = reshape( res_left, [size(res{1},1)*size(res{1},2)*p, size(res{1},3)]);
0301         res_combined =  res_left * unfold(res{2},'right');
0302
0303         if precondition == 2
0304             expBglobal = constr_precond_outer( A, X, mu, mu+1 );
0305             res_size = [size(res{1},1), size(res{1},2), p, size(res{2},2), size(res{2},3)];
0306             res_combined = reshape( res_combined, res_size );
0307             res_combined = permute( res_combined, [1 2 4 5 3]);  % move p to back
0308             res_combined = apply_global_precond( res_combined, res_size([1 2 4 5 3]), expBglobal );
0309             res_combined = ipermute( res_combined, [1 2 4 5 3]);  % move p to back
0310         end
0311
0312         res_combined = reshape( res_combined, [size(res{1},1)*size(res{1},2)*p, size(res{2},2)*size(res{2},3)]);
0313
0314         [uu,ss,vv] = svd( res_combined, 'econ');
0315         s = trunc_singular( diag(ss), opts.tolOP, true, opts.maxrankRes );
0316
0317         res{1} = reshape( uu(:,1:s)*ss(1:s,1:s), [size(res{1},1), size(res{1},2), p, s]);
0318         res{2} = reshape( vv(:,1:s)', [s, size(res{2},2), size(res{2},3)]);
0319
0320         right = cat(1, X.U{mu+1}, res{2} );
0321         res{1} = permute( res{1}, [1 2 4 3]);
0322         U = cat(3, X.U{mu}, res{1});
0323
0324         tmp = permute( U, [1, 2, 4, 3] );
0325         U = reshape( tmp, [size(tmp,1)*size(tmp,2), p*size(tmp,4)] );
0326
0327         [U,S,V] = svd( U, 'econ' );
0328         s = trunc_singular( diag(S), opts.tol, true, opts.maxrank );
0329         U = U(:,1:s);
0330         X.U{mu} = reshape( U, [sz(1), sz(2), s] );
0331         W = S(1:s,1:s)*V(:,1:s)';
0332         W = reshape( W, [s, p, size(tmp,4)]);
0333         W = permute( W, [1, 3, 2]);
0334         for k = 1:p
0335             C{k} = tensorprod_ttemps( right, W(:,:,k), 1);
0336         end
0337         X.U{mu+1} = C{1};
0338
0339         if opts.verbose
0340             disp( ['Augmented system of size (', num2str( [sz(1)*p, sz(2)*sz(3)]), '). Cut-off tol: ', num2str(opts.tol) ])
0341             disp( sprintf( 'Number of SVs: %i. Truncated to: %i => %g %%', length(diag(S)), s, s/length(diag(S))*100))
0342             disp(' ')
0343         end
0344
0345     end
0346     %residuums = [residuums, norm(apply(A,X) - lambda(end)*X)];
0347
0348     fprintf(1, '---------------    finshed sweep.    ---------------\n')
0349     disp(['Current residuum: ', num2str(residuums(:,2*(i-1)+2).')])
0350     disp(' ')
0351     disp(' ')
0352 end
0353
0354 evs = zeros(p,1);
0355 % Evaluate rayleigh quotient for the p TT/MPS tensors
0356 for i=1:p
0357     evec = X;
0358     evec.U{d} = C{i};
0359     evs(i) = innerprod( evec, apply(A, evec));%/ innerprod( evec, evec);
0360 end
0361 evalue = [evalue, evs];
0362 end
0363
0364 function [left, right] = Afun_prepare( A, x, idx )
0365     y = A.apply(x);
0366     if idx == 1
0367         right = innerprod( x, y, 'RL', idx+1 );
0368         left = [];
0369     elseif idx == x.order
0370         left = innerprod( x, y, 'LR', idx-1 );
0371         right = [];
0372     else
0373         left = innerprod( x, y, 'LR', idx-1 );
0374         right = innerprod( x, y, 'RL', idx+1 );
0375     end
0376 end
0377
0378 function res = Afun_block_optim( A, U, sz, left, right, mu )
0379
0380
0381     p = size(U, 2);
0382
0383     V = reshape( U, [sz, p] );
0384     V = A.apply( V, mu );
0385
0386     if mu == 1
0387         tmp = reshape( permute( V, [3 1 2 4] ), [size(V, 3), sz(1)*sz(2)*p]);
0388         tmp = right * tmp;
0389         tmp = reshape( tmp, [size(right, 1), sz(1), sz(2), p]);
0390         tmp = ipermute( tmp, [3 1 2 4] );
0391     elseif mu == A.order
0392         tmp = reshape( V, [size(V,1), sz(2)*sz(3)*p]);
0393         tmp = left * tmp;
0394         tmp = reshape( tmp, [size(left, 1), sz(2), sz(3), p]);
0395     else
0396         tmp = reshape( permute( V, [3 1 2 4] ), [size(V, 3), size(V, 1)*sz(2)*p]);
0397         tmp = right * tmp;
0398         tmp = reshape( tmp, [size(right, 1), size(V, 1), sz(2), p]);
0399         tmp = ipermute( tmp, [3 1 2 4] );
0400
0401         tmp = reshape( tmp, [size(V, 1), sz(2)*sz(3)*p]);
0402         tmp = left * tmp;
0403         tmp = reshape( tmp, [size(left, 1), sz(2), sz(3), p]);
0404
0405     end
0406
0407     res = reshape( tmp, [prod(sz), p] );
0408 end
0409
0410
0411 function res = apply_local_precond( U, sz, expB)
0412
0413     p = size(U, 2);
0414
0415     x = reshape( U, [sz, p] );
0416     res = zeros( [sz, p] );
0417
0418     for i = 1:size( expB, 2)
0419         tmp = reshape( x, [sz(1), sz(2)*sz(3)*p] );
0420         tmp = reshape( expB{1,i}*tmp, [sz(1), sz(2), sz(3), p] );
0421
0422         tmp = reshape( permute( tmp, [2 1 3 4] ), [sz(2), sz(1)*sz(3)*p] );
0423         tmp = ipermute( reshape( expB{2,i}*tmp, [sz(2), sz(1), sz(3), p] ), [2 1 3 4] );
0424
0425         tmp = reshape( permute( tmp, [3 1 2 4] ), [sz(3), sz(1)*sz(2)*p] );
0426         tmp = ipermute( reshape( expB{3,i}*tmp, [sz(3), sz(1), sz(2), p] ), [3 1 2 4] );
0427
0428         res = res + tmp;
0429     end
0430     res = reshape( res, [prod(sz), p] );
0431
0432 end
0433
0434 function res = apply_global_precond( x, res_sz, expB)
0435
0436     res = zeros( res_sz );
0437
0438     for i = 1:size( expB, 2)
0439         tmp = reshape( x, [res_sz(1), prod(res_sz(2:5))] );
0440         tmp = reshape( expB{1,i}*tmp, res_sz );
0441
0442         tmp = reshape( permute( tmp, [2 1 3 4 5] ), [res_sz(2), res_sz(1)*prod(res_sz(3:5))] );
0443         tmp = ipermute( reshape( expB{2,i}*tmp, [res_sz(2), res_sz(1), res_sz(3:5)] ), [2 1 3 4 5] );
0444
0445         tmp = reshape( permute( tmp, [3 1 2 4 5] ), [res_sz(3), prod(res_sz([1:2,4:5]))] );
0446         tmp = ipermute( reshape( expB{3,i}*tmp, [res_sz(3), res_sz([1:2,4:5])] ), [3 1 2 4 5] );
0447
0448         tmp = reshape( permute( tmp, [4 1 2 3 5] ), [res_sz(4), prod(res_sz([1:3,5]))] );
0449         tmp = ipermute( reshape( expB{4,i}*tmp, [res_sz(4), res_sz([1:3,5])] ), [4 1 2 3 5] );
0450
0451         res = res + tmp;
0452     end
0453
0454 end```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005