1 function scm_offline_data =
scm_offline(model, detailed_data, M_train, D_train)
2 %scm_offline_data =
scm_offline(model, detailed_data, M_train, D_train)
4 %
function computing the offline quantities
for the SCM.
9 %- M_train(optional): A training set for the SCM-greedy can be given as Input
10 %- D_train(optional): A set which influences the constraints in the online-Phase.
12 % fields of the model can be set in the model to configurate the
13 % offline-phase of the SCM. If they are not set standard values are taken:
14 %- model.scm_M_alpha : natural number influencing the number of constraints used out of the set C in the online-phase.
15 %- model.scm_M_plus : natural number influencing the number of constraints used out of the set D in the online-phase.
16 %- model.scm_eps_tol : stopping condition for the included greedy-algorithmn. Stops when the error indicator is lower than this value.
17 %- model.scm_size_C : stopping condition for the included greedy-algorithmn. Stops when the the desired size of C is reached.
18 % (in general the size of C should scale up with the complexety of the problem (Q and size or parameter space)
19 %- model.scm_desired_constant: 1 for the coercive case (i.e. SCM produces lower and upper bounds to the coercivity constant)
20 % 2 for the inf-sup stable case (i.e. SCM produces lower and upper bounds to the inf-sup constant)
22 % Output (the structure scm_offline_data containing the following fields):
23 %- sigma_plus : the upper values of the bounding box, required in the online-phase
24 %- sigma_minus : the lower values of the bounding box, required in the online-phase
25 %- D : the set of parameters for the constraints J(mu,y) >= 0, required in the online-phase
26 %- Theta_D : the coefficient vectors of the parameters in D, required in the online-phase
27 %- C : the set of parameters for the constraints J(mu,y) >= alpha(mu), required in the online-phase
28 %- Theta_C : the coefficient vectors of the parameters in C, required in the online-phase
29 %- P_train : the parameter training set (for additional information)
30 %- alpha_C : constains the exact stability constant (coercivity or inf-sup) for every parameter in C, required in the online-phase
31 %- ystern_C : constains the vector y*(mu) for every mu in C, required in the online-phase
32 %- max_err_seq : the error sequence of the greedy algorithmn (for additional information)
33 %- exit_on_tol : 1 if the greedy stopped because of reaching the tolerance, 0 else
34 %- exit_on_emptyTrain : 1 if the greedy stopped because all the parameters in the trainings set are used in C, 0 else
35 %- exit_on_size_C : 1 if the greedy stopped because the desired size of C is reached, 0 else
38 % Dominik Garmatter 20.09 2012
41 disp(
'starting offline-phase of SCM')
43 %% initialise parameterindependant quantities
45 total_time_start = tic;
46 %
if there are no scm-configuration fields in the model -> take standard values
47 if ~isfield(model,
'scm_M_alpha')
48 model.scm_M_alpha = 50;
50 if ~isfield(model, 'scm_eps_tol')
51 model.scm_eps_tol = 1e-10;
53 if ~isfield(model, 'scm_size_C')
54 model.scm_size_C = 100;
56 if ~isfield(model, 'scm_desired_constant')
57 model.scm_desired_constant = 1;
59 if ~isfield(model, 'scm_M_plus')
60 model.scm_M_plus = 200;
62 % write the configuration fields in a simpler form
63 M_alpha = model.scm_M_alpha;
64 eps_tol = model.scm_eps_tol;
65 size_C = model.scm_size_C;
66 desired_constant = model.scm_desired_constant;
67 M_plus = model.scm_M_plus;
68 % save them in the scm_offline_data for later online-phases (i.e. in the rb_simulation)
69 % they also present alot of onformation about the scm_offline_data
70 scm_offline_data.M_alpha = M_alpha;
71 scm_offline_data.M_plus = M_plus;
72 scm_offline_data.eps_tol = eps_tol;
73 scm_offline_data.size_C = size_C;
74 scm_offline_data.desired_constant = desired_constant;
77 % set the options for the linear program in the online-phase, i.e. use the active-set algorithm
78 scm_offline_data.Options = optimset('Display','off','LargeScale','off');
80 if strcmp(model.rb_problem_type, 'lin_evol') == 1
81 % get the components and other quantities out of the detailed_data
82 Kreg = model.get_inner_product_matrix(detailed_data);
83 M = model.nt+1; % number of timesteps
84 model.Q_E = size(detailed_data.L_E_comp(:),1); % number of components of L_E - required for the evaluation in scm_shifted_evaluation
85 components = [{speye(size(detailed_data.L_I_comp{1}))}; detailed_data.L_E_comp; detailed_data.L_I_comp];
86 Q = size(components,1); % number of components
87 %evaluation functions
for the components required in scm_shifted evaluation
88 model.scm_component_evaluation = @(x,q) components{q}
'*Kreg*x;
89 model.scm_transposed_component_evaluation = @(x,q) (x'*components{q}
'*Kreg)';
90 model.scm_shiftmatrix_mult = @(x) Kreg*x; %
for the shifts and innerprodmatrix evaluation
91 if desired_constant == 2
92 Q = Q^2; % in the inf-sup-
case this is the number of components,
else its just Q.
93 model.scm_component_evaluation = @(x,q) components{q}*x;
94 model.scm_transposed_component_evaluation = @(x,q) (x
'*components{q})';
97 components{q+1} = -components{q+1}; % L_E has a -1 in the PG-BLF-Matrix
100 elseif strcmp(model.rb_problem_type,
'lin_stat') == 1
101 %
get the components and other quantities out of the detailed_data
102 Q = size(detailed_data.A_comp(:),1);
103 components = detailed_data.A_comp;
104 if isfield(detailed_data,
'df_info') == 1 % neglect the dirichlet-values in FEM-cases
105 K = detailed_data.df_info.h10_inner_product_matrix;
106 K(:, detailed_data.df_info.dirichlet_gids) = 0;
107 K(detailed_data.df_info.dirichlet_gids, :) = 0;
109 components{q}(detailed_data.df_info.dirichlet_gids, :) = 0;
110 components{q}(:, detailed_data.df_info.dirichlet_gids) = 0;
113 K = model.get_inner_product_matrix(detailed_data);
115 Kreg = model.get_inner_product_matrix(detailed_data);
116 %evaluation functions
for the components required in scm_shifted evaluation
117 model.scm_component_evaluation = @(x,q) components{q}*x;
118 model.scm_transposed_component_evaluation = @(x,q) (x
'*components{q})';
119 model.scm_shiftmatrix_mult = @(x) K*x;
120 if desired_constant == 2
121 Q = Q^2; % number of components in inf-sup
case
122 model.scm_invers_innerprodmatrix_mult = @(x) Kreg\x;
124 M = 1; % only one timestep
126 error(
'rb_problem_type unknown!')
129 H_total = H*M; % total size of the system
131 warningstate_error = warning('error','MATLAB:eigs:NoEigsConverged'); % sets the state of this warning to error. This way you can tackle it with a try catch block and adjust the problem.
133 %% compute the SCM bounding box
134 % always determine the absolute of the largest magnitude. If its positiv it is sigma_plus, if its negativ it is sigma_minus.
135 % The shifted EWP then gives the corresponding smallest or largets value.
138 disp('computing SCM bounding box')
140 sigma_plus = zeros(Q,1);
141 sigma_minus = zeros(Q,1);
142 K_shift = eps; %shift used for eigs in case of all-zero component-matrices (ARPACK info-flag = -9 was the problem)
147 try %first try to solve the EWP with default tolerance eps
148 LM = eigs(@(x)scm_shifted_evaluation(model, x, i, K_shift, Q, M, H), H_total, kron(speye(M),Kreg), 1, 'LM') - K_shift;
150 if (strcmp(err.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err.identifier,'MATLAB:eigs:NoEigsConverged'))...
151 || (strcmp(err.message,'Error with ARPACK routine dnaupd: info = -8'))
152 % in these cases (convergence problems or heavy clustering of eigenvalues) change the options and try again. Can be timeconsuming.
153 if H_total < 600 % opts.p has to be smaller than H_total
154 opts.p = ceil(H_total/2);
162 disp('There was no solution to an eigenvalueproblem in scm offline-phase. I changed the options.');
164 LM = eigs(@(x)scm_shifted_evaluation(model, x, i, K_shift, Q, M, H), H_total, kron(speye(M),Kreg), 1, 'LM', opts) - K_shift;
165 no_solution = 0; %end when there is no error
166 catch err2 %when it still fails change the options and try again!
168 if (strcmp(err2.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err2.identifier,'MATLAB:eigs:NoEigsConverged'))...
169 || (strcmp(err2.message,'Error with ARPACK routine dnaupd: info = -8'))
170 opts.tol = 100 * opts.tol;
171 if H_total < 600 % opts.p has to be smaller than H_total
172 opts.p = opts.p + ceil(H_total/50);
176 opts.maxit = opts.maxit + 100;
178 if opts.tol > 1e-3 %if the options have changed several times and there is still no solution produce an error.
179 error(['I changed the options of an eigenvalueproblem in scm offline-phase up to opts.tol = ' mat2str(opts.tol) ', opts.p = ' mat2str(opts.p) ' and opts.maxit = ' mat2str(opts.maxit) ' but I still could not find a solution. Reconsider your problem']);
183 else %if there is another error rethrow
188 if abs(LM) < 1e-10 %compensating numerics in case of zero eigenvalues
194 try %first try to solve the EWP with default tolerance eps
195 sigma_minus(i) = eigs(@(x)scm_shifted_evaluation(model, x, i, -(LM+K_shift), Q, M, H), H_total, kron(speye(M),Kreg), 1, 'LM') + (LM+K_shift);
197 if (strcmp(err.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err.identifier,'MATLAB:eigs:NoEigsConverged'))...
198 || (strcmp(err.message,'Error with ARPACK routine dnaupd: info = -8'))
199 % in these cases (convergence problems or heavy clustering of eigenvalues) change the options and try again. Can be timeconsuming.
200 if H_total < 600 % opts.p has to be smaller than H_total
201 opts.p = ceil(H_total/2);
209 disp('There was no solution to an eigenvalueproblem in scm offline-phase. I changed the options.');
211 sigma_minus(i) = eigs(@(x)scm_shifted_evaluation(model, x, i, -(LM+K_shift), Q, M, H), H_total, kron(speye(M),Kreg), 1, 'LM', opts) + (LM+K_shift);
212 no_solution = 0; %end when there is no error
213 catch err2 %when it still fails change the options and try again!
215 if (strcmp(err2.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err2.identifier,'MATLAB:eigs:NoEigsConverged'))...
216 || (strcmp(err2.message,'Error with ARPACK routine dnaupd: info = -8'))
217 opts.tol = 100 * opts.tol;
218 if H_total < 600 % opts.p has to be smaller than H_total
219 opts.p = opts.p + ceil(H_total/50);
223 opts.maxit = opts.maxit + 100;
225 if opts.tol > 1e-3 %if the options have changed several times and there is still no solution produce an error.
226 error(['I changed the options of an eigenvalueproblem in scm offline-phase up to opts.tol = ' mat2str(opts.tol) ', opts.p = ' mat2str(opts.p) ' and opts.maxit = ' mat2str(opts.maxit) ' but I still could not find a solution. Reconsider your problem']);
230 else %if there is another error rethrow
237 try %first try to solve the EWP with default options
238 sigma_plus(i) = eigs(@(x)scm_shifted_evaluation(model, x, i, -(LM+K_shift), Q, M, H), H_total, kron(speye(M),Kreg), 1, 'LM') + (LM+K_shift);
240 if (strcmp(err.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err.identifier,'MATLAB:eigs:NoEigsConverged'))...
241 || (strcmp(err.message,'Error with ARPACK routine dnaupd: info = -8'))
242 % in these cases (convergence problems or heavy clustering of eigenvalues) change the options and try again. Can be timeconsuming.
243 if H_total < 600 % opts.p has to be smaller than H_total
244 opts.p = ceil(H_total/2);
252 disp('There was no solution to an eigenvalueproblem in scm offline-phase. I changed the options.');
254 sigma_plus(i) = eigs(@(x)scm_shifted_evaluation(model, x, i, -LM, Q, M, H), H_total, kron(speye(M),Kreg), 1, 'LM', opts) + LM;
255 no_solution = 0; %end when there is no error
256 catch err2 %when it still fails change the options and try again!
258 if (strcmp(err2.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err2.identifier,'MATLAB:eigs:NoEigsConverged'))...
259 || (strcmp(err2.message,'Error with ARPACK routine dnaupd: info = -8'))
260 opts.tol = 100 * opts.tol;
261 if H_total < 600 % opts.p has to be smaller than H_total
262 opts.p = opts.p + ceil(H_total/50);
266 opts.maxit = opts.maxit + 100;
268 if opts.tol > 1e-3 %if the options have changed several times and there is still no solution produce an error.
269 error(['I changed the options of an eigenvalueproblem in scm offline-phase up to opts.tol = ' mat2str(opts.tol) ', opts.p = ' mat2str(opts.p) ' and opts.maxit = ' mat2str(opts.maxit) ' but I still could not find a solution. Reconsider your problem']);
273 else %if there is another error rethrow
278 % again compensating numerics
279 if abs(sigma_plus(i)) < 1e-10
282 if abs(sigma_minus(i)) < 1e-10
290 scm_offline_data.sigma_plus = sigma_plus;
291 scm_offline_data.sigma_minus = sigma_minus;
293 %% initialising parameterdependant stuff
295 % the set D - use a huge randomset if D_train isn't specified as Input
299 build1 = zeros(size(model.mu_ranges,2),1);
300 build2 = zeros(size(model.mu_ranges,2),1);
301 for i = 1:size(model.mu_ranges,2)
302 build1(i,1) = model.mu_ranges{i}(1);
303 build2(i,1) = model.mu_ranges{i}(2);
305 D = unifrnd(repmat(build1,1,200000),repmat(build2,1,200000));
307 % save the coefficient vector
for all elements in D - required
for the constraints in the online-Phase
308 Theta_D = zeros(Q,size(D,2));
309 for l = 1:size(Theta_D,2)
310 model = set_mu(model, D(:,l));
311 Theta_D(:,l) = Theta(model);
313 % the SCM training set - use 2000 random parameters if M_train isn't specified as input
317 build1 = zeros(size(model.mu_ranges,2),1);
318 build2 = zeros(size(model.mu_ranges,2),1);
319 for i = 1:size(model.mu_ranges,2)
320 build1(i,1) = model.mu_ranges{i}(1);
321 build2(i,1) = model.mu_ranges{i}(2);
323 P_train = unifrnd(repmat(build1,1,2000),repmat(build2,1,2000));
325 scm_offline_data.D = D;
326 scm_offline_data.Theta_D = Theta_D;
327 scm_offline_data.P_train = P_train;
329 % use the first mu of the training set as starting mu
330 current_mu = P_train(:,1);
331 model = set_mu(model, current_mu);
333 % initialising loop quantities
341 exit_on_emptyTrain = 0;
345 if model.verbose > 20
346 disp(
'starting extension of the set C')
350 disp(['Detected maximum error indicator for mu = [',...
351 num2str(current_mu'),']']);
354 % preparation for the full EWPs and computation of the new element of Theta_C
355 if strcmp(model.rb_problem_type, 'lin_stat') == 1
356 current_Theta = Theta(model);
357 old_decomp_mode = model.decomp_mode;
358 model.decomp_mode = 2;
359 coeffs = model.operators(model,[]);
360 model.decomp_mode = old_decomp_mode;
361 Matrix = lincomb_sequence(components, coeffs);
362 model.scm_eval_Matrix = @(x) Matrix*x;
363 model.scm_eval_transposed_Matrix = @(x) (x'*Matrix)';
364 else %if strcmp(model.rb_problem_type, 'lin_evol') == 1
365 [current_Theta, coeffs] = Theta(model);
366 L_E = lincomb_sequence(components(2:model.Q_E+1), coeffs(2:model.Q_E+1));
367 L_I = lincomb_sequence(components(model.Q_E+2:end), coeffs(model.Q_E+2:end));
368 if desired_constant == 1
369 % additonal evaluations for the coercive case
370 model.scm_eval_L_E = @(x) L_E'*Kreg*x;
371 model.scm_eval_L_E_transposed = @(x) (x'*L_E'*Kreg)';
372 model.scm_eval_L_I = @(x) L_I'*Kreg*x;
373 model.scm_eval_L_I_transposed = @(x) (x'*L_I'*Kreg)';
374 elseif desired_constant == 2
375 % additional evaluations for the inf-sup case
376 model.scm_eval_L_Et_L_E = @(x) L_E'*Kreg*L_E*x;
377 model.scm_eval_L_It_L_I = @(x) L_I'*Kreg*L_I*x;
378 model.scm_eval_L_It_L_E = @(x) L_I'*Kreg*L_E*x;
379 model.scm_eval_L_Et_L_I = @(x) L_E'*Kreg*L_I*x;
380 model.scm_eval_K_L_E = @(x) Kreg*L_E*x;
381 model.scm_eval_L_Et_K = @(x) (x'*Kreg*L_E)';
384 % again compute the largest eigenvalue LM.
385 % if LM is positiv -> shift with -LM and then the largest eigenvalue
386 % will be the desired smallest eigenvalue (when shifted back with LM).
387 % if LM is negativ it is already the desired smallest eigenvalue.
389 try %first try to solve the EWP with default options
390 [EV_LM, LM] = eigs(@(x)scm_shifted_evaluation(model, x, 0, 0, Q, M, H), H_total, kron(speye(M),Kreg), 1, 'LM');
392 if (strcmp(err.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err.identifier,'MATLAB:eigs:NoEigsConverged'))...
393 || (strcmp(err.message,'Error with ARPACK routine dnaupd: info = -8'))
394 % in these cases (convergence problems or heavy clustering of eigenvalues) change the options and try again. Can be timeconsuming.
395 if H_total < 600 % opts.p has to be smaller than H_total
396 opts.p = ceil(H_total/2);
404 disp('There was no solution to an eigenvalueproblem in scm offline-phase. I changed the options.');
406 [EV_LM, LM] = eigs(@(x)scm_shifted_evaluation(model, x, 0, 0, Q, M, H), H_total, kron(speye(M),Kreg), 1, 'LM', opts);
407 no_solution = 0; %end when there is no error
408 catch err2 %when it still fails change the options and try again!
410 if (strcmp(err2.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err2.identifier,'MATLAB:eigs:NoEigsConverged'))...
411 || (strcmp(err2.message,'Error with ARPACK routine dnaupd: info = -8'))
412 opts.tol = 100 * opts.tol;
413 if H_total < 600 % opts.p has to be smaller than H_total
414 opts.p = opts.p + ceil(H_total/50);
418 opts.maxit = opts.maxit + 100;
420 if opts.tol > 1e-3 %if the options have changed several times and there is still no solution produce an error.
421 error(['I changed the options of an eigenvalueproblem in scm offline-phase up to opts.tol = ' mat2str(opts.tol) ', opts.p = ' mat2str(opts.p) ' and opts.maxit = ' mat2str(opts.maxit) ' but I still could not find a solution. Reconsider your problem']);
425 else %if there is another error rethrow
431 try %first try to solve the EWP with default options
432 [current_EV_in_X, EW] = eigs(@(x)scm_shifted_evaluation(model, x, 0, -LM, Q, M, H), H_total, kron(speye(M),Kreg), 1, 'LM');
434 if (strcmp(err.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err.identifier,'MATLAB:eigs:NoEigsConverged'))...
435 || (strcmp(err.message,'Error with ARPACK routine dnaupd: info = -8'))
436 % in these cases (convergence problems or heavy clustering of eigenvalues) change the options and try again. Can be timeconsuming.
437 if H_total < 600 % opts.p has to be smaller than H_total
438 opts.p = ceil(H_total/2);
446 disp('There was no solution to an eigenvalueproblem in scm offline-phase. I changed the options.');
448 [current_EV_in_X, EW] = eigs(@(x)scm_shifted_evaluation(model, x, 0, -LM, Q, M, H), H_total, kron(speye(M),Kreg), 1, 'LM', opts);
449 no_solution = 0; %end when there is no error
450 catch err2 %when it still fails change the options and try again!
452 if (strcmp(err2.identifier,'MATLAB:eigs:ARPACKroutineErrorMinus14')) || (strcmp(err2.identifier,'MATLAB:eigs:NoEigsConverged'))...
453 || (strcmp(err2.message,'Error with ARPACK routine dnaupd: info = -8'))
454 opts.tol = 100 * opts.tol;
455 if H_total < 600 % opts.p has to be smaller than H_total
456 opts.p = opts.p + ceil(H_total/50);
460 opts.maxit = opts.maxit + 100;
462 if opts.tol > 1e-3 %if the options have changed several times and there is still no solution produce an error.
463 error(['I changed the options of an eigenvalueproblem in scm offline-phase up to opts.tol = ' mat2str(opts.tol) ', opts.p = ' mat2str(opts.p) ' and opts.maxit = ' mat2str(opts.maxit) ' but I still could not find a solution. Reconsider your problem']);
467 else %if there is another error rethrow
471 current_alpha_C = EW + LM; % shifting back the eigenvalue
473 current_alpha_C = LM;
474 current_EV_in_X = EV_LM;
477 % backtransformation of the eigenvector from X to R^Q
478 current_ystern_C = zeros(Q,1);
480 y_q_build = scm_shifted_evaluation(model, current_EV_in_X, q, 0, Q, M, H);
481 current_ystern_C(q,1) = current_EV_in_X'*y_q_build(:);
484 % extend the offline data
485 alpha_C = [alpha_C; current_alpha_C];
486 ystern_C = [ystern_C, current_ystern_C];
488 Theta_C = [Theta_C, current_Theta];
490 scm_offline_data.C = C;
491 scm_offline_data.Theta_C = Theta_C;
492 scm_offline_data.alpha_C = alpha_C;
493 scm_offline_data.ystern_C = ystern_C;
495 % feasability test if the upper bound of the bounding box - scm_offline_data.sigma_plus - is a feasable solution for the linprog in
496 % the upcoming online-phase. If yes, it is chosen as starting vector for the linprog in the online-phase.
497 feasability_C = scm_offline_data.Theta_C'*scm_offline_data.sigma_plus >= scm_offline_data.alpha_C;
498 feasability_D = scm_offline_data.Theta_D'*scm_offline_data.sigma_plus >= 0;
499 if isequal(feasability_C,ones(size(feasability_C))) && isequal(feasability_D,ones(size(feasability_D)))
500 scm_offline_data.isfeasable = 1;
502 scm_offline_data.isfeasable = 0;
503 warning('OFFLINE:infeasable_data', 'offline_data can cause infeasable Problems. Check them!');
507 % compute the error indicator (Upper_bound(mu) - Lower_bound(mu))/Upper_bound(mu) for every mu in the training set
508 Quotient = zeros(size(P_train,2),1);
509 for i = 1 : size(P_train, 2)
513 model = set_mu(model, P_train(:,i));
514 test_Theta = Theta(model);
515 scm_results =
scm_online(P_train(:,i), test_Theta, scm_offline_data, M_alpha, M_plus, desired_constant);
516 Quotient(i,1) = (scm_results.Upper_Bound - scm_results.Lower_Bound)/scm_results.Upper_Bound;
522 [max_Quotient, ind] = max(Quotient); % take the worst error indicator
524 % extend the error sequence
525 max_err = [max_err, max_Quotient];
526 scm_offline_data.max_err_seq = max_err;
528 % check the stopping conditions and stop if one of them hits. Else
529 % choose the mu with the largest error indicator as new mu and delete
530 % it from the training set.
531 if max_Quotient < eps_tol
534 total_time = toc(total_time_start);
536 disp(['Offline_phase finished on e_tol = ' mat2str(eps_tol) ' in ' mat2str(total_time) 'seconds'])
538 elseif isempty(P_train) == 1
539 exit_on_emptyTrain = 1;
541 total_time = toc(total_time_start);
543 disp(['Offline_phase finished on empty Trainingset in ' mat2str(total_time) 'seconds'])
545 elseif size(C,2) > size_C
548 total_time = toc(total_time_start);
550 disp(['Offline_phase finished on desired size of C in ' mat2str(total_time) 'seconds'])
554 current_mu = P_train(:,ind);
555 model = set_mu(model, current_mu);
558 disp(['extending C to ' mat2str(size(C,2)+1) ' Elements!']);
563 warning(warningstate_error) % return the state of the warning to 'on' instead of 'error'
565 scm_offline_data.exit_on_tol = exit_on_tol;
566 scm_offline_data.exit_on_emptyTrain = exit_on_emptyTrain;
567 scm_offline_data.exit_on_size_C = exit_on_size_C;
570 function y = scm_shifted_evaluation(model, x, q, a, Q, M, H)
571 %y = scm_shifted_evaluation(model, x, q, a, Q, M, H)
573 % function can either evaluate a full matrix (q=0) or just the qth
574 % component (q=/=0). Q is the total number of components, M the number of
575 % timesteps (=1 in lin_stat case) and H is the number of DOFs. x is the
576 % vector that will evaluate the desired matrix. Which matrix gets evaluated
577 % is controlled by model.desired_constant and model.rb_problem_type,
578 % because those fields specify if you have a lin_stat or a lin_evol problem
579 % and which constant, i.e. which eigenvalueproblem you want to solve.
580 % It is also possible to shift the evaluated matrix with the inner product
581 % matrix (K or Kreg), i.e. y =A*x + K*x, with A being the matrix to be evaluated.
583 % Dominik Garmatter 03.09 2012
585 if a ~= 0 %if there is a shift, evaluate the shiftmatrix, i.e. the Id-matrix right here
587 y2((i-1)*H+1 : i*H) = model.scm_shiftmatrix_mult(x((i-1)*H+1 : i*H));
595 if model.scm_desired_constant == 1 %evaluation for the coercivity constant
596 if strcmp(model.rb_problem_type, 'lin_stat') == 1
598 y1 = model.scm_eval_Matrix(x);
599 y1_transposed = model.scm_eval_transposed_Matrix(x);
601 y1 = model.scm_component_evaluation(x,q);
602 y1_transposed = model.scm_transposed_component_evaluation(x,q);
605 if q == 0 %eval the combined matrices
606 y1(1:H) = model.scm_shiftmatrix_mult(x(1:H)); %first component is always Id*Kreg = Kreg
607 y1_transposed(1:H) = model.scm_shiftmatrix_mult(x(1:H)) + model.scm_eval_L_E_transposed(x(H+1:2*H));
609 y1((i-1)*H+1 : i*H) = model.scm_eval_L_E(x((i-2)*H+1 : (i-1)*H)) + model.scm_eval_L_I(x((i-1)*H+1 : i*H));
610 y1_transposed((i-1)*H+1 : i*H) = model.scm_eval_L_I_transposed(x((i-1)*H+1 : i*H)) + model.scm_eval_L_E_transposed(x(i*H+1 : (i+1)*H));
612 y1((M-1)*H+1 : M*H) = model.scm_eval_L_E(x((M-2)*H+1 : (M-1)*H)) + model.scm_eval_L_I(x((M-1)*H+1 : M*H));
613 y1_transposed((M-1)*H+1 : M*H) = model.scm_eval_L_I_transposed(x((M-1)*H+1 : M*H));
615 %evaluation of the various components of the symmetric part of the PG-BLF-Matrix
616 elseif q == 1 %Id component
617 y1(1:H) = model.scm_component_evaluation(x(1:H),q);
619 y1_transposed(1:H) = model.scm_transposed_component_evaluation(x(1:H),q);
620 y1_transposed(H+1:M*H) = 0;
622 elseif 2 <= q && q <= 1 + model.Q_E %L_E components with Q_E being the number of L_E components
624 y1_transposed(1:H) = model.scm_transposed_component_evaluation(x(H+1:2*H),q);
626 y1((i-1)*H+1 : i*H) = model.scm_component_evaluation(x((i-2)*H+1 : (i-1)*H),q);
627 y1_transposed((i-1)*H+1 : i*H) = model.scm_transposed_component_evaluation(x(i*H+1 : (i+1)*H),q);
629 y1((M-1)*H+1 : M*H) = model.scm_component_evaluation(x((M-2)*H+1 : (M-1)*H),q);
630 y1_transposed((M-1)*H+1 : M*H) = 0;
634 y1_transposed(1:H) = 0;
636 y1((i-1)*H+1 : i*H) = model.scm_component_evaluation(x((i-1)*H+1 : i*H),q);
637 y1_transposed((i-1)*H+1 : i*H) = model.scm_transposed_component_evaluation(x((i-1)*H+1 : i*H),q);
641 y1_symm = 0.5 * (y1 + y1_transposed);
642 y = y1_symm(:) + a*y2; %shift
645 elseif model.scm_desired_constant == 2 %evaluation for the inf-sup constant
646 % evaluation for the inf-sup-constant
647 Q_comp = sqrt(Q); % number of components
649 if q == 0 %evaluate the combined matrices
650 % no symmetric evaluation here, since B*K^-1*B' is always symmetric
651 if strcmp(model.rb_problem_type, 'lin_stat') == 1 %need an evaluation here, since i dont want to form the components with K^-1 themselves
652 x1 = model.scm_eval_transposed_Matrix(x);
653 x2 = model.scm_invers_innerprodmatrix_mult(x1);
654 y1 = model.scm_eval_Matrix(x2);
656 y1(1:H) = model.scm_shiftmatrix_mult(x(1:H)) + model.scm_eval_K_L_E(x(H+1:2*H));
657 y1(H+1:2*H) = model.scm_eval_L_Et_K(x(1:H)) + model.scm_eval_L_Et_L_E(x(H+1:2*H)) + model.scm_eval_L_It_L_I(x(H+1:2*H)) + ...
658 model.scm_eval_L_It_L_E(x(2*H+1:3*H));
660 y1((i-1)*H+1 : i*H) = model.scm_eval_L_Et_L_I(x((i-2)*H+1 : (i-1)*H)) + model.scm_eval_L_It_L_I(x((i-1)*H+1 : i*H)) + ...
661 model.scm_eval_L_Et_L_E(x((i-1)*H+1 : i*H)) + model.scm_eval_L_It_L_E(x(i*H+1 : (i+1)*H));
663 y1((M-1)*H+1:M*H) = model.scm_eval_L_Et_L_I(x((M-2)*H+1:(M-1)*H)) + model.scm_eval_L_Et_L_E(x((M-1)*H+1:M*H)) + model.scm_eval_L_It_L_I(x((M-1)*H+1:M*H));
665 y = y1(:) + a*y2; %shift
668 else %the various component evaluations
669 % symmetric evaluation in all non-zero non-Id components, since all of them can be non-symmetric
670 %
double for-loop to create the wanted evaluation comp1{i}*K^{-1}*comp2{j}
'
673 if q == (i-1)*Q_comp + j %counts from 1 to Q^2, if it matches the desired q the component evaluates, else it increases by 1
674 if strcmp(model.rb_problem_type, 'lin_stat
') == 1 %need an evaluation here, since i dont want to form the components with K^-1 themselves
675 x1 = model.scm_transposed_component_evaluation(x,j);
676 x2 = model.scm_invers_innerprodmatrix_mult(x1);
677 x3 = model.scm_component_evaluation(x2,i);
678 x1symm = model.scm_transposed_component_evaluation(x,i);
679 x2symm = model.scm_invers_innerprodmatrix_mult(x1symm);
680 x3symm = model.scm_component_evaluation(x2symm,j);
681 y1 = 0.5*(x3 + x3symm);
682 else %if strcmp(model.rb_problem_type, 'lin_evol
') == 1
683 %only IdL_E, L_EId, L_EL_I abd L_IL_E comp are evaluated symmetric -> rest is symm already
684 if j == 1 %if comp2 is the Id-comp
685 if i == 1 %if comp1 is also the Id-comp
686 y1(1:H) = model.scm_shiftmatrix_mult(x(1:H)); % no symmetric eval here since Id is symmetric :)
688 elseif 2 <= i && i <= 1 + model.Q_E %if comp1 is an L_E_comp
690 x1 = model.scm_shiftmatrix_mult(x(1:H));
691 ysymm_1(H+1:2*H) = model.scm_transposed_component_evaluation(x1,i);
692 ysymm_1(2*H+1:M*H) = 0;
693 x1symm = model.scm_component_evaluation(x(H+1:2*H),i);
694 ysymm_2(1:H) = model.scm_shiftmatrix_mult(x1symm);
695 ysymm_2(H+1:M*H) = 0;
696 y1 = 0.5*(ysymm_1 + ysymm_2);
697 else %if comp1 is an L_I_comp
699 end %end of comp1-differentiation
701 elseif 2 <= j && j <= 1 + model.Q_E %if comp2 is an L_E_comp
702 if i == 1 %if comp1 is the Id-comp
703 x1 = model.scm_component_evaluation(x(H+1:2*H),j);
704 ysymm_1(1:H) = model.scm_shiftmatrix_mult(x1);
705 ysymm_1(H+1:M*H) = 0;
707 x1symm = model.scm_shiftmatrix_mult(x(1:H));
708 ysymm_2(H+1:2*H) = model.scm_transposed_component_evaluation(x1symm,j);
709 ysymm_2(2*H+1:M*H) = 0;
710 elseif 2 <= i && i <= 1 + model.Q_E %if comp1 is also an L_E_comp
714 x1 = model.scm_component_evaluation(x((k-1)*H+1:k*H),j);
715 x2 = model.scm_shiftmatrix_mult(x1);
716 ysymm_1((k-1)*H+1:k*H) = model.scm_transposed_component_evaluation(x2,i);
717 x1symm = model.scm_component_evaluation(x((k-1)*H+1:k*H),i);
718 x2symm = model.scm_shiftmatrix_mult(x1symm);
719 ysymm_2((k-1)*H+1:k*H) = model.scm_transposed_component_evaluation(x2symm,j);
721 else %if comp1 is an L_I_comp
725 x1 = model.scm_component_evaluation(x((k+1)*H+1:(k+2)*H),j);
726 x2 = model.scm_shiftmatrix_mult(x1);
727 ysymm_1(k*H+1:(k+1)*H) = model.scm_transposed_component_evaluation(x2,i);
728 c1 = model.scm_component_evaluation(x(k*H+1:(k+1)*H),i);
729 c2 = model.scm_shiftmatrix_mult(c1);
730 ysymm_2((k+1)*H+1:(k+2)*H) = model.scm_transposed_component_evaluation(c2,j);
732 ysymm_1((M-1)*H+1:M*H) = 0;
733 end %end of comp1-differentiation
734 y1 = 0.5*(ysymm_1 + ysymm_2);
736 else %if comp2 is and L_I_comp
737 if i == 1 %if comp1 is the Id-comp
739 elseif 2 <= i && i <= 1 + model.Q_E %if comp1 is an L_E_comp
743 x1 = model.scm_component_evaluation(x(k*H+1:(k+1)*H),j);
744 x2 = model.scm_shiftmatrix_mult(x1);
745 ysymm_1((k+1)*H+1:(k+2)*H) = model.scm_transposed_component_evaluation(x2,i);
746 c1 = model.scm_component_evaluation(x((k+1)*H+1:(k+2)*H),i);
747 c2 = model.scm_shiftmatrix_mult(c1);
748 ysymm_2(k*H+1:(k+1)*H) = model.scm_transposed_component_evaluation(c2,j);
750 ysymm_2((M-1)*H+1:M*H) = 0;
751 y1 = 0.5*(ysymm_1 + ysymm_2);
752 else %if comp1 is also an L_I_comp
756 x1 = model.scm_component_evaluation(x((k-1)*H+1:k*H),j);
757 x2 = model.scm_shiftmatrix_mult(x1);
758 ysymm_1((k-1)*H+1:k*H) = model.scm_transposed_component_evaluation(x2,i);
759 x1symm = model.scm_component_evaluation(x((k-1)*H+1:k*H),i);
760 x2symm = model.scm_shiftmatrix_mult(x1symm);
761 ysymm_2((k-1)*H+1:k*H) = model.scm_transposed_component_evaluation(x2symm,j);
763 y1 = 0.5*(ysymm_1 + ysymm_2);
764 end %end of comp1-differentiation
765 end %end of comp2-differentiation
766 end %end of rb_problem_type-differentiation
767 end %end of check, if the q-th component is matched
768 end %end of inner for-loop
769 end %end of outer for-loop
770 y = y1(:) + a*y2; %shift
773 end %end of if desired_constant
775 function [current_Coefficients, coeffs] = Theta(model)
776 %[current_Coefficients, coeffs] = Theta(model)
778 % function which computes the coefficients of the current mu set in the
779 % model. Includes various, required cases.
781 old_decomp_mode = model.decomp_mode;
782 model.decomp_mode = 2;
783 if strcmp(model.rb_problem_type, 'lin_evol
') == 1
784 [L_I,L_E] = model.operators_ptr(model);
785 coeffs = [1;L_E;L_I];
786 if model.scm_desired_constant == 1
787 current_Coefficients = coeffs(:);
788 elseif model.scm_desired_constant == 2
789 Q = size(coeffs(:),1);
790 current_Coefficients = zeros(Q^2,1);
793 current_Coefficients((i-1)*Q+j) = coeffs(i)*coeffs(j);
798 elseif strcmp(model.rb_problem_type, 'lin_stat
') == 1
799 coeffs = model.operators(model);
800 if model.scm_desired_constant == 1;
801 current_Coefficients = coeffs(:); % we want a column vector
802 elseif model.scm_desired_constant == 2
803 Q = size(coeffs(:),1);
804 current_Coefficients = zeros(Q^2,1);
807 current_Coefficients((i-1)*Q+j) = coeffs(i)*coeffs(j);
813 error('rb_problem_type unknown!
')
815 model.decomp_mode = old_decomp_mode;
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
function scm_results = scm_online(mu, Theta_mu, scm_offline_data, M_alpha, M_plus, desired_constant)
scm_results = scm_online(mu, Theta_mu, scm_offline_data, M_alpha, M_plus, desired_constant) ...
function scm_offline_data = scm_offline(model, detailed_data, M_train, D_train)
scm_offline_data = scm_offline(model, detailed_data, M_train, D_train)