Проблема

Я хочу решить общую систему линейных уравнений A * x = b. Матрица m-by-m является редкой, реальной, квадратной, несколько плохо обусловленной, несимметричной, но она сингулярна (ранг (A) == m-1), поскольку x известен только до аддитивной постоянной,

Я могу создать матрицу A, указав ее ненулевые записи в трех векторах (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 \ b, MATLAB интерпретирует \ как команду mldivide() (), который проводит некоторые тесты на матрице и определяет оптимальный алгоритм для решения этой процедуры (см. ссылку для подробностей).

Эти детали сделаны подробными во время выполнения, установив параметры разрешающей матрицы MATLAB через spparms('spumoni',2)

Когда я вычисляю x и/или y, я замечаю следующее:

MATLAB использует разложение LU через UMFPACK V5.4.0 для квадратного, m-by-m исходного уравнения.
MATLAB использует QR-декомпозицию через SuiteSparseQR 1.1.0 для модифицированного уравнения m-by- (m + 1).

Оба UMFPACK и SuiteSparseQR входят в пакет программного обеспечения SuiteSparse (ссылка).

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

Моя проблема

Теперь, когда я могу решить эту проблему в MATLAB, я хочу сделать это в Fortran. К сожалению, команда MATLAB mldivide() - это черный ящик, в котором я не вижу, как он устанавливает или вызывает процедуры SuiteSparse.

Учитывая, что у меня есть три разреженных вектора в Fortran (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? (Я, вероятно, мог бы извлечь наибольшую выгоду из этой демонстрации, учитывая достаточно времени, но странно, что он вызывает процедуры несколько раз....)

Примечание. Хотя я задаю этот вопрос с конкретной проблемой (кто не является!?), я считаю, что это достаточно общее, чтобы быть полезным и для других.

Как позвонить в UMFPACK как MATLAB?

Проблема

Исходное уравнение

Модифицированное уравнение

Числовые

Моя проблема

Ответ 1