#
#   function [A,X,sigma,num_iter]=als3a0(H,p0,maxit0,A00,epsilon0)
#
#   Normalized scores algorithm. Multiple solutions.
#   H is assumed standardized. Naive orthonormalization of X:
#   X -> X * (X'X)^(-1/2). This algorithm doesn't work. See als3a.
#
#   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
#   maxit0, A00, epsilon0, are optional arguments, to control iteration.
#
#   Output: WARNING! X and A do NOT contain successive intermediate
#           steps, but different orthogonal solutions.
#   A  = weights
#   X  = scores
#   sigma =  loss function
#   num_iter = number of iterations
#
#
function [A,X,sigma,num_iter]=als3a0(H,p0,maxit0,A00,epsilon0)
# ------------------------------------------------------------------
   [n,m]=size(H)
# ------------------------------------------------------------------
#
#   Default values for parameters
#
   default_p = min([n,m]);
   default_epsilon=1.0e-10;
   default_maxit=25;
# ------------------------------------------------------------------
#
#   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
   default_A0=rand(m,p);
# ------------------------------------------------------------------
#
#   Check if maxit0 (=maximum number of iterations) was entered.
#   Otherwise, set a default value.
#
   if nargin>=3,
      maxit=maxit0;
   else
      maxit=default_maxit;
   endif
# ------------------------------------------------------------------
#
#   Check if A0 (=initial weights matrix) was entered. If so,
#   check its dimensions. Otherwise, set a default value.
#
   if nargin>=4
      [m1,m2]=size(A00);
      if m1==m && m2<=p,
         p=min([m2,p]);
         A0=A00(:,1:p);
      else
         disp('*** als3a0: Wrong sized A0. Set default');
         A0=default_A0;
      endif
   else
      A0=default_A0;
   endif
# ------------------------------------------------------------------
#
#   Check if epsilon0 (=tolerance for convergence) was entered.
#   Otherwise, set a default value.
#
   if nargin>=5
      epsilon=epsilon0;
   else
      epsilon=default_epsilon;
   endif
# ------------------------------------------------------------------
#
#   Check if H is centered and standardized.
#   Otherwise correct it and warn user.
#
   s=sum(H);
   if abs(sum(s))> epsilon
      disp('*** als3a0: 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('*** als3a0: Columns of H not standardized. Corrected.');
      H=H*inv(diag(d));
   endif
# ------------------------------------------------------------------
#
#   Initialize A = A0 and X = 0. Other initializations: Sigma is an
#   auxiliary [1,m] row vector, used to compute the loss function.
#   A_conv and sigma_conv are indicators of convergence.
#
   X=zeros(n,p);
   A=A0;
   Sigma=zeros(1,m);
   sigma=0;
   A_conv=0;
   sigma_conv=0;
# ------------------------------------------------------------------
#
#   Iterate
#
    for i=1:maxit
         previous_sigma=sigma;
         previous_A=A;
         X=H*A;
         for j=1:m
            Sigma(j)=ssq(X-H(:,j)*A(j,:));
         endfor
         sigma=sum(Sigma)/m;
         if abs(sigma-previous_sigma)<epsilon
            if !sigma_conv
               disp('*** als3a0: Convergence for sigma attained');
            endif
            sigma_conv=1;
         else
            format long;
            sigma
            format
            input("Press ENTER to continue...");
         endif
# ------------------------------------------------------------------
#
#   Orthonormalization of X.
#
         [U,L,V]=svd(X,0);
         X=U*V';
# ------------------------------------------------------------------
#
#   Check convergence for A
#
         A=H'*X
         format long
         format
         input("Press ENTER to continue...");
         if ssq(A-previous_A)<epsilon
            disp('*** als3a0: Convergence for A attained');
            A_conv=1;
         endif
# ------------------------------------------------------------------
         if A_conv
            break
         endif
    endfor
    num_iter=i;
endfunction