Здравствуйте, Cyberax, Вы писали:

C>>>Дело в том, что у оптимизатора Фортрана намного больше свободы — там нет указателей (и алиасинга) и ряда других трудностей.
A>>есть там и указатели, и динамическая память (но пользуются там этим скажем так несколько меньше чем new в плюсах), и про equivalence не забудь
C>Ну это явная операция, а в С оно неявно может быть.
явно, неявно... хотите помучать фортрановский оптимизатор на предмет алиасинга (я про аргументы функций молчу)? а запросто —

integer, pointer :: P(:)
integer, target :: A(100)
integer, target :: B(100)
integer N
read (*,10) N
if(N .lt. 42) then
  P=>A
else
  P=>B
end if

! а чёрт его знает, есть тут алиасинг или нет
do i=1,100-N
  P(i) = A(i+N)
enddo

10 format (1I5)
end

ну и опять же man <любимый компилятор фортрана> | grep alias ... неспроста оно там, ой неспроста

C>>>Впрочем, часть проблем с алиасингом решается с помощью ключевого слова restrict в С.
A>>это замечательно выглядит для всяких там тройных указателей
C>Ну не хуже const volatile, наверное
точно но ведь это С(99), а не С++?

Erop пишет:

> 1) Это всё сильно зависит от платформы. Всё-таки всякие разные
> веккторные процессоры ещё никто не отменял.

Да ну естественно я на все случаи жизни не смогу расписать.

> 2) Если речь идёт о i86, и если вспомнить о том, что в выч. хадачах
> обычно все поля одного типа -- double или long double, что само по себе
> зависит от обстоятельств. То можно сообразить что в параллельных
> массивах можно использовать предвычисленное смещение в массиве
> Все эмоциональные формулировки не соотвествуют действительному положению
> вещей и приведены мной исключительно "ради красного словца". За
> корректными формулировками и неискажённым изложением идей, следует
> обращаться к их автором или воспользоваться поиском


Это я не понял. Какие--такие предвычисленные смещения можно
вычислить для параллельных массивов ?

Один массив -- это (кстати, независимо от типа и платформы)
-- адрес начала массива
-- размер элемента (возможно, 2 размера -- с выравниванием и без)
-- кол-во элементов (нужно только для контроля выхода за границы)

Далее мы по этим данным на основе индекса в массиве
вычисляем адрес в памяти элемента. Делается это
путём прибавления к базе I размеров элементов.

Так что тут можно ещё предвычислить, и так одна комманда
машинная.

Хотя наверное это уже малоинтересно.

Cyberax пишет:

>> > Дело в том, что у оптимизатора Фортрана намного больше свободы — там нет
>> > указателей (и алиасинга) и ряда других трудностей.
> MZ>Что такое алиасинг ? А указатели (ссылки) в фортране есть.
> http://en.wikipedia.org/wiki/Aliasing_%28computing%29

Тогда он ЕСТЬ и в фортране. EQUIVALENCE, COMMON BLOCK

> У оптимизатора С++ банально меньше возможностей, так как оптимизатор
> Фортрана имеет право делать предположения, которые для С++ часто будут
> неверными.

Что-то я всё больше сомневаюсь. Ну да ладно.

Здравствуйте, MasterZiv, Вы писали:

>> MZ>Что такое алиасинг ? А указатели (ссылки) в фортране есть.
>> http://en.wikipedia.org/wiki/Aliasing_%28computing%29
MZ>Тогда он ЕСТЬ и в фортране. EQUIVALENCE, COMMON BLOCK
Он _явный_.

Здравствуйте, Mephisto666, Вы писали:

M>Здравствуйте, _Ursus_, Вы писали:

TB>>>Например, массив структур хуже по времени доступа, чем параллельные массивы.
_U_>>Хочется услышать обоснование.

M>Вот занимательно расказано
M>http://szotin.livejournal.com/13335.html

Автор отстал от жизни. Современные компиляторы прекрасно паплайнят циклы. Надо только не забывать restrict вставлять в нужных местах.
Оптимизация кода "параллельными" массивами -- вещь ненужная и вредная, т.к. приводит к плохо структурированному коду.

Копай Нео, копай -- летать научишься. © Matrix. Парадоксы

Cyberax пишет:

> MZ>Тогда он ЕСТЬ и в фортране. EQUIVALENCE, COMMON BLOCK
> Он _явный_.

А неявный это как ? ССылки / указатели ? Чем это плохо для
оптимизатора ? что он этот код не имеет ? Имеет. Что не
может разобраться, что это одно и то же ? может.

А в фортране в общем то же самое твориться.

Здравствуйте, Programador, Вы писали:

D14>>А смысл сравнивать РАЗНЫЕ алгоритмы? И с чем? C Intel MKL, AMD ACML?
P>Так я и сравнивал алгоритмы — мой, с обращением матриц из LAPAk. Мой быстрее.

Попробовал потестировать обращение матриц 2000*2000. Ваш алгоритм оказался порядка 20 раз медленнее Matlab. Подумал, что не так померил что. Сравнил с boost ublas. Не так сильно, на снова Ваш алгоритм медленнее.

matrix<double,column_major> m (n, n);
matrix<double,column_major> invm (n, n);
permutation_matrix<double> pm(m.size1());
invm.assign(identity_matrix<double>(m.size1()));
lu_factorize(m,pm);
lu_substitute(m, pm,invm);

Здравствуйте, Programador, Вы писали:

P>Здравствуйте, Sergey Chadov, Вы писали:

SC>>Как бы это помягче, у меня создалось впечатление, что автор намеренно пытался ввести в заблуждение читателей этого кода. Форматирование и названия переменных просто психоделические, константы типа 1e99 вызывают нервный смех Что так или не так в самом алгоритме сказать не могу — мне жалко мозг разбираться.

P>У меня очень понятное название invers, а их DGETRI вводит меня в ступор, как и все остальные имена, кроме одного EXTERNAL ( а их там не мало ) error message related XERBLA — это по нашему .

Это было придумано черт знает когда, когда современных понятий о качестве кода просто не существовало. И в любом случае то, что где-то еще сделано плохо — не повод делать плохо самому. А насчет понятного названия invers — лично мне логика наименования непонятна, либо мы используем распространенное сокращение "inv", либо пишем слово целиком, "inverse". Впрочем, название функции — это еще далеко не самое страшное.

Здравствуйте, Sergey Chadov, Вы писали:

SC>Это было придумано черт знает когда, когда современных понятий о качестве кода просто не существовало. И в любом случае то, что где-то еще сделано плохо — не повод делать плохо самому. А насчет понятного названия invers — лично мне логика наименования непонятна, либо мы используем распространенное сокращение "inv", либо пишем слово целиком, "inverse". Впрочем, название функции — это еще далеко не самое страшное.

Исторически, Фортран различал только первые 6 символов в названиях функций. Так что что INVERS, что INVERSE для него одно и то же. Соответственно, никто и не заморачивался писать "лишние" символы, которые ни на что не влияют.

Здравствуйте, D14, Вы писали:

D14>Здравствуйте, Programador, Вы писали:

D14>>>А смысл сравнивать РАЗНЫЕ алгоритмы? И с чем? C Intel MKL, AMD ACML?
P>>Так я и сравнивал алгоритмы — мой, с обращением матриц из LAPAk. Мой быстрее.

D14>Попробовал потестировать обращение матриц 2000*2000. Ваш алгоритм оказался порядка 20 раз медленнее Matlab. Подумал, что не так померил что. Сравнил с boost ublas. Не так сильно, на снова Ваш алгоритм медленнее.

Там полный выбор ведущего элемента ( с пропуском обработанных строк/столбцов, помеченных -1 ), обычно ограничиваются первой попавшейся строкой, на него половина времени уходит, но я предпочитаю полный.

  for(j=N;j--;)
  {  if(~rn[j])
     for(k=N;k--;)
       if(~cn[k])
       if(vved<fabs(m[j][k]))
          vved=fabs(m[j][k]),jved=j,kved=k;
  }

И если конкретный матричный шаблон DMat, то он тормознутый . На сайте еще тексты есть, где (i,j) вместо [i][j] , както нет стабильных контейнеров, часто тексты допиливаю по месту. А так хороший алгоритм, быстрый.

Здравствуйте, Programador, Вы писали:

P>Там полный выбор ведущего элемента ( с пропуском обработанных строк/столбцов, помеченных -1 ), обычно ограничиваются первой попавшейся строкой, на
Без выбора ведущего элемента там есть версия без permutation_matrix. Выбор в boost тоже есть.
P>И если конкретный матричный шаблон DMat, то он тормознутый . На сайте еще тексты есть, где (i,j) вместо [i][j] , както нет стабильных контейнеров, часто тексты допиливаю по месту.
"(i,j) вместо [i][j]" оптимизаторы давно как орехи щелкают и инлайнят, выбрасывая временные объекты строк/столбцов.

Здравствуйте, Pzz, Вы писали:


Pzz>Исторически, Фортран различал только первые 6 символов в названиях функций. Так что что INVERS, что INVERSE для него одно и то же. Соответственно, никто и не заморачивался писать "лишние" символы, которые ни на что не влияют.

Это понятно, однако ж это не повод так писать сейчас. Поэтому ИМХО нужно либо писать нормально, либо использовать хоть и уродские, но привычные имена из BLAS/LAPACK. Зачем останавливаться на полпути?

Здравствуйте, MasterZiv, Вы писали:

>> MZ>Тогда он ЕСТЬ и в фортране. EQUIVALENCE, COMMON BLOCK
>> Он _явный_.
MZ>А неявный это как ? ССылки / указатели ? Чем это плохо для
MZ>оптимизатора ? что он этот код не имеет ? Имеет. Что не
MZ>может разобраться, что это одно и то же ? может.
нет, со ссылками/указателями далеко не всегда может (я там даже ниже пример — на фортране — привел). если бы алиасинг был бы только явным, оптимизатор точно бы знал что делать. при неявном (а он всегда будет, хотя бы через параметры функций/указатели, что есть в С, и в фортране) делается дофига проверок, и если в результате всё равно "хз", потенциально опасная оптимизация не делается.

а с equivalence/common block таки да, явно знает (в пределах функции), но кому от этого легче? сам факт наличия equivalence уже может напакостить оптимизации, для которой важно выравнивание (то есть практически все simd'ы)

MZ>А в фортране в общем то же самое твориться.
да, там всё так же плохо

Здравствуйте, MasterZiv, Вы писали:

>> MZ>Тогда он ЕСТЬ и в фортране. EQUIVALENCE, COMMON BLOCK
>> Он _явный_.
MZ>А неявный это как ? ССылки / указатели ?
Да.

MZ>Чем это плохо для оптимизатора ? что он этот код не имеет ? Имеет. Что не
MZ>может разобраться, что это одно и то же ? может.
Он не может разобраться в общем случае, что алиасинга нет и делает всё пессимистично.

Antikrot пишет:

> нет, со ссылками/указателями далеко не всегда может (я там даже ниже
> пример — на фортране — привел). если бы алиасинг был бы только явным,
> оптимизатор точно бы знал что делать. при неявном (а он всегда будет,
> хотя бы через параметры функций/указатели, что есть в С, и в фортране)
> делается дофига проверок, и если в результате всё равно "хз",
> потенциально опасная оптимизация не делается.

Родной, common block и EQUIVALENCE такого могут нагородить,
что там чёрт ногу сломит. Оно может произвольно в принципе
накладывать байтики одного массива/структуры/объекта на
другой. О чём ты говоришь тогда ?

Всё, убеждать устал. БОльше не буду.

> а с equivalence/common block таки да, явно знает (в пределах функции),
> но кому от этого легче? сам факт наличия equivalence уже может

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

Хотя кажется, ты мой сторонник тогда извини если что.

Здравствуйте, D14, Вы писали:

D14>Здравствуйте, Programador, Вы писали:

Так значит матрица 2000х2000 весь кэш поела, а размерчик актуальным для меня скоро станет . Правда для отладки можно тренироваться "на кошечках"

А подробней не расскажете? Сколько в минутах выходило, матлаб чем пользуется? AMD ACML мне кажется бесплатно, он к minGW линкуется? Кроме интел, амд еще кто-то может большие матрицы бороть?

Здравствуйте, Programador, Вы писали:
P>А подробней не расскажете? Сколько в минутах выходило, матлаб чем пользуется?
На Celeron 1700:
-boost ublas partial pivoting 151 сек.
-(c) Drobotenko 213 сек.
-Matlab 19 сек.
Насчет Matlaba, оказалось, я держал в голове цифры меньшей для размерности, когда писал про 20х преимущество
P> AMD ACML мне кажется бесплатно, он к minGW линкуется?
Не в курсе.
P> Кроме интел, амд еще кто-то может большие матрицы бороть?
Можно поискать по blas vendors, но по-моему под win32 только MKL и ACML в ходу. Последние версии IMSL дергают тот же MKL внутри. Матлаб использует их обе.

Здравствуйте, D14, Вы писали:

D14>На Celeron 1700:
D14>-boost ublas partial pivoting 151 сек.
D14>-(c) Drobotenko 213 сек.
D14>-Matlab 19 сек.

Странно у меня на Доу-2 2.13гг 233 сек, с полной оптимизацией. Но возможно не стоит осуждать GCC поскольку в DMat стоял мой range-chack с инлайн асм int 03.

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

Переход на указатели во внутреннем цикле + частичный поиск ведущего — 35
Переход на указатели во внутреннем цикле — 63 сек

против исходного при полной оптимизации — 233


Кэш в 10 раз ускоряе http://www.ixbt.com/cpu/rmma/yonah/ynh_mem_bw.png из http://www.ixbt.com/cpu/rmma-yonah.shtml , MKL конечно быстрая, но чудес не показывает. Думал если в 10 раз быстрее, может как-то блочно умеют влазить в кэш. Пока без них буду жить

     inline static  void  invers(DMat m)
     {
       int N= dimA(m);
       int *rn=(int*)alloca((sizeof(int))*N),*cn=(int*)alloca((sizeof(int))*N);
       //int rn[N],cn[N];

       int j,k;

       for(j=N;j--;)
         rn[j]=cn[j]=j;
       gaus_minved=1e99;
       gaus_deter=1;

         //  первый поск ведущего элемента
          double vved=-1,t;
          int next_jved,next_kved;
          for(j=N;j--;)
          for(k=N;k--;)
          if(vved<fabs(m[j][k]))
                  vved=fabs(m[j][k]),next_jved=j,next_kved=k;



       for(gaus_ostatok=N;gaus_ostatok;gaus_ostatok--)
       {  int jved=next_jved,kved=next_kved;
          if(gaus_minved>vved)
             gaus_minved=vved;
          gaus_deter*=m[jved][kved];
          if(vved<1e-99)
          {  for(j=N;j--;)
             {  if(~rn[j]) for(k=N;k--;)
                  m[j][k]=_NaN();
                if(~cn[j]) for(k=N;k--;)
                  m[k][j]=_NaN();
             }
             return;
          }
          int jt=rn[jved],kt=cn[kved];
          for(j=N;j--;)
              t=m[kt][j],m[kt][j]=m[jved][j],m[jved][j]=t;
          for(j=N;j--;)
              t=m[j][jt],m[j][jt]=m[j][kved],m[j][kved]=t;

          rn[jved]=rn[kt];
          cn[kved]=cn[jt];
          rn[kt]=cn[jt]=-1;

          vved=m[kt][jt];   m[kt][jt]=1;
          double next_vved=-1;
          int notvse=5;
          for(j=N;j--;)
          {  if(j==kt)
               continue;
             double mul=m[j][jt]/vved;
             m[j][jt]=0;


             #if 0
             for(k=N;k--;)
             //for(k=0;k<N;k++)
               m[j][k]-=m[kt][k]*mul;
             #else
             //   самый внутренний цикл через указатели
             double *mj=&m[j][0],*mkt=&m[kt][0];
             for(k=0;k<N;k++)
               mj[k]-=mkt[k]*mul;
             #endif

             #if 0
             // ограничение числа строк на поск следующего ведущего элемента
             if(notvse>0)
             #endif
             {
             if(~rn[j])
             { notvse--;
               for(k=N;k--;)
               //for(k=0;k<N;k++)
               if(~cn[k])
               if(next_vved<fabs(m[j][k]))
                  next_vved=fabs(m[j][k]),next_jved=j,next_kved=k;
             }
             }
          }
          for(k=N;k--;)
              m[kt][k]/=vved;
          vved=next_vved;

       }
     }

Здравствуйте, Programador, Вы писали:

P>Здравствуйте, D14, Вы писали:

D14>>На Celeron 1700:
D14>>-boost ublas partial pivoting 151 сек.
D14>>-(c) Drobotenko 213 сек.
D14>>-Matlab 19 сек.

P>Странно у меня на Доу-2 2.13гг 233 сек, с полной оптимизацией. Но возможно не стоит осуждать GCC поскольку в DMat стоял мой range-chack с инлайн асм int 03.
[/ccode]
Дык потому, что Intel C++. С выкинутым прочь ASSERT 130 сек.

Два аспекта:
* код целиком создаете сами;
* интенсивно пользуеетесь сторонними библиотеками.

1) Формально-то можно написать "одинаковый" по скорости. Но сильно сомневаюсь, что обычный програмер готов лезть в доки и изучать особенности архитектуры процессора (редкостные регистры и пр.). А при одинаковом алгоритме фортрановский компилятор лучше оптимизирОВАЛ (последние годы не в теме) индексные выражения. Т.е. фортрановский код чаще будет более быстрым.

2) В фортранАХ мощнейшие библиотеки по вычметодам. И если использовать их — рядового С-ишника "обыграешь" легко.
Приятель годами работает в Интеле, как раз на библиотеках. Видел примеры кода. ЭТО — не фортран и не С++! Скорей, суперпродвинутый ассемблер. Отдельные версии функций для 2/4 ядер... Но, подчеркну, библиотеки ДЛЯ фортрана.


А>Добрый день,
А>Разные люди говорят что вычислительные программы работают быстрее если они были написаны на Fortran-е а не на C++.
А>Сам ответа на этот вопрос не знаю и отношусь к этому скептически.
А>Что можете сказать по этому поводу.
А>С Уважением.