begin

if 1#=рД i=7q then

begin intl: = a[i,p]; mu: = aO,ql;

a[q,i]: = ali.q]: = intl X sint+mu X cost;

a[p,i]: = a[i,p]: = int 1X cost-mu X sint end;

intl: = s[i,p]; mu:=s[i,q]; s[i,q]: = intl Xsint+muXcost; S[i,p]: = intl Xcost-muXsint end i;

mu: = sintf2; omega: = costf2; intl:=sintXcost; a[p,p]:=vlXomega + v3xmu-2Xv2Xintl; a q,q]:=vlXmu+v3Xomega+2Xv2Xintl; alp.q]-==a[q>p]:= (vl-v3) Xintl + v2x (omega-mu) end p;

comment Теперь проверка достижения конечной точности; if ind=l then begin ind:=0; go to mainl end; if thr>norm2 then go to main; fin: end jacobi;

Свидетельство к алгоритму 856

Алгоритм 856 получен из алгоритма 85а путем внесения в него поправки, предложенной В. Ф. Дейкиным в нижеследующем его «Замечании к алгоритму 85а», и модификации, предложенной П. Науром в нижеследующем его «Подтве:рждении к алгоритму 85». Кроме того, было сделано несколько тождественных преобразований для экономии машинного времени.

Свидетельство к алгоритму 85а

Алгоритм 85а получен в результате исправления, некоторого сокращения и ординарной переработки алгоритма 85 (Evans Т. G. «САСМ», 1962, Яо 4).

Более подробные сведения о методе Якоби можно получить, например, в работе Д. К. Фаддеева и В. Н. Фаддеевой [4, § 81].

Алгоритм 85а транслирован с исходными данными п=4, rho= = 10-8 и

Результаты трансляции:

0.579642502 0.459996665 0.433459111 0.514325614 2.32274880

1.00 0.42 0.54 0.66

0.42 1.00 0.32

0.44

0.54 0.32 1.00

0.22

0.66

0.44 0.22 1.00

- 0.0503284495 0.237226458

- 0.812846170 0.529595844

- 0.718845953

- 0.0956989810 0.387435463 0.569206432

0.100051435,0-11 0.479925056,0-14

0.100051435,0-11 0.796706689 0.000000000

0.479925056,0-14 0.000000000 0.242260708

0.000000000 0.179516047,0-14 0.000000000

0.380449881 0.850275473 0.0358896058 0.361941215

0.000000000

0.179516047,0-14

0.000000000

0.638283803

С-оглашо работе Фаддеева Д. ¥ч. -и Фаддеевой Ь. R. \4, с. Ъ1Ъ- исходная матрица должна иметь собственные значения 2.32274884, 0.63828379, 0.79670672, 0.24226073 и третий собственный вектор (-0.71884900, -0.09569928, 0.38743715, 0.5692089).

Замечание к алгоритму 85а

Г. Г. Дядюша, Киев, октябрь 1966

Алгоритм 85а может быть существенно усовершенствован, если учесть следующие обязательства.

1. Для вычисления тригонометрических функций угла элементарного вращения ф достаточно двух операций извлечения квадратного корня.

2. Вместо sin ф можно применять tg ф, который удобно вычисляется по формуле

tg (p-aij/iix+signin+HP)) щ],

0. iO.

где V={и - 0.11)12; i*. = Yv- + оц; 8(1*) =

1, 1л = 0.

3 Диагональные элементы преобразованной матрицы удобно вычисляются по формулам aii, а]з = ац +n±ix,i.

4. Элемент aij должен быть равен нулю.

5. Целесообразно при каждой итерации выбирать преграду чуть ниже среднего квадратичного недиагональных элементов.

6. Достаточно хранить в памяти машины только верхнюю половину симметричной матрицы («ак в алгоритме 67а).

7. Для многих программирующих программ удобно применять только одномерные массивы.

В приводимой ниже процедуре jacobi2 исходная матрица занисы-вается в массиве а[1:пХ{п+1)/2] в следующем порядке: в первых п ячейках - диагональные элементы, в последующих - элементы над-диагональной верхней треугольной матрицы построчно. Матрица s записывается в массиве s[l:nXn] по столбцам, procedure jacobi2(n,rho)dataresult: (a)result: (s);

value n.pho; real rho;integer n; array a,s; begin real norml,norm2,mu,mul,t,cost,d;

Integer i,j,k,p,q,r,il,kl,rl,nl,n2,n3; Boolean g;

p:=0;

for i:= 1 step 1 until n do for j: = 1 step 1 until n do begin p:=:p+l;

s[p]:= if i=j then 1 else 0 end i;

nl: = n-1; n2:=nl-1; n3:-=nx (n-f l)/2; g:= false; thr: norml:=0;

for i: = n+l step 1 until n3 do norml:==norml-l-a[i]f2;

norm 1: = sqrt (2 X norm 1) /n;

if g then g:- false else norm2:=norml Xrho;

if normKnorm2 then norml:=norm2;

k: = n;

for q:=:0 step 1 until n2 do begin d: = a[q-i-l];

for p:=q +1 step 1 until nl do begin k: = k+l; if abs(a[k])>norml then

begin t: = alk]; a[k]:=0; g:= true; mu:= (a[p+l]-d)/2; mul:=sqrt(tf2+muf2); if mu<0 then mul:=-mul; mu:=mu+mul; t:=t/mu; a[p+l]: = d: = d+mu; mul:=2Xmul; d: = d-mul; cost:=sqrt(mu/mul); i:=q+n; r:=p+n; il:=nxq; rl:=nXp; for j:=0 step 1 until nl do begin

if j=7q AJ=5P then begin mu: = costXa[i]; mul:=costXa[r];

a[r]:=muXt+mul; a[i]:=:mu-tXmul end ifj;

il:=il + l; rl:=rl + l; kl: = n2-j; mu:=costXs[il]; mul: = costXs[rl]; s[r 1]: =:mu Xt+mu 1; s[i 1]:=mu-tXmu 1; i:= if j<q then i+kl else i+l; r:= if ]<p then r+kl else r+1 end j; end ifabs end p; a[q+l]: = d end q; if g then go to thr end jacobi2; *

Подтверждение к алгоритму 85a

Г. Г. Дядюша, Киев, октябрь 1967

Алгоритм 85а в разных вариантах неоднократно транслировался с помощью транслятора ТА-2 в виде составной части программ, предназначенных для кнантово-химических расчетов, и всегда работал безотказно.

Кроме ординарной переработки на входной язык транслятора ТА-2 эти варианты отличались от процедуры jacobi2, дриведенной в преды дущем «Замечании», следующими деталями.

1. Конструкции «...:=...:=...» и «if... Д... then) расписывались подробнее из-за того, что транслятор ТА-2 переводит их недостаточно экономно.

2. погт2 не вычислялась, а задавалась вместо rho.

3. norml вычислялась по другим формулам, например,

шгт1 = [1-{-ЯX (« - 1)2]- 2 «г/

»

ТТроцедура jacobi2, приведениая выше, отличается от представленной Г. Дядю-£ией тем, что а ней исправлены некоторые очевидные ошибки. После исправления процедура jacobi2 была транслирована в системе TA-LM и для примера, яриведенного в «Свидетельстве к алгоритму 85а», дала точно такие же результаты, как и процедура jacobi. {Прим. ред.)