function  biga = op_biga( X, Zi, store_a, nonz)
% evaluate system matrix using X and Zi;
% call:  biga = op_biga( X, Zi, store_a, nonz)

[ m, m1] = size( nonz);         % dimension of dual
[ n, m1] = size( X);            % dimension of primal
biga = zeros( m);               % initialise
last = 0;
for i = 1:m;
% first we have to form (X * A_i * inv(Z))
        first = last + 1;
        last  = first + nonz( i) - 1;
        tmp = zeros( n);
        for j = first : last;
                i1 = store_a( 1, j); i2 = store_a( 2, j);
                w  = store_a( 3, j);
                xi1 = X( :, i1); xi2 = X( :, i2);
                zi1 = Zi( :, i1); zi2 = Zi( :, i2);
                if i1 == i2; tmp = tmp + w * xi1 * zi1';
                        else
                tmp = tmp + w * (xi1 * zi2'+ xi2 * zi1');
                        end;
                end;
        biga( :, i) = op_a( tmp, store_a, nonz);
        end;