lufac.m

matlab算法集 matlab算法集

function [L,U,P,Delta] = lufac (A,dm)
%----------------------------------------------------------------
% Usage:      [L,U,P,Delta] = lufac (A,dm)
%
% Descripton: Factor the matrix A into a lower-triangular factor
%             L, and upper-triangular factor U, and a row 
%             permutation matrix P such that:
%
%                LU = PA
%
%             Set the diagonal components of U to one (Crout's
%             method).             
%
% Inputs:     A  = n by n matrix
%             dm = optional display mode. If present,
%                  intermediate results are displayed.   
%
% Outputs:    L     = n by n lower-triangular factor
%             U     = n by n upper-triangular factor with U(k,k) = 1
%             P     = n by n row permutation matrix
%             Delta = determinant of A
%
% Notes:      The function lufac can be used with the function
%             lusub to efficiently solve the system Ax = b
%             multiple times using different right-hand side
%             vectors b. 
%----------------------------------------------------------------

% Initialize 

   Delta = 1;
   [n,n] = size(A);
   U = eye(n);
   P = eye(n); 
   D = [A P];
   display = nargin > 1;
   
% Factor A into L and U 

   for j = 1 : n

% Compute column j of D 

      for k = j : n
         for i = 1 : j-1
            D(k,j) = D(k,j) - D(k,i)*D(i,j);
         end
      end
   
% Select pivot 

      q = j;
      for i = j+1 : n
         if (abs(D(i,j)) > abs(D(q,j)))
            q = i;
         end
      end
      
% Check for singular A 

      if (abs(D(q,j)) < eps) 
         Delta = 0;
         L = zeros(n,n);
         return;
      end

% Swap rows q and j of D, change sign of Delta 

     if (q > j) 
        for i = 1 : 2*n
           d = D(q,i);
           D(q,i) = D(j,i);
           D(j,i) = d;
        end
        Delta = -Delta;
     end
        
% Compute row j of D, update delta 

      for k = j+1 : n
         for i = 1 : j-1
            D(j,k) = D(j,k) - D(j,i)*D(i,k);
         end
         D(j,k) = D(j,k)/D(j,j);
      end
      Delta = D(j,j)*Delta;

% Display intermediate results
      
      if display
         j
         q
         D
         Delta
         wait
      end

   end   
      
% Finalize 

   for k = 1 : n
      for j = 1 : n
         if (k >= j)
            L(k,j) = D(k,j);
         else
            U(k,j) = D(k,j);
         end
         P(k,j) = D(k,j+n);
      end
   end
%----------------------------------------------------------------