ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 06.12.2023
Просмотров: 344
Скачиваний: 6
ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.
Работа sum c dim=2, но без mask:
sum(b,dim=2)=
9
sum(b,dim=2)=
12
sum(b,dim=2)=
177
sum(b,dim=2)=
30
sum(b,dim=2)=
61
sum(b,dim=2)=
159
sum(b,dim=2)=
75
sum(b,dim=2)=
66
Размер левого слоя массива b есть size(b(:,1,:)=
8
Форма левого слоя массива b есть shape(b(:,1,:)=
2 4
Содержимое левого слоя массива b:
1 73 31 19 2
8 43 20
Содержимое cреднего слоя массива b:
3 93 13 33 4
10 63 22
Содержимое правого слоя массива b:
5 11 17 23 6
12 53 24
Работа sum c dim=3, но без mask:
sum(b,dim=3)=
124
sum(b,dim=3)=
73
sum(b,dim=3)=
142
sum(b,dim=3)=
99
sum(b,dim=3)=
56
sum(b,dim=3)=
95
Размер листа массива b есть size(b(:,:,1)=
6
Форма листа массива b есть shape(b)=
2 3
Работа sum c dim=1 и с mask:
sum(b,1,mask)=
0
sum(b,1,mask)=
3
sum(b,1,mask)=
0
sum(b,1,mask)=
73
sum(b,1,mask)=
93
sum(b,1,mask)=
0
sum(b,1,mask)=
43
sum(b,1,mask)=
76
sum(b,1,mask)=
53
sum(b,1,mask)=
0
sum(b,1,mask)=
33
sum(b,1,mask)=
23 255
Работа sum c dim=2 и с mask:
sum(b,2,mask)=
3
sum(b,2,mask)=
0
sum(b,2,mask)=
166
sum(b,2,mask)=
0
sum(b,2,mask)=
13
sum(b,2,mask)=
159
sum(b,2,mask)=
56
sum(b,2,mask)=
0
Размер левого слоя массива b есть size(b(:,1,:)=
8
Форма левого слоя массива b есть shape(b(:,1,:)=
2 4
Содержимое левого слоя массива b:
1 73 31 19 2
8 43 20
Содержимое cреднего слоя массива b:
3 93 13 33 4
10 63 22
Содержимое правого слоя массива b:
5 11 17 23 6
12 53 24
Работа sum c dim=3 и с mask:
sum(b,3,mask)=
73
sum(b,3,mask)=
43
sum(b,3,mask)=
142
sum(b,3,mask)=
63
sum(b,3,mask)=
23
sum(b,3,mask)=
53 256
7.3
Функции умножения векторов и матриц
1. dot_product(vector_a,vector_b) — скалярное произведение векторов:
program tdot; implicit none; integer :: ia(3)=(/1,2,3/), ib(3)=(/4,5,6/)
real
::
a(3)=(/1.0,2.0,3.0/), b(3)=(/4.0,5.0,6.0/)
complex :: ca(3)=(/(1.0,0.0),(2.0,0.0),(3.0,0.0)/)
complex :: cb(3)=(/(4.0,0.0),(5.0,0.0),(6.0,0.0)/)
logical :: la(3)=(/.false.,.true.,.false./)
logical :: lb(3)=(/.false.,.true.,.false./)
logical :: lc(3)=(/.true.,.false.,.true./)
write(*,’(" integer: ia=",3g4.0)’) ia write(*,’("
ib=",3g4.0)’) ib write(*,’("
dot_product(ia,ib)=",i4)’) dot_product(ia,ib)
write(*,’("
real: a=",3g15.7)’) a write(*,’("
b=",3g15.7)’) b write(*,’("
dot_product( a, b)=",e15.7)’) dot_product( a, b)
write(*,’(" complex: ca=",6g7.2)’) ca write(*,’("
cb=",6g7.2)’) cb write(*,’("
dot_product(ca,cb)=",2g15.7)’) dot_product(ca,cb)
write(*,’(" logical: la=",3g7.0,10x,"
la=",3g7.0)’) la,la write(*,’("
lb=",3g7.0,10x,"
lc=",3g7.0)’) lb,lc write(*,’("
dot_product(la,lb)=",g5.0,10x,&
&"
dot_product(la,lc)=",g5.0)’) dot_product(la,lb), dot_product(la,lc)
write(*,’(" complex: ca=",6g7.2)’) ca write(*,’(" integer: ib=",3g7.0)’) ib write(*,’("
dot_product(ca,ib)=",2g15.7)’) dot_product(ca,ib)
end integer: ia=
1 2
3
ib=
4 5
6
dot_product(ia,ib)=
32
real: a=
1.000000 2.000000 3.000000
b=
4.000000 5.000000 6.000000
dot_product( a, b)=
0.3200000E+02
complex: ca=1.0 0.0 2.0 0.0 3.0 0.0
cb=4.0 0.0 5.0 0.0 6.0 0.0
dot_product(ca,cb)=
32.00000 0.000000
logical: la=
F
T
F
la=
F
T
F
lb=
F
T
F
lc=
T
F
T
dot_product(la,lb)=
T
dot_product(la,lc)=
F
complex: ca=1.0 0.0 2.0 0.0 3.0 0.0
integer: ib=
4 5
6
dot_product(ca,ib)=
32.00000 0.000000
Для векторов логического типа dot_product получает результат дизъюнкции поэлементных конъюнкций.
257
2. matmul(matrix_a,matrix_b) — умножение или двух матриц,
или матрицы на вектора, или вектора на матрицу:
rab(i, j) =
k
X
l=1
a(i, l) ∗ b(l, j); i = 1(1)n; j = 1(1)m rac(i) =
m
X
j=1
a(i, j) ∗ c(j);
i = 1(1)n rca(j) =
k
X
i=1
c(i) ∗ a(i, j);
j = 1(1)m program tmatmul; implicit none integer :: a(2,3)=reshape((/1,2,3,4,5,6/), shape=(/2,3/)), i, j integer :: b(3,2)=reshape((/1,2,3,4,5,6/), shape=(/3,2/)), rab(2,2)
integer :: c(2)=(/1,2/), d(3)=(/1,2,3/)
write(*,’(" a: ",3g3.0/4x,3g3.0)’) ((a(i,j),j=1,3),i=1,2)
write(*,’(" b: ",2(2g3.0/4x)2g3.0))’)
((b(i,j),j=1,2),i=1,3)
rab=matmul(a,b)
write(*,’(" matmul(a,b)=",(2g3.0/13x,2g3.0))’) matmul(a,b)
write(*,’("
rab=",2(2g3.0/13x))’) ((rab(i,j),j=1,2),i=1,2)
write(*,’(" a: ",3g3.0/"
",3g3.0)’) ((a(i,j),j=1,3),i=1,2)
write(*,’(" d: ",3g3.0)’) d write(*,’(" matmul(a,d)=",2g3.0)’) matmul(a,d)
write(*,’(/" c: ",2g3.0)’) c write(*,’(" a: ",3g3.0/"
",3g3.0)’) ((a(i,j),j=1,3),i=1,2)
write(*,’(" matmul(c,a)=",3g3.0))’) matmul(c,a)
end a:
1 3
5
|
a:
1 3
5
|
c:
1 2
2 4
6
|
2 4
6
|
|
|
b:
1 4
|
d:
1 2
3
|
a:
1 3
5 2
5
|
|
2 4
6 3
6
|
|
matmul(a,b)= 22 28
|
matmul(a,d)= 22 28
|
1 ... 10 11 12 13 14 15 16 17 ... 23
matmul(c,a)=
5 11 17 49 64
|
|
rab= 22 49
|
|
28 64
|
|
Помним, что умножение матриц определяется лишь для согласован- ных матриц (т.е. число столбцов первой должно равняться числу строк второй). Элемент матрицы-произведения, стоящий на пересе- чении i-ой строки и j-ого столбца, равен сумме произведений эле- ментов i-ой строки первой матрицы на соответствующие элементы j-ого столбца второй.
258
matmul ЧУВСТВУЕТ отсутствие согласованности матриц:
program tmat1; implicit none integer :: a(2,3)=reshape((/1,2,3,4,5,6/), shape=(/2,3/))
integer :: b(2,3)=reshape((/1,2,3,4,5,6/), shape=(/2,3/))
integer, allocatable :: f(:,:), g(:,:), h(:,:), x(:,:)
integer ier, i, j, l, m, n; character(1) sn; character(20) sf
! write(*,*) matmul(a,b)
write(*,’("Работа matmul с динамическими массивами")’)
l=3; m=4; n=5
allocate(f(l,m),stat=ier); if (ier.ne.0) write(*,*) ’Не размещается f’
allocate(g(m,n),stat=ier); if (ier.ne.0) write(*,*) ’Не размещается g’
allocate(h(l,n),stat=ier); if (ier.ne.0) write(*,*) ’Не размещается h’
allocate(x(n,m),stat=ier); if (ier.ne.0) write(*,*) ’Не размещается h’
do i=1,l; f(i,:)=(/(
i+j,j=1,m)/); enddo do i=1,m; g(i,:)=(/(10*i+j,j=1,n)/); enddo;
h=matmul(f,g)
write(*,’(a,i1,a,4i6)’) (’f(’,i,’)=’,f(i,:),i=1,l); write(*,*)
write(*,’(a,i1,a,5i6)’) (’g(’,i,’)=’,g(i,:),i=1,m); write(*,*)
write(*,’(a,i1,a,5i6)’) (’h(’,i,’)=’,h(i,:),i=1,l); write(*,*)
write(sn,’(i1)’) n; sf=’(’//sn//’(i6)’//’)’; write(*,sf)
(h(i,:),i=1,l)
h=matmul(g,f)
end
Так, раскомментировав в tmat1 строку !write(*,*) matmul(a,b),
в случае несогласованных статических массивов a и b получим ещё на шаге компиляции сообщение:
write(*,*) matmul(a,b)
1
ошибка: Different shape on dimension 2 for argument ’matrix_a’
and dimension 1 for argument ’matrix_b’ at (1) for intrinsic matmul
А в случае несогласованных размещаемых массивов при вызо- ве matmul(g,f ) (см. предпоследний оператор программы) получим соответствующее сообщение и на шаге выполнения:
Fortran runtime error: dimension of array B incorrect in MATMUL intrinsic
Полезно ещё раз обратить внимание на возможность использования:
1) сечений для построкового вывода матриц;
2) строковой переменной (в tmat1 переменная sf) для орграниации динамического повторителя (в данном случае для вывода нужного количества столбцов матрицы) при формировке формата вывода.
259
program tmat1; implicit none integer :: a(2,3)=reshape((/1,2,3,4,5,6/), shape=(/2,3/))
integer :: b(2,3)=reshape((/1,2,3,4,5,6/), shape=(/2,3/))
integer, allocatable :: f(:,:), g(:,:), h(:,:), x(:,:)
integer ier, i, j, l, m, n; character(1) sn; character(20) sf
! write(*,*) matmul(a,b)
write(*,’("Работа matmul с динамическими массивами")’)
l=3; m=4; n=5
allocate(f(l,m),stat=ier); if (ier.ne.0) write(*,*) ’Не размещается f’
allocate(g(m,n),stat=ier); if (ier.ne.0) write(*,*) ’Не размещается g’
allocate(h(l,n),stat=ier); if (ier.ne.0) write(*,*) ’Не размещается h’
allocate(x(n,m),stat=ier); if (ier.ne.0) write(*,*) ’Не размещается h’
do i=1,l; f(i,:)=(/(
i+j,j=1,m)/); enddo do i=1,m; g(i,:)=(/(10*i+j,j=1,n)/); enddo;
h=matmul(f,g)
write(*,’(a,i1,a,4i6)’) (’f(’,i,’)=’,f(i,:),i=1,l); write(*,*)
write(*,’(a,i1,a,5i6)’) (’g(’,i,’)=’,g(i,:),i=1,m); write(*,*)
write(*,’(a,i1,a,5i6)’) (’h(’,i,’)=’,h(i,:),i=1,l); write(*,*)
write(sn,’(i1)’) n; sf=’(’//sn//’(i6)’//’)’; write(*,sf)
(h(i,:),i=1,l)
h=matmul(g,f)
end
Так, раскомментировав в tmat1 строку !write(*,*) matmul(a,b),
в случае несогласованных статических массивов a и b получим ещё на шаге компиляции сообщение:
write(*,*) matmul(a,b)
1
ошибка: Different shape on dimension 2 for argument ’matrix_a’
and dimension 1 for argument ’matrix_b’ at (1) for intrinsic matmul
А в случае несогласованных размещаемых массивов при вызо- ве matmul(g,f ) (см. предпоследний оператор программы) получим соответствующее сообщение и на шаге выполнения:
Fortran runtime error: dimension of array B incorrect in MATMUL intrinsic
Полезно ещё раз обратить внимание на возможность использования:
1) сечений для построкового вывода матриц;
2) строковой переменной (в tmat1 переменная sf) для орграниации динамического повторителя (в данном случае для вывода нужного количества столбцов матрицы) при формировке формата вывода.
259
matmul и dot_product осуществляют лишь формальное приведение числовых типов друг к другу, т.е. данное типа real(4), содержащее лишь
7-8 десятичных цифр, после преобразования в тип real(8) будет пере- писано по формату real(8), но из шестнадцати соответствующих цифр результата, верными будут лишь 7-8 старших.
program tmat2; implicit none; real(4), parameter :: c13=1.0/3.0
real(8) :: a(2,2), c(2,2); real(4) :: b(2,2); integer i, j a=1d0; b=c13; c=matmul(a,b)
write(*,’(14x,"a",18x,"b",15x,"matmul(a,b)")’);
write(*,’(d25.16,e15.7,d25.16)’) ((a(i,j),b(i,j),c(i,j),j=1,2),i=1,2)
end a
b matmul(a,b)
0.1000000000000000D+01 0.3333333E+00 0.6666666865348816D+00 0.1000000000000000D+01 0.3333333E+00 0.6666666865348816D+00 0.1000000000000000D+01 0.3333333E+00 0.6666666865348816D+00 0.1000000000000000D+01 0.3333333E+00 0.6666666865348816D+00
В [2] написано, что matmul при работе с матрицами и векторами типа logical выполняет операцию логического умножения. Точнее, ре- зультатом оказывается матрица элементы которой находятся согласно правилу перемножения числовых матриц, в котором сложение заменено дизъюнкцией, а умножение — конъюнкцией.
program tmat3; implicit none logical :: a(2,2)=reshape((/.true.,.true.,.false.,.false./),(/2,2/))
logical :: b(2,2)=reshape((/.true.,.false.,.true.,.false./),(/2,2/))
logical c(2,2), d(2,2),e(2,2); integer i, j write(*,*) ’ a=’,a; write(*,*) ’ b=’,b; d=matmul(a,b); e=a.and.b write(*,’(1x,"i",1x,"j",2x"a",2x,"b",2x,"c",2x,"a.and.b",3x,&
&"matmul(a,b)")’)
do j=1,2; do i=1,2;
c(i,j)= a(i,j).and.b(i,j)
write(*,’(2g2.0,3g3.0,3x,g3.0,g10.0)’)&
&
i,j,a(i,j), b(i,j), c(i,j),e(i,j),d(i,j)
enddo; enddo end a= T T F F
b= T F T F
i j a
b c
a.and.b matmul(a,b)
1 1
T
T
T
T
T
2 1
T
F
F
F
T
1 2
F
T
F
F
T
2 2
F
F
F
F
T
T T T T
260
7-8 десятичных цифр, после преобразования в тип real(8) будет пере- писано по формату real(8), но из шестнадцати соответствующих цифр результата, верными будут лишь 7-8 старших.
program tmat2; implicit none; real(4), parameter :: c13=1.0/3.0
real(8) :: a(2,2), c(2,2); real(4) :: b(2,2); integer i, j a=1d0; b=c13; c=matmul(a,b)
write(*,’(14x,"a",18x,"b",15x,"matmul(a,b)")’);
write(*,’(d25.16,e15.7,d25.16)’) ((a(i,j),b(i,j),c(i,j),j=1,2),i=1,2)
end a
b matmul(a,b)
0.1000000000000000D+01 0.3333333E+00 0.6666666865348816D+00 0.1000000000000000D+01 0.3333333E+00 0.6666666865348816D+00 0.1000000000000000D+01 0.3333333E+00 0.6666666865348816D+00 0.1000000000000000D+01 0.3333333E+00 0.6666666865348816D+00
В [2] написано, что matmul при работе с матрицами и векторами типа logical выполняет операцию логического умножения. Точнее, ре- зультатом оказывается матрица элементы которой находятся согласно правилу перемножения числовых матриц, в котором сложение заменено дизъюнкцией, а умножение — конъюнкцией.
program tmat3; implicit none logical :: a(2,2)=reshape((/.true.,.true.,.false.,.false./),(/2,2/))
logical :: b(2,2)=reshape((/.true.,.false.,.true.,.false./),(/2,2/))
logical c(2,2), d(2,2),e(2,2); integer i, j write(*,*) ’ a=’,a; write(*,*) ’ b=’,b; d=matmul(a,b); e=a.and.b write(*,’(1x,"i",1x,"j",2x"a",2x,"b",2x,"c",2x,"a.and.b",3x,&
&"matmul(a,b)")’)
do j=1,2; do i=1,2;
c(i,j)= a(i,j).and.b(i,j)
write(*,’(2g2.0,3g3.0,3x,g3.0,g10.0)’)&
&
i,j,a(i,j), b(i,j), c(i,j),e(i,j),d(i,j)
enddo; enddo end a= T T F F
b= T F T F
i j a
b c
a.and.b matmul(a,b)
1 1
T
T
T
T
T
2 1
T
F
F
F
T
1 2
F
T
F
F
T
2 2
F
F
F
F
T
T T T T
260
О старом ФОРТРАНе последнее слово.
В старых версиях ФОРТРАНа допускались только статические с це- лочисленными константами в качестве граничных индексов.
Однако во многих задачах размеры матриц естественнее было вводить или вычислять. Возникало противоречие между возможностями описа- ния матриц и требованиями эксплуатации программы. Для его разреше- ния при описании матрицы в главной программе указывали максимально возможные (в разумных пределах) количества строк и столбцов, которые могли потребоваться в задаче, размещая матрицу с нужными размерами в левом верхнем углу описанной. Например:
program tmat4; implicit none; integer, parameter :: ndim=5
integer a(ndim,ndim), b(ndim,ndim), c(ndim,ndim), d(ndim,ndim)
integer l, m, n, i, j l=2; m=3; n=4
! Для краткости моделируем расчёт l,m и n a=-1; b=-2; c=-3; d=-4
! Для уяснения результата засеем массивы do i=1,l; do j=1,m; a(i,j)=i+j; enddo; enddo do i=1,m; do j=1,n; b(i,j)=i*j; enddo; enddo do i=1,ndim; write(*,’(" a(",i2,",:)=",50i3)’) i,(a(i,j),j=1,ndim); enddo do i=1,ndim; write(*,’(" b(",i2,",:)=",50i3)’) i,(b(i,j),j=1,ndim); enddo call matmul1(a,b,c,ndim,l,m,n)
do i=1,ndim; write(*,’(" c(",i2,",:)=",50i5)’) i,(c(i,j),j=1,ndim); enddo call matmul2(a,b,d,l,m,n)
do i=1,ndim; write(*,’(" d(",i2,",:)=",50i5)’) i,(d(i,j),j=1,ndim); enddo end
В tmat4 принято, что могут потребоваться матрицы размером не боль- ше чем из пяти строк и пяти столбцов. Процедура умножения матриц
(назовём её matmul1) могла иметь вид:
subroutine matmul1(a,b,c,ndim,l,m,n); implicit none integer ndim, l, m, n, i, j, k, s; integer a(ndim,m),b(ndim,n),c(ndim,n)
do i=1,l; do k=1,n; s=0; do j=1,m; s=s+a(i,j)*b(j,k); enddo;
c(i,k)=s enddo; enddo end
В ней число строк матриц, служащих формальными аргументами рав- но заявленному ndim, что принципиально важно. Часто встречалась ситуация, когда программист при описании формального аргумен- та (матрицы) указывал в качестве её строковой размерности значение переменной, хранящей не заявленное в главной программе, а нужное по проблемному разумению число строк, т.е. так
261
subroutine matmul2(a,b,c,l,m,n); implicit none integer ndim, l, m, n, i, j, k, s;
integer a(l,m),b(m,n), c(l,n)
do i=1,l; do k=1,n; s=0; do j=1,m; s=s+a(i,j)*b(j,k); enddo; c(i,k)=s enddo; enddo end забывая, что число строк матрицы, описанной в главной программе, рав- но ПЯТИ. В итоге программист обманывал и подпрограмму, и себя. Уяс- ним результат работы test4 при l=2, m=3 и n=4:
a( 1,:)=
2 3
4 -1 -1 При работе matmul1 подпрограмма "считает",
a( 2,:)=
3 4
5 -1 -1 что a(1,2), b(1,2), c(1,2) следуют a( 3,:)= -1 -1 -1 -1 -1 за a(5,1), b(5,1) и с(5,1) соответственно,т.е a( 4,:)= -1 -1 -1 -1 -1 так как определено главной программой.
a( 5,:)= -1 -1 -1 -1 -1 Поэтому результат расчёта, например, c(1,1)
по формуле:
b( 1,:)=
1 2
3 4 -2
a(1,1)*b(1,1)+a(1,2)*b(2,1)+a(1,3)*b(3,1)=
b( 2,:)=
2 4
6 8 -2 2*1 + 3*2 + 4*3 = 2+6+12= 20
ВЕРЕН.
b( 3,:)=
3 6
9 12 -2
b( 4,:)= -2 -2 -2 -2 -2
MATMUL2 по описаниям a(l,m), b(m,n), c(l,n),
b( 5,:)= -2 -2 -2 -2 -2
что при вызове даёт a(2,3), b(3,4), c(2,4)
"полагает", что a(1,2), b(1,2), c(1,2)
c( 1,:)=
20 40 60 80
-3
следуют за a(2,1), b(3,1), c(2,1),
c( 2,:)=
26 52 78 104
-3
что противоречит определению главной c( 3,:)=
-3
-3
-3
-3
-3
программы. Поэтому подпрограмма,
c( 4,:)=
-3
-3
-3
-3
-3
беря элемент a(1,2), c точки зрения c( 5,:)=
-3
-3
-3
-3
-3
главной программы берёт элемент a(3,1)=-1
Cоответственно элемент, именуемый d( 1,:)=
-3 0
-4
-4
-4
подпрограммой как a(1,3), оказывается d( 2,:)=
10 -13
-4
-4
-4
для главной a(4,1). Так что для расчёта d( 3,:)=
-4 9
-4
-4
-4
c(1,1) берутся вовсе не те элементы:
d( 4,:)=
2
-4
-4
-4
-4
d( 5,:)=
4
-4
-4
-4
-4
MATMUL2 "думает" об элементах a(1,1)*b(1,1)+a(1,2)*b(2,1)+a(1,3)*b(3,1)=
а на самом деле берёт a(1,1)*b(1,1)+a(3,1)*b(2,1)+a(4,1)*b(3,1)=
= 2*1 + (-1*2)+ (-1*3)=2-2-3=2-5=-3
В ФОРТРАНе аргументы процедур по умолчанию передаются по ад- ресу. Поэтому в качесве адреса начального элемента матрицы a(1:2,1:3)
будет передан адрес начального элемента матрицы a(1:5,1:5), описанной в главной программе. Так как ФОРТРАН располагает элементы матри- цы в оперативной памяти непосредственно друг за другом по столбцам,
то третий и четвёртый элементы первого столбца матрицы a(1:5,1:5) с точки зрения MATMUL2 ,полагающей согласно описанию формальных аргументов, что она имеет дело с матрицей a(1:2,1:3), окажутся двумя последними элементами её первой строки.
262
integer a(l,m),b(m,n), c(l,n)
do i=1,l; do k=1,n; s=0; do j=1,m; s=s+a(i,j)*b(j,k); enddo; c(i,k)=s enddo; enddo end забывая, что число строк матрицы, описанной в главной программе, рав- но ПЯТИ. В итоге программист обманывал и подпрограмму, и себя. Уяс- ним результат работы test4 при l=2, m=3 и n=4:
a( 1,:)=
2 3
4 -1 -1 При работе matmul1 подпрограмма "считает",
a( 2,:)=
3 4
5 -1 -1 что a(1,2), b(1,2), c(1,2) следуют a( 3,:)= -1 -1 -1 -1 -1 за a(5,1), b(5,1) и с(5,1) соответственно,т.е a( 4,:)= -1 -1 -1 -1 -1 так как определено главной программой.
a( 5,:)= -1 -1 -1 -1 -1 Поэтому результат расчёта, например, c(1,1)
по формуле:
b( 1,:)=
1 2
3 4 -2
a(1,1)*b(1,1)+a(1,2)*b(2,1)+a(1,3)*b(3,1)=
b( 2,:)=
2 4
6 8 -2 2*1 + 3*2 + 4*3 = 2+6+12= 20
ВЕРЕН.
b( 3,:)=
3 6
9 12 -2
b( 4,:)= -2 -2 -2 -2 -2
MATMUL2 по описаниям a(l,m), b(m,n), c(l,n),
b( 5,:)= -2 -2 -2 -2 -2
что при вызове даёт a(2,3), b(3,4), c(2,4)
"полагает", что a(1,2), b(1,2), c(1,2)
c( 1,:)=
20 40 60 80
-3
следуют за a(2,1), b(3,1), c(2,1),
c( 2,:)=
26 52 78 104
-3
что противоречит определению главной c( 3,:)=
-3
-3
-3
-3
-3
программы. Поэтому подпрограмма,
c( 4,:)=
-3
-3
-3
-3
-3
беря элемент a(1,2), c точки зрения c( 5,:)=
-3
-3
-3
-3
-3
главной программы берёт элемент a(3,1)=-1
Cоответственно элемент, именуемый d( 1,:)=
-3 0
-4
-4
-4
подпрограммой как a(1,3), оказывается d( 2,:)=
10 -13
-4
-4
-4
для главной a(4,1). Так что для расчёта d( 3,:)=
-4 9
-4
-4
-4
c(1,1) берутся вовсе не те элементы:
d( 4,:)=
2
-4
-4
-4
-4
d( 5,:)=
4
-4
-4
-4
-4
MATMUL2 "думает" об элементах a(1,1)*b(1,1)+a(1,2)*b(2,1)+a(1,3)*b(3,1)=
а на самом деле берёт a(1,1)*b(1,1)+a(3,1)*b(2,1)+a(4,1)*b(3,1)=
= 2*1 + (-1*2)+ (-1*3)=2-2-3=2-5=-3
В ФОРТРАНе аргументы процедур по умолчанию передаются по ад- ресу. Поэтому в качесве адреса начального элемента матрицы a(1:2,1:3)
будет передан адрес начального элемента матрицы a(1:5,1:5), описанной в главной программе. Так как ФОРТРАН располагает элементы матри- цы в оперативной памяти непосредственно друг за другом по столбцам,
то третий и четвёртый элементы первого столбца матрицы a(1:5,1:5) с точки зрения MATMUL2 ,полагающей согласно описанию формальных аргументов, что она имеет дело с матрицей a(1:2,1:3), окажутся двумя последними элементами её первой строки.
262
Выводы:
1. Размещаемые массивы современного ФОРТРАНа позволяют просто решать те задачи, для которых на старом ФОРТРАНе приходилось использовать некие искусственные приёмы.
2. При незначительной модификации старых программ (для сохране- ния единого стиля исходного текста) приходится в качестве допол- нительного аргумента процедуры передавать наряду с матрицей и её заявленный строковый размер.
3. Современный ФОРТРАН позволяет обойтись без упомянутой пере- дачи строкового размера за счёт указания о заимствовании фор- мальным аргументом, являющимся матрицей, формы фактическо- го, т.е. список формальных аргументов процедуры можно умень- шить (меньше аргументов — меньше вероятность ошибки):
program tmat5; implicit none interface; subroutine matmul3(a,b,c,l,m,n);
integer l,m,n,a(:,:),b(:,:),c(:,:); end subroutine matmul3
end interface integer, parameter :: ndim=5
integer a(ndim,ndim), b(ndim,ndim), c(ndim,ndim), d(ndim,ndim)
integer l, m, n, i, j; l=2; m=3; n=4
! Моделируем расчёт l, m и n.
a=-1; b=-2; c=-3; d=-4
! Инициализация массивов.
do i=1,l; do j=1,m; a(i,j)=i+j; enddo; enddo do i=1,m; do j=1,n; b(i,j)=i*j; enddo; enddo do i=1,l; write(*,’(" a(",i2,",:)=",5i3)’) i,(a(i,j),j=1,m); enddo do i=1,m; write(*,’(" b(",i2,",:)=",5i3)’) i,(b(i,j),j=1,n); enddo call matmul3(a,b,c,l,m,n)
do i=1,l; write(*,’(" c(",i2,",:)=",5i5)’) i,(c(i,j),j=1,n); enddo end subroutine matmul3(a,b,c,l,m,n); implicit none; integer ndim,l,m,n,i,j,k,s integer a(:,:), b(:,:), c(:,:)
do i=1,l; do k=1,n; s=0; do j=1,m; s=s+a(i,j)*b(j,k); enddo; c(i,k)=s enddo; enddo end a( 1,:)=
2 3
4
a( 2,:)=
3 4
5
b( 1,:)=
1 2
3 4
b( 2,:)=
2 4
6 8
b( 3,:)=
3 6
9 12
c( 1,:)=
20 40 60 80
c( 2,:)=
26 52 78 104 263
7.4
Транспонирование матриц
Функция transpose(matrix) транспонирует матрицу-аргумент, возвра- щая транспонированную матрицу через имя transpose.
! Результат работы trans:
program trans
!
a:
implicit none
!
1 3 5
integer a(2,3), b(3,2), i, j
!
2 4 6
a=reshape((/ 1,2,3,4,5,6 /),(/2,3/))
!
b:
b=transpose(a)
!
1 2
write(*,*) ’ a:’; write(*,’(3i2)’) (a(i,:),i=1,2) !
3 4
write(*,*) ’ b:’; write(*,’(2i2)’) (b(i,:),i=1,3) !
5 6
end
7.5
Функция слияния массивов
Функция merge(tsource,fsource,mask) — поэлементная функция трёх аргументов — массивов одинаковой формы, из которых первые два оди- накового типа, а третий — логического. merge возвращает через своё
имя массив той же формы, в котором элементы первого аргумента заме- нены элементами второго аргумента в том и только в том случае, если соответствующие элементы маски имеют значение .false..
Имя merge иногда именует процедуру слияния двух yпорядоченных массивов в один упорядоченный. К упомянутой сортировке встроенная merge ФОРТРАНа никакого отношения не имеет. Правильнее её назвать функцией совмещения элементов двух массивов.
program tmerge
! tmerg посредством вызова implicit none
! merge получает матрицу C,
integer a(2,3), b(2,3), c(2,3),
i
! элементы которой равны logical l(2,3)
! элементам матриц A или B,
a=reshape((/ 1, 2,
3, 4, 5,
6/),(/2,3/));
! в зависимости от того b=reshape((/ 0, 8, -1, 4, 7, -5/),(/2,3/));
! больше ли элементы матрицы l=a>b
! A соответствующих элементов c=merge(a,b,l)
! матрицы B или меньше.
write(*,’(5x,"a",11x,"b",11x,"l",11x,"c")’)
write(*,’(3i3,3x,3i3,3x,3l3,3x,3i3)’) (a(i,:),b(i,:),l(i,:),c(i,:),i=1,2)
end a
b l
c
1 3
5 0 -1 7
T
T
F
1 3
7 2
4 6
8 4 -5
F
F
T
8 4
6 264
7.6
Функции упаковки и распаковки массивов
Функция pack(array,mask,[vector]) возвращает вектор, собранный из тех и только тех элементов массива array, для которых значения маски mask есть true.
array — массив любой формы и типа. mask — скаляр или логиче- ский массив той же формы, что и array. vector — необязательный ар- гумент (вектор, элементы которого доопределяют хвостовые элементы массива-результата). Размер vector не меньше числа элементов mask со значением .true..
program tpack1; implicit none; integer a(3,3),r1(2),r2(7)
a=0; a(2,2)=1; a(3,3)=2;
r1=pack(a,a/=0); write(*,’(a,2i3)’) ’ r1=’,r1
r2=pack(a,a/=0,(/10,20,30,40,50,60,70/)); write(*,’(a,7i3)’) ’ r2=’,r2
end r1=
1 2
r2=
1 2 30 40 50 60 70
Функция unpack(vector,mask,field) возвращает массив, получен- ный расфасовкой вектора vector в соответствии с содержимым логи- ческого массива mask той же формы, что и unpack.
Именно, содержимое очередного элемента вектора vector присваива- ется тому и только тому элементу unpack, если соответствующий эле- мент mask имеет значение .true.. Если элемент mask имеет значение
.false., то соответствующему элементу unpack присваиваивается либо значение скаляра field, либо значение соответствующего элемента мас- сива field той же формы, что и mask.
program tunpack; implicit none logical mask(2,3); integer :: matr(2,3), vec(3)=((/3,4,5/))
mask=reshape((/.true.,.false.,.false., .true., .true., .false./),(/2,3/))
write(*,’(" mask(1,:)=",3l3/" mask(2,:)=",3l3)’) mask(1,:), mask(2,:)
matr=unpack(vec,mask,77)
write(*,’(" matr(1,:)=",3i3/" matr(2,:)=",3i3)’) matr(1,:), matr(2,:)
end mask(1,:)=
T
F
T
! Под термином "упаковка" часто понимают сокращение mask(2,:)=
F
T
F
! объёма информации без её потери. PACK "трактует"
matr(1,:)=
3 77 5
! термин "упаковка" более широко, допуская априори matr(2,:)= 77 4 77
! и возможную потерю. По сути выполняемой работы
! к функциям PACK и UNPACK на русском языке более
! подходят термины "фильтрация" и "расфасовка".
265
7.7
Сборка массива через добавление измерения
Функция spread( source, dim, ncopies) добавляет ncopies копий ис- ходного массива source вдоль заданного измерения dim массива, полу- чаемого spread.
program tspread; implicit none integer m(3) /1,2,3/, m2 (4,3), m3(4,3,2), i, j, k m2=spread(m,1,4)
write(*,*) ’ m=’,m write(*,’(" m2:")’); write(*,’(3i3)’) ((m2(i,j),j=1,3),i=1,4)
write(*,’("m3: к матрице m2 добавились ещё две такие же,")’)
write(*,’("
сформировав второй и третий листы книги m3:")’)
m3=spread(m2,3,2)
m3(:,:,2)=10*m3(:,:,2) ! меняем содержимое второго и третьего листа,
m3(:,:,3)=20*m3(:,:,3) ! чтобы убедиться, что выводим действительно
! разные листы, а не 3 копии первого.
do k=1,3; write(*,’(i3,a$)’) k,"-й лист"; enddo; write(*,*)
write(*,’(3i3," ",3i3," ",3i3$)’) (((m3(i,j,k),j=1,3),k=1,3),i=1,4)
write(*,*)
end m=
1 2
3
m2:
1 2
3 1
2 3
1 2
3 1
2 3
m3: к матрице m2 добавились ещё две такие же,
сформировав второй и третий листы книги m3:
1-й лист
2-й лист
3-й лист
1 2
3 10 20 30 20 40 60 1
2 3
10 20 30 20 40 60 1
2 3
10 20 30 20 40 60 1
2 3
10 20 30 20 40 60 1. m2=spread(m,1,4) формирует матрицу из четырёх строк и трех столбцов, копируя в неё построчно исходный m четыре раза.
2. m3=spread(m2,3,4) формирует из трёхмерный массив (книгу, со- стоящую из трёх листов), копируя в неё полистно матрицу m2.
3. При указании номера измерения, не согласованного с формой копи- руемого исходного массива, получаем сообщение об ошибке:
m3=spread(m2,1,2)
1
Error: different shape for Array assignment at (1) on dimension 1 (4/2)
266