%
%   version 29.07.03
%
% matrix_sequence determines the index of a DAE described in the routine daefile
%
% function index=matrix_sequence(daefile,xfunction,t,s_A,s_D,s_B,threshold,fid,argx,varargin);
%
% daefile of structure f((d(x,t))',x,t)=f(y,x,t)
% xfunction - x(t), where its parameters are the first argx parameters in varargin
% t - timepoint
% 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
% given_index - reduces the computation if you give the known index
%               -1  - check of the index (default)
%                k  - integer value of the given index
% threshold - all singular values < threshold are set to zero
%             (rank determination)
% fid - file identifier
% argx - number of parameters of xfunction in varargin
% varargin - additional arguments of xfunction and daefile
%
function [ok,index,A,R,D,rankD,B0,Dm,Q0,P0,W0,DP1,rankDP1,W1,rankW1,U]=...
    matrix_sequence(daefile,xfunction,t,s_A,s_D,s_B,given_index,threshold,fid,argx,varargin);
ok=1;
fprintf(fid,'%s \n','in matrix_sequence');
y=cal_dxtprime(daefile,xfunction,t,argx,varargin{:});
x=feval(xfunction,t,varargin{1:argx});
[Ui_tilde,diagSi,Sim,Faci,B0,A,D,Gi,Gim,Am,Dm,dim,m_i]=...
    chain_start(daefile,y,x,t,s_A,s_D,s_B,threshold,fid,varargin{argx+1:end});
%
%   initialization of vector with the ranks of the different levels
%
rank_vec=[dim];
%
mS=m_i(:,1:dim)*diagSi;
%
%   storage of V_D for later use (computation of rank DP_1)
%
V_D=Faci;
%
%   how to compute the matrices
%
text=' (computed ';
num='numerically)';
ana='analytically)';
if s_A
    textA=['A',text,ana];
else
    textA=['A',text,num];
end;
if s_D
    textD=['D',text,ana];
else
    textD=['D',text,num];
end;
if s_B
    textB=['B',text,ana];
else
    textB=['B',text,num];
end;
print_mat(fid,textA,A);
print_mat(fid,textD,D);
print_mat(fid,textB,B0);
%
%   calucation of R
%
R=Am*A;
rankD=dim;
% check of the correct relations between the different matrices
% all the matrices diff* have to have the value zero
%diffP=norm(Gim*Gi-Dm*D,inf);
%diffR=norm(Am*A-D*Dm,inf);
%diffAD=norm(A*D-Gi,inf);
%fprintf(fid,'\n\n%s \n',['different representations of projector P: G0m*G0-Dm*D = ',num2str(diffP)]);
%fprintf(fid,'%s \n',['different representations of projector R: Am*A-D*Dm   = ',num2str(diffR)]);
%fprintf(fid,'%s \n',['different representations of projector A*D: A*D-G0    = ',num2str(diffAD)]);
%clear diffP diffR diffAD;
%%%%%%%%%%%%%%%%%%%%%%
index=-1;
n=size(B0,1);
%
%   print of rank
%
fprintf(fid,'\n%s\n',['Rank of 0th level (G0):',num2str(dim)]);
%
if dim==n | given_index == 0
    %disp('index = 0');
    fprintf(fid,'\n%s\n','index = 0');
    index=0;
else
  Pi=zeros(n,n,n);
  Qi=zeros(n,n,n);
  Wi=zeros(n,n,n);
  %
  % only for representation
  %
  %disp(['norm(G',num2str(0),'*Q',num2str(0),')    = ',num2str(norm(Gi*Qi(:,:,1),inf))]);
  %
  % calculation of the projectors P0, Q0 and W0
  %
  [Pi(:,:,1),Qi(:,:,1),Wi(:,:,1)]=P_Q_and_W(Gi,Gim);
  print_mat(fid,' P0 =',Pi(:,:,1));
  print_mat(fid,' Q0 =',Qi(:,:,1));
  print_mat(fid,' W0 =',Wi(:,:,1));
  P0=Pi(:,:,1);
  Q0=Qi(:,:,1);
  W0=Wi(:,:,1);
  if given_index == 1
      index = given_index;
  else
      %
      % next level
      %
      B0q=Ui_tilde'*B0*Faci;
      dimp1=dim+1;
      [Ui,Sigma,Vi] = svd(B0q(dimp1:n,dimp1:n));
      [ri,next_sing_value]=rank_of_Sigma(Sigma,threshold,fid);
      Si=Sigma(1:ri,1:ri);
      bb=Faci(:,dimp1:n);
      Zni=[];
      DPi=D;
      if given_index == -1
          ende=n;
      else
          ende=given_index-1;
      end;
      % the loop variable has n as its maximal value, because that
      % is the maximal possible index for an n-dimensional dae
      for k=1:ende
        %
        %   updating of ranks vector of the different levels
            rank_vec=[rank_vec,ri];
        %
        %   print of rank
        %
            fprintf(fid,'\n%s\n',['rank of ',num2str(k),'th level: ',num2str(ri)]);
        % 
        %
        % check of index
        %
        if dim+ri == n
          %disp(['index = ',num2str(k)]);
          fprintf(fid,'\n%s\n',[' index = ',num2str(k)]);
          index=k;
          break;
        else
          %--------------------------------------------------------
          if k==1
              Bi=B0;
              prodPi=speye(n);
          else
              %disp('sequence:mit ableitung');
              DPiDm=DPi*Dm;
              [Bi,ok_Bi]=Bi_calc(k,Bi,Gi,D,Dm,prodPi,Pi(:,:,k-1),DPiDm,daefile,xfunction,t,s_A,s_D,s_B,...
                  threshold,fid,rank_vec,argx,varargin{:});
              if ~ok_Bi==1
                  fprintf(fid,'\n%s\n',['matrix_sequence: Error computing B',num2str(k-1),...
                          ' ! It looks like a singularity. Check also your parameters.']);
                  ok=0;
                  return;
              end;
              prodPi=prodPi*Pi(:,:,k-1);
          end;
          BQi=Bi*Qi(:,:,k);
          Gi=Gi+BQi;
          Biq=Ui_tilde'*Bi*Faci;
          %
          %  G(i+1)=Gi+Bi*Qi=(Gi+Wi*Bi*Qi)*F=(Gi+Wi*Bi*Qi)(I+Gim*Bi*Qi)
          %  F^(-1)=Fm=I-Gim*Bi*Qi
          %
          aux=Gim*BQi;
          Fm=speye(n)-aux;
          %
          if k == 1
            %
            %     preparation for the calculation of U = F P1 Fm for k=1
            %
            F0=speye(n)+aux;
            F0m=Fm;
            %
            %   rank of W1
            %
            rankW1=n-dim-ri;
          end;
          %
          % F(i+1)=Fm*Fi*Ti*Vqi
          %
          Faci=Fm*Faci;
          Faci=trans_mat_mult(Faci,mS);
          Faci=block_mult(Faci,Vi);
          %
          % updating of Zni
          %
          Biq(dim+1:n,dim+1:n)=Ui'*Biq(dim+1:n,dim+1:n);
          Zni=[Zni,[Biq(:,dim+1:n)]];
          if ~ri==0
             diagSi=[diagSi,sparse(dim,ri);sparse(ri,dim),Si];
             Sim=[Sim,sparse(dim,ri);sparse(ri,dim),diag(1./diag(Si))];
          end;
          dim=dim+ri;
          wwi=Faci\bb;
          wwi22q=wwi(dim+1:n,:);
          [m_i,sing]=WZm(wwi22q,Zni,threshold);
          if sing
              fprintf(fid,'%s \n','matrix_sequence: Matrix pencil singular!');
              ok=0;
              %
              % only to assign something
              %
              rankDP1=rankD;
              W1=sparse(n);
              U=sparse(n);
              DP1=D;
              rankW1=0;
              return;
              %error('pencil singular');
          end;
          mS=m_i(:,1:dim)*diagSi;
          %
          % computation of rank DP_1
          %
          if k==1
              aux=V_D'*Faci*[speye(dim);mS];
              diags=svd(aux(1:rankD,:));
              rankDP1=rank_of_Sigma(diag(diags),threshold,fid);
          end;
          %
          % updating of bb
          %
          bb=[bb,Faci(:,dim+1:n)];
          %
          %                     ( I     )
          % Ui_tilde = Ui_tilde*(       )
          %                     ( 0  Ui )
          %
          Ui_tilde=block_mult(Ui_tilde,Ui);
          %
          % Gim - reflexive generalized inverse of Gi
          %
          Gim=Faci*[Sim,sparse(dim,n-dim);m_i]*Ui_tilde';
          %
          % calculation of the projectors Pk, Qk and Wk
          %
          [Pi(:,:,k+1),Qi(:,:,k+1),Wi(:,:,k+1)]=P_Q_and_W(Gi,Gim);
          print_mat(fid,[' P',num2str(k),' ='],Pi(:,:,k+1));
          print_mat(fid,[' Q',num2str(k),' ='],Qi(:,:,k+1));
          print_mat(fid,[' W',num2str(k),' ='],Wi(:,:,k+1));
          %
          % only for representation
          %
          %disp(['norm(G',num2str(k),'*Q',num2str(k),')    = ',num2str(norm(Gi*Qi(:,:,k+1),inf))]);
          %
          % DPi = D*P1*P2...Pk
          %
          DPi=DPi*Pi(:,:,k+1);
          if given_index == k+1
              index=given_index;
          else
              %
              % next level
              % 
              B0q=Ui_tilde'*B0*Faci;
              [Ui,Sigma,Vi] = svd(full(B0q(dim+1:n,dim+1:n)));
              [ri,next_sing_value]=rank_of_Sigma(Sigma,threshold,fid);  
              Si=Sigma(1:ri,1:ri);
          end;   
       end;
     end;
   end;
   % test output
   %for k=2:index
   %  text=['Q',num2str(k-2)];
   %  disp(['norm(',text,'*',text,'-',text,') = ',num2str(norm(Qi(:,:,k-1)*Qi(:,:,k-1)-Qi(:,:,k-1),inf))]);
   %  for j=1:k-1
   %     disp(['norm(Q',num2str(k-1),'*Q',num2str(j-1),')    = ',num2str(norm(Qi(:,:,k)*Qi(:,:,j),inf))]);
   %  end;
   %end;
   %text=['Q',num2str(index-1)];
   %disp(['norm(',text,'*',text,'-',text,') = ',num2str(norm(Qi(:,:,index)*Qi(:,:,index)-Qi(:,:,index),inf))]);
end;
switch index
case 2
    %
    %   projectors for index 2 calculations
    %
    W1=sparse(Wi(:,:,2));
    DP1=D*Pi(:,:,2);
    %save 'DP1neu' DP1
    U=F0*Pi(:,:,2)*F0m;
case 1
    rankDP1=rankD;
    W1=sparse(n);
    U=sparse(n);
    DP1=D;
    rankW1=0;
case 0
    P0=speye(n);
    Q0=sparse(n);
    W0=speye(n);
    rankDP1=n;
    DP1=speye(n);
    W1=sparse(n);
    U=sparse(n);
    rankW1=0;
otherwise
    W1=sparse(Wi(:,:,2));
    DP1=D*Pi(:,:,2);
    U=F0*Pi(:,:,2)*F0m;
end;
fprintf(fid,'%s \n','out matrix_sequence');
%disp('out matrix_sequence');