Узкое место программы в решении СЛАУ. Несколько вариантов:
а) использовать имеющийся алгоритм (уже портировано и работает)
б) найти другой алгоритм, который был бы более оптимизирован для CUDA (меньше обращений к памяти, использование меньшей памяти -> возможность использования общей памяти блока)
в) использовать CUBLAS - что кардинально изменит программу.

Вопросы:
1. Какие алгоритмы вообще существуют - какими вы уже пользовались.
2. Стоит ли переносить все вычисления кроме решения СЛАУ на хост - иначе говоря, каким быстродействием обладает CUBLAS? Скорее всего не стоит, потому что решение системы делается в цикле с количеством итераций 10000 и больше. Но все же, может кто пробовал?