Как позвонить в UMFPACK, как это делает MATLAB?
Эта проблема
Я хочу решить общую систему линейных уравнений A*x=b. Матрица m-by-m является разреженной, вещественной, квадратной, несколько плохо обусловленной, несимметричной, но она является единственной (rank(A)==m-1), поскольку x известен только с точностью до аддитивной константы.
Я могу создать матрицу А, указав ее ненулевые записи в трех векторах (i, j, а также v) такой что A(i(k),j(k)) = v(k):
A = sparse( i, j, v, m, m );
Исходное уравнение
Я могу решить это оригинальное уравнение следующим образом:
x = A \ b;
Если мне нужно уникальное решение, я могу наложить ограничение (скажем, x(4) == 3.14159) после вычисления неуникального решения:
x = x - x(4) + 3.14159;
Модифицированное уравнение
Я могу сделать новую, полноценную матрицу C добавив дополнительное ограничение уникальности, следующим образом:
% Add the constraint x(4) == 3.14159
extraRow = zeros(1,m);
extraRow(4) = 1.0;
C = [A; extraRow]; % Add to the matrix A
d = [b; 3.14159]; % Add to the RHS vector, b
% Solve C*y = d for y
y = C \ d;
Числовые
Я понимаю, что когда я решаю эти уравнения с помощью x = A \ b или же y = C \ bMATLAB интерпретирует \ как mldivide() команда ( ссылка), которая запускает некоторые тесты на матрице и определяет оптимальный алгоритм для решения этой процедуры (подробности см. в ссылке).
Эти подробности становятся подробными во время выполнения, устанавливая параметры разреженной матрицы MATLAB через spparms('spumoni',2)
Когда я вычисляю x и / или yЯ заметил следующее:
- MATLAB использует декомпозицию LU через UMFPACK V5.4.0 для квадратного, м-м-м исходного уравнения.
- MATLAB использует QR-разложение через SuiteSparseQR 1.1.0 для модифицированного уравнения m-by-(m+1).
UMFPACK и SuiteSparseQR являются частью программного пакета SuiteSparse ( ссылка).
(Несколько неожиданно, решение модифицированного уравнения дает больше ошибок, чем исходное уравнение. Несмотря на значимость, эта ошибка все еще находится на приемлемо низком пороге.)
Моя проблема
Теперь, когда я могу решить эту проблему в MATLAB, я хочу сделать это в Fortran. К сожалению, MATLAB's mldivide() Команда - это черный ящик, в котором я не вижу, как она устанавливает или вызывает процедуры SuiteSparse.
Учитывая, что у меня есть три разреженных вектора в Фортране (90+), как показано ниже, как я могу решить проблему с помощью SuiteSparse?
Кроме того, кто-нибудь знает какие-либо процедуры-оболочки F90, которые взаимодействуют с UMFPACK, чтобы сделать это проще?
Я более чем счастлив, если кто-нибудь может помочь с этим - либо с исходным уравнением, либо с измененным уравнением. (Если вы поможете с одним, я мог бы получить другой.)
subroutine solveSparseMatrixEqnViaSuiteSparse( m, n, nnz, i, j, v, x )
implicit none
integer, intent(in) :: m ! sparse matrix rows
integer, intent(in) :: n ! sparse matrix columns
integer, intent(in) :: nnz ! number of nonzero entries
integer, dimension(1:nnz), intent(in) :: i ! row indices of nonzero entries
integer, dimension(1:nnz), intent(in) :: j ! column indices of nonzero entries
real*8, dimension(1:nnz), intent(in) :: v ! values of nonzero entries
real*8, dimension(1:n), intent(out) :: x ! solution vector
! I have no idea what to do next!
end function solveSparseMatrixEqnViaSuiteSparse
Что меня смущает, так это следующее:
- Что делает MATLAB за кулисами, чтобы настроить вызов SuiteSparse? (Кажется, это не задокументировано....)
- Что нужно сделать, чтобы вызвать SuiteSparse из Fortran? (Я мог бы, вероятно, понять больше всего из этой демонстрации, если бы у меня было достаточно времени, но странно, что он вызывает процедуры несколько раз....)
Примечание. Хотя я задаю этот вопрос с учетом конкретной проблемы (а кто нет!?), я считаю, что она достаточно общая, чтобы быть полезной и для других.
2 ответа
Я могу ответить на часть вашего вопроса. Я аспирант и работаю над UMFPACK под руководством доктора Тима Дэвиса, и мы не используем фортран.
MATLAB использует свой собственный шаг, чтобы увидеть, использовать ли UMFPACK или нет. На самом деле, он начинается с факторизации Холецкого, если матрица имеет симметричную структуру, затем он выбирает между простым левым взглядом или UMFPACK.
Я не знаю, предоставил ли доктор Дэвис интерфейс Fortran для самого UMFPACK. Мы также не используем фортран в нашей группе. Однако я знаю, что вы можете использовать интерфейс SuiteSparse matlab. Перейдите в /SuiteSparse/UMFPACK/MATLAB, и вы можете использовать различные интерфейсы, и вам не нужно быть знакомым с его базовым кодом C. Вы можете вызвать любую функцию напрямую. Реализована функция MATLAB mex, и вы можете отлаживать и следовать коду, если хотите.
Здесь есть оболочка для Фортрана для UMFPACK [ http://geo.mff.cuni.cz/~lh/Fortran/UMFPACK/README.html ], но она старая и я не думаю, что она поддерживает последнюю версию UMFPACK
Здесь две вещи. Во-первых, для примеров вызова UMFPACK в fortran/C ознакомьтесь с исходным кодом в одном из моих проектов, например:
https://github.com/fangq/blit/blob/master/src/blit_ilupcond.f90#L156-L162
здесь я вызываю umfpack из fortran90, используя модуль umf4_f77wrapper.c
https://github.com/fangq/blit/blob/master/src/umf4_f77wrapper.c
Во-вторых, что касается оператора обратной косой черты, то он, по сути, решает задачу псевдообращения Мура-Пенроуза, чтобы получить решение методом наименьших квадратов, см.
https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_inverse
по сути, любое матричное уравнение A*x=b (включая недоопределенные - size(A,1)<size(A,2) - или переопределенные - size(A,1)>size(A,2) - проблемы), решениеx=A\bэто то же самое, что решить
(A'A) x = A'b=>x=inv(A'A)*A'b
однако вы не решаете правое уравнение путем инверсии, а применяете линейный решатель, решая левое уравнение.(A'A)x=A'b
вы также можете применить формулу Шермана-Моррисона-Вудбери для построения эквивалентного, но более компактного решения, когда исходное уравнение недоопределено, чтоx=A'*inv(A'A)*b