/ Регистрация

Разложение методом Холецкого

Тема в разделе "В помощь учащимся", создана пользователем SunEye, 12.07.05.

  1. 12.07.05 #1
    SunEye

    SunEye Новичок

    5
    0
    Нужно написать программу на паскале....
    Дана матрица А (симметричная и положительно определенная). Нужно написать программу, которая находит треугольную матрицу В, чтобы А=ВхВ', где В' транспонированная матрица.
    Помогите... дайте идею, как это можно реализовать, я нашел что можно разложением метода Холецкого, но сам алгоритм ни где не нашел...:bawling:
     
  2. 19.07.05 #2
    Гость

    Гость Гость

    Методом Холецкого.
    function [A,p] = chol(A)
    % Embedded MATLAB Library function.
    %
    % Limitations:
    % 1) Does not support CHOL with two output arguments
    % 2) May factor matrices that MATLAB will not.
    % 3) May fail to factor matrices that MATLAB will factor.

    % $INCLUDE(DOC) toolbox\matlab\matfun\chol.m$
    % Copyright 2002-2004 The MathWorks, Inc.
    % $Revision: 1.1.4.2 $ $Date: 2004/07/01 04:50:25 $

    % Translated from LAPACK functions zpotf2.f/cpotf2.f/dpotf2.f/spotf2.f

    eml_assert(isfloat(A),['Function ''chol'' is not defined for values of class ''' class(A) '''.']);
    eml_assert(nargin == 1, 'Not enough input arguments.');
    eml_assert(size(A,1) == size(A,2), 'Matrix must be square.');
    eml_assert(nargout<2,'Embedded MATLAB''s CHOL function does not support the two output version.');

    if isempty(A)
    return;
    end

    n = size(A,1);
    A = triu(A);
    p = 0;
    for j = 1 : n

    %s = real(A(j,j)) - A(1:j-1,j)' * A(1:j-1,j);
    t = cast(0,class(A));
    if j>1
    for k = 1 : j-1
    akj = A(k,j);
    % t = t + real(conj(A(k,j)) * A(k,j));
    t = t + real(akj)*real(akj) + imag(akj)*imag(akj);
    end
    end
    ajj = real(A(j,j)) - t;

    if ajj <= 0
    A(j,j) = ajj;
    p=j;
    if nargout <= 1
    error(['Not positive definite. For a matrix with eigenvalues ' ...
    'close to EPS, MATLAB and Embedded MATLAB may ' ...
    'disagree on whether or not it is positive definite ' ...
    'due to rounding errors.']);
    end
    return;
    end

    ajj = sqrt(ajj);
    A(j,j) = ajj;

    if j < n
    %A(j,j+1:n) = A(j,j+1:n) - A(1:j-1,j)'*A(1:j-1,j+1:n);
    for k = j+1:n
    if isreal(A)
    temp = cast(0,class(A));
    if j > 1
    for w = 1 : j-1
    temp = temp + A(w,j) * A(w,k);
    end
    end
    A(j,k) = (A(j,k) - temp)/ajj;
    else
    tre = cast(0,class(A));
    tim = cast(0,class(A));
    tim2 = cast(0,class(A));
    if j > 1
    for w = 1 : j-1
    awj = A(w,j);
    awk = A(w,k);
    tre = tre + real(awj) * real(awk) + imag(awj) * imag(awk);
    tim = tim + real(awj) * imag(awk);
    tim2 = tim2 + imag(awj) * real(awk);
    end
    end
    A(j,k) = (A(j,k) - complex(tre,tim - tim2))/ajj;
    end
    end
    end
    end

    function B = cast(A,cls)

    switch cls
    case 'double',
    B = double(A);
    case 'single',
    B = single(A);
    otherwise,
    eml_assert(false,'error','Internal error: unexpected class.');
    end
    Только это не Паскаль!!!!!! Это стандартная функция из MATLAB.