gram_schmidt_parallel.m

Go to the documentation of this file.

function Xon = gram_schmidt_parallel(X, K, epsilon, n_on)
% function Xon = gram_schmidt_parallel(X, k, epsilon, n_on)
% performs classical gram schmidt orthonormalization algorithm with
% re-orthogonalization. loss of orthogonality is a multiple of machine
% precision.
% - X: vectors to be orthonormalized
% - K: inner product matrix
% - epsilon: norm threshold
% - n_on: if given, first n_on vectors are assumed to be already orthonormal

if nargin < 4
    n_on = 0;
end

pool = gcp;

nworkers = pool.NumWorkers;
nvecs = size(X, 2);
ndofs = ceil(size(X, 1) / nworkers);

X_p = Composite(nworkers);
K_p = Composite(nworkers);

for i = 1:nworkers
    X_p{i} = X(1+(i-1)*ndofs : min(i*ndofs, size(X, 1)), :);
    K_p{i} = K(1+(i-1)*ndofs : min(i*ndofs, size(K, 1)), :);
end

clear('X');
clear('K');

spmd
    
    % storing matrix avoids multiple computation of inner products
    KX_p = K_p * gcat(X_p, 1);
    
    for i = (n_on + 1):nvecs
        
        v_p = X_p(:, i);
        
        norm = sqrt(gplus(v_p' * KX_p(:, i)));
        
        niter = 0;
        
        % two orthogonalization iterations suffice
        while norm > epsilon && niter < 2
            
            KXv = gplus(KX_p' * v_p);
            KXv(n_on+1 : end) = 0;
            
            v_p = v_p - X_p * KXv;
            
            w = K_p * gcat(v_p, 1);
            
            norm = sqrt(gplus(v_p' * w));
            
            niter = niter + 1;
            
        end
        
        % include new vector
        if norm > epsilon
            n_on = n_on + 1;
            
            X_p(:, n_on) = v_p / norm;
            KX_p(:, n_on) = w / norm;
        end
    end
end

Xon = vertcat(X_p{:});
Xon = Xon(:, 1:n_on{1});

end