#
#   function [A1,X1,A2,X2,sigma]=svdcomp(H,p0)
#
#   Singular value decomposition of H (to compare with ALS solutions).
#
#   Input:
#   H       Data matrix. Columns are assumed to be normalized to
#           zero mean and unit length. If H does not verify this
#           condition, a warning is issued and it is corrected.
#   p0      Requested number of solutions
#
#   Output:
#   A1      weights (for standardized scores)
#   X1      scores  (standardized)
#   A2      weights (standardized)
#   X2      scores  (for standardized weights)
#   sigma   loss function (minimum value)
#
function [A1,X1,A2,X2,sigma]=svdcomp(H,p0)
#
# ------------------------------------------------------------------
   if nargin < 1
      help svdcomp;
      return;
   endif
# ------------------------------------------------------------------
   [n,m]=size(H);
# ------------------------------------------------------------------
#
#   Default values for parameters
#
   default_p = min([n,m]);
   epsilon=1.0e-8;
# ------------------------------------------------------------------
#
#   Check if p0 (= requested number of solutions) was entered.
#   Otherwise, set a default value.
#
   if nargin>=2
     p=min([p0,default_p]);
   else
     p=default_p;
   endif
# ------------------------------------------------------------------
#
#   Check if H is centered and standardized.
#   Otherwise correct it and warn user.
#
   s=sum(H);
   if abs(sum(s))> epsilon
      disp('*** svdcomp: Columns of H not centered. Corrected');
      H=H-ones(n,1)*s/n;
   endif
   d=diag(H'*H);
   if abs(sum(d)-m) > epsilon
      disp('*** svdcomp: Columns of H not standardized. Corrected');
      H=H*inv(diag(d));
   endif
# ------------------------------------------------------------------
   [U,L,V]=svd(H,0);
   L=L(1:p,1:p);
   U=U(:,1:p);
   V=V(:,1:p);
# ------------------------------------------------------------------
#
#   Standardized scores solution
#
   X1=U;
   A1=p*V*L;
# ------------------------------------------------------------------
#
#   Standardized weights solution
#
   X2=U*L/m;
   A2=V;
# ------------------------------------------------------------------
#
#   sigma (not yet)
#
   sigma=0;
endfunction