home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / cod_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  1.8 KB  |  64 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Complete orthogonal decomposition.
  4.  
  5. // Syntax:    C = cod ( A , TOL )
  6.  
  7. // Description:
  8.  
  9. //    Cod computes a decomposition A = U*T*V, where U and V are
  10. //    unitary, T = [R, 0; 0, 0] has the same dimensions as  A, and R
  11. //    is upper triangular and nonsingular of dimension rank(A). Rank
  12. //    decisions are made using TOL, which defaults to approximately
  13. //    MAX(SIZE(A))*NORM(A)*EPS. 
  14.  
  15. //    Cod returns a list with elements `u', `r', and `v', such that
  16. //        A = u*r*v
  17.  
  18. //      Reference:
  19. //      G.H. Golub and C.F. Van Loan, Matrix Computations, Second
  20. //      Edition, Johns Hopkins University Press, Baltimore, Maryland,
  21. //      1989, sec. 5.4.2.
  22.  
  23. //    This file is a translation of cod.m from version 2.0 of
  24. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  25. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  26.  
  27. // Dependencies
  28.    require trap2tri
  29.  
  30. //-------------------------------------------------------------------//
  31.  
  32. cod = function ( A , tol )
  33. {
  34.   global (eps)
  35.  
  36.   m = A.nr; n = A.nc;
  37.  
  38.   // QR decomposition.
  39.   qra = qr(A, "p");        // A*qra.p = qra.u*qra.r (A*P = U*R)
  40.   V = qra.p';            // A = URV
  41.   if (!exist (tol)) 
  42.   {
  43.     // |R[1;1]| approx NORM(A).
  44.     tol = max(m,n)*eps*abs(qra.r[1;1]);
  45.   }
  46.  
  47.   // Determine r = effective rank.
  48.   r = sum(abs(diag(qra.r)) > tol);
  49.   r = r[1];        // Fix for case where R is vector.
  50.   qra.r = qra.r[1:r;];    // Throw away negligible rows (incl. all zero rows, m>n).
  51.  
  52.   if (r != n)
  53.   {
  54.     // Reduce nxr R' =  r  [L]  to lower triangular form: QR' = [Lbar].
  55.     //                 n-r [M]                                  [0]
  56.  
  57.     tr = trap2tri(qra.r');    //    [Q, R] = trap2tri(qra.r');
  58.     V = tr.q*V;
  59.     qra.r = tr.r';
  60.   }
  61.  
  62.   return << u = qra.q; r = qra.r; v = V >>;
  63. };
  64.