[Xenia 38]

Вам надо на каждое собственное значение по несколько таких факторизаций выполнить, то есть по памяти Вы 2*N*N*K раз обратитесь, где N - размерность матрицы, K - среднее число итераций на каждое собственное значение.

Отчего же, к моменту вычисления собственных значений матрица уже сведена к тридиагональному виду, и существует виде всего двух векторов - главной диагонали и ее поддиагонали (вторую поддиагональ не хранят из-за ее полной симметричности первой). Т.е. итерации здесь идут на "двух палочках" :), а съезжающий вниз боковой выступ динамически храниться в отдельной переменной (а у меня в регистре FPU87). Там в среднем не более 2-х итераций на каждое собственное значение (к хвосту быстрее).

 

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

 

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

 

Если бы Вы одновременно итерировали бы несколько ( собственных значений, то обращений по памяти у Вас было бы всего-то 2*N*N*K/B при том же количестве арифметических операций (7*N*N*K), но, так как доступ к памяти занимает много тактов при больших матрицах, это время будет определяющим в этом алгоритме, Вы бы получили решение несравненно быстрее.

Так не получится, т.к. QR/QL - это все-таки метод исчерпывания, в котором собственные значения вычисляются в монотонном порядке, причем обычно сначала ловятся бОльшие по абсолютной величине. Не "исчерпав" большее, трудно найти меньшее, т.к. нет уверенности, что итерационный процесс, начатый с разными стартовыми условиями (начальными сдвигами) сойдется к разным собственным значениям, а не к большему среди них.

[iiv 19]

Отчего же, к моменту вычисления собственных значений матрица уже сведена к тридиагональному виду, и существует виде всего двух векторов - главной диагонали и ее поддиагонали (вторую поддиагональ не хранят из-за ее полной симметричности первой). Т.е. итерации здесь идут на "двух палочках" :), а съезжающий вниз боковой выступ динамически храниться в отдельной переменной (а у меня в регистре FPU87). Там в среднем не более 2-х итераций на каждое собственное значение (к хвосту быстрее).

 

начнем здесь. Вы говорили, что поместили все в фпустек, а именно факторизацию этой трехдиагональной матрицы, заданной главной и поддиагоналями. На каждое собственное значение Вам надо однажды прогнать факторизацию по этой матрице - N*2 обращений по памяти, и около 6*N арифметических операций. Далее так надо сделать для каждого собственного значения. Вот у Вас уже есть N*N*2 операций доступа к памяти. Итак, тут алгоритм на стек и только чтение по памяти.

 

Так не получится, т.к. QR/QL - это все-таки метод исчерпывания, в котором собственные значения вычисляются в монотонном порядке, причем обычно сначала ловятся бОльшие по абсолютной величине. Не "исчерпав" большее, трудно найти меньшее, т.к. нет уверенности, что итерационный процесс, начатый с разными стартовыми условиями (начальными сдвигами) сойдется к разным собственным значениям, а не к большему среди них.

 

Правильно ли я понимаю, что мы с Вами говорим о методе DSYEVD, а не о DSYEV и DSYEVX? Первый эффективен именно в блочном виде, а второй в блочном почти не реализуем если не считать приведение матрицы к трехдиагональному виду. Я все-таки склоняюсь, что мы говорим именно о DSYEVD...

 

Почитайте пожалуйста, статьи James Demmel и Juan Molera где-то в конце 90ых и начале 2000 они оба очень много об этом на конференциях говорили как можно выполнять эту факторизицию блочно с последующим вычислением собственных векторов. Так как я многое тогда прочитал и сейчас лениво перечитывать, найти на раз не смогу, да и немного лениво, но я точно помню, что они хвалились, что именно так надо делать и в лапак они это встроили. Еще им какой-то аспирант из Вупперталя помогал, но уже не помню его фамилию.

 

Врямя же тянет не итерация, а перестройка собственных векторов, которая вынуждена сопровождать перестроение.

 

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

 

И еще, можно поинтересоваться попросить Вас сделать такой примерчик:

 

возьмите Ваш хорошо оптимизированный алгоритм, полностью сидящий в ФПУстеке, прогоните, пожалуйста, на Вашем компьютере, и посчитайте, пожалуйста, достигнутые гигафлопсы.

 

Теперь возьмите две большие (под 10к) матрицы и вызовите DGEMM из MKL с включенной поддержкой многоядерности (она у него по умолчанию), и сравните гигафлопсы. У меня разница на шестиядерниках зачастую в 10 раз зашкаливает, а ведь DGEMM-то в память обращается, вроде должен быть медленнее.

[Xenia 38]

Правильно ли я понимаю, что мы с Вами говорим о методе DSYEVD, а не о DSYEV и DSYEVX? Первый эффективен именно в блочном виде, а второй в блочном почти не реализуем если не считать приведение матрицы к трехдиагональному виду. Я все-таки склоняюсь, что мы говорим именно о DSYEVD...

Честно говоря, я понятия не имею о том, как работает DSYEVD. :) Я его выбрала по списку из "Reference manual", где про него сказано только то, как под эту функция параметры подсовывать. А про существование DSYEV и DSYEVX, и их отличие от DSYEVD, я даже понятия не имею :). Просветите меня пожалуйста, откуда брать информацию об алгоритмах процедур MKL? В каком документе они описаны?

 

Почитайте пожалуйста, статьи James Demmel и Juan Molera где-то в конце 90ых и начале 2000 они оба очень много об этом на конференциях говорили как можно выполнять эту факторизицию блочно с последующим вычислением собственных векторов. Так как я многое тогда прочитал и сейчас лениво перечитывать, найти на раз не смогу, да и немного лениво, но я точно помню, что они хвалились, что именно так надо делать и в лапак они это встроили. Еще им какой-то аспирант из Вупперталя помогал, но уже не помню его фамилию.

От этих чтениев голова пухнёт. :) В конце концов, мой интерес к библиотеке MKL можно расценивать, как белый флаг капитуляции. Я поняла, что вникнуть во все эти премудрости не смогу. Пусть я один алгоритм победила, другой, но когда-то допущу ошибку. А в случае библиотеки есть все-таки надежда на то, что функция тестировалась умными дяденьками :), а то, что эти дяденьки не заметили, то могли заметить многочисленные пользователи, которые не пройдут мимо ошибки, если за продукт свои деньги платили. Отчасти поэтому я склоняюсь к коммерческим продуктам, а не к опенсорсу, когда непонятно с кого за ошибки спрашивать.

 

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

Нет, не оптимизировала. Ибо не понимаю, чем я могу в эту сторону оптимизировать алгоритм. Я сама вообще брала алгоритмы из ... "Справочника алгоритмов на языке Алгол, линейная алгебра" Райнша и Уилкинсона. Потому как алгоритмы там понятные, а переводить с Алгола на С одно удовольствие, чего про Фортран никак не скажешь. Понимаю, что это очень старая версия алгоритмов, но все предельно ясно написано, а в современной версии LAPACK черт ногу сломит. :)

 

Теперь возьмите две большие (под 10к) матрицы и вызовите DGEMM из MKL с включенной поддержкой многоядерности (она у него по умолчанию), и сравните гигафлопсы. У меня разница на шестиядерниках зачастую в 10 раз зашкаливает, а ведь DGEMM-то в память обращается, вроде должен быть медленнее.

Нема у меня многоядерности, я на Pentium-4 живу. :) Однако сравнивала скорость с Core Dio 2 той же частоты. Впечатление такое, что SSE2 стали на последнем чуть быстрее, но скорость FPU87 упала даже по сравнению с Pentium-4. Поэтому мой (ассемблерный) вариант выигрыша на Core Dio 2 не дает, а MKL-ный вариант становится заметно быстрее.

 

Экстра вопрос: не знаете ли вы, на каком языке программирования эта MKL-библиотека компилировалась? И как она хавает матрицы: уложенные в памяти столбцами (на фортрановский манер) или строками (на С-шный манер)?

[iiv 19]

Честно говоря, я понятия не имею о том, как работает DSYEVD. :) Я его выбрала по списку из "Reference manual", где про него сказано только то, как под эту функция параметры подсовывать. А про существование DSYEV и DSYEVX, и их отличие от DSYEVD, я даже понятия не имею :). Просветите меня пожалуйста, откуда брать информацию об алгоритмах процедур MKL? В каком документе они описаны?

 

На сайте http://www.netlib.org/lapack довольно много и понятно написано о структуре лапака. Именно его функциональность поддерживает MKL. Но в MKL есть еще некоторые функции, которые находятся за пределами лапака.

 

Кстати, во времена Уилкинсона об алгоритме, который запрограммирован в DSYEVD еще ничего не было известно - как я помню, DSYEVD только в середине 80-х придумали какие-то американцы с незапоминающимися китайскими фамилиями. Различие DSYEVD и DSYEV в том, что первый требует дополнительно две копии матрицы в памяти, но имеет меньшую (в разы) арифметическую сложность.

 

Нема у меня многоядерности, я на Pentium-4 живу. :) Однако сравнивала скорость с Core Dio 2 той же частоты. Впечатление такое, что SSE2 стали на последнем чуть быстрее, но скорость FPU87 упала даже по сравнению с Pentium-4. Поэтому мой (ассемблерный) вариант выигрыша на Core Dio 2 не дает, а MKL-ный вариант становится заметно быстрее.

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

 

И как она хавает матрицы: уложенные в памяти столбцами (на фортрановский манер) или строками (на С-шный манер)?

фортрановский манер.

 

Экстра вопрос: не знаете ли вы, на каком языке программирования эта MKL-библиотека компилировалась?

MKL состоит из нескольких частей - BLAS, LAPACK, FFT, Sparse, все остальное.

 

Про все остальное - не знаю как сделано, в основном, из-за того, что мне это не нужно.

 

LAPACK написан на фортране, вернее взят именно тот вариант, который лежит на лапаковском официальном сайте и перекомпилирован. Для С интерфейса написаны враперы.

 

Сам лапак очень интенсивно использует BLAS - вот они-то написаны частично на ассемблере, частично на С, и, как мне известно, фортрана там уже нет. В них-то зарыта вся оптимизация по кешу, конвейеру и многоядерности, более того, МКЛ имеет разные ветки одной и той же функции в зависимости от типа процессора.

 

FFT написан на С, но имеет оптимизации по кешу, конвейеру и многоядерности, частично писанные на ассемблере.

 

Sparse - ядро алгоритма разрабатывалось изначально на фортране Саадом (Миннеаполис) и Болхофером (тогда еще Берлин), потом в 2002 это перенял Шенк (Базель, сейчас Лузана), являющийся ответственным за это в настоящее время, все это недавно было переведено на С, интенсивно дергает BLAS, но имеет несколько специализированных кеш оптимизированных алгоритмов для работы с целочисленными векторами.

 

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

 

Бласы с нетлиба - самые неоптимальные по скорости, вот тут надо найти что-то максимально дешевое или бесплатное, но и максимально быстрое. Как я говорил, я использую ACML - как очень хорошую и бесплатную альтернативу MKL.