В Eigen существует несколько методов для решения линейных систем, когда матрица коэффициентов разреженная. Из-за специального представления этого класса матриц необходимо проявлять особую внимательность, чтобы обеспечить хорошую производительность. Подробное введение в разреженные матрицы в Eigen можно найти в разделе Работа с разреженными матрицами. На этой странице перечислены доступные в Eigen разреженные решатели. Также представлены основные шаги, общие для всех этих линейных решателей. В зависимости от свойств матрицы, желаемой точности пользователь может настроить эти шаги, чтобы улучшить производительность своего кода. Обратите внимание, что глубокое понимание того, что скрывается за этими шагами, не требуется: в последнем разделе представлена процедура бенчмарка, которую можно легко использовать, чтобы получить представление о производительности всех доступных решателей.

Список разреженных решателей

Eigen в настоящее время предоставляет широкий набор встроенных решателей, а также обёртки для внешних библиотек решателей. Они обобщены в следующих таблицах:

Встроенные прямые решатели

Класс	Тип решателя	Тип матрицы Matrix	Особенности, связанные с производительностью	Лицензия	Примечания
SimplicialLLT #include<Eigen/SparseCholesky>	Прямая факторизация LLt	SPD	Снижение заполнения	LGPL	SimplicialLDLT зачастую предпочтительнее
SimplicialLDLT #include<Eigen/SparseCholesky>	Прямая факторизация LDLt	SPD	Снижение заполнения	LGPL	Рекомендуется для очень разреженных и не очень больших задач (например, уравнение Пуассона в 2D)
SparseLU #include<Eigen/SparseLU>	Факторизация LU	Квадратная	Снижение заполнения, использование быстрой алгебры плотных матриц	MPL2	Оптимизировано для небольших и больших задач с нерегулярными шаблонами
SparseQR #include<Eigen/SparseQR>	Факторизация QR	Любая, прямоугольная	Снижение заполнения	MPL2	рекомендуется для задач наименьших квадратов, имеет базовый признак определения ранга

Встроенные итеративные решатели

Класс	Тип решателя	Тип матрицы Matrix	Поддерживаемые прекондиционеры, [по умолчанию]	Лицензия	Примечания
ConjugateGradient #include<Eigen/IterativeLinearSolvers>	Классический итеративный CG	SPD	IdentityPreconditioner, [DiagonalPreconditioner], IncompleteCholesky	MPL2	Рекомендуется для больших симметричных задач (например, уравнение Пуассона в 3D)
LeastSquaresConjugateGradient #include<Eigen/IterativeLinearSolvers>	CG для прямоугольной задачи наименьших квадратов	Прямоугольная	IdentityPreconditioner, [LeastSquareDiagonalPreconditioner]	MPL2	Solve для min |A'Ax-b|^2 без формирования A'A
BiCGSTAB #include<Eigen/IterativeLinearSolvers>	Итеративный стабилизированный бисопряжённый градиент	Квадратная	IdentityPreconditioner, [DiagonalPreconditioner], IncompleteLUT	MPL2	Для ускорения сходимости попробуйте использовать прекондиционер IncompleteLUT.

Обёртки для внешних решателей

Класс	Модуль	Тип решателя	Тип матрицы Matrix	Особенности, связанные с производительностью	Зависимости, лицензия	Примечания
PastixLLT PastixLDLT PastixLU	PaStiXSupport	Прямые факторизации LLt, LDLt, LU	SPD SPD Квадратная	Снижение заполнения, использование быстрой алгебры плотных матриц, многопоточность	Требуется пакет PaStiX, CeCILL-C	Оптимизировано для сложных задач и симметричных шаблонов
CholmodSupernodalLLT	CholmodSupport	Прямая факторизация LLt	SPD	Снижение заполнения, использование быстрой алгебры плотных матриц	Требуется пакет SuiteSparse, GPL
UmfPackLU	UmfPackSupport	Прямая факторизация LU	Квадратная	Снижение заполнения, использование быстрой алгебры плотных матриц	Требуется пакет SuiteSparse, GPL
KLU	KLUSupport	Прямая факторизация LU	Квадратная	Снижение заполнения, подходит для моделирования схем	Требуется пакет SuiteSparse, GPL
SuperLU	SuperLUSupport	Прямая факторизация LU	Квадратная	Снижение заполнения, использование быстрой алгебры плотных матриц	Требуется библиотека SuperLU, (BSD-подобная)
SPQR	SPQRSupport	Факторизация QR	Любая, прямоугольная	снижение заполнения, многопоточность, быстрая алгебра плотных матриц	требуется пакет SuiteSparse, GPL	рекомендуется для задач линейных наименьших квадратов, имеет функцию определения ранга
PardisoLLT PardisoLDLT PardisoLU	PardisoSupport	Прямые факторизации LLt, LDLt, LU	SPD SPD Квадратная	Снижение заполнения, использование быстрой алгебры плотных матриц, многопоточность	Требуется пакет Intel MKL, Проприетарная	Оптимизировано для сложных задач, шаблонов, см. также использование MKL с Eigen

Здесь SPD означает симметрично-положительно определённую.

Концепция разреженного решателя

Все эти решатели следуют одной общей концепции. Вот типичный и общий пример:

#include <Eigen/RequiredModuleName>
// ...
SparseMatrix<double> A;
// fill A
VectorXd b, x;
// fill b
// solve Ax = b
SolverClassName<SparseMatrix<double> > solver;
solver.compute(A);
if(solver.info()!=Success) {
  // decomposition failed
  return;
}
x = solver.solve(b);
if(solver.info()!=Success) {
  // solving failed
  return;
}
// solve for another right hand side:
x1 = solver.solve(b1);

Для SPD решателей второй необязательный шаблонный аргумент позволяет указать, какая треугольная часть должна быть использована, например:

#include <Eigen/IterativeLinearSolvers>
 
ConjugateGradient<SparseMatrix<double>, Eigen::Upper> solver;
x = solver.compute(A).solve(b);

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

В случае, когда необходимо решить несколько задач с одинаковым шаблоном разреженности, шаг "вычисления" можно разложить следующим образом:

SolverClassName<SparseMatrix<double> > solver;
solver.analyzePattern(A);   // for this step the numerical values of A are not used
solver.factorize(A);
x1 = solver.solve(b1);
x2 = solver.solve(b2);
...
A = ...;                    // modify the values of the nonzeros of A, the nonzeros pattern must stay unchanged
solver.factorize(A);
x1 = solver.solve(b1);
x2 = solver.solve(b2);
...

Метод compute() эквивалентен вызову как analyzePattern(), так и factorize().

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

Наконец, большинство итеративных решателей также могут быть использованы в контексте без матрицы, см. пример по этой ссылке .

Шаг вычисления

В функции compute() матрица обычно разлагается: LLT для самосопряжённых матриц, LDLT для общих эрмитовых матриц, LU для неэрмитовых матриц и QR для прямоугольных матриц. Это результаты использования прямых решателей. Для этого класса решателей этап compute дополнительно подразделяется на analyzePattern() и factorize().

Цель analyzePattern() — переупорядочить ненулевые элементы матрицы таким образом, чтобы этап факторизации создавал меньше заполнения. Этот этап использует только структуру матрицы. Следовательно, результаты этого этапа могут быть использованы для других систем линейных уравнений, где матрица имеет ту же структуру. Однако обратите внимание, что иногда некоторые внешние решатели (например, SuperLU) требуют, чтобы значения матрицы были заданы на этом этапе, например, для выравнивания строк и столбцов матрицы. В этой ситуации результаты этого этапа не должны использоваться с другими матрицами.

Eigen предоставляет ограниченный набор методов для переупорядочения матрицы на этом этапе, либо встроенные (COLAMD, AMD), либо внешние (METIS). Эти методы задаются в списке параметров шаблона решателя:

DirectSolverClassName<SparseMatrix<double>, OrderingMethod<IndexType> > solver;

См. модуль OrderingMethods для списка доступных методов и связанных опций.

В factorize() вычисляются факторы матрицы коэффициентов. Этот этап должен вызываться каждый раз, когда значения матрицы меняются. Однако структурная структура матрицы не должна меняться между несколькими вызовами.

Для итеративных решателей этап compute используется для создания предварительного условия. Например, с предварительным условием ILUT неполные факторы L и U вычисляются на этом этапе. Помните, что основная цель предварительного условия — ускорить сходимость итерационного метода путём решения модифицированной системы линейных уравнений, где матрица коэффициентов имеет более сгруппированные собственные значения. Для реальных задач итерационный решатель всегда должен использоваться с предварительным условием. В Eigen предварительное условие выбирается просто добавлением его как параметра шаблона к объекту итерационного решателя.

IterativeSolverClassName<SparseMatrix<double>, PreconditionerName<SparseMatrix<double> > solver; 

Член-функция preconditioner() возвращает ссылку на предварительное условие для прямого взаимодействия с ним. См. модуль итеративных решателей и документацию каждого класса для списка доступных методов.

Этап решения

Функция solve() вычисляет решение систем линейных уравнений с одним или несколькими правыми частями.

X = solver.solve(B);

Здесь B может быть вектором или матрицей, где столбцы образуют различные правые части. Функция solve() также может вызываться несколько раз, например, когда все правые части недоступны сразу.

x1 = solver.solve(b1);
// Get the second right hand side b2
x2 = solver.solve(b2); 
//  ...

Для прямых методов решение вычисляется с точностью до машинной. Иногда решение не должно быть слишком точным. В этом случае итерационные методы более подходят, и желаемая точность может быть установлена перед шагом решения с помощью setTolerance(). Для всех доступных функций см. документацию модуля итеративных решателей.

Процедура тестирования производительности

Большую часть времени всё, что вам нужно знать, это сколько времени займёт решение вашей системы и, надеюсь, какой решатель является наиболее подходящим. В Eigen мы предоставляем процедуру тестирования производительности, которую можно использовать для этой цели. Её очень легко использовать. В каталоге сборки перейдите в bench/spbench и скомпилируйте процедуру, набрав make spbenchsolver. Запустите её с опцией –help, чтобы получить список всех доступных опций. В принципе, тестируемые матрицы должны быть в формате MatrixMarket Coordinate, а процедура возвращает статистику от всех доступных решателей в Eigen.

Для экспорта ваших матриц и векторов правой части в формате MatrixMarket можно использовать неподдерживаемый модуль SparseExtra:

#include <unsupported/Eigen/SparseExtra>
...
Eigen::saveMarket(A, "filename.mtx");
Eigen::saveMarket(A, "filename_SPD.mtx", Eigen::Symmetric); // if A is symmetric-positive-definite
Eigen::saveMarketVector(B, "filename_b.mtx");

Следующая таблица даёт пример XML-статистики нескольких встроенных и внешних решателей Eigen.

Matrix	N	NNZ		UMFPACK	SUPERLU	PASTIX LU	BiCGSTAB	BiCGSTAB+ILUT	GMRES+ILUT	LDLT	CHOLMOD LDLT	PASTIX LDLT	LLT	CHOLMOD SP LLT	CHOLMOD LLT	PASTIX LLT	CG
vector_graphics	12855	72069	Время вычисления	0.0254549	0.0215677	0.0701827	0.000153388	0.0140107	0.0153709	0.0101601	0.00930502	0.0649689
			Solve Time	0.00337835	0.000951826	0.00484373	0.0374886	0.0046445	0.00847754	0.000541813	0.000293696	0.00485376
			Общее время	0.0288333	0.0225195	0.0750265	0.037642	0.0186552	0.0238484	0.0107019	0.00959871	0.0698227
			Ошибка(Итерации)	1.299e-16	2.04207e-16	4.83393e-15	3.94856e-11 (80)	1.03861e-12 (3)	5.81088e-14 (6)	1.97578e-16	1.83927e-16	4.24115e-15
poisson_SPD	19788	308232	Время вычисления	0.425026	1.82378	0.617367	0.000478921	1.34001	1.33471	0.796419	0.857573	0.473007	0.814826	0.184719	0.861555	0.470559	0.000458188
			Solve Time	0.0280053	0.0194402	0.0268747	0.249437	0.0548444	0.0926991	0.00850204	0.0053171	0.0258932	0.00874603	0.00578155	0.00530361	0.0248942	0.239093
			Общее время	0.453031	1.84322	0.644241	0.249916	1.39486	1.42741	0.804921	0.862891	0.4989	0.823572	0.190501	0.866859	0.495453	0.239551
			Ошибка(Итерации)	4.67146e-16	1.068e-15	1.3397e-15	6.29233e-11 (201)	3.68527e-11 (6)	3.3168e-15 (16)	1.86376e-15	1.31518e-16	1.42593e-15	3.45361e-15	3.14575e-16	2.21723e-15	7.21058e-16	9.06435e-12 (261)
sherman2	1080	23094	Время вычисления	0.00631754	0.015052	0.0247514	-	0.0214425	0.0217988
			Solve Time	0.000478424	0.000337998	0.0010291	-	0.00243152	0.00246152
			Общее время	0.00679597	0.01539	0.0257805	-	0.023874	0.0242603
			Ошибка(Итерации)	1.83099e-15	8.19351e-15	2.625e-14	1.3678e+69 (1080)	4.1911e-12 (7)	5.0299e-13 (12)
bcsstk01_SPD	48	400	Время вычисления	0.000169079	0.00010789	0.000572538	1.425e-06	9.1612e-05	8.3985e-05	5.6489e-05	7.0913e-05	0.000468251	5.7389e-05	8.0212e-05	5.8394e-05	0.000463017	1.333e-06
			Solve Time	1.2288e-05	1.1124e-05	0.000286387	8.5896e-05	1.6381e-05	1.6984e-05	3.095e-06	4.115e-06	0.000325438	3.504e-06	7.369e-06	3.454e-06	0.000294095	6.0516e-05
			Общее время	0.000181367	0.000119014	0.000858925	8.7321e-05	0.000107993	0.000100969	5.9584e-05	7.5028e-05	0.000793689	6.0893e-05	8.7581e-05	6.1848e-05	0.000757112	6.1849e-05
			Ошибка(Итерации)	1.03474e-16	2.23046e-16	2.01273e-16	4.87455e-07 (48)	1.03553e-16 (2)	3.55965e-16 (2)	2.48189e-16	1.88808e-16	1.97976e-16	2.37248e-16	1.82701e-16	2.71474e-16	2.11322e-16	3.547e-09 (48)
sherman1	1000	3750	Время вычисления	0.00228805	0.00209231	0.00528268	9.846e-06	0.00163522	0.00162155	0.000789259	0.000804495	0.00438269
			Solve Time	0.000213788	9.7983e-05	0.000938831	0.00629835	0.000361764	0.00078794	4.3989e-05	2.5331e-05	0.000917166
			Общее время	0.00250184	0.00219029	0.00622151	0.0063082	0.00199698	0.00240949	0.000833248	0.000829826	0.00529986
			Ошибка(Итерации)	1.16839e-16	2.25968e-16	2.59116e-16	3.76779e-11 (248)	4.13343e-11 (4)	2.22347e-14 (10)	2.05861e-16	1.83555e-16	1.02917e-15
young1c	841	4089	Время вычисления	0.00235843	0.00217228	0.00568075	1.2735e-05	0.00264866	0.00258236
			Solve Time	0.000329599	0.000168634	0.00080118	0.0534738	0.00187193	0.00450211
			Общее время	0.00268803	0.00234091	0.00648193	0.0534865	0.00452059	0.00708447
			Ошибка(Итерации)	1.27029e-16	2.81321e-16	5.0492e-15	8.0507e-11 (706)	3.00447e-12 (8)	1.46532e-12 (16)
mhd1280b	1280	22778	Время вычисления	0.00234898	0.00207079	0.00570918	2.5976e-05	0.00302563	0.00298036	0.00144525	0.000919922	0.00426444
			Solve Time	0.00103392	0.000211911	0.00105	0.0110432	0.000628287	0.00392089	0.000138303	6.2446e-05	0.00097564
			Общее время	0.0033829	0.0022827	0.00675918	0.0110692	0.00365392	0.00690124	0.00158355	0.000982368	0.00524008
			Ошибка(Итерации)	1.32953e-16	3.08646e-16	6.734e-16	8.83132e-11 (40)	1.51153e-16 (1)	6.08556e-16 (8)	1.89264e-16	1.97477e-16	6.68126e-09
crashbasis	160000	1750416	Время вычисления	3.2019	5.7892	15.7573	0.00383515	3.1006	3.09921
			Solve Time	0.261915	0.106225	0.402141	1.49089	0.24888	0.443673
			Общее время	3.46381	5.89542	16.1594	1.49473	3.34948	3.54288
			Ошибка(Итерации)	1.76348e-16	4.58395e-16	1.67982e-14	8.64144e-11 (61)	8.5996e-12 (2)	6.04042e-14 (5)