home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / gfpp_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  1.7 KB  |  63 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Matrix giving maximal growth factor for Gaussian 
  4. //              elim. with pivoting.
  5.  
  6. // Syntax:      G = gfpp ( T )
  7.  
  8. // Description:
  9.  
  10. //      G is a matrix of order N for which Gaussian elimination with
  11. //      partial pivoting yields a growth factor 2^(N-1). T is an
  12. //      arbitrary nonsingular upper triangular matrix of order
  13. //      N-1. gfpp(T, C) sets all the multipliers to C  (0 <= C <= 1)
  14. //      and gives growth factor (1+C)^(N-1). gfpp(N, C) (a special
  15. //      case) is the same as GFPP(EYE(N-1), C) and generates the
  16. //      well-known example of Wilkinson. 
  17.  
  18. //      Reference:
  19. //       N.J. Higham and D.J. Higham, Large growth factors in
  20. //       Gaussian elimination with pivoting, SIAM J. Matrix Analysis and
  21. //       Appl., 10 (1989), pp. 155-164.
  22.  
  23. //    This file is a translation of gfpp.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. //-------------------------------------------------------------------//
  28.  
  29. gfpp = function ( T , c )
  30. {
  31.   local (T, c)
  32.  
  33.   if (any(any(T != triu(T))) || any(!diag(T)))
  34.   {
  35.     error("First argument must be a nonsingular upper triangular matrix.");
  36.   }
  37.  
  38.   if (!exist (c)) { c = 1; }
  39.   if (c < 0 || c > 1)
  40.   {
  41.     error("Second parameter must be a scalar between 0 and 1 inclusive.");
  42.   }
  43.  
  44.   m = T.nr; n = T.nc;
  45.   if (m == 1)    // Handle the special case T = scalar
  46.   {
  47.     n = T;
  48.     m = n-1;
  49.     T = eye(n-1,n-1);
  50.   else
  51.     n = m+1;
  52.   }
  53.  
  54.   d = 1+c;
  55.   L =  eye(n,n) - c*tril(ones(n,n), -1);
  56.   U = [T,  (d.^[0:n-2])'; zeros(1,m), d^(n-1)];
  57.   A = L*U;
  58.   theta = max(abs(A[:]));
  59.   A[;n] = (theta/norm(A[;n],inf())) * A[;n];
  60.  
  61.   return A;
  62. };
  63.