Версия для печати темы Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > Fortran > [Pro] Этюды. Сюжет 2. Симметричные матрицы:

Автор: Cr@$h 24.9.2006, 00:55

Как известно, массивов массивов в Fortran нет. Например, в С++ и Java есть массивы указателей на массивы, но нет истинно многомерных массивов. Зато это позволяет создавать в них рваные массивы (jagged arrays), содержащие разное число элементов в каждой строке. Рваные массивы могут оказаться полезными при формировании симметричных матриц. Для того, чтобы в Fortran организовать рваные массивы, нужно создать производный тип-строку, содержащий всего один компонент -- размещаемый массив, а после объявить массив таких производных-типов-строк. Это будет массив строк, каждая из которых может размещаться со своей длиной, а значит, этот "двумерный массив" будет рваным. При выводе симметричных матриц воспользуемся расширением над стандартом языка: в дескрипторах преобразования везде вместо буквальных констант можно использовать целочисленные выражения заключённые в угловые скобки. Код print '(f<n>.<m>)' Здесь n и m -- целочисленные выражения. Это расширение в ряде случаев позволяет не прибегать к программированию спецификации формата. Рассматриваемое расширение, хоть и является расширением над стандартом F03 (на сколько я знаю), позволяет часто подобным образом сократить код по выводу чисел. Использовать его следует на своё усмотрение (например, не использовать для переносимости кода). Мы же воспользуемся им при выводе симметричных матриц на экран, чтобы получить элегантный код. Рассмотрим пример реализации, заполнения и вывода симметричной матрицы согласно изложенному подходу. Код program Simmetric     implicit none     ! Строка.     type Row         integer, allocatable :: Row(:)     end type Row     integer :: N = 24, i, j     ! Симметричная матрица -- массив строк.     type(Row), allocatable :: SimmetricMatrix(:)     ! Размещаем симметричную матрицу.     allocate( SimmetricMatrix(N) )     ! Заполняем симметричную матрицу.     ! Будем её строить как нижнюю треугольную матрицу.     do i = 1, N         ! Размещаем i-ую строку.         allocate( SimmetricMatrix(i)%Row(i) )               ! Заполняем её слева направо натуральными числами.         ! Используем неявый цикл.         SimmetricMatrix(i)%Row = [ (j, j = 1, i) ]     end do     ! Выводим симметричную матрицу на экран, используя неявный цикл.     ! Прелесть в том, что переменная неявного цикла участвует также и в     ! дескрипторе формата, задавая каждый раз своё число целых элементов на строке.     print '(<i>i3)', (SimmetricMatrix(i)%Row(1:i), i = 1, N)     read * end program Simmetric Вот, что будет на экране: Код   1   1  2   1  2  3   1  2  3  4   1  2  3  4  5   1  2  3  4  5  6   1  2  3  4  5  6  7   1  2  3  4  5  6  7  8   1  2  3  4  5  6  7  8  9   1  2  3  4  5  6  7  8  9 10   1  2  3  4  5  6  7  8  9 10 11   1  2  3  4  5  6  7  8  9 10 11 12   1  2  3  4  5  6  7  8  9 10 11 12 13   1  2  3  4  5  6  7  8  9 10 11 12 13 14   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Прошу обратить внимание и изучить строку Код print '(<i>i3)', (SimmetricMatrix(i)%Row(1:i), i = 1, N) Здесь переменная неявного цикла присутствует и в дескрипторе формата, т.е. для каждой строки вывода дескриптор формата будет свой -- он будет содержать разное число выводимых значений, i. Если немного доработать оператор вывода Код print '(<N>i3)', (SimmetricMatrix(i)%Row(1:i), (SimmetricMatrix(j)%Row(i), j = i + 1, N), i = 1, N) то он будет выводить целиком симметричную матрицу: Код   1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1   1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2   1  2  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3   1  2  3  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4   1  2  3  4  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5  5   1  2  3  4  5  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6  6   1  2  3  4  5  6  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7   1  2  3  4  5  6  7  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8   1  2  3  4  5  6  7  8  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9  9   1  2  3  4  5  6  7  8  9 10 10 10 10 10 10 10 10 10 10 10 10 10 10 10   1  2  3  4  5  6  7  8  9 10 11 11 11 11 11 11 11 11 11 11 11 11 11 11   1  2  3  4  5  6  7  8  9 10 11 12 12 12 12 12 12 12 12 12 12 12 12 12   1  2  3  4  5  6  7  8  9 10 11 12 13 13 13 13 13 13 13 13 13 13 13 13   1  2  3  4  5  6  7  8  9 10 11 12 13 14 14 14 14 14 14 14 14 14 14 14   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 15 15 15 15 15 15 15 15 15   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 16 16 16 16 16 16 16 16   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 17 17 17 17 17 17 17   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 18 18 18 18 18 18   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 19 19 19 19 19   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 20 20 20 20   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 21 21 21   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 22 22   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 23   1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Тут полезно сделать одно замечание. Чтобы выводить матрицу полностью, мы добавили в неявном цикле ещё один, вложенный цикл. При этом его границы и выводимые значения зависят от переменной i внешнего по отношению к нему неявного цикла: вложенный цикл добавляет для вывода i-ые элементы N - i + 1 последних строк. Первое, что хочется сделать после этого, это избавиться от вложенного цикла, добавив его выводимые значения к внешнему циклу, выражая j через i и используя сечения массива. Попробуем так преобразовать оператор полного вывода симметричной матрицы: Код print '(<N>i3)', (SimmetricMatrix(i)%Row(1:i), (SimmetricMatrix(j)%Row(i), j = i + 1, N), i = 1, N) в оператор без вложенного неявного цикла: Код print '(<N>i3)', (SimmetricMatrix(i)%Row(1:i), SimmetricMatrix(i+1:N)%Row(i), i = 1, N) Кажется, так просто, но компилятор тут же заругается: компонент справа от массива не может быть ссылкой или быть размещаемым. Думаю, это связано с тем, что не известно заранее, ссылается ли на что-то или размещен ли этот компонент. Полагаю, ни один язык в мире не способен одним оператором вывести симметричную матрицу на экран, даже.. Fortran 2003. Нам же удалось это сделать. Здесь мы рвали строки матрицы -- они имели каждая свою длину: от 1 до 24. Мы делали так. потому что так привычнее -- первый индекс отвечает за номер строки, второй -- за номр столбца. Искушённый Fortran-программист подозревает, что большая часть операций с матрицами, например, сингулярное разложение (SVD), а не какое-нибудь там исключение методом Гаусса, работает со столбцами матрицы. По этой причине можно использовать сфомированный тип не как массив строк разной длины, а как массив столбцов разной длины (достаточно заменить имя типа Row на Column). Так, в операциях над столбцами будут участвовать элементы столбцов, лежащие в памяти в одном массиве-столбце, а значит будет возможно поддержать векторизацию, задавая в операциях целиком массив-столбец. Так или иначе, что делать рваными, строки или столбцы матрицы, определяется задачей, а точнее, методом её решения. Так что feel free. Интересно посмотреть, есть ли другие удобные реализации симметричных матриц, да и обсудить, в общем.

Автор: Spike 26.9.2006, 12:44

Насколько я помню, в памяти фортран массивы размещает именно по столбцам

Автор: Cr@$h 27.9.2006, 01:30

Цитата(Spike @  26.9.2006,  13:44 ) Насколько я помню, в памяти фортран массивы размещает именно по столбцам Ну, да. Для векторизации важно, чтобы элементы лежали друг за другом в памяти, тогда будет возможна выборка сразу нескольких чисел, например 4. Это будет первый операнд. Также выбирается 4 значения второго операнда. И те, и другие помещаются в 4 АЛУ по работе с плавающей точкой практически так, как если бы помещались в одно (параллельно это делается). Затем параллельно все АЛУ выполняют одну и ту же операцию над своими операндами. Получается, скорость увеличивается примерно в 4 раза. Для осуществления этого в коде должны содержаться специальные операции: SSE (SIMD ... Extension). Просто тут выполняется одна SSE-операция. SIMD -- абревиатура параллельной архитектуры Single Instruction Multiple Data. То и означает -- одна операция обрабтывает сразу целый набор данных. Некоторые компиляторы поддерживают автоматическую векторизацию (Automatic Vectorizer), вставляя, где надо и где можно, такие инструкции. Например, это делает Intel Fortran Compiler. Цитата Automatic Vectorizer parallelizes code and aligns data, including loop peeling to generate aligned loads and loop unrolling to match the prefetch of a full cache line. Вот более подробная информация об этом. Цитата Automatic Vectorizer Vectorization automatically parallelizes code to maximize underlying processor capabilities. This advanced optimization analyzes loops and determines when it is safe and effective to execute several iterations of the loop in parallel by utilizing MMX™, SSE, SSE2, and SSE3 instructions. Figure is a graphical representation of a vectorized loop that shows four iterations computed with one SSE2 operation. Use vectorization to optimize your application code and take advantage of these new extensions when running on Intel processors. Features include support for advanced, dynamic data alignment strategies, including loop peeling to generate aligned loads and loop unrolling to match the prefetch of a full cache line. Figure. The Vectorizer in action. Давайте вернёмся к нашим баранам. Итак, в первый раз мы сформировали тип симметричной матрицы как Код ! Строка. type Row     integer, allocatable :: Row(:) end type Row ! Симметричная матрица -- массив строк. type(Row), allocatable :: SimmetricMatrix(:) Для безопасности и надёжности мы используем размещаемые массивы. Идея создания симметричной матрицы проста как апельсин. Задаём производный тип Row, содержащий всего один компонент -- размещаемый массив, Row. Его имя подсказывает, что он позиционируется нами как строка. Затем, объявляем размещаемый массив, каждый элемент которого имеет производный тип Row. При работе с таким массивом мы не сможем делать сечения строк, например, присваивать первым элементам строк с помощью сечения: Код SimmetricMatrix(1:4)%Row(1) = Vector(1:4) Любой уважающий себя компилятор сразу заругается (IVFC): Цитата Error: The part-name to the right of a part-ref with nonzero rank has the ALLOCATABLE attribute (6.1.2).[ROW] Не знаю, на что ссылается компилятор  , но именно в п. 6.1.2 Fortran 2003 Draft International Standard сказано: Цитата There shall not be more than one part-ref with nonzero rank. A part-name to the right of a part-ref with nonzero rank shall not have the ALLOCATABLE or POINTER attribute. Вот, как в нём определены компоненты структуры: Цитата data-ref    is part-ref [ % part-ref ] ... part-ref     is part-name [ ( section-subscript-list ) ] Теперь можно нормально осознать это. В обращении должно быть не более одного.. хм.. данного с ненулевым рангом. Компонент справа от обращения с ненулевым рангом не должен иметь атрибут allocatable или pointer. Иными словами, в сечениях массивов структур нельзя обращаться к размещаемым или ссылочным компонентам этих структур. Для нашей симметричной матрицы это означает, что мы не можем делать сечение строк. Но это и не важно с точки зрения векторизации, ведь, даже если бы запись SimmetricMatrix(1:4)%Row(1) = Vector(1:4) была ликвидна, элементы SimmetricMatrix(1)%Row(1), ..., SimmetricMatrix(4)%Row(1) располагаются совсем не рядом в памяти, а для SSE-операции это необходимо. Теперь мы видим, что построенный тип симметричной матрицы поддерживает векторизацию только для сечений строк матрицы: Код SimmetricMatrix(1)%Row(1:4) = Vector(1:4) Если в применяемом алгоритме происходит работа преимущественно со стобцами матрицами, то для поддержки векторизации нужно поменять роли: сделать массив рваных столбцов. Код ! Строка. type Column     integer, allocatable :: Column(:) end type Column ! Симметричная матрица -- массив строк. type(Column), allocatable :: SimmetricMatrix(:) Тогда код Код SimmetricMatrix(1)%Column(1:4) = Vector(1:4) будет эквивалентен по назначению коду (который не компилируется) для предыдущей реализации симметричной матрицы Код SimmetricMatrix(1:4)%Row(1) = Vector(1:4) И при этом, естественно, векторизация будет возможна: теперь первые 4 элемента первого столбца расположены друг за дружкой в памяти. В зависимости от алгоритма стоит применять ту или иную релизацию симметричной матрицы. На сколько это важно? Ну, например, при правильном выборе векторизация будет использоваться чаще при работе метода.

Автор: Cr@$h 3.10.2006, 04:47

M Cr@$h Эта тема продолжила, цикл сюжетных тем "Этюды для Fortran-пограммистов". Тема выделилас из первой сюжетной темы "Вывод вещественных чисел без лишних пробелов и нулей". Конечно, эта тема открыта для обсуждений, вопросов и новых реализаций, тем более, после такого материала. Наверное, со временем такие темы в переработанном и аккумулированном виде хорошо было бы помещать на Fortran Vingrad Wiki (Fortran Viki). Создавайте или предлагайте новые интересные темы, которые можно будет так же хорошо и занимательно разобрать и включить в цикл этюдов для Fortran-программистов. С длинными названиями тем и их описаний помогу;-)

Powered by Invision Power Board (http://www.invisionboard.com)
© Invision Power Services (http://www.invisionpower.com)