4.7. Умножение уплотненной симметрической матрицы на прямоугольную

Произведение матрицы

A=||a_ij|| (i = 1, ..., m; j = 1, ..., n)

на матрицу

B=||b_ij || (i = 1, ..., n; j = 1, ..., p)

есть матрица

C=||c_ij|| (i = 1, ..., m; j = 1, ...,p),

эле­мен­ты которой определяются фор­му­лой

, i=1, ..., m; j = 1,..., p.(1.50)

Из равенства (1.50) следует, что операция пе­­­ре­мно­жения матриц определена только тогда, ког­да чис­­­ло столбцов матрицы Aравно количеству строк мат­ри­цыB. В случае если одна из пе­ре­мно­жа­е­мых мат­­риц, на­при­мер матрицаA,является сим­мет­­ри­ч­­еской, т.е.a_ij = a_ji, то она мо­­жет быть уплот­не­на пу­тем сведения к почти тре­у­голь­ной мат­рице (верх­­няя треугольная часть матрицыA,вклю­чая глав­­ную ди­а­гональ), которую из прак­ти­чес­ких вы­­чис­­ли­­­тель­ных со­ображений (экономия ма­шин­ной па­мя­ти, уско­­­ре­ние поиска нужного элемента) удоб­но по­строч­но пе­­­ре­­пи­сать в виде одномерного мас­си­ваd [1:n×(n+1)/2]. Заметим, чтов этом случае в про­цес­се вы­­пол­нения со­от­вет­ст­вующей про­це­ду­рызна­че­­ния эле­­ментов исходной мат­ри­цыa[i, j] могут при­­­­­сва­и­ваться элементамd[m_i+j], гдеm_i = =(2n -i) (i - 1) / 2. От­­сю­да сле­ду­ет, чтоm_i+1=m_i +n -i, т.е. зна­­че­нияm_i мо­гут вы­чис­лять­ся рекурсивно, если по­ло­житьm₀ = 0. Та­­кое уп­лот­не­ние симметрической мат­­ри­цы мо­жет быть по­лез­­но, на­при­мер, в том слу­чае, ког­да она целиком не по­­­мещается в сво­бод­ную часть опе­ра­тив­ной памяти ЭВМ. Раз­ме­ще­ние ре­зуль­тирующей мат­рицыCна мес­те вво­ди­мой мат­­ри­цыBтакже по­зво­ляет оп­ти­ми­зи­ро­вать зат­­раты опе­­ративной па­мя­ти.

Изложенный алгоритм реализован в про­це­ду­ре multcram.

Формальные параметры про­цедуры.Входные:n(типinteger) -порядок квадратной мат­рицы A и количество строк пря­моугольной мат­ри­цыB;d[1:n×(n+1)/2] (типreal) - вво­ди­мая в виде од­но­мерного массиваdуплотненная мат­ри­ца A;r (типinteger) - количество столб­цов вводимой пря­мо­­у­голь­ной матрицыB;b[1:n, 1:r] (типreal) - матрицаB.Выходные:b[1:n, 1:r] (типreal) - вы­во­димая на месте вхо­дной матрицы Bматрица-про­из­ведение C = A*B.

procedure multcram (d: mas1; n,r: integer;

var b:mas2);

var s:real; i,j,k,m : integer; v array [1..n] of real;

begin

for j:=1 to r do

begin

for i:=1 to n do

begin

s:=0; m:=i-n;

for k:=1 to n do

begin

if k <= i then m := m+n+1-k

else m:=m+1;

s:=s+d[m]*b[k,j]

end {k}; v[i]:=s

end {i};

for i:=1 to n do

b[i,j]:=v[i]

end {j};

end {multcram}.

Процедура multcram была получена из про­це­ду­­рыcram[Агеев, 1976]¹путем ее перевода с язы­каALGOLна язы­кPASCAL. Тестирование про­це­ду­ры про­­водилось на IBM PC/AT-386 для тех же ис­ход­­ных данных, что и в работе Каффрей (1961):

.

При этом был получен верный результат:

.

4.8. Корректировка обратной матрицы после изменения одного элемента в прямой матрице

Пусть дана матрица A = M ^-1размерностью[n´n] и мат­ри­цаM', полученная в результате уве­­ли­чения эле­мен­таm[i,j] матрицыMна величинуd. Тогда скор­рек­ти­рованная мат­ри­цаB = (M')^-1мо­жет быть эффективно вы­­числена по эле­мен­там мат­рицыA без обращения мат­рицыMпри по­мо­щи пред­­став­лен­ной процедурыadjust.

Формальные параметры про­цедуры.Входные:a[n,n] - (типreal) - исходная матрицаA;n(типin­te­ger) - по­ря­док мат­рицыA;i, j(типin­­te­ger) - ин­дек­сы измененного эле­мен­таm[i, j] мат­ри­цы M;d(типreal) - величина из­ме­нения эле­­мен­таm[i, j].Вы­ходные:b[n, n] (типreal) - скор­рек­ти­ро­ванная мат­рица B.

procedure adjust (a:mas1;n,i,j: integer;

d:real; var b:mas1);

var t : real; r,s : integer;

begin t:=d/(a[j,i]*d+1);

for r:=1 to n do

for s:=1 to n do

b[r,s]:=a[r,s]-t*a[r,i]*a[j,s]

end {adjust};

Процедура adjustпредставляет собой перевод на языкPASCALопубликованной в работе [Агеев и др., 1976], ис­правленной, сокращенной и ор­ди­нар­­но пе­ре­ра­бо­­танной про­цедуры, опубликованной в работе [Herndon, 1961]. Тес­ти­ро­ва­ние про­цедурыadjustбыло вы­пол­не­но наIBM PC/AT-286 для следующих ис­ход­ных дан­ных:

,n = 2,i = 1,j = 2,d= 3.

Результат корректировки

с точностью до ошибок округления совпал с точ­ным решением

.

Замечание.Интересно, что исходной ал­го­ритм[Hern­don, 1961] проверялся дополнительно Р. Георгом (1962) для мат­риц порядкаn = 2, 3, ..., 10 при про­из­воль­­­ных зна­че­ни­яхi, jислучайных при­­ращенияхэле­мен­­­таm[i, j] прямой мат­рицыM: -1.0 <d<< 1.0. По ре­зультатам корректировки рас­счи­ты­ва­лось про­из­ве­де­ниеC=B*M ' ивы­чис­ля­лась по­греш­­ность каждого опы­­тапо фор­му­ле

,

при этом во всех произведенных опытах зна­че­ние eне превышало 2e-8.