Выполнение нескольких умножений матрицы на матрицу за одну операцию

Я реализую алгоритм, который, по сути, представляет собой серию умножений матрицы на матрицу, например:

Res = M1.M2.M3. ... .Mn

Мои матрицы действительно маленькие поплавки 100x100, но последовательность действительно длинная, порядка миллиардов.

Я пробовал использовать CUBLAS для умножения матриц, но это было медленно, однако я заметил кое-что интересное.

умножение матрицы 100x100 на матрицу 100x100 было медленным, но умножение 1.000.000x100 на 100x100 было относительно быстрым, это заставило меня задуматься. Если бы я вместо того, чтобы сканировать слева направо, имел 10.000 сканирований параллельно. Это должно быть довольно быстро, и если я умножу свои матрицы, когда закончу с этим, я получу тот же результат - только быстрее.

Res1 = M1.M2.M3. ... .Mn/1000-1
Res1 = M1+n/1000.M2+n/1000.M3+n/1000. ... .M2(n/1000)-1
...
Res1  = M1+999*n/1000.M2+999*n/1000.M3+999*n/1000. ... .M1000*(n/1000)-1
Res = Res1*Res2* ... *Res999 

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

Теперь вот моя проблема. Я реализовал матрицу-матрицу (sgemm), вдохновленную той, которую nvidia демонстрирует в своей документации, но она примерно в 4 раза медленнее, чем cublas. Кто-нибудь знает, как работает CUBLAS? А если где-то доступен код?