function [err] = laplearn_rmm(V,D,m,tind,y,testy,C,B)

  %Inputs:
  % V - eigenvectors of the Laplacian
  % D - eigenvalues of the Laplacian
  % m - the number of eigenvectors to be considered (q in the paper)
  % tind - ndex of the training examples 
  % y - training labels ( + or - 1 )
  % testy - test labels ( + or - 1 )
  % C,B - parameters of the RMM

  %Outputs
  % err - error rate on the unlabeled examples

  param.MSK_DPAR_INTPNT_CO_TOL_PFEAS = 1.0e-10;

  % Dual feasibility tolerance for the dual solution
  param.MSK_DPAR_INTPNT_CO_TOL_DFEAS = 1.0e-10;

  % Relative primal-dual gap tolerance.
  param.MSK_DPAR_INTPNT_CO_TOL_REL_GAP = 1.0e-10;
  param.MSK_IPAR_INTPNT_NUM_THREADS = 3;
  eps = 1e-6;
  [temp,ind] =  sortrows(diag(D));
  V = V(:,ind(m:-1:1));
  temp1 = diag(D);
  orig = temp1(ind(m:-1:1));
  p=m-1;
  l = length(y);
  n = size(V,1);
  dy = diag(y);
  %variables [ alpha(1..l), eta(1..n), etas(1..n) tau, delta(1...m) gamma(1..n)] 
  prob.c = [ ones(l,1); -B*ones(2*l,1) ; -1 ;  eps*ones(m,1); zeros(l,1);  ];
  prob.blx = [  zeros(l,1); zeros(2*l,1) ; -inf; zeros(m,1); -inf(l,1);  ];
  prob.bux = [ C*ones(l,1); inf*ones(2*l,1); +inf; inf(m,1); inf(l,1); ];

  % alpha'*d(y)*V(tind,j:n)*V(tind,j:n)'*d(y)*alpha <= tau - delta_(j)
  prob.qcsubk = [];
  prob.qcsubi = [];
  prob.qcsubj = [];
  prob.qcval = [];
  prob.a = [];
  for i=1:p
    temp = [ sparse(l+2*l+1+m ,l+3*l+1+m);  sparse(l,l+2*l+1+m)  V(tind,i:p)*V(tind,i:p)' ];
    [tempi,tempj,tempval] = find(sparse(tril(temp)));
    prob.qcsubi = [ prob.qcsubi; tempi ];
    prob.qcsubj = [ prob.qcsubj; tempj ];
    prob.qcval = [ prob.qcval; tempval ];
    prob.qcsubk = [ prob.qcsubk; i*ones(size(tempi))];
    prob.a   = [ prob.a; sparse(1,l+2*l) -(p-i+1)  sparse(1,i-1) 1 sparse(1,m-i) sparse(1,l)];
  end
  for i=p+1:m
    temp = [ sparse(l+2*l+1+m ,l+3*l+1+m);  sparse(l,l+2*l+1+m)  V(tind,i)*V(tind,i)' ];
    [tempi,tempj,tempval] = find(sparse(tril(temp)));
    prob.qcsubi = [ prob.qcsubi; tempi ];
    prob.qcsubj = [ prob.qcsubj; tempj ];
    prob.qcval = [ prob.qcval; tempval ];
    prob.qcsubk = [ prob.qcsubk; i*ones(size(tempi))];
    prob.a   = [ prob.a; sparse(1,l+2*l) -1  sparse(1,i-1) 1 sparse(1,m-i) sparse(1,l)];
  end
  prob.blc = [ -inf*ones(m,1) ];
  prob.buc = [  zeros(m,1) ];

  prob.a = sparse([ prob.a; y' -ones(1,l) ones(1,l) zeros(1,1+m+l)]);
  prob.blc = [ prob.blc; zeros(l+1,1) ];
  prob.buc = [ prob.buc; zeros(l+1,1) ];
  prob.a = sparse([prob.a; diag(y) -eye(l) eye(l) zeros(l,m+1) -eye(l)]);
  prob.blc = [ prob.blc;  ];
  prob.buc = [ prob.buc;  ];
  [r,res]=mosekopt('maximize echo(0)',prob,param);
  gamma = res.sol.itr.xx(2*l+l+1+m+1:2*l+l+1+m+l);
  alpha = res.sol.itr.xx(1:l);
  eta   = res.sol.itr.xx(l+1:l+l);
  etas   = res.sol.itr.xx(2*l+1:2*l+l);

%  tau = res.sol.itr.xx(l+1);
%  delta = res.sol.itr.xx(l+2:l+1+m)
  beta = - res.sol.itr.suc(1:m) +  res.sol.itr.slc(1:m);
  bias = - res.sol.itr.suc(m+1) + res.sol.itr.slc(m+1);
  lambda = zeros(size(beta));
  for i=1:p 
     lambda(i) = sum(beta(1:i));        
  end
  for i=p+1:m
     lambda(i) = beta(i);
  end
  pred = (gamma'*V(tind,1:m)*diag(lambda)*V(:,1:m)' + bias)';
  teind = setdiff([1:n]',tind);
  err = sum(pred(teind).*testy <= 0)/length(testy);