Используя gmres с Матричным Входным параметром

A = gallery('wilk',21);  
b = sum(A,2);
tol = 1e-12;  
maxit = 15; 
M1 = diag([10:-1:1 1 1:10]);

x = gmres(A,b,10,tol,maxit,M1);

отображения следующее сообщение:

gmres(10) converged at outer iteration 2 (inner iteration 9) to 
a solution with relative residual 3.3e-013

Используя gmres с Указателем на функцию

Этот пример заменяет матричный A в предыдущем примере с указателем на функцию матричного векторного произведения afun и предварительный формирователь M1 с указателем на функцию backsolve mfun. Пример содержится в функциональном run_gmres это

Следующее показывает код для run_gmres:

function x1 = run_gmres
n = 21;
b = afun(ones(n,1));
tol = 1e-12;  maxit = 15; 
x1 = gmres(@afun,b,10,tol,maxit,@mfun);
 
    function y = afun(x)
        y = [0; x(1:n-1)] + ...
              [((n-1)/2:-1:0)'; (1:(n-1)/2)'].*x + ...
              [x(2:n); 0];
    end
 
    function y = mfun(r)
        y = r ./ [((n-1)/2:-1:1)'; 1; (1:(n-1)/2)'];
    end
end

Когда вы входите

x1 = run_gmres;

^MATLAB отображает сообщение

gmres(10) converged at outer iteration 2 (inner iteration 10) 
to a solution with relative residual 1.1e-013.

Используя предварительный формирователь без перезапуска

Этот пример демонстрирует использование предварительного формирователя, не перезапуская gmres.

Загрузите west0479, действительное 479 479 несимметричная разреженная матрица.

load west0479;
A = west0479;

Установите допуск и максимальное количество итераций.

tol = 1e-12;
maxit = 20;

Задайте b так, чтобы истинное решение было вектором из всех единиц.

b = full(sum(A,2));
[x0,fl0,rr0,it0,rv0] = gmres(A,b,[],tol,maxit);

fl0 равняется 1, потому что gmres не сходится к требуемому допуску 1e-12 в требуемых 20 итерациях. Лучшее приближенное решение, что возвраты gmres являются последним (как обозначено it0(2) = 20). MATLAB хранит остаточную историю в rv0.

Постройте поведение gmres.

semilogy(0:maxit,rv0/norm(b),'-o');
xlabel('Iteration number');
ylabel('Relative residual');

График показывает, что решение медленно сходится. Предварительный формирователь может улучшить результат.

Используйте ilu, чтобы сформировать предварительный формирователь, поскольку A несимметричен.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-5));

Error using ilu
There is a pivot equal to zero. Consider decreasing 
the drop tolerance or consider using the 'udiag' option.

Обратите внимание, что MATLAB не может создать неполный LU, когда это привело бы к сингулярному фактору, который бесполезен как предварительный формирователь.

Как обозначено сообщением об ошибке, попробуйте еще раз с уменьшаемым допуском отбрасывания.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
[x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);

fl1 0, потому что gmres управляет относительной невязкой к 9.5436e-14 (значение rr1). Относительная невязка является меньше, чем предписанный допуск 1e-12 в шестой итерации (значение it1(2)), когда предобусловленный неполной LU-факторизацией с допуском отбрасывания 1e-6. Выводом, rv1(1), является norm(M\b), где M = L*U. Выводом, rv1(7), является norm(U\(L\(b-A*x1))).

Следуйте прогресс gmres путем графического вывода относительных невязок в каждой итерации, начинающей с первоначальной оценки (выполните итерации номера 0).

semilogy(0:it1(2),rv1/norm(b),'-o');
xlabel('Iteration number');
ylabel('Relative residual');

Используя предварительный формирователь с перезапуском

Этот пример демонстрирует использование предварительного формирователя с перезапущенным gmres.

Загрузите west0479, действительное 479 479 несимметричная разреженная матрица.

load west0479;
A = west0479;

Задайте b так, чтобы истинное решение было вектором из всех единиц.

b = full(sum(A,2));

Создайте неполный предварительный формирователь LU как в предыдущем примере.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));

Преимущество к использованию перезапущенного gmres должно ограничить объем памяти, требуемый выполнить метод. Без перезапуска gmres требует, чтобы векторы maxit устройства хранения данных сохранили основание подпространства Крылова. Кроме того, gmres должен ортогонализировать против всех предыдущих векторов на каждом шаге. Перезапуск ограничивает количество используемой рабочей области и объем работы, сделанный на внешнюю итерацию. Обратите внимание на то, что даже при том, что предобусловленный gmres сходился в шести итерациях выше, алгоритм допускал целых двадцать базисных векторов и поэтому, выделил все то место впереди.

Выполните gmres(3), gmres(4) и gmres(5)

tol = 1e-12;
maxit = 20;
re3 = 3;
[x3,fl3,rr3,it3,rv3] = gmres(A,b,re3,tol,maxit,L,U);
re4 = 4;
[x4,fl4,rr4,it4,rv4] = gmres(A,b,re4,tol,maxit,L,U);
re5 = 5;
[x5,fl5,rr5,it5,rv5] = gmres(A,b,re5,tol,maxit,L,U);

fl3, fl4 и fl5 являются всем 0, потому что в каждом случае, перезапущенном, gmres управляет относительной невязкой к меньше, чем предписанный допуск 1e-12.

Следующие графики показывают, что истории сходимости каждого перезапустили метод gmres. gmres(3) сходится во внешней итерации 5, внутренняя итерация 3 (it3 = [5, 3]), который совпал бы с внешней итерацией 6, внутренняя итерация 0, следовательно маркировка 6 на итоговой отметке деления.

figure
semilogy(1:1/3:6,rv3/norm(b),'-o');
h1 = gca;
h1.XTick = [1:1/3:6];
h1.XTickLabel = ['1';' ';' ';'2';' ';' ';'3';' ';' ';'4';' ';' ';'5';' ';' ';'6';];
title('gmres(3)')
xlabel('Iteration number');
ylabel('Relative residual');

figure
semilogy(1:1/4:3,rv4/norm(b),'-o');
h2 = gca;
h2.XTick = [1:1/4:3];
h2.XTickLabel = ['1';' ';' ';' ';'2';' ';' ';' ';'3'];
title('gmres(4)')
xlabel('Iteration number');
ylabel('Relative residual');

figure
semilogy(1:1/5:2.8,rv5/norm(b),'-o');
h3 = gca;
h3.XTick = [1:1/5:2.8];
h3.XTickLabel = ['1';' ';' ';' ';' ';'2';' ';' ';' ';' '];
title('gmres(5)')
xlabel('Iteration number');
ylabel('Relative residual');