The Datafile PD-CD 5/DATAFILE_PDCD5.iso/utilities/r/rlab/CTB/tzreduce

Text File | 1995-11-15 | 5.5 KB | 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 >>
};