Как известно, массивов массивов в Fortran нет. Например, в С++ и Java есть массивы указателей на массивы, но нет истинно многомерных массивов. Зато это позволяет создавать в них рваные массивы (jagged arrays), содержащие разное число элементов в каждой строке. Рваные массивы могут оказаться полезными при формировании симметричных матриц.

Для того, чтобы в Fortran организовать рваные массивы, нужно создать производный тип-строку, содержащий всего один компонент -- размещаемый массив, а после объявить массив таких производных-типов-строк. Это будет массив строк, каждая из которых может размещаться со своей длиной, а значит, этот "двумерный массив" будет рваным.

При выводе симметричных матриц воспользуемся расширением над стандартом языка: в дескрипторах преобразования везде вместо буквальных констант можно использовать целочисленные выражения заключённые в угловые скобки.

Здесь n и m -- целочисленные выражения. Это расширение в ряде случаев позволяет не прибегать к программированию спецификации формата.

Рассматриваемое расширение, хоть и является расширением над стандартом F03 (на сколько я знаю), позволяет часто подобным образом сократить код по выводу чисел. Использовать его следует на своё усмотрение (например, не использовать для переносимости кода). Мы же воспользуемся им при выводе симметричных матриц на экран, чтобы получить элегантный код.

Рассмотрим пример реализации, заполнения и вывода симметричной матрицы согласно изложенному подходу.

Вот, что будет на экране:

Прошу обратить внимание и изучить строку

Здесь переменная неявного цикла присутствует и в дескрипторе формата, т.е. для каждой строки вывода дескриптор формата будет свой -- он будет содержать разное число выводимых значений, i.
Если немного доработать оператор вывода

то он будет выводить целиком симметричную матрицу:

Тут полезно сделать одно замечание. Чтобы выводить матрицу полностью, мы добавили в неявном цикле ещё один, вложенный цикл. При этом его границы и выводимые значения зависят от переменной i внешнего по отношению к нему неявного цикла: вложенный цикл добавляет для вывода i-ые элементы N - i + 1 последних строк. Первое, что хочется сделать после этого, это избавиться от вложенного цикла, добавив его выводимые значения к внешнему циклу, выражая j через i и используя сечения массива. Попробуем так преобразовать оператор полного вывода симметричной матрицы:

в оператор без вложенного неявного цикла:

Кажется, так просто, но компилятор тут же заругается: компонент справа от массива не может быть ссылкой или быть размещаемым. Думаю, это связано с тем, что не известно заранее, ссылается ли на что-то или размещен ли этот компонент.

Полагаю, ни один язык в мире не способен одним оператором вывести симметричную матрицу на экран, даже.. Fortran 2003.
Нам же удалось это сделать.

Здесь мы рвали строки матрицы -- они имели каждая свою длину: от 1 до 24. Мы делали так. потому что так привычнее -- первый индекс отвечает за номер строки, второй -- за номр столбца.

Искушённый Fortran-программист подозревает, что большая часть операций с матрицами, например, сингулярное разложение (SVD), а не какое-нибудь там исключение методом Гаусса, работает со столбцами матрицы. По этой причине можно использовать сфомированный тип не как массив строк разной длины, а как массив столбцов разной длины (достаточно заменить имя типа Row на Column). Так, в операциях над столбцами будут участвовать элементы столбцов, лежащие в памяти в одном массиве-столбце, а значит будет возможно поддержать векторизацию, задавая в операциях целиком массив-столбец.

Так или иначе, что делать рваными, строки или столбцы матрицы, определяется задачей, а точнее, методом её решения. Так что feel free.

Интересно посмотреть, есть ли другие удобные реализации симметричных матриц, да и обсудить, в общем.

Это сообщение отредактировал(а) Cr@$h - 25.9.2006, 01:14

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

Ну, да. Для векторизации важно, чтобы элементы лежали друг за другом в памяти, тогда будет возможна выборка сразу нескольких чисел, например 4. Это будет первый операнд. Также выбирается 4 значения второго операнда. И те, и другие помещаются в 4 АЛУ по работе с плавающей точкой практически так, как если бы помещались в одно (параллельно это делается). Затем параллельно все АЛУ выполняют одну и ту же операцию над своими операндами. Получается, скорость увеличивается примерно в 4 раза. Для осуществления этого в коде должны содержаться специальные операции: SSE (SIMD ... Extension). Просто тут выполняется одна SSE-операция. SIMD -- абревиатура параллельной архитектуры Single Instruction Multiple Data. То и означает -- одна операция обрабтывает сразу целый набор данных. Некоторые компиляторы поддерживают автоматическую векторизацию (Automatic Vectorizer), вставляя, где надо и где можно, такие инструкции. Например, это делает Intel Fortran Compiler.

Вот более подробная информация об этом.

Давайте вернёмся к нашим баранам. Итак, в первый раз мы сформировали тип симметричной матрицы как

Для безопасности и надёжности мы используем размещаемые массивы. Идея создания симметричной матрицы проста как апельсин. Задаём производный тип Row, содержащий всего один компонент -- размещаемый массив, Row. Его имя подсказывает, что он позиционируется нами как строка. Затем, объявляем размещаемый массив, каждый элемент которого имеет производный тип Row. При работе с таким массивом мы не сможем делать сечения строк, например, присваивать первым элементам строк с помощью сечения:

Любой уважающий себя компилятор сразу заругается (IVFC):

Не знаю, на что ссылается компилятор  , но именно в п. 6.1.2 Fortran 2003 Draft International Standard сказано:

Вот, как в нём определены компоненты структуры:

Теперь можно нормально осознать это. В обращении должно быть не более одного.. хм.. данного с ненулевым рангом. Компонент справа от обращения с ненулевым рангом не должен иметь атрибут allocatable или pointer. Иными словами, в сечениях массивов структур нельзя обращаться к размещаемым или ссылочным компонентам этих структур.

Для нашей симметричной матрицы это означает, что мы не можем делать сечение строк. Но это и не важно с точки зрения векторизации, ведь, даже если бы запись SimmetricMatrix(1:4)%Row(1) = Vector(1:4) была ликвидна, элементы SimmetricMatrix(1)%Row(1), ..., SimmetricMatrix(4)%Row(1) располагаются совсем не рядом в памяти, а для SSE-операции это необходимо. Теперь мы видим, что построенный тип симметричной матрицы поддерживает векторизацию только для сечений строк матрицы:

Если в применяемом алгоритме происходит работа преимущественно со стобцами матрицами, то для поддержки векторизации нужно поменять роли: сделать массив рваных столбцов.

Тогда код

будет эквивалентен по назначению коду (который не компилируется) для предыдущей реализации симметричной матрицы

И при этом, естественно, векторизация будет возможна: теперь первые 4 элемента первого столбца расположены друг за дружкой в памяти.

В зависимости от алгоритма стоит применять ту или иную релизацию симметричной матрицы. На сколько это важно? Ну, например, при правильном выборе векторизация будет использоваться чаще при работе метода.

M Cr@$h

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