#
#   function [A,X,sigma,num_iter]=als3b(H,p0,maxit0,X00,epsilon0)
#
#   Normalized weights algorithm. Multiple solutions.
#   H is assumed standardized. Gram-Schmidt orthonormalization of A.
#
#   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 (optional).
#   maxit0, X00, 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]=als3b(H,p0,maxit0,X00,epsilon0)
#
   [n,m]=size(H);
# ------------------------------------------------------------------
#
#   Default values for parameters
#
   default_p = min([n,m]);
   default_epsilon=1.0e-10;
   default_maxit=50;
# ------------------------------------------------------------------
#
#   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_X0=rand(n,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 X0 (=initial scores matrix) was entered. If so,
#   check its dimensions. Otherwise, set a default value.
#
   if nargin>=4
      [n1,n2]=size(X00);
      if n1==n && n2<=p,
         p=min([n2,p]);
         X0=X00(:,1:p);
      else
         disp('*** als3b: Wrong sized X0. Set default');
         X0=default_X0;
      endif
   else
      X0=default_X0;
   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('*** als3b: 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('*** als3b: Columns of H not standardized. Corrected.');
      H=H*inv(diag(d));
   endif
# ------------------------------------------------------------------
#
#   Initialize X = X0 and A = 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=X0;
   A=zeros(m,p);
   Sigma=zeros(1,m);
   sigma=0;
   A_conv=0;
   sigma_conv=0;
# ------------------------------------------------------------------
#
#   Iterate
#
    for i=1:maxit
         previous_sigma=sigma;
         previous_A=A;
         A=H'*X;
# ------------------------------------------------------------------
#
#   Orthonormalization of A
#
         U=grschm(A);
         A=U;
# ------------------------------------------------------------------
#
#   Compute next X.
#
         X=H*A/m;
# ------------------------------------------------------------------
#
#   Check convergence for A. Display A if sigma already converged.
#
         if sigma_conv
            format long
            A
            format
            input("Press ENTER to continue...");
         endif
         if ssq(A-previous_A)<epsilon
            disp('*** als3b: Convergence for A attained');
            A_conv=1;
         endif
# ------------------------------------------------------------------
#
#   Check convergence for sigma.
#
         for j=1:m
            Sigma(j)=ssq(X-H(:,j)*A(j,:));
         endfor
         sigma=sum(Sigma)/m;
         if !sigma_conv
            format long
            sigma
            format
         endif
         if abs(sigma-previous_sigma)<epsilon
            if !sigma_conv
               disp('*** als3b: Convergence for sigma attained');
            endif
            sigma_conv=1;
         endif
# ------------------------------------------------------------------
         if A_conv
            break
         endif
    endfor
    num_iter=i;
endfunction