Форум программистов, компьютерный форум, киберфорум
Наши страницы
palva
Войти
Регистрация
Восстановить пароль
Оценить эту запись

Работа с массивами в современном фортране

Запись от palva размещена 03.01.2014 в 23:55
Обновил(-а) palva 03.01.2014 в 23:57

В новом фортране появились операции над массивами. Их можно поэлементно складывать, умножать и т. д. Матрицы и вектора умножаются не в математическом, а в поэлементном смысле. Для матричного умножения все же приходится использовать встроенные функции. Сечения массива позволяют работать с подматрицами – клетками матрицы. Клетки можно не только извлекать из матрицы, клеткам, как и любым матрицам можно что-то присвоить. Мы хотим продемонстрировать этот механизм на примере одного алгоритма обращения матриц. Процитируем книгу Ахо А., Хопкрофт Дж., Ульман Дж. Построение и анализ эффективных алгоритмов. (М.: 1979) с. 262.

Цитата:
Лемма 6.5. Пусть матрица A разбита так:

http://www.cyberforum.ru/cgi-bin/latex.cgi?<br />
\begin{bmatrix}<br />
A_{11} & A_{12}\\ <br />
A_{21} & A_{22}<br />
\end{bmatrix}

Предположим, что существует http://www.cyberforum.ru/cgi-bin/latex.cgi?A_{11}^{-1}. Обозначим http://www.cyberforum.ru/cgi-bin/latex.cgi?A_{22}-A_{21}A_{11}^{-1}A_{12} через http://www.cyberforum.ru/cgi-bin/latex.cgi?\Delta и допустим, что http://www.cyberforum.ru/cgi-bin/latex.cgi?\Delta^{-1} существует. Тогда

http://www.cyberforum.ru/cgi-bin/latex.cgi?<br />
A^{-1}=\begin{bmatrix}<br />
A_{11}^{-1}+A_{11}^{-1}A_{12}\Delta^{-1}A_{21}A_{11}^{-1} & -A_{11}^{-1}A_{12}\Delta^{-1}\\ <br />
-\Delta^{-1}A_{21}A_{11}^{-1} & \Delta^{-1}<br />
\end{bmatrix}
Далее в книге идет доказательство, которое мы здесь приводить не будем.

Предположим, нам надо обратить матрицу четвертого порядка. По изложенной Лемме такую задачу можно свести к двум обращениям матриц второго порядка и нескольким умножениям матриц. А уж для обращения матриц второго порядка нетрудно написать очевидную подпрограмму. Вопрос о том, будут ли матрицы второго порядка обратимыми, мы оставим открытым. Скорее всего будут. В крайнем случае программа, которую мы ниже предлагаем, при специально подобранных плохих данных упадет из-за деления на нуль.

Fortran
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
! И. В. Проскуряков.
! Сборник задач по линейной алгебре. 2005.
! № 857. Найти обратную матрицу к матрице
! |3  3 -4 -3|
! |0  6  1  1|
! |5  4  2  1|
! |2  3  3  2|
! В задачнике имеется следующий ответ
! | -7   5  12  -19|
! |  3  -2  -5    8|
! | 41 -30 -69  111|
! |-59  43  99 -159|
real*8 :: a(4,4) = &
! (Обратите внимание, что матрица в фортране записывается по столбцам.)
    (/ 3, 0, 5, 2, &
       3, 6, 4, 3, &
      -4, 1, 2, 3, &
      -3, 1, 1, 2 /)
interface
    function inv2(a)
    real(8), dimension(2,2) :: a, inv2
    end function inv2
end interface
integer i,j
real*8 :: a11_inv(2,2), d_inv(2,2)
 
a11_inv=inv2(a(1:2,1:2))
d_inv=inv2(a(3:4,3:4)-matmul(a(3:4,1:2),matmul(a11_inv,a(1:2,3:4))))
a(1:2,1:2)=a11_inv+matmul(a11_inv,matmul(a(1:2,3:4), &
                            matmul(d_inv,matmul(a(3:4,1:2),a11_inv))))
a(1:2,3:4)=-matmul(a11_inv,matmul(a(1:2,3:4),d_inv))
a(3:4,1:2)=-matmul(d_inv,matmul(a(3:4,1:2),a11_inv))
a(3:4,3:4)=d_inv
! Результат
print '(4F7.1)', ((a(i,j), j=1,4), i=1,4)
end
    
function inv2(a)
! Обращение матрицы порядка 2
real(8) :: a(2,2), inv2(2,2), det
inv2(1,1)=a(2,2)
inv2(1,2)=-a(1,2)
inv2(2,1)=-a(2,1)
inv2(2,2)=a(1,1)
det=1.0/(a(1,1)*a(2,2)-a(1,2)*a(2,1))
inv2=inv2*det
end function inv2
Вряд ли приведенный алгоритм имеет практическую ценность. Книга, из которой он взят, посвящена теоретической эффективности алгоритмов. Для обращения матриц даже такого небольшого порядка как 4 практичнее использовать готовые библиотеки. Хотя заняться замерами эффективности было бы интересно.

Если же говорить об эффективности, то сама реализация, пусть и прозрачная, вызывает сомнения с точки зрения эффективности. Про эффективность оптимизатора Intel-Fortran сложены легенды. Но в данном случае человек может заметить такие свойства задачи, которых не заметит никакой оптимизатор. Вот здесь http://www.delphimaster.net/view/15-1181196883/0-40 в 2007 году я обсуждал этот алгоритм и пытался реализовать его наиболее эффективно.
Размещено в Без категории
Просмотров 1058 Комментарии 0
Всего комментариев 0
Комментарии
 
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Рейтинг@Mail.ru