итеративный линейный решатель для положительно определенной и плохо обусловленной матрицы

Я собираюсь предположить, что ваши проблемы возникают из-за неточного вычисления произведения матрица-вектор. (Я никогда не видел, чтобы сопряженный градиент полностью не уменьшал остаток, если только произведение матрица-вектор не было плохим. Первая итерация после перезапуска просто выполняет самый крутой спуск.)

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

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

Глядя на загруженную вами матрицу, некоторые вещи кажутся немного странными:

1. Ваша матрица K является относительно небольшой (400 x 400) плотной матрицей.
2. Ваша матрица K содержит значительное количество записей, близких к нулю, с (abs(K(i,j)) < 1.E-16*max(abs(K))).

Для матриц такого размера непосредственное вычисление факторизации Холецкого должно быть наиболее эффективным подходом. Я не уверен, почему вы говорите, что не можете этого сделать?

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

При решении таких систем разреженных линейных уравнений, как эта, важно не количество строк / столбцов в матрице, а шаблон разреженности самой матрицы.

Если, например, ваша матрица A очень разреженная, возможно, удастся вычислить разреженную факторизацию Холецкого A = L*L' напрямую. Однако будьте осторожны, порядок уравнений определяет образец разреженности результирующих факторов, и выбор плохой стратегии упорядочивания для A может привести к катастрофическому заполнению для L*L' и низкой производительности.

Существует ряд стратегий, таких как приблизительная минимальная степень и Многоуровневое вложенное рассечение, которое следует использовать для изменения порядка A для получения псевдооптимальной разреженности для L*L'.

Существует ряд хороших пакетов, которые реализуют высокопроизводительную разреженную факторизацию, включая реализации схем переупорядочения, описанных выше. Я бы порекомендовал изучить пакет CHOLMOD Дэвиса.

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

Вы должны всегда использовать хотя бы простой диагональный предварительный кондиционер Jacobi, хотя обычно гораздо лучшая производительность может быть достигнута с использованием более сложных подходов, таких как неполная факторизация Холецкого или, возможно, алгебраические многосеточные или многоуровневые методы. В этом отношении вам может быть полезна библиотека PETSc, поскольку она содержит высокопроизводительные реализации количество итерационных решателей и схем предварительной обработки.

Надеюсь это поможет.

Я видел (но не особо усмотрел) недавние работы по этому поводу, такие как http://www.cs.yale.edu/homes/spielman/precon/precon.html. Связывая сказанное с Википедией, вы можете посмотреть http://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method, который является частным случаем http://en.wikipedia.org/wiki/Successive_Over-relaxation.

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