rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
scm_offline.m
Go to the documentation of this file.
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)
3 %
4 % function computing the offline quantities for the SCM.
5 %
6 % Input:
7 %- model : the model
8 %- detailed_data :
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.
11 %
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)
21 %
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
36 
37 
38 % Dominik Garmatter 20.09 2012
39 
40 if model.verbose > 20
41  disp('starting offline-phase of SCM')
42 end
43 %% initialise parameterindependant quantities
44 
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;
49 end
50 if ~isfield(model, 'scm_eps_tol')
51  model.scm_eps_tol = 1e-10;
52 end
53 if ~isfield(model, 'scm_size_C')
54  model.scm_size_C = 100;
55 end
56 if ~isfield(model, 'scm_desired_constant')
57  model.scm_desired_constant = 1;
58 end
59 if ~isfield(model, 'scm_M_plus')
60  model.scm_M_plus = 200;
61 end
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;
75 
76 
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');
79 
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})';
95  end
96  for q = 1:model.Q_E
97  components{q+1} = -components{q+1}; % L_E has a -1 in the PG-BLF-Matrix
98  end
99 
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;
108  for q = 1:Q
109  components{q}(detailed_data.df_info.dirichlet_gids, :) = 0;
110  components{q}(:, detailed_data.df_info.dirichlet_gids) = 0;
111  end
112  else
113  K = model.get_inner_product_matrix(detailed_data);
114  end
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;
123  end
124  M = 1; % only one timestep
125 else
126  error('rb_problem_type unknown!')
127 end
128 H = size(Kreg,1);
129 H_total = H*M; % total size of the system
130 
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.
132 
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.
136 
137 if model.verbose > 20
138  disp('computing SCM bounding box')
139 end
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)
143 for i = 1 : Q
144  if model.verbose > 20
145  fprintf('.');
146  end
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;
149  catch err
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);
155  else
156  opts.p = 300;
157  end
158  opts.maxit = 1500;
159  opts.tol = 1e-14;
160  no_solution = 1;
161  while no_solution
162  disp('There was no solution to an eigenvalueproblem in scm offline-phase. I changed the options.');
163  try
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!
167  no_solution = 1;
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);
173  else
174  opts.p = opts.p +50;
175  end
176  opts.maxit = opts.maxit + 100;
177  end
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']);
180  end
181  end
182  end
183  else %if there is another error rethrow
184  rethrow(err)
185  end
186  end
187 
188  if abs(LM) < 1e-10 %compensating numerics in case of zero eigenvalues
189  LM = 0;
190  end
191 
192  if LM >= 0
193  sigma_plus(i) = LM;
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);
196  catch err
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);
202  else
203  opts.p = 300;
204  end
205  opts.maxit = 1500;
206  opts.tol = 1e-14;
207  no_solution = 1;
208  while no_solution
209  disp('There was no solution to an eigenvalueproblem in scm offline-phase. I changed the options.');
210  try
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!
214  no_solution = 1;
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);
220  else
221  opts.p = opts.p +50;
222  end
223  opts.maxit = opts.maxit + 100;
224  end
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']);
227  end
228  end
229  end
230  else %if there is another error rethrow
231  rethrow(err)
232  end
233  end
234 
235  else %LM < 0
236  sigma_minus(i) = LM;
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);
239  catch err
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);
245  else
246  opts.p = 300;
247  end
248  opts.maxit = 1500;
249  opts.tol = 1e-14;
250  no_solution = 1;
251  while no_solution
252  disp('There was no solution to an eigenvalueproblem in scm offline-phase. I changed the options.');
253  try
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!
257  no_solution = 1;
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);
263  else
264  opts.p = opts.p +50;
265  end
266  opts.maxit = opts.maxit + 100;
267  end
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']);
270  end
271  end
272  end
273  else %if there is another error rethrow
274  rethrow(err)
275  end
276  end
277  end
278  % again compensating numerics
279  if abs(sigma_plus(i)) < 1e-10
280  sigma_plus(i) = 0;
281  end
282  if abs(sigma_minus(i)) < 1e-10
283  sigma_minus(i) = 0;
284  end
285 
286 end
287 if model.verbose > 20
288  fprintf('\n');
289 end
290 scm_offline_data.sigma_plus = sigma_plus;
291 scm_offline_data.sigma_minus = sigma_minus;
292 
293 %% initialising parameterdependant stuff
294 
295 % the set D - use a huge randomset if D_train isn't specified as Input
296 if nargin > 3
297  D = D_train;
298 else
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);
304  end
305  D = unifrnd(repmat(build1,1,200000),repmat(build2,1,200000));
306 end
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);
312 end
313 % the SCM training set - use 2000 random parameters if M_train isn't specified as input
314 if nargin > 2
315  P_train = M_train;
316 else
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);
322  end
323  P_train = unifrnd(repmat(build1,1,2000),repmat(build2,1,2000));
324 end
325 scm_offline_data.D = D;
326 scm_offline_data.Theta_D = Theta_D;
327 scm_offline_data.P_train = P_train;
328 
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);
332 
333 % initialising loop quantities
334 continue_loop = 1;
335 alpha_C = [];
336 ystern_C = [];
337 C = [];
338 Theta_C = [];
339 max_err = [];
340 exit_on_tol = 0;
341 exit_on_emptyTrain = 0;
342 exit_on_size_C = 0;
343 
344 %% while-loop
345 if model.verbose > 20
346  disp('starting extension of the set C')
347 end
348 while continue_loop
349  if model.verbose > 20
350  disp(['Detected maximum error indicator for mu = [',...
351  num2str(current_mu'),']']);
352  end
353 
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)';
382  end
383  end
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.
388 
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');
391  catch err
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);
397  else
398  opts.p = 300;
399  end
400  opts.maxit = 1500;
401  opts.tol = 1e-14;
402  no_solution = 1;
403  while no_solution
404  disp('There was no solution to an eigenvalueproblem in scm offline-phase. I changed the options.');
405  try
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!
409  no_solution = 1;
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);
415  else
416  opts.p = opts.p +50;
417  end
418  opts.maxit = opts.maxit + 100;
419  end
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']);
422  end
423  end
424  end
425  else %if there is another error rethrow
426  rethrow(err)
427  end
428  end
429 
430  if LM >= 0
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');
433  catch err
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);
439  else
440  opts.p = 300;
441  end
442  opts.maxit = 1500;
443  opts.tol = 1e-14;
444  no_solution = 1;
445  while no_solution
446  disp('There was no solution to an eigenvalueproblem in scm offline-phase. I changed the options.');
447  try
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!
451  no_solution = 1;
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);
457  else
458  opts.p = opts.p +50;
459  end
460  opts.maxit = opts.maxit + 100;
461  end
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']);
464  end
465  end
466  end
467  else %if there is another error rethrow
468  rethrow(err)
469  end
470  end
471  current_alpha_C = EW + LM; % shifting back the eigenvalue
472  else
473  current_alpha_C = LM;
474  current_EV_in_X = EV_LM;
475  end
476 
477 % backtransformation of the eigenvector from X to R^Q
478  current_ystern_C = zeros(Q,1);
479  for q = 1:Q
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(:);
482  end
483 
484  % extend the offline data
485  alpha_C = [alpha_C; current_alpha_C];
486  ystern_C = [ystern_C, current_ystern_C];
487  C = [C, current_mu];
488  Theta_C = [Theta_C, current_Theta];
489 
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;
494 
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;
501  else
502  scm_offline_data.isfeasable = 0;
503  warning('OFFLINE:infeasable_data', 'offline_data can cause infeasable Problems. Check them!');
504  keyboard;
505  end
506 
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)
510  if model.verbose > 20
511  fprintf('.');
512  end
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;
517  end
518  if model.verbose > 20
519  fprintf('\n');
520  end
521 
522  [max_Quotient, ind] = max(Quotient); % take the worst error indicator
523 
524  % extend the error sequence
525  max_err = [max_err, max_Quotient];
526  scm_offline_data.max_err_seq = max_err;
527 
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
532  exit_on_tol = 1;
533  continue_loop = 0;
534  total_time = toc(total_time_start);
535  if model.verbose > 20
536  disp(['Offline_phase finished on e_tol = ' mat2str(eps_tol) ' in ' mat2str(total_time) 'seconds'])
537  end
538  elseif isempty(P_train) == 1
539  exit_on_emptyTrain = 1;
540  continue_loop = 0;
541  total_time = toc(total_time_start);
542  if model.verbose > 20
543  disp(['Offline_phase finished on empty Trainingset in ' mat2str(total_time) 'seconds'])
544  end
545  elseif size(C,2) > size_C
546  exit_on_size_C = 1;
547  continue_loop = 0;
548  total_time = toc(total_time_start);
549  if model.verbose > 20
550  disp(['Offline_phase finished on desired size of C in ' mat2str(total_time) 'seconds'])
551  end
552  else
553  continue_loop = 1;
554  current_mu = P_train(:,ind);
555  model = set_mu(model, current_mu);
556  P_train(:,ind) = [];
557  if model.verbose > 20
558  disp(['extending C to ' mat2str(size(C,2)+1) ' Elements!']);
559  end
560  end
561 end % of while-loop
562 
563 warning(warningstate_error) % return the state of the warning to 'on' instead of 'error'
564 
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;
568 
569 
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)
572 %
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.
582 
583 % Dominik Garmatter 03.09 2012
584 
585 if a ~= 0 %if there is a shift, evaluate the shiftmatrix, i.e. the Id-matrix right here
586  for i = 1:M
587  y2((i-1)*H+1 : i*H) = model.scm_shiftmatrix_mult(x((i-1)*H+1 : i*H));
588  end
589  y2 = y2(:);
590 else
591  y2 = zeros(size(x));
592  y2 = y2(:);
593 end
594 
595 if model.scm_desired_constant == 1 %evaluation for the coercivity constant
596  if strcmp(model.rb_problem_type, 'lin_stat') == 1
597  if q == 0
598  y1 = model.scm_eval_Matrix(x);
599  y1_transposed = model.scm_eval_transposed_Matrix(x);
600  else
601  y1 = model.scm_component_evaluation(x,q);
602  y1_transposed = model.scm_transposed_component_evaluation(x,q);
603  end
604  else %lin_evol case
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));
608  for i = 2:M-1
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));
611  end
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));
614 
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);
618  y1(H+1:M*H) = 0;
619  y1_transposed(1:H) = model.scm_transposed_component_evaluation(x(1:H),q);
620  y1_transposed(H+1:M*H) = 0;
621 
622  elseif 2 <= q && q <= 1 + model.Q_E %L_E components with Q_E being the number of L_E components
623  y1(1:H) = 0;
624  y1_transposed(1:H) = model.scm_transposed_component_evaluation(x(H+1:2*H),q);
625  for i = 2:M-1
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);
628  end
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;
631 
632  else %L_I components
633  y1(1:H) = 0;
634  y1_transposed(1:H) = 0;
635  for i = 2:M
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);
638  end
639  end
640  end
641  y1_symm = 0.5 * (y1 + y1_transposed);
642  y = y1_symm(:) + a*y2; %shift
643  y = y(:);
644 
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
648 
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);
655  else %lin_evol case
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));
659  for i = 3:M-1
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));
662  end
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));
664  end
665  y = y1(:) + a*y2; %shift
666  y = y(:);
667 
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}'
671  for i = 1:Q_comp
672  for j = 1:Q_comp
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 :)
687  y1(H+1:M*H) = 0;
688  elseif 2 <= i && i <= 1 + model.Q_E %if comp1 is an L_E_comp
689  ysymm_1(1:H) = 0;
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
698  y1(1:M*H) = 0;
699  end %end of comp1-differentiation
700 
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;
706  ysymm_2(1: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
711  ysymm_1(1:H) = 0;
712  ysymm_2(1:H) = 0;
713  for k = 2:M
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);
720  end
721  else %if comp1 is an L_I_comp
722  ysymm_1(1:H) = 0;
723  ysymm_2(1:2*H) = 0;
724  for k = 1:M-2
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);
731  end
732  ysymm_1((M-1)*H+1:M*H) = 0;
733  end %end of comp1-differentiation
734  y1 = 0.5*(ysymm_1 + ysymm_2);
735 
736  else %if comp2 is and L_I_comp
737  if i == 1 %if comp1 is the Id-comp
738  y1(1:M*H) = 0;
739  elseif 2 <= i && i <= 1 + model.Q_E %if comp1 is an L_E_comp
740  ysymm_1(1:2*H) = 0;
741  ysymm_2(1:H) = 0;
742  for k = 1:M-2
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);
749  end
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
753  ysymm_1(1:H) = 0;
754  ysymm_2(1:H) = 0;
755  for k = 2:M
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);
762  end
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
771  y = y(:);
772  end %end of q
773 end %end of if desired_constant
774 
775 function [current_Coefficients, coeffs] = Theta(model)
776 %[current_Coefficients, coeffs] = Theta(model)
777 %
778 % function which computes the coefficients of the current mu set in the
779 % model. Includes various, required cases.
780 
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);
791  for i = 1:Q
792  for j = 1:Q
793  current_Coefficients((i-1)*Q+j) = coeffs(i)*coeffs(j);
794  end
795  end
796  end
797 
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);
805  for i = 1:Q
806  for j = 1:Q
807  current_Coefficients((i-1)*Q+j) = coeffs(i)*coeffs(j);
808  end
809  end
810  end
811 
812 else
813  error('rb_problem_type unknown!')
814 end
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...
Definition: verbose.m:17
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) ...
Definition: scm_online.m:17
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)
Definition: scm_offline.m:17