Комментарии 33
Почему не хотите использовать С++ библиотеки вроде Eigen (на замену собственной библиотеки на .NET)? Какую библиотеку используете для переупорядочивания вершин?
1. От нативных библиотек было решено отказаться на этапе концептуального проектирования ПВК. Eigen мы используем в другом проекте, где расчетный модуль разрабатывается на С++.
2. Какие вершины вы имеете ввиду?
2. Какие вершины вы имеете ввиду?
Переменные и уравнения системы (они у вас скорее всего соответствуют узлам сетки, на которой решаются дифференциальные уравнения в частных производных).
При решении дифференциальных уравнений в частных производных у нас задача переименования вершин не возникает. Но возникает при моделировании потокораспределения в системах газоснабжения (СГ). На основе графической схемы СГ мы формируем объектный граф. У каждого узла имеется признак, который показывает известны ли в нем параметры газового потока. Далее узлы графа разбиваются на два множества в соответствии с этим признаком обычным linq запросом. И уже для каждого множества узлов строятся два блока матрицы инциденций. На одном из них осуществляется поиск решения. Если будет интересно пользователям, то можно подготовить отдельную публикацию.
Переупорядочивание вершин нужно для уменьшения количества ненулевых элементов в LU разложении. Если только у вас сетка не статическая и вы заранее знаете оптимальный порядок, оно нужно всегда. Ну либо если переупорядочивание приводит к неустойчивости и экспоненциальному увеличению погрешности при увеличении размера матрицы.
Насколько я знаю, это обычно делается внутри LU разложения и называется по-английски pivoting (row и\или column). Например, в BLAS\LAPACK реализации Intel MKL. Не уверен, что есть смысл это делать до передачи матрицы к «решателю». Или вы о чем-то другом?
Pivoting делается для уменьшения погрешности решения, а переупорядочивание — для уменьшения времени работы. Вот, к примеру, замеры времени при использовании библиотеки METIS: matrixprogramming.com/2008/05/metis.
Дифуры решаем на статической сетке, т.к. решаем систему уравнений параболического типа. LU разложение в cuSPARSE является неполным, поэтому и количество ненулевых элементов матриц LU регулируется исходя из параметра tolerance. Так я понимаю. Хотя здесь есть над чем задуматься. Спасибо за замечание, с этим стоит разобраться.
На всякий случай, было бы не плохо проверить, с какой скоростью решают С-шные библиотеки.
Согласен. Однако на это время у нас не было. В диссертацию придется обязательно включать сравнение Eigen.
Кроме того, насколько я помню, Math.Net не имеет вообще ни какого распараллеливания. Вы пробовали использовать AsParallel? Я понимаю, что умножение матриц вот так просто не распараллелить, но вам скорее всего нужно решать задачу не раз. Если у вас есть доступ к Micrsoft Visual Studio Ultimate, то можно поймать какое же действие наиболее прожорливое.
Сразу советую сравнивать не с Eigen, а с приличной реализацией BLAS/LAPACK. Или как минимум дать возможность Eigen использовать эту реализацию — точно есть опция через Intel MKL.
Или я чего-то не понимаю. Или времена счета какие-то странные. Сейчас просто открыл Matlab и написал в нем следующее:
n = 1000; A = rand(n,n); b = rand(n,1); tic; x = A\b; toc
Выход: Elapsed time is 0.118833 seconds.
Макбук, i7, 2.6 ГГц
Одна десятая секунды. На плотной матрице.
P.S. Прошу прощения за неоформленный листинг, видимо кармы моей недостаточно.
n = 1000; A = rand(n,n); b = rand(n,1); tic; x = A\b; toc
Выход: Elapsed time is 0.118833 seconds.
Макбук, i7, 2.6 ГГц
Одна десятая секунды. На плотной матрице.
P.S. Прошу прощения за неоформленный листинг, видимо кармы моей недостаточно.
500x500 — одна сотая секунды. Еще раз: плотная матрица, счет безо всяких GPU.
Еще пара вопросов: что параллелится? LU-факторизация? Вы используете double. Пробовали float? Если да, насколько быстрей, насколько деградирует точность?
Еще пара вопросов: что параллелится? LU-факторизация? Вы используете double. Пробовали float? Если да, насколько быстрей, насколько деградирует точность?
Согласен, 500x500 – вообще смешная размерность.
numpy плотную матрицу 500x500 решает за 3 мс на достаточно старом i5. 5000x5000 – за 2 секунды.
numpy плотную матрицу 500x500 решает за 3 мс на достаточно старом i5. 5000x5000 – за 2 секунды.
In [1]: n = 5000
In [2]: a = np.random.rand(n, n); b = np.random.rand(n, 1)
In [3]: t = time.time(); c = np.linalg.solve(a, b); print time.time() - t
2.17847013474
np.linalg.solve использует процедуру gesv из LAPACK. производительность которой сильно зависит от конкретной матрицы (адаптивном используется одинарная и двойная точность арифметики — как в реализации от Intel). Поэтому замеряя производительность лучше пользоваться другими способами — например прямой вызов (или через интерфейс) пары getrf, getrs из LAPACK, которые вычисляют, соответственно, LU-декомпозицию и решение в той арифметике, которую вы укажете.
Хотя я с вами абсолютно согласен, что ваш пример явно демонстрирует слишком малый размер примера 500 х 500. Нужно подниматься как минимум до 5-10 тысяч.
Хотя я с вами абсолютно согласен, что ваш пример явно демонстрирует слишком малый размер примера 500 х 500. Нужно подниматься как минимум до 5-10 тысяч.
в пример вставил 500x500, т.к. при больших размерностях уже Math.Net неприлично долго считает и графики Mani.Net и ManagedCuda визуально неразличимы.
На моем ноутбуке матлабу потребовалось также около одной сотой секунды. В принципе, ваш результат согласуется с ManagedCuda. В расчетах использую double. К сожалению, float пока не пробовал. Это может дать приличный прирост?
Здесь очень серьезно нужно разбираться в конкретной программе.
Сравнивается прямой метод LU-декомпозиции (сложности порядка O(N^3)) и, если я правильно понимаю, итеративная реализация алгоритма — так как была указана точность решения ( сложность O(N_it * N^2), N_it — число итераций). Количество итераций может меняться в зависимости от обусловленности матрицы. Кроме того, делать стандартную декомпозицию LU для разреженных матриц попросту неэффективно — так как декомпозиция будет плотной матрицей (или как минимум banded). Распараллеливать LU и распараллеливать умножение матрицы на вектор (что используется в итеративных методах) — большая разница. Плюс, к итеративный методам можно и нужно добавлять предобуславливатели (preconditioners), которые уменьшают количество итераций значительно.
Все это сочетает в себе как проблемы решения СЛАУ для конкретной проблемы, так и общую теорию. Хотелось бы больше характеристик того, что именно решалось и как именно решалось. Иначе получается сравнение слонов с носорогами.
Сравнивается прямой метод LU-декомпозиции (сложности порядка O(N^3)) и, если я правильно понимаю, итеративная реализация алгоритма — так как была указана точность решения ( сложность O(N_it * N^2), N_it — число итераций). Количество итераций может меняться в зависимости от обусловленности матрицы. Кроме того, делать стандартную декомпозицию LU для разреженных матриц попросту неэффективно — так как декомпозиция будет плотной матрицей (или как минимум banded). Распараллеливать LU и распараллеливать умножение матрицы на вектор (что используется в итеративных методах) — большая разница. Плюс, к итеративный методам можно и нужно добавлять предобуславливатели (preconditioners), которые уменьшают количество итераций значительно.
Все это сочетает в себе как проблемы решения СЛАУ для конкретной проблемы, так и общую теорию. Хотелось бы больше характеристик того, что именно решалось и как именно решалось. Иначе получается сравнение слонов с носорогами.
В своем посте я допустил огрешность, не указав, что в cuSPARSE используется incomplete-LU factorization.
Все-таки я не ставил себе целью провести сравнение различных технологий, а предлагаю решение задачи решения СЛАУ с помощью GPU под .Net. Приведенное сравнение носит иллюстративный характер. Исходил из того, что если бы мне такая статья попалась сразу, то я сэкономил бы много усилий и времени.
Я подготовлю статью о решаемой мной проблеме. Однако беспокоит, что проблема является весьма специфичной — широкому кругу может показаться неинтересной.
Все-таки я не ставил себе целью провести сравнение различных технологий, а предлагаю решение задачи решения СЛАУ с помощью GPU под .Net. Приведенное сравнение носит иллюстративный характер. Исходил из того, что если бы мне такая статья попалась сразу, то я сэкономил бы много усилий и времени.
Я подготовлю статью о решаемой мной проблеме. Однако беспокоит, что проблема является весьма специфичной — широкому кругу может показаться неинтересной.
А почему просто не использовать итерационные методы рашения? Тот же GMRES.
Нужно узнавать у автора, каковы характеристики его матрицы (и, соответственно, самой проблемы):
Обусловленность (Conditioning) Использование предобуславливателя (PreConditioning) Количество «правых» частей уравнения Необходимая точность решения
Делать слепой общий «решатель» СЛАУ итеративным, не зная характеристик матрицы, проблемы и необходимой точности, довольно опрометчиво.
Обусловленность (Conditioning) Использование предобуславливателя (PreConditioning) Количество «правых» частей уравнения Необходимая точность решения
Делать слепой общий «решатель» СЛАУ итеративным, не зная характеристик матрицы, проблемы и необходимой точности, довольно опрометчиво.
Это я как раз к тому, рассматривались ли другие методы. И есть ли вообще представление о том, какие методы решения СЛАУ бывают, какие у них особенности.
У автора задача оптимизации вроде бы, обычно там высока точность не требуется.
Если размерность системы велика(500 это маленькая система, прямо скажем), то PC(ilu или еще что, может можно придумать «физический» предобуславливатель) + gmres вполне можно и расмотреть как вариант.
У автора задача оптимизации вроде бы, обычно там высока точность не требуется.
Если размерность системы велика(500 это маленькая система, прямо скажем), то PC(ilu или еще что, может можно придумать «физический» предобуславливатель) + gmres вполне можно и расмотреть как вариант.
Да, я что-то не сразу соотнес оптимизацию и требования к точности. Согласен. Хотя для 500 может и никакого PC не нужно. Маааленькое все очень. Просто нужен цивилизованный метод. Вопрос, насколько они собираются расти в размере системы — и тогда асимптотичекая сложность прямого LU в лоб — может быть фатальна.
Я бы, конечно, мог предложить быстрые прямые методы типа иерархических матриц — но здесь это, скорее, overkill.
Вы абсолютно правы, скорее нужно смотреть хорошую реализацию итеративного алгоритма GMRES, BiCGstab (если матрица симметрична или имеет определенные свойства — можно сэкономить на упрощенных реализациях), TFMQR с потенциалом внедрения предобуславливателя. И, разумеется, очень желательно использование библиотеки LAPACK для низкоуровневых операций с векторами и матрицами. Иначе — изобретение велосипеда может получиться.
Я бы, конечно, мог предложить быстрые прямые методы типа иерархических матриц — но здесь это, скорее, overkill.
Вы абсолютно правы, скорее нужно смотреть хорошую реализацию итеративного алгоритма GMRES, BiCGstab (если матрица симметрична или имеет определенные свойства — можно сэкономить на упрощенных реализациях), TFMQR с потенциалом внедрения предобуславливателя. И, разумеется, очень желательно использование библиотеки LAPACK для низкоуровневых операций с векторами и матрицами. Иначе — изобретение велосипеда может получиться.
Нет, в статье я указывал, что оптимизационная задача решается методом динамического программирования, но при переборе вариантов приходится решать системы нелинейных уравнений.
О методах решения СЛАУ автор осведомлен в достаточной степени. Свои изыскания я начинал именно с итерационных методов (метод Гаусса-Зейделя). Однако было принято решение использовать прямые методы.
О методах решения СЛАУ автор осведомлен в достаточной степени. Свои изыскания я начинал именно с итерационных методов (метод Гаусса-Зейделя). Однако было принято решение использовать прямые методы.
Ресь про то, что для опимизационных задач обычно не нужно очень тоное решение СЛАУ, обычно 1E-6 вполне хватает. Естественно желательно обазразмеривать коэффициенты СЛУ чтобы они были, по возможности, одного порядка — для улучшения точности.
Тогда поясните пожалуйста, почему именно было решено использовать именно полное LU разложение.
Метод Гаусса-Зейделя не самый распространенный на сегодняшний день, очень редко когда подходит.
Тогда поясните пожалуйста, почему именно было решено использовать именно полное LU разложение.
Метод Гаусса-Зейделя не самый распространенный на сегодняшний день, очень редко когда подходит.
Автор указал, насколько я понимаю, что для вычислительных экспериментов использовалась случайно заполненная на 10% матрица. О каком диагональном преобладании, кстати, тут можно говорить, мне непонятно. Может быть, конечно, заполнение не совсем случайное, а, скажем, случайное оно только вне главной диагонали и числами из диапазона [-1, 0], а на диагональ идет сумма по строке. Тем не менее, об этом автор умолчал.
Это конечно замечательно, что вы применили ГПУ и воспользовались уже написанной библиотекой, но
1) почему для расчетов используется шарп — ведь по моему ясно, что Си для этих целей куда лучше — и портируемость и поддержка со стороны CUDA/OpenCL/ MPI
2) ваш Mani.net вполне возможно лучше бы работал, если бы был на Си, так как используя CUDA вы выходите на хорошо оптимизированные библиотеки, написанные на Си, так еще и на GPU.
Также может я что то упустил, но не заметил — параллельная ли ваша версия или последовательная? именно то, что сравнивалось на графике. По вашим цифрам получается где то в 16 раз ускорилось.
PS: еще как пожелание — график было бы лучше смотреть, если бы на нем была только ваша версия и GPU, то есть увеличенный масштаб. А то он лишь показывает, что то, с чем вы сравниваетесь — хуже, а конкретно на сколько лучше — не ясно.
1) почему для расчетов используется шарп — ведь по моему ясно, что Си для этих целей куда лучше — и портируемость и поддержка со стороны CUDA/OpenCL/ MPI
2) ваш Mani.net вполне возможно лучше бы работал, если бы был на Си, так как используя CUDA вы выходите на хорошо оптимизированные библиотеки, написанные на Си, так еще и на GPU.
Также может я что то упустил, но не заметил — параллельная ли ваша версия или последовательная? именно то, что сравнивалось на графике. По вашим цифрам получается где то в 16 раз ускорилось.
PS: еще как пожелание — график было бы лучше смотреть, если бы на нем была только ваша версия и GPU, то есть увеличенный масштаб. А то он лишь показывает, что то, с чем вы сравниваетесь — хуже, а конкретно на сколько лучше — не ясно.
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Решение разреженных СЛАУ больших размерностей средствами ManagedCuda в .NET