function [Ck,Gk,Bk,Lk]=romlin(n,C,G,B,L,k,s0) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Purpose
% =======
%
%    Given a linear SISO system: 
%
%          C*x'(t) + G*x(t) = B*u(t) 
%                      y(t) = L^T*x(t) 
%  
%    This function computes the reduced-order model (ROM)
%
%          Ck*z'(t) + Gk*z(t) = Bk*u(t) 
%                       y1(t) = Lk^T*z(t) 
%
%    where Ck and Gk are k x k matrices, Bk and Lk are k x 1 vectors.    
%
% Input
% =====
%
%    n         state-space dimension  
%    C, G      system matrices, n x n
%    B, L      input and output influence vectors 
%    k         <= n, the specified order for the ROM
%              (in a future version, it should be replaced by 
%               a user-specified accuracy)
%    s0        expansion point in the frequency domain
%              (related to the frequency interval of interest)
%    
% Output
% ======
%
%    Ck, Gk    system matrices, k x k, for the ROM
%    Bk,Lk     input and output influence vectors for the ROM    
%
%
% Notes  
% =====
%
%    *  current version only deals with SISO system, i.e., 
%       u(t) and y(t) are scalar functions. In the future, the ROM of
%       a MIMO system may be developed. 
%
% 
% Written by Z. Bai, bai@cs.ucdavis.edu, 4/24/00
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[L0,U0]=lu(s0*C + G);

R1 = U0\(L0\B);

coef = L'*R1; 
R1norm = norm(R1);  
L1norm = norm(L);  

[alpha,beta,delta,eta,gamma,rho,W,V] = lanczos2(L0,U0,C,k,R1,L,n);
Tk = tridiag(k,rho(2:k),alpha(1:k),beta(2:k));

e1 = zeros(k,1); e1(1) = 1; 
Ck = -Tk;
Gk = speye(k) + s0*Tk; 
Bk = R1norm*e1;   
Lk = (coef/R1norm)*e1;