rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
model_orthonormalize_gram_schmidt.m
1 function onvec = model_orthonormalize_gram_schmidt(model, model_data, vec,epsilon)
2 % function onvec = model_orthonormalize_gram_schmidt(model, model_data, vec[,epsilon])
3 %
4 % Gram-Schmidt orthonormalization of vectors in columns of
5 % matrix vec to onvec. Almost zero vectors are deleted.
6 %
7 % note: a threshold epsilon is involved in this routine for detecting
8 % zero columns. This is a quite sensitive quantity! Change with care, as
9 % many functions build on this!! Default is 1e-7
10 
11 % Bernard Haasdonk 13.6.2002
12 % Martin Drohmann 14.5.2009
13 warning('gram-schmidt orthonormalization might be inaccurate!');
14 
15 if nargin < 4
16  epsilon = 1e-7; % => for nonlinear evolution this value is required.
17 end;
18 %epsilon = 1e-10; => for linear evolution this value was OK.
19 %epsilon = 1e-12; -> This causes error by returning non-orthogonal
20 % vectors (e.g. two identical vectors)
21 % So take epsilon larger than 1e-12!
22 
23 if isempty(vec)
24  onvec = zeros(size(vec));
25  return;
26 end;
27 
28 A = model.get_inner_product_matrix(model_data);
29 
30 % check on identity of vectors (can happen, that numerics does not detect
31 % this afterwards due to rounding errors!!)
32 
33 for i=1:(size(vec,2)-1)
34  for j=(i+1):size(vec,2)
35  if isequal(vec(:,i),vec(:,j))
36  vec(:,j) = 0;
37  end;
38  end;
39 end;
40 
41 onvec = vec;
42 
43 for i = 1:size(vec,2);
44  % orthogonalize next vector i and assume, that it is already
45  % orthogonal to previous ones
46  n = sqrt(onvec(:,i)'*A* onvec(:,i));
47  if (n<epsilon)
48  onvec(:,i) = 0;
49  else
50  onvec(:,i) = onvec(:,i)/n;
51  end;
52 
53  % orthogonalize remaining vectors wrt this one:
54 
55  A_mult_onvec_i = model.get_inner_product_matrix(model_data)*onvec(:,i);
56 
57  % create row-vector with projections on the created on-vector
58  coeffs= A_mult_onvec_i' * onvec(:,i+1:end);
59 
60  % create matrix of scaled versions of actual orthonormalized vector
61  coeffmat = onvec(:,i) * coeffs;
62 
63  % perform orthogonalization
64  onvec(:,i+1:end) = onvec(:,i+1:end) - coeffmat;
65 
66 end
67 
68 % eliminate zero-columns
69 nsqr = sum(onvec.^2);
70 i = nsqr > 0.1;
71 onvec = onvec(:,i);
72 
73 %resort zero-columns to back
74 %nsqr = sum(onvec.^2);
75 %[s,i] = sort(-nsqr);
76 %onvec = onvec(:,i);
77 
78 K = onvec'*A*onvec;
79 if max(max(abs(eye(size(onvec,2))-K)))>1e-3
80  K
81  error(['check orthonormalization!! Non orthogonal vectors' ...
82  ' generated'])
83 end;