function [p,q,w,t,u,b,ssqdif] = pls(x,y,lv)


[mx,nx] = size(x);
[my,ny] = size(y);
if nx < lv
  error('No. of LVs must be <= no. of x-block variables')
end 
p = zeros(nx,lv);
q = zeros(ny,lv);
w = zeros(nx,lv);
t = zeros(mx,lv);
u = zeros(my,lv);
b = zeros(1,lv);
ssq = zeros(lv,2);
ssqx = 0;
for i = 1:nx
  ssqx = ssqx + sum(x(:,i).^2);
end
ssqy = sum(sum(y.^2)');
for i = 1:lv
  [pp,qq,ww,tt,uu] = plsnipal(x,y);
  b(1,i) = uu'*tt/(tt'*tt);
  x = x - tt*pp';
  y = y - b(1,i)*tt*qq';
  ssq(i,1) = (sum(sum(x.^2)'))*100/ssqx;
  ssq(i,2) = (sum(sum(y.^2)'))*100/ssqy;
  t(:,i) = tt(:,1);
  u(:,i) = uu(:,1);
  p(:,i) = pp(:,1);
  w(:,i) = ww(:,1);
  q(:,i) = qq(:,1);
end
ssqdif = zeros(lv,2);
ssqdif(1,1) = 100 - ssq(1,1);
ssqdif(1,2) = 100 - ssq(1,2);
for i = 2:lv
  for j = 1:2
    ssqdif(i,j) = -ssq(i,j) + ssq(i-1,j);
  end
end