Используя 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');