Поправка к ранее опубликованным замечаниям к алгоритмам 42, 107, 120 и алгоритму gir

П. Hayip (Naur P. «САСМ», 1963, № 8)

Г. Форсайт (Стенфордский университет) в личной пе-)еписке информировал меня о двух главных недостатках моих замечаниях к вышеуказанным алгоритмам.

1. Полученные результаты обращения округленных датриц Гильберта сравнивались с точными результатами обращения неокругленных матриц Гильберта вместо сравнения с очень точными результатами обращения округленных матриц Гильберта.

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

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

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

W Таблица 18

I. Этот критерий соответствует критике Форсайта. Действительно, на основании этого критерия нельзя отдавать предпочтение ни процедуре INVERSION, ни процедуре gjn .

Вычисления .производились в системе GIER-ALGOL, которая .имеет плавающую запятую для чисел с 29 значащими разрядами.

АЛГОРИТМ 436

Метод Краута для решения системы линейных алгебраических ура1знений [F4]

Процедура croui является модификацией алгоритма 16 (Horsythe G. Е. «САСМ», 1960, № 9). Для решения системы уравнений ау = Ь тл цреобраэования расширенной матрицы \[аЬ\\ к ее треугольной составляющей L со всеми элементами L[k; /г] = 1 здесь использован .метод Краута с вращением [161, с. 44-100]. Если матрица а вырожденная, то происходит выход к нелокализованной -MemQ-signaL43, Переменнаg-4?it;o/[fe] принимает значение текущего индекса строки главного элемента в k-ы столбце. Таким образом, сохраняется достаточно информации, чтобы .процедура могла быть выполнена с новыми правыми частями без повторного вычисления треугольной матрицы, если логической переменной repeat задано значение true. Точность, обеспечиваемая процедурой crout, могла бы значительно новыситься при использова-кии более точной процедуры скалярного произведения, чем процедура scalar, приведенная ниже. Глобальные идентификаторы процедуры crout: 1) метка signal 43, к которой происходит выход при det=0, 2) процедура-функция scalar, приведенная ниже.

procedure crout (a,b,n,repeat) result: (y,pivot); value n; integer n; Boolean repeat; array a,b,y; integer array pivot; begin real t; integer i,k,imax,p; if repeat then go to labl; for k:=: 1 step 1 until n do begin t: = 0 for i:= к step 1 until n do begin

a[i,k]: = a[i,k]-scalar (a[i,p],a[p,k],p,l, k-1); if abs(a[i,k])>t then begin t: = absa[i, k]); imax: = i end

end; pivot)Ik]: = imax;

comment Мы нашли, что a [imax, к] является

наибольшим элементом в й-м столбце. Теперь переставим строки k и imax; if \т&хфк tlien begin

for i:= 1 step 1 until n do begin t: = a[k,i]; a[k,i]: = ai[imax,i]; a[imax,i]:=t end;

t:=:b[k]; bi[k]: = b[imax]; b[imax]: = t end перестановки строк; comment Перестановка сделана. Переходим к исключению;

if а{к,к]=0 then go to signal43; t: = 1.0/atk,k];

for i:= k+1 step 1 until n do a[i,k]:=tXa[i,k]; for i; = k+l step 1 until n do

atk.i]: = atk.i]-scalar (a[k,p],a [p.i] ,p, 1 ,k- 1); b [k]: = b [k] -scalar (а(к,р],Ыр],р, 1 ,k-l) end k; go to lab2;

comment Триангуляризация закончена, переходим к обратной подстановке. Нижеследующий оператор цикла используется, когда переменная repeat=true, т. е. дано указание .на то, что матрица а уже приведена к треугольному виду с помощью процедуры crout и нужно решать линейную систему уравнений с -новыми правыми частями без повторения триангуляризации; labl: for к:= 1 step 1 until n do

begin t: = b[pivot{k]]; b[pivot[k]]: = b [k];-

b[k]: = t;

b k]: = b[k]-scalar (a[k, p], b[p], p, 1, k-1) end;

lab2: for k:= n step -1 until 1 do

уМ: = b[k]/atk,k]-scalar (a[k,p]/a [k,k] ,y[p] ,p, -k+l,n) end crout;

Нижеследующая процедура scalar вычисляет сумму Произведений uVh для k=s, s+1,..., f. Бели s>.f, то значение указателя функции scalar равно нулю. Эта процедура может быть описана в начале любого блока, включающего блок, в котором описана процедура crout. Процедура scalar может использоваться и самостоятельно для получения скалярного произведения векторов, умножения матриц и т. д. Эта процедура не пригодна для