function [x,Lce,ErrorQ,rvalues,U]=thqrbp1(F1,F2,Ic,Ntot,LceFlag,Nskip,OrthFlag,p);
%  COMPUTATION OF THE FIRST p LCE'S using HQRBp1 using the method.
%
%  This method is efficient when n is close to p. If p is small, hqrbp2 is
%  more efficient than hqrbp1 (for this case use hqrbp2 instead of hqrbp1).
%  Here n denotes the number of states of the system and p the number of LCEs sought.
%
%  thqrbp1 Computes the trajectory and the Lyapunov Characteristic       
%  Exponents for discrete systems via The Householder QR Based Method.
%
% Information on the method can be found in: "A note on the computation of the
% largest p LCEs of discrete dynamical systems" by Firdaus  E. Udwadia, 
% Hubertus F. von Bremen and Wlodek Proskuroswki, to appear in Applied Mathematics and
% Computation.
%
% Additional information on computing all LCEs, see "An efficient QR based method for the 
% computation of Lyapunov Exponents", by Hubertus F. von Bremen, 
% Firdaus  E. Udwadia and Wlodek Proskuroswki, Physica D, vol 101 pp. 1-16, 1997.
%
% Program written by H. von Bremen 4/17/2000 
% ************************************************************************
%                                                                    
%[x,Lce,ErrorQ]=thqrbp1('xfile','tanmfile',Ic,Ntot,LceFlag,Nskip,OrthFlag,p)
%  computes the trajectory of a discrete dynamical system which is 	
%  specified in the M-file xfile.m .  It also computes the first s Lyapunov	
%  characteristic exponents (Lce's), provided the tangent maps are	
%  specified in the M-file tanmfile.m.	The trajectory is computed for 
%  up to Ntot iterations.  The LCE's are computed starting at the Nskip+1
%  iteration of the trajectory, and are computed for Ntot-Nskip iterations.				
%  									
%  INPUT:								
%  xfile	string containing the user-supplied iteration scheme for 	
%		the trajectory.						
%			CAll: xmap=fun(x,iter)  where fun = 'xfile'	
%			x  	- state vector				
%			iter	- iteration number			
%			xmap	- returned next iteration 		
%				{x(iter+1)=xmap=fun(x(iter),iter)}	
%  tanmfile	string containing the user-supplied file that defines the	
%  		tangent map to the trajectory.				
%			CAll: tanmap=fun(x,iter)  where fun = 'tanmfile'
%			x  	- state vector				
%			iter	- iteration number			
%			tanmap	- returned next iteration 		
%			{Tangent_Map(iter+1)=tanmap=fun(x(iter),iter)}	
%  Ic:		Initial conditions for the Trajectory, given as a column
%		vector									
%  Ntot:		Number of iterations (time steps) to be used to calculate
%		the trajectory of the system.				
%  LceFlag:	flag used to determine if Lce's are computed.		
%		-If LceFlag = 1, then Lce's are computed.		
%		-If LceFlag = 0, then LCE's are not computed and only 	
%		 the trajectory is computed.			
%  Nskip:	Number of iterations to be skipped before the LCE's are	
% 		to be computed.					
%  OrthFlag:	flag used to determine if error in orthogonality is to  
%		be computed, and if computed what type.
%		-If OrthFlag = 0, then the error in orthogonality is not 
%		 computed. (Recommended value)
%		-If OrthFlag = 1, then the determinant error in orthogonality is 	
%		 is computed.  Error in Orthogonality = abs(1-abs(det(Q)))
%                When s is not equal to n, OrthFlag is then set to 2.
%		-If OrthFlag = 2, then the Two-Norm error in orthogonality is 	
%		 is computed.  Error in Orthogonality = Norm(Q'*Q-I)
%     NOTE: computing the error in orthogonality will increase the computation time.
%
%  p		Number of first Lces' to be computed.
%
%  OUTPUT:								
%	x		Trajectory.  Given in matrix form, the k-th 
%			column are the trajectory coordinates for the 
%			k-th iteration (the first column is Ic ). size(x)=[n,Ntot+1].
%	Lce		Computed Lyapunov Exponents.  Given in matrix
%			form, the k-th column are the LCE's for the k-th
%			iteration of LCE's.  size(Lce)=[n,Ntot-Nskip].
%	ErrorQ		Error in othogonality	
%			-If OrthFlag = 1, ErrorQ = abs(1-abs(det(Q)))
%			-If OrthFlag = 2, ErrorQ = Norm(Q'*Q-I)			
%			Given in vector form, the k-th entry is the error
%			in orthogonality of the k-th Lce iteration.
%			size(NormQ)=[Ntot-Nskip,1].
%
%**EXAMPLE:
%  --------------------------------
%  Sample Call:	
%  Ic=[0.2  0.3]'; Ntot= 10000; LceFlag=1; Nskip=1000; OrthFlag=1;
%  [x,Lce,ErrorQ]=thqrbp1('cplgst','tglgst',Ic,Ntot,LceFlag,Nskip,OrthFlag,s);
%
%  --------------------------------
%
%  ---------------------------
%  Sample for xfile='cplgst'; coupled logistic maps (stored in m-file cplgst.m)
%  function xmap=cplgst(x,iter)
%  d=.2; r=3.6;
%  xm(1)=d*r*x(1)*(1-x(1))     +    (1-d)*r*x(2)*(1-x(2));
%  xm(2)=(1-d)*r*x(1)*(1-x(1)) +    d*r*x(2)*(1-x(2));
%  xmap=[xm(1) xm(2)]';
%  ---------------------------
%  Sample for tanmfile='tglgst'; Jacobian for coupled logistic maps (stored in m-file tglgst.m)
%  function [T]=tglgst(x)
%  d=.2; r=3.6;
%  T=  [  d*r*(1-2*x(1))        (1-d)*r*(1-2*x(2))
%          (1-d)*r*(1-2*x(1))       d*r*(1-2*x(2)) ];
%  ---------------------------
%  --------------------------------
%  Sample to plot output:
%  To plot the k-th LCE for each iteration:  plot(Lce(k,:))
%  To plot the k-tk coordinate of the trajectory:  plot(x(k,:))
%  To plot the error in orthogonality:  plot(ErrorQ)
%  --------------------------------
%
% NOTE: As an additional example, the n-dimensional generalization of coupled logistic maps 
% is provided as follows.
%  ---------------------------
%  Sample for xfile='cplgstn'; the generalized coupled logistic maps (stored in m-file cplgstn.m)
%  function xmap=cplgstn(x,iter)
%  This is an nd-map generalization of the 2d-logistic map . . .
%  d=2.01;r=1.8;
%  n=length(x);
%  y=(x-x.^2);
%			%	y=[x(1)*(1-x(1)) x(2)*(1-x(2)) x(3)*(1-x(3))...
%			%	x(4)*(1-x(4)) x(5)*(1-x(5)) x(6)*(1-x(6))...
%			%	x(7)*(1-x(7)) x(8)*(1-x(8)) x(9)*(1-x(9))...
%			%  x(10)*(1-x(10))]';
%  A=d*r*diag(ones(n,1))+(-1+d/2)*r*diag(ones(n-1,1),1)+...
%     (-1+d/2)*r*diag(ones(n-1,1),-1);
%  xmap=A*y;
%  ---------------------------
%  Sample for tanmfile='tglgstn'; Jacobian for coupled logistic maps (stored in m-file tglgstn.m)
%  function [T]=tglgstn(x)
%  % this ia an nd-map generalization of the 2d-logistic map . . .
%  d=2.01;r=1.8;
%  n=length(x);
%  y=(1-2*x)';
%  A=d*r*diag(ones(n,1))+(-1+d/2)*r*diag(ones(n-1,1),1)+...
%     (-1+d/2)*r*diag(ones(n-1,1),-1);
%  for i=1:n
%    B(i,:)=A(i,:).*y;
%  end;
%  T=B;
%  --------------------------------


%  INITIALIZATION PART
  m=Ntot-Nskip;s=p;
  n=max(size(Ic));
  r=ones(1,s);
  q=eye(n);sl=zeros(s,1);Lce=zeros(s,m);rvalues=zeros(s,m);
  gama=zeros(s,1);
  if OrthFlag==1|OrthFlag==2
     ErrorQ=zeros(m,1);
  end  % (if)
  x=zeros(n,Ntot+1);
  
temp=2*n^2*p+(5/6)*p^3+5.5*n*p-4.5*n*p^2-2*n^2;
if (temp>0)
   disp('For the number of LCEs you are computing it may take')
   disp('fewer flops if you use HQRBp2 instead of HQRBp1')
end

   
if (OrthFlag==1) & (s~=n)
   disp('Not a square orthogonal matrix (s not eq. n) ')
   disp('OrthFlag will be set to 2')
   OrthFlag=2;
end
% END INITIALIZATION

x(:,1)=Ic;
for j=2:Nskip+1
  Ic=feval(F1,Ic);
  x(:,j)=Ic;
end  % (j)

if LceFlag==0
   for j=1:m
     Ic=feval(F1,Ic);
   x(:,j+Nskip+1)=Ic;
   end  % (j)
else 

  A=feval(F2,Ic);
  if OrthFlag==0
     a=A;
  end;
  N2=Nskip+1;
  for i=1:m
     if OrthFlag~=0
        a=A*q;
        Ic=feval(F1,Ic);
        x(:,i+N2)=Ic;
        A=feval(F2,Ic);
        B=eye(n);
     elseif OrthFlag == 0 
       Ic=feval(F1,Ic);
       x(:,i+N2)=Ic;
       B=feval(F2,Ic);
     end % (if OrthFlag)
% Computation of the Factorization
        if s==n
           w=n-1;
        else
           w=s;
        end
	for k=1:w
% computation of the reflector
%******************************************************
  	if a(k,k)<0
  	  b=-1;
	  else
	    b=1;
	  end
	  sig=sqrt(a(k:n,k)'*a(k:n,k));
	  gama(k)=sig*(sig+abs(a(k,k)));
	  r(k)=-b*sig;
	  a(k,k)=a(k,k)-r(k);
% end computation of the reflector
     p=k+1;
	  a(p:n,p:s)=a(p:n,p:s)-repmat((a(k:n,k)'*a(k:n,p:s))/gama(k),[n-p+1,1]).*repmat(a(p:n,k),[1,s-p+1]);
% Computation of B*Q

if k~=w
			tem=a(k:n,k)/sqrt(gama(k)); 
         B(:,k:n)=B(:,k:n)-repmat(tem',[n,1]).*repmat((B(:,k:n)*tem),[1,n-k+1]);
         
     else % k==w
         for j=1:n  
	         b=B(j,k:n)*a(k:n,k)/gama(k);
	         B(j,k:s)=B(j,k:s)-a(k:s,k)'*b;
	      end
      end % if k
	end  % end for k
        if s==n
           r(n)=a(n,n);
        end
%*********************************************
% End Computation of the factorization 
     if OrthFlag~=0
        q=B(:,1:s);
        if OrthFlag==1
           ErrorQ(i)=abs(1-abs(det(q)));
        elseif OrthFlag==2
           ErrorQ(i)=norm(q'*q-eye(s));
        end  % (if)
     elseif OrthFlag == 0 % (if OrthFlag)
        a=B(:,1:s);
     end  % (if OrthFlag)
     sl=sl+log(abs(r'));
rvalues(:,i)=r';
     Lce(:,i)=sl/i;
  end  % (i)
end % (else)