% chain_start calculates matrix chain components using the matrices A, D and B 
% 
% function [UA,Sigma,Sigma_m,VD,B,A,D,G_0,G_0m,Am,Dm,r,m_0]...
%    =chain_start(daefile,y,x,t,s_A,s_D,s_B,threshold,fid,varargin);
% daefile of structure f((d(x,t))',x,t)=f(y,x,t)
% y, x, t are the parameters of f
% s_A  - .true.  analytic representation of df/dy=A,
%        .false. numerical calculation of A
% s_D -  .true.  analytic representation of d d(x,t)/dx=D in daefile
%        .false. numerical calculation of D
% s_B  - .true.  analytic representation of df/dx=B in daefile
%        .false. numerical calculation of B
% threshold - all singular values < threshold are set to zero
%             (rank determination)
% fid - file identifier 
% varargin - additional parameters of daefile
% Output:
%
% UA
% Sigma       = diag(sing. values of G_0)
% Sigma_m     reflexive generalized inverse of Sigma
% VD
% B = df/dx
% A = df/dy = UA SA VA'
% D = d d(x,t)/dx    D = UD SD VD'
% G_0 = A*D = UA * Sigma * VD'
% G_0m   reflexive generalized inverse of G_0
% Am  reflexive generalized inverse of A
% Dm  reflexive generalized inverse of D
% r - rank of Sigma
% m_0 - free parameter in G_0m ( set to zero )
%
function [UA,Sigma,Sigma_m,VD,B,A,D,G_0,G_0m,Am,Dm,r,m_0]...
    =chain_start(daefile,y,x,t,s_A,s_D,s_B,threshold,fid,varargin);
fprintf(fid,'%s \n','in chain_start');
%disp('in chain_start');
[A,D,B]= calc_A_D_B(daefile,y,x,t,s_A,s_D,s_B,varargin{:});
[m,n]=size(D);
%  D = m x n matrix
%  A = n x m matrix
[UA,SA,VA] = svd(full(A));
[UD,SD,VD] = svd(full(D));
% UA*SA*VA'-A
[rA,next_sing_value]=rank_of_Sigma(SA,threshold,fid);
[rD,next_sing_value]=rank_of_Sigma(SD,threshold,fid);
SAsigma=SA(1:rA,1:rA);
SDsigma=SD(1:rD,1:rD);
%test rA = rD
if rA == rD
  VAt=VA';
  H1=VAt(1:rA,:)*UD(:,1:rD);
  H2=VAt(1:rA,:)*UD(:,rD+1:m);
  H3=VAt(rA+1:m,:)*UD(:,1:rD);
  S=SAsigma*H1*SDsigma;
  [U,Sigma,V] = svd(S);
  [r,next_sing_value]=rank_of_Sigma(Sigma,threshold,fid);   
  if rA == r
    % The reflexive inverse of G_0 has a free parameter m_0 = (n-r,r)
    m_0=sparse(n-r,n);
    %m_0(1,1)=2.3;
    SAsigma_m=diag(1./diag(SAsigma));
    if m==r
      Am=VA*[SAsigma_m,sparse(r,n-r)]*UA';
    else
      Am=VA*[SAsigma_m,sparse(r,n-r);H3*(H1\SAsigma_m),sparse(m-r,n-r)]*UA';
    end;
    SDsigma_m=diag(1./diag(SDsigma));
    m_D=m_0(:,1:r)*Sigma*V'*SDsigma_m;
    if n==r
      Dm=VD*[SDsigma_m,SDsigma_m*(H1\H2)]*UD';
    else
      %Dm=VD*[SDsigma_m,SDsigma_m*(H1\H2);m_D,zeros(n-r,m-r)]*UD';
      H1H2=H1\H2;
      Dm=VD*[SDsigma_m,SDsigma_m*H1H2;m_D,m_D*H1H2]*UD';
    end;
    % 
    UA(:,1:rA)=UA(:,1:rA)*U;
    VD(:,1:rD)=VD(:,1:rD)*V;
    G_0=UA*[Sigma,sparse(rD,n-rD);sparse(n-rD,n)]*VD';
    Sigma_m=diag(1./diag(Sigma));
    if n==r
      G_0m=VD*[Sigma_m]*UA';
    else
      G_0m=VD*[Sigma_m,sparse(r,n-r);m_0]*UA';
    end;
  else
    error('A and D not proper - do not match');
  end;
else
  rA
  rD
  SA
  SD
  error('A and D not proper - different rank');
end;
%disp('out chain_start');
fprintf(fid,'%s \n','out chain_start');