Вычислительные методы линейной алгебры. Горбаченко В.И - 37 стр.

  [N, N]=size(A);
  P=eye(N);
  L=zeros(N,N);
  U=eye(N); % формирование единичной диагонали матрицы U
  flag=1; % Признак невырожденности матрицы
  epsilon=1e-15; % Малая величина для определения вырождения матрицы
for r=1:N
     % вычисление элементов столбца r матрицы L (на месте элементов A)
     for i=r:N
         s=0;
         for k=1:r-1
              s=s+A(i,k)*A(k,r);
         end
         A(i,r)=A(i,r)-s;
     end
     % Выбор главного элемента для столбца r
     % Y - максимальный элемент одномерного массива abs(A(r:N,r))
     % j - номер максимаьного элемента в массиве abs(A(r:N,r))
     [Y,j]=max(abs(A(r:N,r)));
     % Перестановка строк r и j матрицы A
     C=A(r,:); % Запоминание строки r
     A(r,:)=A(r+j-1,:);%r-номер строки в A, j-номер строки в A(r:N,r)
     A(r+j-1,:)=C;
     % Перестановка строк в матрице перестановок P
     C=P(r,:); % Запоминание строки r
     P(r,:)=P(r+j-1,:);%i-номер строки в A, j-номер строки в A(r:N,r)
     P(r+j-1,:)=C;
     % Проверка вырожденности матрицы
     if abs(A(r,r))<=epsilon
         flag=0;
         return % выход из функции
     end
    % вычисление элементов строки r матрицы U (на месте элементов A)
     if r~=N % диагональ матрицы U равна 1
         for i=r+1:N
              s=0;
              for k=1:r-1
                   s=s+A(r,k)*A(k,i);
              end
               A(r,i)=(A(r,i)-s)/A(r,r);
         end
     end
  end
for i=1:N
     L(i,1:i)=A(i,1:i);
     U(i,i+1:N)=A(i,i+1:N);
end

     В функции реализован алгоритм с выбором главного элемента по
столбцу. Элементы матриц L и U строятся на месте исходной матрицы A
(единичные диагональные элементы матрицы U не хранятся). На каждом
шаге r вычисляются элементы lir столбца r матрицы L и размещаются на

местах   элементов air .   Если max lir = l jr , то строки r и j меняются
                                 i≥r

                                                                         37