  function [Ku,lowinds]=Fpmuup(V,H,n1)
% evaluate the upper/constant part of the operator F'mu
%  Also find the matrix of 0,1's for calculating lower part
%  called from main.m
%   Input:   
%            A    -    given approx. distance matrix
%            H    -    matrix of weights
V=sparse(V);
H=sparse(H);
tn1=(n1+1)*(n1)/2;   % for dim n-1
n2=n1^2;
indn1=find(triu(ones(n1))==1);   % indices for upper triang. part
indn2=find(triu(ones(n1),1)==1);   % indices for strict upper triang. part
I=speye(n1);
colcnt=0;
colcnt2=0;
lowinds=speye(n2,tn1);
Ku=sparse(ones(tn1));
Ku2=speye(tn1);
%  block column one
for j=1:n1,
  for i=1:(j-1),
    eij=I(:,i)*I(j,:);
    eij=sparse(eij+eij');
    temp=sparse(2*KWs(KW(eij,V,H),V,H));
    temp(indn2)=sqrt(2)*temp(indn2);   % strict upper triang part redundant
    colcnt=colcnt+1;
    lowinds(:,colcnt)=eij(:);
    Ku(:,colcnt) =temp(indn1);    % use only upper triangular part
    colcnt2=colcnt2+1;
     %  improve the following using formula for position of the 1!!
    Ku2(:,colcnt2) =-sqrt(2)*eij(indn1);    % use only upper triangular part
  end
  ejj=sparse(I(:,j)*I(j,:));
  temp=sparse(2*KWs(KW(ejj,V,H),V,H));
  temp(indn2)=sqrt(2)*temp(indn2);   % strict upper triang part redundant
  colcnt=colcnt+1;
  lowinds(:,colcnt)=ejj(:);
  Ku(:,colcnt) =temp(indn1);
  colcnt2=colcnt2+1;
  Ku2(:,colcnt2) =-ejj(indn1);   %  change since Ku2 is now simple -I!!!
end
%  add block column two
Ku = sparse([Ku Ku2]);
