Оптимизация кернела CUDA

У меня вот проблема.
Алгоритм решает систему линейных уравнений, в цикле. То есть каждый поток делает n итераций, в каждой - решает систему. Система небольшая, каких то 16*16. Матрица хранится в локальной памяти (читай - глобальной). Решает методом Гаусса (приводит к верхней треугольной матрице, и затем обратной прогонкой). В процессе приведения к верхней треуг матрице производится очень много чтений/записей в матрицу (три вложенных цикла). Каждое обращение стоит 400-600 тиков. Программа работает медленно!
При этом в решении нам надо только последние 4 из 16 неизвестных.
Вот код решателя - то где тратится большая часть времени программы:

  1. __device__ void solverKernel(double a[16][16], double b[16], double x[16]){
  2.         double c;
  3.         // k- номер уравнения, вычитаемого из остальных, и номер исключаемого неизвестного
  4.         // i- номер уравнения, из которого исключается неизвестное
  5.         // j- номер столбца
  6.  
  7.     for (int i = 0; i < 16; i++){
  8.                 x[i] = 0;
  9.                 b[i] = 0.0;
  10.         }
  11.         b[0] = -1.0;
  12.  
  13.         for (int k = 0; k < 15; k++){
  14.         for(int i = (k + 1); i < 16; i++){
  15.                         c = a[i][k] / a[k][k];
  16.                         a[i][k] = 0.0;
  17.                         for (int j = (k + 1); j < 16; j++) {
  18.                                 a[i][j] = a[i][j] - c * a[k][j];
  19.                         }
  20.                         b[i] = b[i] - c * b[k];
  21.             }
  22.         }
  23.  
  24.         //   обратная подстановка
  25.         x[15] = b[15] / a[15][15];
  26.        
  27.         for(int i = (14); i >= 0; i--){
  28.                 c = 0.0;
  29.             for(int j = i + 1; j < 16; j++) {
  30.                         c = c + a[i][j] * x[j];
  31.                 }
  32.                 x[i] = (b[i] - c) / a[i][i];
  33.         }
  34. }

Здесь вектор x - вектор неизвестных, из него нам нужны элементы [11...14]
Оптимально решается в 64 потока на блок, при увеличении числа потоков программа работает медленнее (так как дырок всего 8, насколько мне известно, отсюда partition camping).

Задача: изменить решатель СЛАУ так чтобы было меньше обращений к памяти, либо использовался более быстрый тип памяти. Либо изменить алгоритм решения СЛАУ вообще.
Какие предложения?

Я уже пробовал использовать shared память для хранения части матрицы - но она ограничена, не запустишь больше чем 32 потока.
Пробовал использовать shared память как буфер в решателе, но опять таки - конфликты в банках, прироста практически никакого.

Forums: 

пока параллелизм не на уровне

пока параллелизм не на уровне приведенного алгоритма а выше. то есть приведенный код исполняется последовательно во многих потоках много раз.