1 function Xon = orthonormalize_qr(X,A,epsilon)
2 %
function Xon = orthonormalize_qr(X[,A,epsilon])
4 % orthonormalization of vectors in columns of
5 % matrix X to Xon by qr-decomposition of vectors and
6 % the inner product matrix A can be specified additionally, i.e.
7 % <x,x> = x' * A * x. Then the Cholesky factoriziation of A is
8 % performed A = R_M' * R_M
10 % Then a QR-decomposition of R_M * X is performed which yields the
13 % Bernard Haasdonk 22.8.2008
16 epsilon = 1e-7; % =>
for nonlinear evolution
this value is required.
28 % check on identity of vectors (can happen, that numerics does not detect
29 % this afterwards due to rounding errors!!)
33 if isequal(X(:,i),X(:,j))
41 % incomplete cholesky of inner-product matrix:
44 R_M = ichol(sparse(A), opts);
46 % qr decomposition of R_M * X with permutation indices E
47 [Q, R, E]= qr(R_M * Xon,0 );
49 % sort such that linear independent first and Q * R = R_M * Xon
52 % search nonvanishing diagonal entries of R
53 ind = find(abs(diag(R))>epsilon);
56 Xon = Xon(:,ind) / Rind;
58 % resort vectors such that original order is recovered most:
60 % search in each row maximum entry
63 [max_val, max_i] = max(K(i,:));
65 error('found multiple matching vectors!');
70 K(:,max_i) = 0; % "delete" from further search
72 all_ind = ones(1,size(Xon,2));
74 rest_ind = find(all_ind);
77 Xon = [Xon(:,ind),Xon(:,rest_ind)];
82 %disp('check orthonormality!!');
85 %if max(max(abs(eye(size(Xon,2))-K)))>1e-3
87 % error(['check orthonormalization!! Non orthogonal vectors' ...%
91 % eliminate zero-columns
93 %i = find(nsqr > 0.1);
96 %resort zero-columns to back