home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / comp_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  1.6 KB  |  68 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Comparison matrices.
  4.  
  5. // Syntax:      comp ( A )
  6.  
  7. // Description:
  8.  
  9. //      COMP(A) is diag(B) - tril(B,-1) - triu(B,1), where B =
  10. //      abs(A).  COMP(A, 1) is A with each diagonal element replaced by
  11. //      its absolute value, and each off-diagonal element replaced by
  12. //      minus the absolute value of the largest element in absolute
  13. //      value in its row.  However, if A is triangular COMP(A, 1) is
  14. //      too.  COMP(A, 0) is the same as COMP(A).  COMP(A) is often
  15. //      denoted by M(A) in the literature.
  16.  
  17. //      Reference (e.g.):
  18. //      N.J. Higham, A survey of condition number estimation for
  19. //      triangular matrices, SIAM Review, 29 (1987), pp. 575-596.
  20.  
  21. //    This file is a translation of comp.m from version 2.0 of
  22. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  23. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  24.  
  25. //-------------------------------------------------------------------//
  26.  
  27. comp = function (A, k)
  28. {
  29.   local (A, k)
  30.  
  31.   if (!exist (k)) { k = 0; }
  32.   
  33.   m = A.nr; n = A.nc;
  34.   p = min(m, n);
  35.  
  36.   if (k == 0)
  37.   {
  38.     // This code uses less temporary storage than 
  39.     // the `high level' definition above.
  40.  
  41.     C = -abs(A);
  42.     for (j in 1:p)
  43.     {
  44.       C[j;j] = abs(A[j;j]);
  45.     }
  46.  
  47.   else if (k == 1) {
  48.  
  49.     C = A';
  50.     for (j in 1:p)
  51.     {
  52.       C[k;k] = 0;
  53.     }
  54.     mx = max (abs (C));
  55.     C = -mx'*ones(1,n);
  56.     for (j in 1:p)
  57.     {
  58.       C[j;j] = abs(A[j;j]);
  59.     }
  60.     if (all (A == tril(A))) { C = tril(C); }
  61.     if (all (A == triu(A))) { C = triu(C); }
  62.  
  63.   } }
  64.  
  65.   return C;
  66. };
  67.  
  68.