contract

CONTRACT Contraction of two TT/MPS tensors with inner TT/MPS operator.
   Z = CONTRACT(A,X,IDX) contracts all coress of the two TT/MPS tensors X and Y=apply(A*X) except 
   core IDX. Result Z is a matrix of size 
       [ X.rank(IDX)*A.size_col(IDX)*X.rank(IDX+1), X.rank(IDX)*A.size_row(IDX),X.rank(IDX+1) ].

   See also APPLY.

0001 function res = contract( A, X, mu )
0002     %CONTRACT Contraction of two TT/MPS tensors with inner TT/MPS operator.
0003     %   Z = CONTRACT(A,X,IDX) contracts all coress of the two TT/MPS tensors X and Y=apply(A*X) except
0004     %   core IDX. Result Z is a matrix of size
0005     %       [ X.rank(IDX)*A.size_col(IDX)*X.rank(IDX+1), X.rank(IDX)*A.size_row(IDX),X.rank(IDX+1) ].
0006     %
0007     %   See also APPLY.
0008 
0009     %   TTeMPS Toolbox.
0010     %   Michael Steinlechner, 2013-2016
0011     %   Questions and contact: michael.steinlechner@epfl.ch
0012     %   BSD 2-clause license, see LICENSE.txt
0013 
0014     V = cell(1, A.order);
0015 
0016     Y = A.apply(X);
0017     if mu > 1
0018         left = innerprod( X, Y, 'LR', mu-1 );
0019         left = reshape( left, [X.rank(mu), A.rank(mu), X.rank(mu)] );
0020         left = reshape( permute(left, [1 3 2]), [X.rank(mu)*X.rank(mu), A.rank(mu)] );
0021         res = left * reshape(A.U{mu}, [A.rank(mu), A.size_col(mu)*A.size_row(mu)*A.rank(mu+1)] );
0022     else
0023         res = A.U{1};
0024     end
0025     if mu < X.order
0026         right = innerprod( X, Y, 'RL', mu+1 );
0027         right = reshape( right, [X.rank(mu+1), A.rank(mu+1), X.rank(mu+1)] );
0028         right = matricize( right, 2);
0029         res = reshape( res, [ X.rank(mu)*X.rank(mu)*A.size_col(mu)*A.size_row(mu), A.rank(mu+1) ]);
0030         res = res * right; 
0031     end
0032 
0033     res = reshape( res, [X.rank(mu), X.rank(mu), A.size_col(mu), A.size_row(mu), X.rank(mu+1), X.rank(mu+1) ] );
0034     res = permute( res, [1 3 5 2 4 6] );
0035     res = reshape( res, [X.rank(mu)*A.size_col(mu)*X.rank(mu+1), X.rank(mu)*A.size_row(mu)*X.rank(mu+1) ] );
0036 
0037 end
0038