home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / bandred_r < prev    next >
Encoding:
Text File  |  1995-01-23  |  2.1 KB  |  82 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Band reduction by two-sided unitary transformations.
  4.  
  5. // Syntax:    bandred ( A , KL, KU )
  6.  
  7. // Description:
  8.  
  9. //    bandred(A, KL, KU) is a matrix unitarily equivalent to A with
  10. //    lower bandwidth KL and upper bandwidth KU (i.e. B(i,j) = 0 if
  11. //    i > j+KL or j > i+KU).  The reduction is performed using
  12. //    Householder transformations. If KU is omitted it defaults to
  13. //    KL. 
  14.  
  15. //    Called by randsvd.
  16. //    This is a `standard' reduction.  Cf. reduction to bidiagonal
  17. //    form prior to computing the SVD.  This code is a little
  18. //    wasteful in that it computes certain elements which are
  19. //    immediately set to zero! 
  20. //
  21. //      Reference:
  22. //      G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  23. //      Johns Hopkins University Press, Baltimore, Maryland, 1989.
  24. //      Section 5.4.3.
  25.  
  26. //    This file is a translation of bandred.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. // Dependencies
  31.    require house
  32.  
  33. //-------------------------------------------------------------------//
  34.  
  35. bandred = function ( A , kl , ku )
  36. {
  37.   local (A, kl, ku)
  38.  
  39.   if (!exist (ku)) { ku = kl; else ku = ku; }
  40.  
  41.   if (kl == 0 && ku == 0) {
  42.     error ("You''ve asked for a diagonal matrix.  In that case use the SVD!");
  43.   }
  44.  
  45.   // Check for special case where order of left/right transformations matters.
  46.   // Easiest approach is to work on the transpose, flipping back at the end.
  47.  
  48.   flip = 0;
  49.   if (ku == 0)
  50.   {
  51.     A = A';
  52.     temp = kl; kl = ku; ku = temp; flip = 1;
  53.   }
  54.  
  55.   m = A.nr; n = A.nc; 
  56.  
  57.   for (j in 1 : min( min(m, n), max(m-kl-1, n-ku-1) ))
  58.   {
  59.     if (j+kl+1 <= m)
  60.     {
  61.        <beta; v> = house(A[j+kl:m;j]);
  62.        temp = A[j+kl:m;j:n];
  63.        A[j+kl:m;j:n] = temp - beta*v*(v'*temp);
  64.        A[j+kl+1:m;j] = zeros(m-j-kl,1);
  65.     }
  66.  
  67.     if (j+ku+1 <= n)
  68.     {
  69.        <beta; v> = house(A[j;j+ku:n]');
  70.        temp = A[j:m;j+ku:n];
  71.        A[j:m;j+ku:n] = temp - beta*(temp*v)*v';
  72.        A[j;j+ku+1:n] = zeros(1,n-j-ku);
  73.     }
  74.   }
  75.  
  76.   if (flip) {
  77.     A = A';
  78.   }
  79.  
  80.   return A;
  81. };
  82.