1 SUBROUTINE SOLVE(N,A,JPVT,B)

3 С SOLVES LINEAR SYSTEM, А*Х = В

4 С ASSUMES А IS FACTORISED INTO L.U FORM (BY FACT)

5 с RETURNS SOLUTION, X, IN В

7 DIMENSION A(50,50),JPVT(50),B(50)

9 С FORWARD ELIMINATION

10 С

11 NM = N - 1

12 DO 2 К s 1,NM

13 KP К + 1

14 L JPVT(K)

15 S B(L)

16 B(L) B(K)

17 B(K) S

18 DO 1 I = KP,N

19 B(I) = B(I) + A(I,K)*S

20 1 CONTINUE

21 2 CONTINUE

22 С

23 С BACK SUBSTITUTION

24 С

25 DO 4 KA = 1,NM

26 KM N - KA

27 К = KM + 1

28 B(K) B(K)/A(K,K)

29 S = - B(K)

30 DO 3 I = 1,KM

31 B(I) B(I) + A(I,K)*S

32 3 CONTINUE

33 4 CONTINUE

34 B(l) = B(l)/A(l,l)

35 RETURN

36 END

Рис. 6.16. Распечатка подпрограммы SOLVE.

наибольший ведущий элемент ai,k (FACT, строки 17-24). Если наибольший ведущий элемент равен нулю, то программа выдает сигнал ошибки [JPVT(N)=-1] и осуществляется возврат к основной программе.

В предположении о том, что наибольший ведущий элемент не равен нулю, осуществляется перестановка рядов (FACT строки 39-41), чтобы переместить наибольший ведущий элемент ai,k в положение (й, k). В каждом столбце / (/ = Л+ 1 N) коэффициенты t-й строки (t = fe+l, модифицируются (FACT, строки 38-46) в соответствии с формулой

Для экономии памяти решение V записывается на место использованного массива и выдается уже как массив В.

Реализация исключения по Гауссу в два этапа, соответствующих подпрограммам FACT и SOLVE, полезна в тех случаях, когда требуется строить решения для многих вариантов величины В в правой части уравнения (6.23). Это связано с тем, что для выполнения подпрограммы FACT требуется O(N) операций, тогда как SOLVE требует O(N) операций. Обе эти подпрограммы, представленные в несколько расширенной форме, подвергаются более тщательному обсуждению в книге [Forsythe et al., 1977, разд. 3.3].

Как правило, при использовании спектральных методов (§ 5.6) и панельных методов (§ 14.1) получаются матрицы А с плотной структурой. В противоположность этому методы конечных элементов и конечных разностей приводят к матрицам А разреженной структуры и, если используется методика расщепления (факторизации, § 8.2), матрицам узколенточной структуры.

Если матрица А разреженная, но не обязательно ленточная, то исключение по Гауссу часто осуществляется при введении А в память в форме одного массива, содержащего значения а-/ вместе со связанным с ним индикаторным массивом IA, обеспечивающим информацию о положении точки (f, /). Основная трудность, возникающая при программировании процесса гаус-совского исключения с разреженными матрицами, связана с вероятностью заполнения матрицы А. Иначе говоря, в процессе

Массивы хранятся в памяти вычислительных машин по столбцам (индекс i в а/,/). Поэтому наиболее целесообразно сделать так, чтобы самый внутренний цикл подпрограммы FACT (строки 43-45) реализовался по первому индексу.

Подпрограмма SOLVE прежде всего преобразует элементы вектора В (строки 12-21) с целью взаимодействия с ведущей строкой и модификации элемента bi, i=fe-f 1, Л, путем вычитания из него ведущей строки с соответствующим множителем, т. е.

&. = &/-()&ь (6.26)

что эквивалентно операции (6. 25) в подпрограмме FACT. Окончательное решение получается с помощью последовательных подстановок, реализуемых посредством подпрограммы SOLVE (строки 25-34). Это дает

исключения элементы матрицы А, первоначально равные нулю, становятся отличными от нуля. Различные методы надлежащего представления А и IA обсуждаются в книгах [Jennings, 1977а, гл. 5; Duff, 1981], последняя из которых предлагает более высокий уровень анализа. Существуют пакеты программ для реализации гауссовского исключения при разреженных матрицах; можно рекомендовать разработанные в Харуэлле программы МА28 и т. д., рассмотренные в вышеназванной книге Даффа.

6.2,2, Трехдиагональные системы: алгоритмы Томаса

Использование трехточечных конечно-разностных формул или конечных элементов с линейной интерполяцией приводит после расщепления (§ 8.2) к появлению матриц А в уравнении (6.23), имеющих трехдиагональную структуру. Использование конечно-разностных схем более высокого порядка или конечных элементов более высокого порядка приводит к ленточной структуре А, где ширина области ненулевых элементов будет больше. Алгоритм Томаса пригоден для решения уравнения (6.23) в случае трехдиагональной матрицы А. Обобщение алгоритма Томаса на случай пятидиагональной матрицы А описывается в п. 6.2.4. Если речь идет о системах уравнений, то матрица А имеет, как правило, блочную (трехдиагональную) структуру. Порядок действий в этом случае рассматривается в п. 6.2.5.

Если ненулевые элементы располагаются вблизи главной диагонали, то полезно рассмотреть те варианты исключения по Гауссу, которые используют преимущества ленточной структуры А. Один из примеров такого варианта соответствует рассмотренной в § 9.3 задаче о конвекции - диффузии. Если воспользоваться формулами с центральными разностями, то при обозначениях данного параграфа получается следующий алгоритм:

- (1 + 0.5?cell) Vi + 2Vi - (1 - 0.5?сеи) = О, (6.27)

при повторении которого по отношению к каждому узлу получим

(h 2 2

at bi Ci

(6.28)