%============================================================================
%here starts script.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Script matlab file for Example 1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' MATLAB - Script 1                    ')
      disp(' Problem from BDG - Example 1 + QR factorization       ')
      disp('      QR fact. of an 8 x 8 Hamiltonian matrix, H = QR  ')
      disp('      a0  : R                       ')
      disp('      b0  : Q^T                     ')
%
%     Set the environment parameters here!!
%
%%    clear all
      format short e
%
%     Set the parameters for the solvers here!!
%
      maxiter   = 50;             % Maximum number of iterations
      toliter   = 10*sqrt( eps ); % Tolerance for convergence
      tolrank   = 1000000*eps;    % Tolerance for rank
%
%
%     Set the parameters for the matrix here!!
%
      k = 4;
      n = 2*k;
%%    Nu      = 1.0;            % Parameter
      setseed = 1563561233;     % Seed for random matrix
      randn('seed',setseed)     % Set the seed
%
%     Matrix generation
%
      A = [ -Nu    1   0   0;
             -1  -Nu   0   0;
              0    0  Nu   1;
              0    0  -1  Nu ];
      S = [ 1 1 1 1 ]'*[ 1 1 1 1 ];
      Q = S;
      H = [ A S; Q -A' ];
      [ q, r ] = qr( randn( n ) );
      H = q'*H*q;
%
%     Compute sep, dif and Delta(A)
%
      [v,d] = eig( H );
      [sorted,index] = sort( real( diag( d ) ) );
      v = v( :, index );

      [q,r]= qr(v);
      r = q'*H*q;

      A11 = r(1:k,1:k);
      A22 = r(k+1:n,k+1:n);
      B11 = eye( k );
      B22 = eye( k );
%
%     sep
%
      sep = min( svd( kron(eye(k),A11)-kron(A22',eye(k))))
%
%     Delta(A)
%
      DoA = min( abs( real(eig(H))));
%
%     dif
%
      Adif = [ kron( A11', eye( k ) ), kron( eye( k ), A22 ) ;...
               kron( B11', eye( k ) ), kron( eye( k ), B22  ) ];
      diff = min( svd( Adif ) )
%
%     Matrices of the problem
%
      [ q, r ] = qr( H );
      a0 = r;
      b0 = q';
%
%     Start running the program
%
      format long e
      disp('****************************************************')
      disp('Parameters of the problem ')
      disp(' n - seed ')
      [ n setseed ]
      disp(' Nu ')
      [ Nu ]
      disp(' maxiter - toliter - tolrank ')
      [ maxiter toliter tolrank ]
      disp(' cond(A) cond(B) - || [ A, B ] ||_F ')
      [ cond(a0) cond(b0) norm([a0, b0],'fro') ]
      disp(' D(A) sep   dif ')
      [ DoA sep diff ]
      disp('****************************************************')
      format short e
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Script matlab file for Example 2
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' MATLAB - Script 2                    ')
      disp(' Problem from BDG - Example 2 + QR factorization       ')
      disp('      QR fact. of an 2k x 2k Hamiltonian matrix, H = QR  ')
      disp('      mat1: R                       ')
      disp('      mat2: Q^T                     ')
%
%     Set the environment parameters here!!
%
%     clear all
%
%     Set the parameters for the solvers here!!
%
      maxiter   = 50;             % Maximum number of iterations
      toliter   = 10*sqrt( eps ); % Tolerance for convergence
      tolrank   = 1000000*eps;    % Tolerance for rank
%
%     Set the parameters for the matrix here!!
%
      k       = 10;             % Size of the matrix
      n       = 2*k;            % Size of the matrix
%%    alpha   = 0.4995;         % Parameter
      setseed = 1563561233;     % Seed for random matrix
      randn('seed',setseed)     % Set the seed
%
%     Matrix generation
%
      A11 = diag( ones( k, 1 )*( 1 - alpha ) ) + ...
            diag( ones( k-1, 1 )* alpha, -1 );
      A11( 1, k ) = alpha;
      A12 = randn( k );
      A22 = -A11';
      [ q, r ] = qr( randn( 2*k ) );
      H = q' * [ A11 A12; zeros( k ) A22 ] * q;
%
%     Compute sep and Delta(A)
%
      sep = min( svd( kron(eye(k),A11)-kron(A22',eye(k))) );
      DoA = min( abs( real(eig(H))) );
      Adif = [ kron( A11', eye( k ) ), kron( eye( k ), A22 ) ;...
               kron( eye( k ), eye( k ) ), kron( eye( k ), eye( k )  ) ];
      diff = min( svd( Adif ) );
%
%     Matrices of the problem
%
      [ q, r ] = qr( H );
      a0 = r;
      b0 = q';
%
%     Start running the program
%
      format long e
      disp('****************************************************')
      disp('Parameters of the problem ')
      disp(' n - seed ')
      [ n setseed ]
      disp(' alpha k ')
      [ alpha k ]
      disp(' maxiter - toliter - tolrank ')
      [ maxiter toliter tolrank ]
      disp(' cond(A) cond(B) - || [ A, B ] ||_F ')
      [ cond(a0) cond(b0) norm([a0, b0],'fro') ]
      disp(' D(A) sep   dif ')
      [ DoA sep diff ]
      disp('****************************************************')
      format short e
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Script matlab file for Example 3
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp(' MATLAB - Script 3                    ')
      disp(' Problem from BDG - Example 3 + QR factorization       ')
      disp('      QR fact. of an n x n special diagonal matrix, H = QR  ')
      disp('      mat1: R                       ')
      disp('      mat2: Q^T                     ')
%
%     Set the environment parameters here!!
%
%%    clear all
%
%     Set the parameters for the solvers here!!
%
      maxiter   = 50;             % Maximum number of iterations
      toliter   = 100*sqrt( eps ); % Tolerance for convergence
      tolrank   = 1000000*eps;    % Tolerance for rank
%
%     Set the parameters for the matrix here!!
%
      n       = 5;              % Size of the matrix
%%    d       = 1.0;            % Parameter
      setseed = 1563561231;     % Seed for random matrix 
      randn('seed',setseed)     % Set the seed
%
%     Matrix generation
%
      A11 = triu( randn( n ), +1 );
      A22 = triu( randn( n ), +1 );
      A12 = randn( n );
%     A11 = A11 + diag( abs( randn( n, 1 ) ) )*d;
%     A22 = A22 - diag( abs( randn( n, 1 ) ) )*d;
      vector = abs( randn( n, 1 ) )*d;
      A11 = A11 + diag( vector );
      A22 = A22 - diag( vector );

      [ q, r ] = qr( randn( 2*n ) );
      H = q' * [ A11 A12; zeros( n ) A22 ] * q;
      sep = min( svd( kron(eye(n),A11)-kron(A22',eye(n))) );
      DoA = min(abs(real(eig(H))));
      Adif = [ kron( A11', eye( n ) ), kron( eye( n ), A22 ) ;...
               kron( eye( n ), eye( n ) ), kron( eye( n ), eye( n )  ) ];
      diff = min( svd( Adif ) );
%
%     Matrices of the problem
%
      [ q, r ] = qr( H );
      a0 = r;
      b0 = q';
%
%     Start running the program
%
      format long e
      disp('****************************************************')
      disp('Parameters of the problem ')
      disp(' n - seed ')
      [ n setseed ]
      disp(' d ')
      [ d ]
      disp(' maxiter - toliter - tolrank ')
      [ maxiter toliter tolrank ]
      disp(' cond(A) cond(B) - || [ A, B ] ||_F ')
      [ cond(a0) cond(b0) norm([a0, b0],'fro') ]
      disp(' D(A) sep   dif ')
      [ DoA sep diff ]
      disp('****************************************************')
      format short e
%

%here ends script.m
%============================================================================