home *** CD-ROM | disk | FTP | other *** search
/ CD Concept 6 / CD Concept 06.iso / mac / UTILITAIRE / RLaB / testmatrix / bandred.r < prev    next >
Encoding:
Text File  |  1995-01-28  |  2.1 KB  |  82 lines  |  [TEXT/RLAB]
  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.