parameterdependent

0001 classdef parameterdependent
0002     % Operator from a parameter dependent heat equation,
0003     % the so-called cookie problem.
0004     % Supported number of cookies: 4 or 9.
0005     %
0006     % See the following paper, Section 5.3, for a description:
0007     %
0008     % D. Kressner, M. Steinlechner, and B. Vandereycken.
0009     % Preconditioned low-rank Riemannian optimization for
0010     % linear systems with tensor product structure.
0011     % Technical report, July 2015. Revised February 2016.
0012     % To appear in SIAM J. Sci. Comput.
0013     %
0014 
0015     %   TTeMPS Toolbox.
0016     %   Michael Steinlechner, 2013-2016
0017     %   Questions and contact: michael.steinlechner@epfl.ch
0018     %   BSD 2-clause license, see LICENSE.txt
0019 
0020     
0021     properties( SetAccess = protected, GetAccess = public )
0022 
0023         A           % cell storing the galerkin matrices A0,...Ad
0024         param
0025         rank
0026         order
0027         Lapl
0028         size_row
0029         size_col
0030     end
0031 
0032     methods( Access = public )
0033 
0034         function A = parameterdependent( n, numcookie )
0035             
0036             % only one matrix passed
0037             A.param = transpose(linspace(1,10,n));
0038             
0039             if numcookie == 4
0040                 l = load('four_cookies.mat');
0041                 A.A{1} = l.A0;
0042                 A.A{2} = l.A1;
0043                 A.A{3} = l.A2;
0044                 A.A{4} = l.A3;
0045                 A.A{5} = l.A4;
0046             
0047                 A.rank = 5;
0048                 A.order = 5;
0049                 A.size_row = [1169,n,n,n,n];
0050                 A.size_col = [1169,n,n,n,n];
0051                 
0052                 A.Lapl = parameter_to_TTeMPS_op_laplace( A );
0053                 
0054             elseif numcookie == 9
0055                 l = load('nine_cookies.mat');
0056                 A.A{1} = l.A0;
0057                 A.A{2} = l.A1;
0058                 A.A{3} = l.A2;
0059                 A.A{4} = l.A3;
0060                 A.A{5} = l.A4;
0061                 A.A{6} = l.A5;
0062                 A.A{7} = l.A6;
0063                 A.A{8} = l.A7;
0064                 A.A{9} = l.A8;
0065                 A.A{10} = l.A9;
0066             
0067                 A.rank = 10;
0068                 A.order = 10;
0069                 A.size_row = [2796,n,n,n,n];
0070                 A.size_col = [2796,n,n,n,n];
0071                 
0072                 A.Lapl = parameter_to_TTeMPS_op_laplace( A );
0073             else
0074                 error('Only 4 or 9 cookies supported atm!')
0075             end
0076     
0077         end
0078         
0079         % Apply member function copied over from TTeMPS_op_laplace
0080         function y = apply( A, x, idx )
0081             %APPLY Application of TT/MPS parameter-dependent operator to a TT/MPS tensor
0082             %   Y = APPLY(A, X) applies the TT/MPS Laplace operator A to the TT/MPS tensor X.
0083             %
0084             %   Y = APPLY(A, X, idx) is the application of A but only in mode idx.
0085             %       note that in this case, X is assumed to be a standard matlab array and
0086             %       not a TTeMPS tensor.
0087             %
0088             %   In both cases, X can come from a block-TT format, that is, with a four-dimensional core instead.
0089             %
0090             if ~exist( 'idx', 'var' )
0091                 y = apply(A.Lapl, x );
0092             else
0093                 y = apply(A.Lapl, x, idx);
0094             end
0095         end
0096         
0097         function B = parameter_to_TTeMPS_op_laplace( A )
0098             if A.order == 5
0099                 a_0 = sparse( 1, 1, 1, 5, 1 );
0100                 a_1 = sparse( 2, 1, 1, 5, 1 );
0101                 a_2 = sparse( 3, 1, 1, 5, 1 );
0102                 a_3 = sparse( 4, 1, 1, 5, 1 );
0103                 a_4 = sparse( 5, 1, 1, 5, 1 );
0104                 
0105                 U = cell( 1, 5 );
0106                 U{1} = kron( A.A{1}, a_0 ) + kron( A.A{2}, a_1 ) ...
0107                          + kron( A.A{3}, a_2 ) + kron( A.A{4}, a_3 ) ...
0108                          + kron( A.A{5}, a_4 ); 
0109                    
0110                 n = length(A.param)      
0111                 D = spdiags(A.param,0,n,n);
0112                 E = speye(n);
0113                 
0114                 b_0 = sparse( 1 , 1, 1, 25, 1 );
0115                 b_1 = sparse( 7 , 1, 1, 25, 1 );
0116                 b_2 = sparse( 13, 1, 1, 25, 1 );
0117                 b_3 = sparse( 19, 1, 1, 25, 1 );
0118                 b_4 = sparse( 25, 1, 1, 25, 1 );
0119                 
0120                 U{2} = kron( E, b_0 ) + kron( D, b_1 ) ...
0121                        + kron( E, b_2 ) + kron( E, b_3 ) ...
0122                        + kron( E, b_4 );
0123                 U{3} = kron( E, b_0 ) + kron( E, b_1 ) ...
0124                        + kron( D, b_2 ) + kron( E, b_3 ) ...
0125                        + kron( E, b_4 );
0126                 U{4} = kron( E, b_0 ) + kron( E, b_1 ) ...
0127                        + kron( E, b_2 ) + kron( D, b_3 ) ...
0128                        + kron( E, b_4 );
0129                 U{5} = kron( E, a_0 ) + kron( E, a_1 ) ...
0130                        + kron( E, a_2 ) + kron( E, a_3 ) ...
0131                        + kron( D, a_4 );      
0132                        
0133                 B = TTeMPS_op_laplace( U );
0134                 B.rank = [1 5 5 5 5 1];
0135                 B.size_col = B.size_row;
0136             elseif A.order == 10
0137                 a_0 = sparse( 1,  1, 1, 10, 1 );
0138                 a_1 = sparse( 2,  1, 1, 10, 1 );
0139                 a_2 = sparse( 3,  1, 1, 10, 1 );
0140                 a_3 = sparse( 4,  1, 1, 10, 1 );
0141                 a_4 = sparse( 5,  1, 1, 10, 1 );
0142                 a_5 = sparse( 6,  1, 1, 10, 1 );
0143                 a_6 = sparse( 7,  1, 1, 10, 1 );
0144                 a_7 = sparse( 8,  1, 1, 10, 1 );
0145                 a_8 = sparse( 9,  1, 1, 10, 1 );
0146                 a_9 = sparse( 10, 1, 1, 10, 1 );
0147                 
0148                 U = cell( 1, 10 );
0149                 U{1} = kron( A.A{1}, a_0 ) + kron( A.A{2}, a_1 ) ...
0150                          + kron( A.A{3}, a_2 ) + kron( A.A{4}, a_3 ) ...
0151                          + kron( A.A{5}, a_4 ) + kron( A.A{6}, a_5 ) ...
0152                          + kron( A.A{7}, a_6 ) + kron( A.A{8}, a_7 ) ...
0153                          + kron( A.A{9}, a_8 ) + kron( A.A{10}, a_9 );
0154                    
0155                 n = length(A.param)      
0156                 D = spdiags(A.param,0,n,n);
0157                 E = speye(n);
0158                 
0159                 b_0 = sparse(   1, 1, 1, 100, 1 );
0160                 b_1 = sparse(  12, 1, 1, 100, 1 );
0161                 b_2 = sparse(  23, 1, 1, 100, 1 );
0162                 b_3 = sparse(  34, 1, 1, 100, 1 );
0163                 b_4 = sparse(  45, 1, 1, 100, 1 );
0164                 b_5 = sparse(  56, 1, 1, 100, 1 );
0165                 b_6 = sparse(  67, 1, 1, 100, 1 );
0166                 b_7 = sparse(  78, 1, 1, 100, 1 );
0167                 b_8 = sparse(  89, 1, 1, 100, 1 );
0168                 b_9 = sparse( 100, 1, 1, 100, 1 );
0169 
0170                 U{2} = kron( E, b_0 ) + kron( D, b_1 ) ...
0171                          + kron( E, b_2 ) + kron( E, b_3 ) ...
0172                          + kron( E, b_4 ) + kron( E, b_5 ) ...
0173                          + kron( E, b_6 ) + kron( E, b_7 ) ...
0174                          + kron( E, b_8 ) + kron( E, b_9 );
0175                 U{3} = kron( E, b_0 ) + kron( E, b_1 ) ...
0176                          + kron( D, b_2 ) + kron( E, b_3 ) ...
0177                          + kron( E, b_4 ) + kron( E, b_5 ) ...
0178                          + kron( E, b_6 ) + kron( E, b_7 ) ...
0179                          + kron( E, b_8 ) + kron( E, b_9 );
0180                 U{4} = kron( E, b_0 ) + kron( E, b_1 ) ...
0181                          + kron( E, b_2 ) + kron( D, b_3 ) ...
0182                          + kron( E, b_4 ) + kron( E, b_5 ) ...
0183                          + kron( E, b_6 ) + kron( E, b_7 ) ...
0184                          + kron( E, b_8 ) + kron( E, b_9 );
0185                 U{5} = kron( E, b_0 ) + kron( E, b_1 ) ...
0186                          + kron( E, b_2 ) + kron( E, b_3 ) ...
0187                          + kron( D, b_4 ) + kron( E, b_5 ) ...
0188                          + kron( E, b_6 ) + kron( E, b_7 ) ...
0189                          + kron( E, b_8 ) + kron( E, b_9 );
0190                 U{6} = kron( E, b_0 ) + kron( E, b_1 ) ...
0191                          + kron( E, b_2 ) + kron( E, b_3 ) ...
0192                          + kron( E, b_4 ) + kron( D, b_5 ) ...
0193                          + kron( E, b_6 ) + kron( E, b_7 ) ...
0194                          + kron( E, b_8 ) + kron( E, b_9 );
0195                 U{7} = kron( E, b_0 ) + kron( E, b_1 ) ...
0196                          + kron( E, b_2 ) + kron( E, b_3 ) ...
0197                          + kron( E, b_4 ) + kron( E, b_5 ) ...
0198                          + kron( D, b_6 ) + kron( E, b_7 ) ...
0199                          + kron( E, b_8 ) + kron( E, b_9 );
0200                 U{8} = kron( E, b_0 ) + kron( E, b_1 ) ...
0201                          + kron( E, b_2 ) + kron( E, b_3 ) ...
0202                          + kron( E, b_4 ) + kron( E, b_5 ) ...
0203                          + kron( E, b_6 ) + kron( D, b_7 ) ...
0204                          + kron( E, b_8 ) + kron( E, b_9 );
0205                 U{9} = kron( E, b_0 ) + kron( E, b_1 ) ...
0206                          + kron( E, b_2 ) + kron( E, b_3 ) ...
0207                          + kron( E, b_4 ) + kron( E, b_5 ) ...
0208                          + kron( E, b_6 ) + kron( E, b_7 ) ...
0209                          + kron( D, b_8 ) + kron( E, b_9 );
0210                 U{10} = kron( E, a_0 ) + kron( E, a_1 ) ...
0211                          + kron( E, a_2 ) + kron( E, a_3 ) ...
0212                          + kron( E, a_4 ) + kron( E, a_5 ) ...
0213                          + kron( E, a_6 ) + kron( E, a_7 ) ...
0214                          + kron( E, a_8 ) + kron( D, a_9 );
0215                 B = TTeMPS_op_laplace( U );
0216                 B.rank = [1 10*ones(1,9) 1];
0217                 B.size_col = B.size_row;
0218             else
0219                 error('Wrong number of cookies!');
0220             end
0221 
0222         end
0223     end         
0224         
0225 end