Файл: А. Б. Шнейвайс основы программирования.pdf

program vec1; implicit none; real p(10); integer :: vi(3)
integer i, j vi=(/7,1,4/)
p=1;
p(vi)=(/61,21,31/)
write(*,’(19x,"Массив p после модификации по векторному индексу vi=(/7,1,4/)")’)
write(*,’(" Индекс массива p:",
10i5)’)(i,i=1,10)
write(*,’(" Содержимое p: ",10f5.0)’) (p(i),i=1,10)
write(*,’(" Индексы секции
:",3i5)’) vi write(*,’(" Содержимое секции p(vi): ",3f5.0)’) p(vi)
vi=(/7,1,7/);
p(vi)=(/61,21,77/)
write(*,’(19x,"Массив p после модификации по векторному индексу vi=(/7,1,7/)")’)
write(*,’(" Индекс массива p:",
10i5)’)(i,i=1,10)
write(*,’(" Содержимое p: ",10f5.0)’) (p(i),i=1,10)
write(*,’(" Индексы секции
:",3i5)’) vi write(*,’(" Содержимое секции p(vi): ",3f5.0)’) p(vi)
p((/7,1,7/))=(/35,6,25/)
write(*,’(19x,"Массив p после модификации по векторному индексу (/7,1,7/)")’)
write(*,’(" Индекс массива p:",
10i5)’)(i,i=1,10)
write(*,’(" Содержимое p: ",10f5.0)’) (p(i),i=1,10)
write(*,’(" Индексы секции
:",3i5)’) vi write(*,’(" Содержимое секции p(vi): ",3f5.0)’) p(vi)
end
Массив p после модификации по векторному индексу vi=(/7,1,4/)
Индекс массива p:
1 2
3 4
5 6
7 8
9 10
Содержимое p:
21.
1.
1.
31.
1.
1.
61.
1.
1.
1.
Индексы секции
:
7 1
4
Содержимое секции p(vi):
61.
21.
31.
Массив p после модификации по векторному индексу vi=(/7,1,7/)
Индекс массива p:
1 2
3 4
5 6
7 8
9 10
Содержимое p:
21.
1.
1.
31.
1.
1.
77.
1.
1.
1.
Индексы секции
:
7 1
7
Содержимое секции p(vi):
77.
21.
77.
Массив p после модификации по векторному индексу (/7,1,7/)
Индекс массива p:
1 2
3 4
5 6
7 8
9 10
Содержимое p:
6.
1.
1.
31.
1.
1.
25.
1.
1.
1.
Индексы секции
:
7 1
7
Содержимое секции p(vi):
25.
6.
25.
Выше в 2b) отмечалось (см.также [2, 3]), что секция с повторяющими- ся значениями индексов не может встречаться в левой [2] (правой [3]) ча- сти оператора присваивания. Как видим, используемая версия gfortranа разрешает подобные присваивания, беря в качестве окончательного зна- чения, то, которое соответствует последнему правому из повторяющихся в списке индексов. Другими словами, контроль ситуации возлагается на человека, а не на компилятор.
213

6.3
Индексный триплет (временные оценки)
1. Шаг индексного триплета равен 1:
program testop1; implicit none
! Программа kk раз выполняет integer, parameter :: n=1000, kk=1000000
! заполнение массива из n real(8) a(n)
! элементов integer k,i; real t(2),
r0, r1, r2, r3, r4, r5, r10
call etime(t,r0); do k=1,kk do i=1,n; a(i)=1;enddo;
! обычным циклом enddo call etime(t,r1); do k=1,kk;
a=2;
enddo ! простым присваиванием call etime(t,r2); do k=1,kk; a(1:n:1)=3;enddo ! на базе индексного триплета call etime(t,r3); do k=1,kk;
a(1:n)=4; enddo !
call etime(t,r4); do k=1,kk; a=(/ (i,i=1,n) /)! через конструктор массива enddo call etime(t,r5); write(*,1000)
write(*,’(a,i7,5(f10.3))’) ’n=’,n,r1-r0,r2-r1,r3-r2,r4-r3,r5-r4 1000 format(1x,11x,’do i=1,n’,5x,’a=1’,3x,’a(1:n:1)=2’,2x,&
&
’a(1:n)=3’,3x,’a=(/ (i=1,n) /)’)
end program testop1
do i=1,n a=1
a(1:n:1)=2
a(1:n)=3
a=(/ (i=1,n) /)
n=
1000 6.993 4.928 4.936 4.978 7.737
-O0
n=
1000 0.597 0.596 0.597 0.597 1.386
-O1
n=
4000 28.980 19.952 20.069 19.957 31.148
-O0
n=
4000 2.354 2.354 2.347 2.340 5.459
-O1 2. Шаг индексного триплета равен 2:
Заполнение нечетных элементы одномерного массива единицами, а четных нулями в современном ФОРТРАНе можно осуществить так a(1:n:2)=1;
a(2:n:2)=0
a(:n:2)=1;
a(2::2)=0
a(::2)=1;
a(:n:-2)=0
Ниже приведена программа, позволяющая оценить быстродействие подобного присваивания по сравнению с обычным оператором цик- ла и с вызовом соответствующей подпрограммы. После исходного текста приведены результаты двух её пропусков. Первый — испол- нимый файл получался при ключе оптимизации -O0; второй — при
-O1.
214

program testop2
! Программа kk раз заполняет два сечения (с чётными и implicit none
! нечетными индексами) массива из n элементов значениями integer, parameter :: n=20000, kk=100000
! 3.7d0 и -8.4d0 соответственно.
real(8) a(n)
integer k,i; real t(2),
r0, r1, r2, r3, r4, r5
call etime(t,r0); do k=1,kk
! Обычный цикл:
do i=1,n,2; a(i)=-8.4d0; enddo;
! нечётные индексы do i=2,n,2; a(i)= 3.7d0; enddo
! чётные индексы enddo call etime(t,r1); do k=1,kk; call init(a,n,1)
! Обычный цикл call init(a,n,2)
! в подпрограмме.
enddo call etime(t,r2); do k=1,kk; a(1:n:2)=-8.4d0 ! Индекс.триплет: нечётные a(2:n:2)= 3.7d0 !
и чётные индексы enddo call etime(t,r3); do k=1,kk; a(:n:2)= -8.4d0
! нечетные индексы a(2::2)= 3.7d0
! чётные индексы enddo call etime(t,r4); do k=1,kk; a(::2)= -8.4d0
! нечетные индексы a(2::2)= 3.7d0
! чётные индексы enddo call etime(t,r5); write(*,1000); write(*,2000)
write(*,’(a,i7,6(1x,f10.2,1x))’) ’n=’,n,r1-r0,r2-r1,r3-r2,r4-r3,r5-r4 1000 format(12x,’do i=1,n,2’,6x,’init’,5x,’a(1:n:2)’,5x,’a(:n:2)’,5x,&
&
’a(::2)’)
2000 format(12x,’do i=2,n,2’,6x,’init’,5x,’a(2:n:2)’,5x,’a(2::2)’,5x,&
&
’a(2::2)’)
end program testop2
subroutine init(a,n,k) ! Подпрограмма init(a,n,k) заполняет при k=1
implicit none
! элементы вектора a c нечётными индексами real(8) a(n)
! 1(2)n константой -8.4,
integer i,n,k
! а при k=2 элементы с чётными индексами real(8) t(2)
! 2(2)n константой 3.7
data t / -8.4, 3.7 /
do i=k,n,2 ; a(i)=t(k); enddo end
Результат при ключе компиляции -O0:
do i=1,n,2
init a(1:n:2)
a(:n:2)
a(::2)
do i=2,n,2
init a(2:n:2)
a(2::2)
a(2::2)
n=
20000 9.27 13.55 11.04 11.10 11.09
Результат при ключе компиляции -O1:
do i=1,n,2
init a(1:n:2)
a(:n:2)
a(::2)
do i=2,n,2
init a(2:n:2)
a(2::2)
a(2::2)
n=
20000 3.55 3.75 3.51 3.54 3.56 215

6.4
Примеры использования сечений
6.4.1
Изменение порядка следования элементов массива program revers; implicit none; integer, parameter :: kmax=10000, nn=100001
integer :: k, w, n; real t1, t2, t3, t4, t5, t6
integer :: p(kmax)
do k=1,kmax; p(k)=k; enddo write(*,’("
k =",(10i5))’) (k,k=1,kmax,1000)
write(*,’("p(k)=",(10i5),"
Время работы ")’) (p(k),k=1,kmax,1000)
do n=1,nn do k=1,kmax/2
w=p(k);
p(k)=p(kmax-k+1); p(kmax-k+1)=w enddo enddo call cpu_time(t2)
write(*,’("p(k)=",(10i5)$)’) (p(k),k=1,kmax,1000); write(*,’(f7.2)’) t2-t1
call cpu_time(t3)
do n=1,nn; call rev(p,kmax); enddo call cpu_time(t4)
write(*,’("p(k)=",(10i5)$)’) (p(k),k=1,kmax,1000); write(*,’(f7.2)’) t4-t3
call cpu_time(t5)
do n=1,nn; p=p(kmax:1:-1);
enddo call cpu_time(t6)
write(*,’("p(k)=",(10i5)$)’) (p(k),k=1,kmax,1000); write(*,’(f7.2)’) t6-t5
end subroutine rev(a,kk)
implicit none integer kk,k,a(kk),w do k=1,kk/2; w=a(k); a(k)=a(kk-k+1); a(kk-k+1)=w; enddo end
Результат при ключе оптимизации -O0:
k =
1 1001 2001 3001 4001 5001 6001 7001 8001 9001
p(k)=
1 1001 2001 3001 4001 5001 6001 7001 8001 9001
Время работы p(k)=10000 9000 8000 7000 6000 5000 4000 3000 2000 1000 3.32
p(k)=
1 1001 2001 3001 4001 5001 6001 7001 8001 9001 6.48
p(k)=10000 9000 8000 7000 6000 5000 4000 3000 2000 1000 2.77
Результат при ключе оптимизации -O1:
k =
1 1001 2001 3001 4001 5001 6001 7001 8001 9001
p(k)=
1 1001 2001 3001 4001 5001 6001 7001 8001 9001
Время работы p(k)=10000 9000 8000 7000 6000 5000 4000 3000 2000 1000 0.76
p(k)=
1 1001 2001 3001 4001 5001 6001 7001 8001 9001 0.76
p(k)=10000 9000 8000 7000 6000 5000 4000 3000 2000 1000 2.04

6.4.2
Заполнение матриц
Здесь демонстрируется, как просто, используя сечения можно в первые три столбца матрицы x(3,6) занести содержимое матрицы p(3,3), а в оставшиеся три содержимое матрицы q program form0
implicit none integer x(3,6); integer p(3,3) /
1,
2,
3,
4,
5,
6,
7,
8,
9/
integer i, j;
integer q(3,3) / 10, 20, 30, 40, 50, 60, 70, 80, 90/
write(*,’(" p=",3i3/2("
",3i3/))’) p write(*,’(" q=",3i3/2("
",3i3/))’) q;
x(1:3,1:3)=p;
x(1:3,4:6)=q write(*,’(" x=",3i3/5("
",3i3/)//)’) x write(*,’(" p=",3i3/2("
",3i3/))’) ((p(i,j),j=1,3),i=1,3)
write(*,’(" q=",3i3/2("
",3i3/))’) ((q(i,j),j=1,3),i=1,3)
write(*,’(" x=",6i3/5("
",6i3/))’) ((x(i,j),j=1,6),i=1,3)
end
Результат:
p=
1 2
3
На первый взгляд интуитивно кажется, что матрицы p и q
4 5
6
выведены верно. И это, действительно, так.
7 8
9
ОДНАКО, ПОМНИМ, что ФОРТРАН при указании в операторах инициализации, ввода и вывода только имени матрицы q= 10 20 30
заполняет её, вводит и выводит её элементы по СТОЛБЦАМ.
40 50 60
Поэтому все выведенные строки есть содержимое СТОЛБЦОВ
70 80 90
МАТРИЦЫ!
Сказанное, в частности, относится и матрице x, в которой x=
1 2
3
согласно описанию ТРИ строки и ШЕСТЬ столбцов.
4 5
6
Так что, в контексте данной программы содержимое каждой
7 8
9
строки любой из матриц находится в соответствующем
10 20 30
выведенном СТОЛБЦЕ. Если же хотим вывести содержимое
40 50 60
матрицы так, чтобы СТРОКИ ВЫВОДА оказались действительно
70 80 90
СТРОКАМИ МАТРИЦЫ, то об этом необходимо позаботиться особо, например, указывая очерёдность вывода элементов в циклоподобном списке вывода:
p=
1 4
7
write(*,*) ((p(i,j),j=1,3),i=1,3)
2 5
8
При такой форме его записи индекс j будет меняться
3 6
9
чаще индекса i, и поэтому в строку ВЫВОДА будет выводиться действительно строка МАССИВА.
q= 10 40 70 20 50 80
write(*,*) ((q(i,j),j=1,3),i=1,3)
30 60 90
x=
1 4
7 10 40 70
write(*,*) ((x(i,j),j=1,6),i=1,3)
2 5
8 20 50 80 3
6 9 30 60 90 217

В примере из программы form0 заполнение массивов p и q выполня- лось при описании их типа посредством добавления списка констант, за- ключенного между двух прямых слешей. При таком задании начальных значений элементов матрицы данные из списка, трактуются по умолча- нию, как расположенные последовательно друг за другом значения эле- ментов очередного столбца матрицы. Современный ФОРТРАН предо- ставляет программисту и альтернативный способ посредством встроен- ной функции reshape.
1. reshape — размещает содержимое элементов массива одной формы нужным образом по элементам массива другой формы.
2. Вызов функции reshape в общем случае имеет вид:
reshape(source, shape [, pad] [,order])
reshape — не подпрограмма, а именно функция, которая возвра- щает через своё имя массив. Формальные аргументы:
source — исходный массив, значения элементов которого должны попасть в массив требуемой параметром shape формы;
shape — вектор, задающий форму массива-результата (очередной элемент shape должен хранить число элементов по очередному из- мерению массива-результата).
pad — Если source по размеру меньше результата, то незатронутые элементы последнего заполняются содержимым pad.
order — указывает последовательность заполнения элементов мас- сива (результата) данными из массива source.
program form1; implicit none integer a(9) / 1, 2, 3, 4, 5, 6, 7, 8, 9 /,
b(3,3), c(9,4), i, j write(*,’(" a=",9i3)’) a b=reshape(a,(/3,3/));
write(*,’(" b="3i3/2("
",3i3/))’) b b=reshape(a,(/3,3/),order=(/1,2/)); write(*,’(" b="3i3/2("
",3i3/))’) b b=reshape(a,(/3,3/),order=(/2,1/)); write(*,’(" b="3i3/2("
",3i3/))’) b c=reshape(a,(/9,4/), pad=(/ -1 /)); write(*,’(" c="9i3/3("
",9i3/))’) c c=reshape(a,(/9,4/),(/ -1,-2/));
write(*,’(" c="9i3/3("
",9i3/))’) c c=reshape(pad=(/ -1,-2,-3,(i,i=-90,-80,1)/),source=a,shape=(/9,4/))
write(*,’(" c="9i3/3("
",9i3/))’) c c=reshape(a,(/9,4/),(/ -1,-2,-3,(i,i=-90,-80,1)/),(/2,1/))
write(*,’(" c="9i3/3("
",9i3/))’) c end
218

a=
1 2
3 4
5 6
7 8
9
Как видим, применение reshape в виде reshape(a,(/3,3/)) привело к тому же b=
1 2
3 4
5 6
результату, что и в программе form0 7
8 9
integer p(3,3) /1,2,3,4,5,6,7,8,9/
b=
1 2
3 4
5 6
Тот же результат и при reshape(a,(/3,3/),order=(/1,2/))
7 8
9
Однако,
reshape(a,(/3,3/),order=(/2,1/))
b=
1 4
7
<---= даёт. Здесь содержимое order
2 5
8
указывает, что при заполнении b быстрее всего должен
3 6
9
меняться индекс её второго измерения, т.е. номер столбца,
так что 2 из А должна попасть во 2-й столбец 1-ой строки.
c=
1 2
3 4
5 6
7 8
9
В матрице С на 27 элементов больше чем в A.
-1 -1 -1 -1 -1 -1 -1 -1 -1
Параметр pad=(/-1/) при вызове
-1 -1 -1 -1 -1 -1 -1 -1 -1
c=reshape(a,(/4,9/),pad=(/ -1/))
-1 -1 -1 -1 -1 -1 -1 -1 -1
требует, чтобы элементы C, в которые ничего не попало из A, заполнились -1.
c=
1 2
3 4
5 6
7 8
9
Если хотим, чтобы при их заполнении
-1 -2 -1 -2 -1 -2 -1 -2 -1
чередовались -1 и -2, то используем
-2 -1 -2 -1 -2 -1 -2 -1 -2
-1 -2 -1 -2 -1 -2 -1 -2 -1
c=reshape(a,(/4,9/),pad=(/-1,-2/))
c=
1 2
3 4
5 6
7 8
9
Если же при заполнении оставшихся нужна
-1 -2 -3-90-89-88-87-86-85
некоторая особая цикличность, то можно
-84-83-82-81-80 -1 -2 -3-90
использовать соответствующий циклоподобный
-89-88-87-86-85-84-83-82-81
элемент в конструкторе массива.
c=
1 5
9-90-86-82 -2-88-84
Результат подключения к предыдущему вызову
2 6 -1-89-85-81 -3-87-83
reshape ключевого параметра:
3 7 -2-88-84-80-90-86-82
order=(/2,1/).
4 8 -3-87-83 -1-89-85-81 3. Как видно из текста программы form1 при вызове reshape воз- можны различные способы задания её аргументов:
1) b=reshape(a,(/3,3/)) — заданы лишь два обязательных: ис- ходный массив, поставщик данных и форма массива-результата.
2) b=reshape(a,(/3,3/),order=(/2,1/)). Ключевой аргумент order указывает, что второй индекс при заполнении массива- результата должен меняться быстрее, а первый — медленнее.
3) c=reshape(a,(/9,4/), pad=(/ -1 /)) — через ключевой аргу- мент pad передаётся “наполнитель” ;
219

4) c=reshape(a,(/9,4/),(/ -1,-2/)) — все аргументы позицион- ные (трактовка каждого определяется его позицией).
5) При вызове c=reshape(pad=(/-1,-2,-3,(i,i=-90,-80,1)/),source=a,shape=(/9,4/))
очерёдность ключевых аргументов не имеет значения. При задании “наполнителя” использован циклоподобный список.
6)
c=reshape(a,(/9,4/),(/-1,-2,-3, (i,i=-90,-80,1)/),(/2,1/))
при вызове нет ключевых слов: все аргументы — позиционные.
4. При вызове любой процедуры можно использовать ключевые име- на аргументов. Ключевое имя — имя формального параметра, со- единённое с именем фактического знаком операции присваивания.
program key; implicit none; integer a, b, r interface function f(x,y); integer y,f; integer, optional :: x; end function f;
end interface
! Результат a=85; b=17; r=f( a , b ); write(*,*) r
!
5
r=f(x=a,y=b); write(*,*) r
!
5
r=f(y=b,x=a); write(*,*) r
!
5
r=f(y=a,x=b); write(*,*) r
!
0
r=f( a, y=a); write(*,*) r
!
1
r=f(y=0);
write(*,*) r
!
333
!
r=f(x=a, b )
Позиционный аргумент недопустим после ключевого end function f(x,y); implicit none; integer f, y; integer, optional :: x if (y==0) then; f=333
!
optional else; f=x/y
! сообщает компилятору, что аргумент x --- endif
!
необязательный end
Подробнее см. [2, 3]
5. Ешё раз: reshape — именно функция, которая возвращает через своё имя массив. Современный ФОРТРАН предоставляет програм- мисту возможность описывать функции, обладающие таким свой- ством. Ниже приводится программа testfunarr с тремя вариантами реализации подобного описания:
1) внутренней функцией vecpow0 (описана в главной программе);
2) внешней функцией vecpow1 (описанной в отдельном файле);
3) модульной функцией vecpow2 (описана в module myvec).
220

и соответствующие результаты её работы. Каждая функция решает одну и ту же задачу — по заданному целому n вычисляет вектор,
первые n элементов которого содержат соответствующие степени двойки.
program testfunarr; use myvec; implicit none interface function vecpow1(n); integer n;
real(8) vecpow1(n);
end function vecpow1;
end interface integer, parameter :: n=5; integer i real(8) a(n), b(n), c(n)
a=vecpow0(n); b=vecpow1(n); c=vecpow2(n)
write(*,’(1x," i",13x,"a",14x,"b",13x,"c")’)
write(*,’(i3,3f15.2)’) (i,a(i),b(i),c(i),i=1,n)
contains function vecpow0(n); implicit none integer n, i real(8) vecpow0(n); do i=1,n;
vecpow0(i)=2.0**i; enddo end function vecpow0
end function vecpow1(n)
implicit none integer n, i real(8) vecpow1(n); do i=1,n; vecpow1(i)=2.0**i; enddo end module myvec; implicit none; contains function vecpow2(n)
integer n, i real(8) vecpow2(n); do i=1,n; vecpow2(i)=2.0**i; enddo end function vecpow2
end module myvec i
a b
c
! Результат работы
1 2.00 2.00 2.00
! программы testfunarr
2 4.00 4.00 4.00 3
8.00 8.00 8.00 4
16.00 16.00 16.00 5
32.00 32.00 32.00 221

6.4.3
Сечение в качестве параметра процедуры
Процедура поиска в каждой строке матрицы первого отрицательного элемента (пример из [3]).
1. Реализация через внутреннюю функцию (т.е. функцию, описание которой находится в теле вызывающей программы; не путать со встроенной функцией).
program negel0; implicit none integer, parameter :: m=10, n=5; integer i integer b(m,n)
read (*,’(10i4)’) b; write(*,’(10i4)’) b do i=1,m write(*,’(i3,a,10i3)’) i,’...’,fneg1(b(i,:))
enddo contains integer function fneg1(a); implicit none integer a(:), k fneg1=0
do k=1,size(a)
if (a(k)<0) then fneg1=a(k)
return endif enddo end function fneg1
end
Результат:
1 -11 21 31 41 51 61 71 -81 91
! Обратите внимание, что при
2 12 22 -32 42 52 62 -72 82 92
! выводе матрицы b eё строки
3 13 23 33 43 -53 63 -73 83 93
! оказываются в столбцах
4 14 24 34 -44 -54 64 74 84 -94
! выведенного, т.к. элементы
5 15 25 35 45 55 65 75 -85 -95
! ФОРТРАН-матрицы в
1...
0
! оперативной памяти
2...-11
! расположены по столбцам.
3...
0 4...-32 5...-44 6...-53 7...
0 8...-72 9...-81 10...-94 222

2. Реализация через внешнюю функцию с использованием сечений.
program negel1; implicit none
! Функция с перенимающим interface
! форму массивом должна в integer function fneg1(a); integer a(:); ! вызывающей её программе end function fneg1
! иметь ЯВНО заданный end interface
!
интерфейс.
integer, parameter :: m=10, n=5; integer i integer b(m,n)
! Программе negel0 он не read (*,’(10i4)’) b; write(*,’(10i4)’) b ! нужен, ибо fneg1
была do i=1,m
! внутренней, т.е. имела
! имела ЯВНЫЙ интерфейс.
write(*,’(i3,a,10i3)’) i,’...’,fneg1(b(i,:))
enddo end integer function fneg1(a); implicit none integer a(:), k fneg1=0
do k=1,size(a)
if (a(k)<0) then fneg1=a(k); return endif enddo end
3. Реализация через внешнюю функцию без использования сечений.
program negel2; implicit none integer, parameter :: m=10, n=5; integer i, j, fneg2
integer a(m), b(m,n)
read (*,’(10i4)’) b; write(*,’(10i4)’) b do i=1,m do j=1,m a(j)=b(i,j)
enddo write(*,’(i3,a,10i3)’) i,’...’,fneg2(a,m)
enddo end function fneg2(a,m); implicit none
! Забавный факт. Если имя fneg2
integer fneg2, m, a(m), k
! убрать из списка integer и fneg2=0
! integer поместить перед do k=1,m
! function, то вывод для j=1
if (a(k)<0) then
! осуществится *** а не нулём.
fneg2=a(k)
! Почему? Проверить!
return endif enddo end
223