%Implement the joint Lanczos/Arnoldi method 3 for reducing to Hessenberg
%            V'*A*V = H_{m x m}  where A=A_{n x n}
%Usage:  [V H]=hess_as(A,m)       for a n x n matrix A (id=0)
%        [V H]=hess_as(A,m,id)    for a n x n matrix A (id=0/1)
%        [V H]=hess_as(A,m,q1)    for a n x n matrix A (id=0)
%        [V H]=hess_as(A,m,q1,id) for a n x n matrix A (id=0/1)
%% e.g.                   [If id=1, some intermediate results are shown]
%  A=rand(6),  [V H]=hess_as(A,3)
%  A=A'+A, q=rand(6,1);q=q/norm(q),  [V T]=hess_as(A,3,q)
function [P, H] = hess_Lanczos_Arnoldi(A, m, q1, id)
if nargin<1, help hess_as, return; end
[n1 n]= size(A);
if n ~= n1, disp('Error in Hess_AS.m: Matrix dimensions mismatch...'),
          Size_Input_Matrix = size(A),return,end
if nargin<2, m=n; end  %% Full Hessenberg Form Requested
     if m>n, m=n, disp('WARNING:  m is over the limit'); end
if nargin<3
   disp('                           Warning: q1 is not given... USE q1=e1 !')
   q1 = zeros(n,1); q1(1)=1;
end
   isymm = 0; if norm(A-A',1)<eps, isymm=1; end  % check symmetry (only Info)
if nargin<4, id=0; end
if nargin>=3 & max(size(q1))==1
   id=q1; disp('Warning: q1 is not given... USE q1=e1 !')
   q1 = zeros(n,1); q1(1)=1;
end
if id==1
   echo on
   disp('  ')
end

q_now = q1;  Current_A = A;  H=[]; Q=[]; Q(:,1)=q1;
 
for j=1:m

if id==1
fprintf('============ Column Step %d of the hess_as method =============== (((Start)))\n',j)
end

k_start=1; if isymm==1,k_start=j; end %% vector products (save symm)
for k=k_start:j
   H(k,j) = Q(:,k)'*A*q_now;
end
k_start=1; if isymm==1,k_start=j-1; end
              k_start=max(k_start,1);
   vec_sum = 0; 
for k=k_start:j
   vec_sum=vec_sum + Q(:,k)*H(k,j);
end
   r = A*q_now - vec_sum;
   beta = norm(r);
if beta<1.0E-12 & j<n
   disp('Beta is too small --- Invariance subspace found...')
   break
end
if j<n
   q_now = r / beta;  % new q is found
   Q(:,j+1) = q_now;
   H(j+1,j) = beta ; 
  if j<m
   if isymm==1, H(j,j+1) = beta; end %% beta - Upper (col)
  end
else
   q_now = '         No more q`s'; q_now=q_now';
end


if id==1
if isymm==0
fprintf('r_now = r_%d = Aq_%d - SUM_{k=1}^{%d}',j,j,j)
fprintf(' q_k h_{k,%d}\n', j)
else
if j<2
fprintf('r_now = r_%d = Aq_%d - alpha_%d q_%d\n',j,j,j,j)
else
fprintf('r_now = r_%d = Aq_%d - beta_%d q_%d - alpha_%d q_%d\n',j,j,j-1,j-1,j,j)
end % j
end % isymm
r_now = r'
r_norm = beta
v_sum = vec_sum'
next_q=q_now'
lab_q = sprintf('Q_%d',j);
lab_h = sprintf('H_%d',j);
mat_prt(Q,H,' | ', lab_q, lab_h);
end % info id

end % column j work end  %--------------------------Main Loop----------------

a=Q'*Current_A*Q;
Check = norm(a(:,1:m) - H(:,1:m));
if Check > 1.E-7, 
 disp('Panic: Q`AQ decomp is inaccurate'), 
 Check, return, 
else
 if isymm==0
 fprintf('%s: successful for V''*A*V=H_{%dx%d} GS-Arnoldi n=%d (m=%d)\n',...
    mfilename, m+1,m, n,m), else
 fprintf('%s: successful for V''*A*V=T_{%dx%d} GS-Lanczos n=%d (m=%d)\n',...
    mfilename, m+1,m, n,m)
 end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% (1)
P=Q;  
if id==1
echo off
disp('===================hess_as================END')
end
