home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / housh < prev    next >
Encoding:
Text File  |  1995-11-15  |  1.2 KB  |  53 lines
  1. //--------------------------------------------------------------------------------
  2. //
  3. // housh
  4. //
  5. // Syntax: A=housh(u,j,eps)
  6. //
  7. // This routine constructs a Householder Transformation H=I-s*UU'
  8. // that "mirrors" a vector u to the Jth unit vector. If Norm(U) < eps
  9. // then zero = 1 (True) else zero = 0 (False).
  10. //
  11. // Used in tzero and tzreduce.
  12. //
  13. // REFERENCE: Ported to RLaB from the FORTRAN Code contained in:
  14. // Emami-Naeini, A., and Van Dooren, A., "Computation of Zeros of Linear
  15. // Multivaraible Systems," Automatica, Vol. 18, No. 4, April 1982, pp. 415-430.
  16. //
  17. // The routine returns a list.
  18. //
  19. // A.u     - Transformation Vector
  20. // A.s     - S Vector
  21. // A.zero  - Zero Flag
  22. //
  23. // Copyright (C), by Jeffrey B. Layton, 1994
  24. // Version JBL 931018
  25. //--------------------------------------------------------------------------------
  26.  
  27. housh = function(u,j,eps)
  28. {
  29.    local(s,alfa,dum1,u,s,zero)
  30.  
  31.    s=sum(u.*u);
  32.    alfa=sqrt(s);
  33.    if (alfa <= eps) {
  34.        zero=1;
  35.        return << u=u; s=s; zero=zero >>
  36.    }
  37.  
  38.    zero=0;
  39.    dum1=u[j];
  40.    if (dum1 > 0) {
  41.        alfa=-alfa;
  42.    }
  43.    u[j]=dum1-alfa;
  44.    s=1 ./(s-alfa*dum1);
  45.  
  46. // make sure u is column vector
  47.    if (u.nr == 1) {
  48.        u=u';
  49.    }
  50.  
  51.    return << u=u; s=s; zero=zero >>
  52. };
  53.