home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / ohess_r < prev    next >
Encoding:
Text File  |  1994-09-24  |  1.8 KB  |  65 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Random, orthogonal upper Hessenberg matrix.
  4.  
  5. // Syntax:      H = ohess ( N )
  6.  
  7. // Description:
  8.  
  9. //      H is an N-by-N real, random, orthogonal upper Hessenberg
  10. //      matrix. Alternatively, H = OHESS(X), where X is an arbitrary
  11. //      real N-vector (N > 1) constructs H non-randomly using the
  12. //      elements of X as parameters. In both cases H is constructed
  13. //      via a product of N-1 Givens rotations. 
  14.  
  15. //      Note: See Gragg (1986) for how to represent an N-by-N
  16. //      (complex) unitary Hessenberg matrix with positive subdiagonal
  17. //      elements in terms of 2N-1 real parameters (the Schur
  18. //      parametrization). This M-file handles the real case only and
  19. //      is intended simply as a convenient way to generate random or
  20. //      non-random orthogonal Hessenberg matrices.
  21.  
  22. //      Reference:
  23. //      W.B. Gragg, The QR algorithm for unitary Hessenberg matrices,
  24. //      J. Comp. Appl. Math., 16 (1986), pp. 1-8.
  25.  
  26. //    This file is a translation of ohess.m from version 2.0 of
  27. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  28. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  29.  
  30. //-------------------------------------------------------------------//
  31.  
  32. ohess = function ( x )
  33. {
  34.   local (x)
  35.   global (pi)
  36.  
  37.   if (any (imag (x))) { error("Parameter must be real."); }
  38.  
  39.   n = max(size(x));
  40.  
  41.   if (n == 1)
  42.   {
  43.     //  Handle scalar x.
  44.     n = x;
  45.     x = rand(n-1, 1)*2*pi;
  46.     H = eye(n,n);
  47.     H[n;n] = sign(rand());
  48.   else
  49.     H = eye(n,n);
  50.     H[n;n] = sign(x[n]) + (x[n]==0);   // Second term ensures H[n;n] nonzero.
  51.   }
  52.  
  53.   for (i in n:2:-1)
  54.   {
  55.     // Apply Givens rotation through angle x(i-1).
  56.     theta = x[i-1];
  57.     c = cos(theta);
  58.     s = sin(theta);
  59.     H[ [i-1, i] ;] = [ c*H[i-1;]+s*H[i;] ;
  60.                        -s*H[i-1;]+c*H[i;] ];
  61.   }
  62.  
  63.   return H;
  64. };
  65.