Задайте свою Систему Линейных уравнений с помощью Разреженной матрицы

Когда вы используете distributed функция, MATLAB автоматически запускает параллельный пул с помощью кластерных настроек по умолчанию. Этот пример использует матрицу Wathen от gallery MATLAB функция. Эта матрица является разреженной, симметричной, и случайной матрицей с габаритным размером $$N=3n^2+4n+1$$.

n = 400;
A = distributed(gallery('wathen',n,n));

Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 6).

N = 3*n^2+4*n+1

N = 481601

Можно теперь задать правый вектор $$b$$. В этом примере, $$b$$ задан как сумма строки $$A$$, который приводит к точному решению $$Ax=b$$ из формы $$x_{\text{exact}}=[1,...,1{]}^T$$.

b = sum(A,2);

Начиная с sum действия на распределенном массиве, b также распределяется и его данные хранятся в памяти о рабочих вашего параллельного пула. Наконец, можно задать точное решение для сравнения с решениями, полученными с помощью итерационных методов.

xExact = ones(N,1,'distributed');

Решите свою Систему Линейных уравнений с Методом сопряженных градиентов

Функция MATLAB pcg предоставляет метод метода сопряженных градиентов (CG), который итеративно генерирует серию приближенных решений для $\mathit{x}$, улучшение решения с каждым шагом.

[xCG_1,flagCG_1,relres_CG1,iterCG_1,resvecCG_1] = pcg(A,b);

Когда система решена, можно проверять ошибку между каждым элементом полученного результата xCG_1 и ожидаемые значения xExact. Ошибка в вычисленном результате относительно высока.

errCG_1 = abs(xExact-xCG_1);

figure(1)
hold off
semilogy(errCG_1,'o');
title('System of Linear Equations with Sparse Matrix');
ylabel('Absolute Error');
xlabel('Element in x');

Итеративный расчет заканчивается, когда серия приближенных решений сходится к определенному допуску или после максимального количества шагов итерации. И для распределенных и для массивов на клиенте, $$pcg$$ использование те же настройки по умолчанию:

Как второй выходной аргумент, pcg функционируйте также возвращает флаг сходимости, который дает вам больше информации о полученном результате, включая то, сходилось ли вычисленное решение к желаемому допуску. Например, значение 0 указывает, что решение правильно сходилось.

flagCG_1

flagCG_1 = 1

В этом примере решение не сходится в максимальном количестве по умолчанию итераций, которое приводит к высокой ошибке.

Чтобы увеличить вероятность сходимости, можно настроить настройки для допуска и максимального количества шагов итерации.

tolerance = 1e-12;
maxit = N;

tCG = tic;
    [xCG_2,flagCG_2,relresCG_2,iterCG_2,resvecCG_2] = pcg(A,b,tolerance,maxit);
tCG = toc(tCG);

flagCG_2

flagCG_2 = 0

С пользовательскими настройками сходится решение. Это решение имеет улучшенную абсолютную погрешность по сравнению с предыдущим решением.

errCG_2 = abs(xExact-xCG_2);
figure(2)
hold off
semilogy(errCG_1,'o');
hold on
semilogy(errCG_2,'d');
title('Comparison of Absolute Error');
ylabel('Absolute Error');
xlabel('Element in x');
legend('Default tolerance and iterations','Improved tolerance and iterations');
hold off

pcg метод также возвращает вектор из нормы невязки на каждом шаге итерации, norm(b-A*x)/norm(b). Относительная норма невязки показывает отношение точности между последовательными шагами итерации. Эволюция остаточных значений во время итеративного процесса может помочь вам изучить, почему решение не сходилось без пользовательских настроек.

figure(3)
f=semilogy(resvecCG_2./resvecCG_2(1));
hold on
semilogy(f.Parent.XLim,[1e-6 1e-6],'--')
semilogy([20 20], f.Parent.YLim,'--')
semilogy(f.Parent.XLim,[1e-12 1e-12],'--')
title('Evolution of Relative Residual');
ylabel('Relative Residual');
xlabel('Iteration Step');
legend('Residuals of CG','Default Tolerance','Default Number of Steps','Custom Tolerance')
hold off

Ясно, что количества по умолчанию шагов недостаточно, чтобы достигнуть хорошего решения для этой системы.

Решите свою Систему Линейных уравнений с Предобусловленным Методом сопряженных градиентов

Можно повысить эффективность решения системы с помощью метода предобусловленного метода сопряженных градиентов (PCG). Во-первых, предусловие ваша система линейных уравнений с помощью матрицы перед формирователем $\mathit{M}$. Затем решите свою предобусловленную систему с помощью метода CG. Метод PCG может взять намного меньше итераций, чем метод CG.

Функция MATLAB pcg также используется для метода PCG. Можно предоставить подходящую матрицу перед формирователем $\mathit{M}$как дополнительный вход.

Идеальная матрица перед формирователем является матрицей чья инверсия ${\mathit{M}}^{-1}$ близкое приближение к инверсии матрицы коэффициентов, ${\mathit{A}}^{-1}$, но легче вычислить. Этот пример использует диагональ $$A$$ к предусловию система линейных уравнений.

M = spdiags(spdiags(A,0),0,N,N);
tPCG = tic;
    [xPCG,flagPCG,relresPCG,iterPCG,resvecPCG]=pcg(A,b,tolerance,maxit,M);
tPCG = toc(tPCG);

figure(4)
hold off;
semilogy(resvecCG_2./resvecCG_2(1))
hold on;
semilogy(resvecPCG./resvecPCG(1))
title('Evolution of Relative Residual');
ylabel('Relative Residual');
xlabel('Iteration Step');
legend('Residuals of CG','Residuals of PCG with M \approx diag(A)')

Предыдущий рисунок показывает, что для метода PCG нужно решительно меньше шагов, чтобы сходиться по сравнению с непредобусловленной системой. Этот результат также отражается во времена выполнения.

fprintf([...
    '\nTime to solve system with CG:  %d s', ...
    '\nTime to solve system with PCG: %d s'],tCG,tPCG);

Time to solve system with CG:  1.244593e+01 s
Time to solve system with PCG: 7.657432e-01 s

А также решая эту систему в качестве примера на меньшем количестве шагов итерации, метод PCG также возвращает более точное решение.

errPCG = abs(xExact-xPCG);
figure(5)
hold off
semilogy(errCG_1,'o');
hold on
semilogy(errCG_2,'d');
semilogy(errPCG,'x');
title('Comparison of absolute error');
ylabel('Absolute error');
xlabel('Element in x');
legend('CG default','CG custom','PCG');

После того, как вы будете сделаны с вашими расчетами, можно удалить параллельный пул. gcp функция возвращает текущий параллельный объект пула, таким образом, можно удалить текущий пул.

delete(gcp('nocreate'))

Матрица Wathen, используемая в этом примере, является хорошей демонстрацией того, как хороший предварительный формирователь может существенно повысить эффективность решения. Матрица Wathen имеет относительно маленькие недиагональные компоненты, таким образом выбирая$$M=\text{diag}(A)$$ дает подходящий предварительный формирователь. Для произвольной матрицы $\mathit{A}$, нахождение предварительного формирователя не может быть таким образом прямо.

Для примера того, как аппроксимировать дифференциальное уравнение линейной системой и решить его с помощью распределенного итеративного решателя с многосеточным предварительным формирователем, смотрите, Решают Дифференциальное уравнение Используя Многосеточный Предварительный формирователь на Распределенной Дискретизации.