Вычисление определителя методом QR-разложения

Этот способ основывается на QR-разложении матрицы и является более устойчивым к накоплению вычислительных ошибок. Матрица A приводится к виду A = QR, где определитель матрицы Q равен 1, а матрица R верхнетреугольная. После такого разложения определитель матрицы A равен определителю R и равен произведению диагональных элементов.

Если нашли ошибку в алгоритме - сообщите!

Реализации алгоритма на различных языках:

Реализация алгоритма на AlgoPascal:

unit DetQRUnit;

interface DetQR;

implementation

(*******************************************************************
Функция вычисления определителя матрицы A методом QR-разложения.

Входные данные:
    A:  array [1..N, 1..N] - разлагаемая матрица
    N:  - размер матрицы
*******************************************************************)
function DetQR(
    var     A:  array of array of Real;
    const   N:  Integer):Real;
var
    Col     :   Integer;
    Row     :   Integer;
    I       :   Integer;
    J       :   Integer;
    X1      :   Real;
    X2      :   Real;
    T       :   Real;
    T1      :   Real;
    T2      :   Real;
    C       :   Real;
    S       :   Real;
    Denomin :   Real;
begin
    //вращение
    for Col:=1 to N do
        for Row:=Col+1 to N do
        begin
            //выбор угла вращения
            X1:=A[Col,Col];
            X2:=A[Row,Col];
            if (X1=0) and (X2=0) then
                Continue;
            if AbsReal(X1)>AbsReal(X2) then
            begin
                T:=X2/X1;
                Denomin:=AbsReal(X1)*Sqrt(1+T*T)
            end
            else
            begin
                T:=X1/X2;
                Denomin:=AbsReal(X2)*Sqrt(1+T*T);
            end;
            C:=X1/Denomin;
            S:=-X2/Denomin;

            //поворот матрицы коэффициентов и обратной матрицы
            for I:=1 to N do
            begin
                T1:=A[Col,I];
                T2:=A[Row,I];
                A[Col,I]:=C*T1-S*T2;
                A[Row,I]:=S*T1+C*T2;
            end;
        end;
        
    //подсчет определителя
    Result:=1;
    for I:=1 to N do
        Result:=Result*A[I,I];
end;

end.

Бочканов Сергей, Быстрицкий Владимир
Copyright © 1999-2004
При поддержке проекта MANUAL.RU