> В подпрограммах NBTRIP и PBTRIP используются при работе другие подпрограммы, текст которых также приведен в приложении. Укажем на их основные функции, приведенные в комментариях к программе.

Подпрограмма LUDECO представляет матрицу А в виде произведения верхней (U) и нижней (L) треугольных матриц. Результат записывается на место матрицы А. Главный элемент при этом не выделяется.

Подпрограмма N[lJLPVT(AyByC,ORDER) производит умножение вектора В на матрицу А и вычитает полученное произведение из другого вектора С. Результат расчета запоминается в С, поэтому массив С при работе программы изменяется.

Подпрограмма LVSOLV{А,В,С,ORDER) решает систему алгебраических уравнений /*С = Б и записывает решение в С. Матрица А должна быть заранее представлена в виде произведения верхней (и) и нижней (L) треугольных матриц.- /7риж. перев.

В подпрограммах NBTRIP и PBTRIP выделение главного элемента при решении систем уравнений методом исключения не производится. Отметим, что для каждого значения N можно создать специальную подпрограмму решения системы уравнений с блочной трехдиагональной матрицей, которая будет работать быстрее, чем общая подпрограмма, приведенная ниже>.

Текст подпрограммы NBTRIP

с...subroutine то solve non-periodic block tridiagonal с...system of equations witkout pivoting strategy c...with the dimensions of the block matrices being c...n x n (n is any number greater than 1). c. . .

subroutine nbtrip(a,в,с,d,il,iu,order)

c...

dimension a(1),b(1),c(1),d(1) integer order,ordsq

c. . .

c...a = sub diagonal matrix c...b = diagonal matrix c...c = sup diagonal matrix c...d = right hand side vector

C...IL = LOWER VALUE OK INDEX FOR WHICH MATRICES ARE DEFINED C...IU .= UPPER VALUE OF INDEX FOR WHICH MATRICES ARE DEFINED C... (SOLUTION IS SOUGHT FOR BTRI(A,В,C)*X = D

C... FOR INDICES OF X BETWEEN IL AND IU (INCLUSIVE).

C... SOLUTION WRITTEN IN D VECTOR (ORIGINAL CONTENTS

C... ARE OVERWRITTEN)).

C...ORDER = ORDER OF A,B,C MATRICES AND LENGTH OF D VECTOR

C... AT EACH POINT DENOTED BY INDEX I

C... (ORDER CAN BE ANY INTEGER GREATER THAN 1).

C. . .

C...THE MATRICES AND VECTORS ARE STORED IN SINGLE SUBSCRIPT FORM C. . .

ORDSQ = ORDER--a

C. . .

C...FORWARD ELIMINATION C...

I = IL

lOMAT = 1+(I-1)*0RDSQ lOVEC = 1+(I-1)*0RDER CALL LUDECO(B(I0MAT),ORDER)

CALL LUSOLV(B(I0MAT),D(IOVEC),D(IOVEC),ORDER)

DO 100 J=l,ORDER

lOMATJ = IOMAT+(J-1)*ORDER

CALL LUSOLV(B(IOMAT),C(IOMATJ),C(IOMATJ),ORDER) 100 CONTINUE 200 CONTINUE

I = I+l

lOMAT = 1+(I-1)*0RDSQ lOVEC = 1+(I-1)*0RDER IIMAT = lOMAT-ORDSQ IIVEC = lOVEC-ORDER

CALL MULPUT(A(IOMAT),D(I1VEC),D(IOVEC),ORDER)

DO 300 J=l,ORDER

lOMATJ = IOMAT+(J-1)*ORDER

IIMATJ = I1MAT+(J-1)*0RDER

CALL MULPUT(A(IOMAT),C(I1MATJ),B(IOMATJ),ORDER) 300 CONTINUE

CALL LUDECO(B(I0MAT),ORDER)

CALL LUSOLV(B(I0MAT),D(IOVEC),D(I0VEC),ORDER)

IF(I.EQ.IU) GO TO 500

DO 400 J=l,ORDER

lOMATJ = IOMAT+(J-1)*ORDER

CALL LUSOLV(B(IOMAT),С(IOMATJ),С(IOMATJ),ORDER) 400 CONTINUE

GO TO 200 500 CONTINUE C...

C...BACK SUBSTITUTION C...

I = IU 600 CONTINUE I = I-l

lOMAT = 1+(I-1)*0RDSQ lOVEC = 1+(I-1)*0RDER IIVEC" = lOVEC+ORDER

CALL MULPUT(C(IOMAT),D(I1VPC),D(IOVEC),ORDER) IF (I.GT.IL) GO TO 600

Текст подпрограммы PBTRIP

subroutine то solve periodic block tridiagonal system of equations without pivoting strategy. each block matrix may be of dimension n with n any number greater than 1.

subroutine pbtrip(a,b,c,d,il,iu,order) dimension a(1),b(1),c(1),d(1) dimension ad(25),cd(25) integer order,ordsq

a = sub diagonal matrix в = diagonal matrix с = sup diagonal matrix d = right hand side vector

il = lower value of index for which matrices are defined

iu = upp::r value of index for which matrices are defined

(solution is sought for btri(a,b,c)*x = d for indices of x betweeen il and iu (inclusive). solution written in d vector (original contents are overwritten)). order = order of a,b,C matrices and length of d vector at each point denoted by index I (order can be any integer greater than 1) (arrays ad and cd must be at least of length 0rder*2) (current length of 25 anticipates maximum order of 5).

is = il+1 ie = iu-1 ordsq = 0rder**2 lUMAT « 1+(iu-1)*drdsq lUVEC = 1+(iu-1)*0rder lEMAT = 1+(ie-1)*0rdsq lEVEC = 1+(ie-1)*0rder

..forward elimination

i = il

lOMAT = 1+(i-1)*0rdsq lOVEC » 1+(i-1)*0rder call ludeco(b(iomat),order)

call lusolv(b(i0mat),d(iovec),d(iovec).order)

do 10 J=l,order

lOMATJ - i0mat+(J-1)*0rder

call lus0lv(b(i0mat),c(iomatj),c(iomatj).order) call lusolv(b(i0mat),a(i0matj),a(i0matj).order) 10 continue

DO 200 i e IS,ie lOMAT « 1+(i-1)*0rdsq lOVEC = 1+(i-1)*0rder iimat = lOMAT-ORDSQ iivec = lOVEC-order DO 20 J=l.ORDSQ