Программа метода сопряженных градиентов с предобуславливанием

Программа  учебная на Фортране 77.  https://www.mediafire.com/file/yw7p6zsqlimrj5e/example+(1).zip      фаол исходных данных  *.dat

Решается система линейных уравнений Ах=b. Где A — матрица 3x3,   b — матрица столбец правой части.

В программе использyется матрица А (4х3). В ней 4-й столбец это матрица b.В итоге для исходных данных

А=  0.40000000E+01  0.30000000E+01  0.00000000E+00  0.24000000E+02  

0.30000000E+01  0.40000000E+01 -0.10000000E+01  0.30000000E+02  

0.00000000E+00 -0.10000000E+01  0.40000000E+01 -0.24000000E+02

Исходная система линейных уравнений преобразуется в систему C*Аx=C*b. Матрица С (3х3) называется предобуславливателем. Для этого примера она не выбирается как  C=А-1. Матрица С — это не обратная матрица А, перед решением задачи должна быть подсчитана.  C = 

0.5   0   0

0   0.5  0

0   0  0.5

Матрица С — всегда выбирается как только диагональная. Элементы ее диагонали — это соответствующий элемент матрицы А, в нашем случае (3Х3). Он поделен на единицу, и затем из него взят квадратный корень. Все недиагональные элементы матрицы С равны нулю.

1/sqrt(4)=1/2=0.5

 Можно попробовать еще уменьшить величины невязок. Поделите каждыи элемент матрицы А  на 1000000 и запишите в качестве исходной матрицы А(3х4), а элементы матрицы С умножьте на  300.

Получим первый элемент матрицы А равен не 4, а 4е-6. А первый элемент матрицы С равен 0,5*300=150. Все остальные элементы диагональной матрицы С равны этой величине.

Пoсмотрите на невязку, если С(1)  = 50000. И при этой величине увеличьте элементы матрицы А до 4е-4. Надо пробовать разные варианты.

Задание начального значения вектора Х вместо (0,0,0 )  на (2,3,-4) (они на 1 отличаются от  точного ответа)  не меняет сильно число итераций. Иногда величина невязки становится больше. 

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

Три начальных значений вектора х для первой итерации задано как:

0.00000000E+00  0.00000000E+00  0.00000000E+00

В итоге, через 3 итерации получаем ответ. 0.29999981E+01  0.39999995E+01 -0.49999990E+01

Точный ответ (3; 4; -5)

Значение невязок для каждого вектора х1, х2, х3 равно на 3-й итерации 

RESIDUAL IS   0.80093741E-05  0.81062317E-05 -0.16540289E-05

Текст программы на фортране 77, файл исходных данных, файл ответа  CON1.txt скачайте по ссылке.

https://ru.files.fm/u/56hyj7mfr9

или тут

https://fex.net/s/vvdtt90

Файл readme.txt -  дана команда в Linux для компиляции, файл  а.out  -  исполняемый файл программы в Линуксе. Если есть вопросы, тo задавайте.

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

Иногда метод решает слабо несимметричные матрицы А. Если число обyсловленности большое, метод Гаусса не работает. Метод сопряженных градиентов решает точно и быстро.

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

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