home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
The Datafile PD-CD 4
/
DATAFILE_PDCD4.iso
/
languages
/
rlab1_23a
/
CTB
/
tzreduce
< prev
next >
Wrap
Text File
|
1995-11-14
|
6KB
|
201 lines
//-------------------------------------------------------------------------------
//
// tzreduce
//
// Syntax: A=tzreduce(abcd,m,n,p,zeps,ro,sigma)
//
// This routine extracts the reduced system (a,b,c,d) with d having full row
// rank from the original system (a,b,c,d).
//
// Calling the routine as,
// [ABCD,MU,NU] = TZREDUCE(ABCD,M,N,P,Zeps,RO,SIGMA)
//
// Extracts from the (N+P)x(M+N) system [ B A ], a (NU+MU)x(M+NU) 'reduced'
// [ B' A' ] [ D C ]
// system [ D' C' ] having the same transmission zeros but with D' Full Row
// Rank. The system [ A' B' C' D' ] overwrites the old system.
//
// REFERENCE: Ported to RLaB from the FORTRAN Code contained in:
// Emami-Naeini, A., and Van Dooren, A., "Computation of Zeros of Linear
// Multivaraible Systems," Automatica, Vol. 18, No. 4, April 1982, pp. 415-430.
//
// The results are returned in a list:
//
// A.abcd = abcd;
// A.mu = mu;
// A.nu = nu;
//
// Copyright (C), by Jeffrey B. Layton, 1994
// Version JBL 931024
//-------------------------------------------------------------------------------
rfile housh
tzreduce = function(abcd,m,n,p,Zeps,ro,sigma)
{
local(Sum2,mu,nu,ro1,mnu,numu,dummy,A,rows,icol,irow,s,i,...
zero,mu,nu,dum,ibar,Rows,Adum,tau,i1,mm1,n1,i2,mn1,cols,m1)
Sum2=zeros(1,max(p,m));
mu=p;
nu=n;
while (mu != 0) {
ro1=ro;
mnu=m+nu;
numu=nu+mu;
if (m != 0 ) {
ro1=ro1+1;
irow=nu;
// *** Compress rows of D. First exploit triangular shape
m1=sigma-1;
for (icol in 1:m1) {
rows=irow+[1:ro1];
dummy=abcd[rows,icol];
A=housh(dummy,1,Zeps);
dummy=A.u;
s=A.s;
zero=A.zero;
// The following performs the same as the subroutine TR1 in the ref.
abcd[rows;icol:mnu]=(eye(ro1,ro1)-s*dummy*dummy')*abcd[rows;icol:mnu];
irow=irow+1;
}
// *** Continue with Householder with Pivoting ***
if (sigma == 0) {
sigma=1;
ro1=ro1-1;
}
if (sigma != m) {
// Special case for ro1=1
if (ro1 == 1) {
Sum2[sigma:m]=abcd[irow+1;sigma:m].*abcd[irow+1;sigma:m];
else
Sum2[sigma:m]=sum(abcd[irow+1:irow+ro1;sigma:m].*abcd[irow+1:irow+ro1;sigma:m]);
}
}
for (icol in sigma:m) {
// *** Pivot if necessary ***
if (icol != m) {
// The following 4 lines perform the same as the subroutine PIVOT in the ref.
Rows=1:numu;
dum=max(Sum2[icol:m]);
ibar=maxi(Sum2[icol:m]);
ibar=ibar+icol-1;
if (ibar != icol) {
Sum2[ibar]=Sum2[icol];
Sum2[icol]=dum;
dum=abcd[Rows;icol];
abcd[Rows;icol]=abcd[Rows;ibar];
abcd[Rows;ibar]=dum;
}
}
// *** Perform Householder Transformation ***
dummy=abcd[irow+1:irow+ro1;icol];
Adum=housh(dummy,1,Zeps);
dummy=Adum.u;
s=Adum.s;
if (Adum.zero) {
break;
}
if (ro1 == 1) {
return << abcd=abcd; mu=mu; nu=nu >>
}
// This next line replaces routine TR1 in the ref.
abcd[irow+1:irow+ro1;icol:mnu]=(eye(ro1,ro1)-s*dummy*dummy')*abcd[irow+1:irow+ro1;icol:mnu];
irow=irow+1;
ro1=ro1-1;
Sum2[icol:m]=Sum2[icol:m]-abcd[irow;icol:m] .* abcd[irow;icol:m];
}
}
tau=ro1;
sigma=mu-tau;
// *** Compress the columns of C ***
if (nu <= 0) {
mu=sigma;
nu=0;
return << abcd=abcd; mu=mu; nu=nu >>
}
i1=nu+sigma;
mm1=m+1;
n1=nu;
if (tau != 1) {
// Special case of mm1=1
if (mm1 == mnu) {
Sum2[1:tau]=(abcd[i1+1:i1+tau;mm1].*abcd[i1+1:i1+tau;mm1])';
else
// General case of mm1 != 1
Sum2[1:tau]=sum([abcd[i1+1:i1+tau;mm1:mnu].*abcd[i1+1:i1+tau;mm1:mnu]]');
}
}
for (ro1 in 1:tau) {
ro=ro1-1;
i=tau-ro;
i2=i+i1;
// *** Pivot if necessary ***
if (i != 1) {
// The next 2 lines replace the routine PIVOT in the ref.
dum=max(Sum2[1:i]);
ibar=maxi(Sum2[1:i]);
if (ibar != i) {
Sum2[ibar]=Sum2[i];
Sum2[i]=dum;
dum = abcd[i2;mm1:mnu];
abcd[i2;mm1:mnu] = abcd[ibar+i1;mm1:mnu];
abcd[ibar+i1;mm1:mnu] = dum;
}
}
// *** Perform Householder Transformation ***
cols=m+[1:n1];
dummy=abcd[i2;cols];
A=housh(dummy,n1,Zeps);
dummy=A.u;
s=A.s;
if (A.zero) {
break;
}
if (n1 != 1) {
// The following line replaces routine TR2 in the ref.
abcd[1:i2;cols] = abcd[1:i2;cols]*(eye(n1,n1)-s*dummy*dummy');
mn1=m+n1;
abcd[1:n1;1:mn1]=(eye(n1,n1)-s*dummy*dummy')*abcd[1:n1;1:mn1];
Sum2[1:i]=Sum2[1:i]-abcd[i1+1:i1+i;mn1]' .* abcd[i1+1:i1+i;mn1]';
mnu=mnu-1;
}
n1=n1-1;
}
if (!A.zero) {
ro=tau;
}
nu=nu-ro;
mu=sigma+ro;
if (ro == 0) {
return << abcd=abcd; mu=mu; nu=nu >>
}
}
return << abcd=abcd; mu=mu; nu=nu >>
};